Introduction to ERGMs

Let’s dig into exponential random graph models in R.

Recall Robins and Lusher (2013:13): “Exponential random graph models (ERGMs) are statistical models for network structure permitting inferences about how network ties are patterned. Put it another way, ERGMs are tie-based models for understanding how and why social network ties arise.”

ERGMs are the most popular ways of thinking about the factors that led to an observed network coming about for several reasons: 1. computationally OK 2. allow for theory-testing (often highlight competing theories) 3. easy to implement in R with good support materials.

For example, you may wonder what factors lead to romantic encounters between fictional medical professionals. What is the effect of monogamy in this network (or at least the effect of someone who is only willing to date one colleague)? What is the effect of racial homophily?

We will begin with a network based on Grey’s Anatomy Romantic Encounters. This case was developed by Gary Weissman and Benjamin Lind years ago here. As the network is dated, I updated by tracing this updated graph here. Here, we include one-time hook ups to marriages/divorces.

In addition to the code developed by Lind, I also follow along with Luke (2015), Goordeau et al. (2008), among others.

Luke, D. A. (2015). A user’s guide to network analysis in R. London: Springer.

Goodreau, S. M., Handcock, M. S., Hunter, D. R., Butts, C. T., & Morris, M. (2008). A statnet tutorial. Journal of statistical software, 24(9), 1.

The first thing that we have to do as always is bring in our packages for the workshop. We use igraph to read the graph and plot it in R. And we use intergraph to move into statnet and statnet for making exponential random graph models.

library(igraph)
library(intergraph)
library(statnet)

Next, let’s load the data and explore in igraph.

greys_graph <- readRDS("~/Documents/GitHub/soc613spr2023/lectures/009_ergms/data/greys_graph.RDS")

Let’s plot in igraph in plot.igraph.

plot.igraph(greys_graph, layout=layout_with_kk, edge.arrow.size=.1,     
            vertex.color=as.factor(V(greys_graph)$gender), vertex.label.cex=.75, 
            vertex.size=igraph::degree(greys_graph))

Now we can bring into statnet using asNetwork in intergraph.

greys <- asNetwork(greys_graph)

And we can use plot in statnet to plot it.

la <- plot(greys)

Let’s begin to observe differences in gender by looking at degree distributions.

We can use the degree function in statnet. As the network is undirected, we use degree(g)/2.

degree(greys)
##  [1] 18  8  8 12  6 12 26  6 20 10 14  4 10  6 12  6 10  8  4 10 10 16  4  4  6
## [26]  4 10  6  4  6  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  4  2
## [51]  2  2  2  2  2  2  2  2  2  2  4  2  2  2  2  2  2  2  2  2  2  2  2  2  2
## [76]  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2
gry.degreedist <- table(degree(greys)/2)

gry.degreedist 
## 
##  1  2  3  4  5  6  7  8  9 10 13 
## 58  8  7  3  6  3  1  1  1  1  1
#Compare degree distribution by sex

knitr::kable(summary(greys ~ degree(5:8, "gender")))
x
deg5.genderF 6
deg6.genderF 2
deg7.genderF 1
deg8.genderF 0
deg5.genderM 0
deg6.genderM 1
deg7.genderM 0
deg8.genderM 1
deg5.genderNB 0
deg6.genderNB 0
deg7.genderNB 0
deg8.genderNB 0

We can count triangles using summary function and triangle in statnet.

summary(greys~triangle)
## triangle 
##        0

We can look at the mixing matrices by sex and race using the mixingmatrix function in statnet.

mixingmatrix(greys, "gender")
##     F  M NB
## F   9 90  1
## M  90  1  0
## NB  1  0  0
mixingmatrix(greys, "race")
##       black other white
## black     8     6    10
## other     6     0    14
## white    10    14    63

Fitting an ERGM: Null Model

Now we can start fitting an ERGM. Let’s start with the null model or the equivalent of running a regression with just its intercept. The null model for ergms is called Bernoulli model where the is an equal probability for all edges in a network. We are basically adding the number of edges in addition to the number of nodes as a factor in the model.

There are numerous network terms that we use in these equations from statnet. The term for this baseline model is edges.

So the function will look like this:

ergm(g ~ edges)

model1 <- ergm(greys ~ edges)
## Starting maximum pseudolikelihood estimation (MPLE):
## Evaluating the predictor and response matrix.
## Maximizing the pseudolikelihood.
## Finished MPLE.
## Stopping at the initial estimate.
## Evaluating log-likelihood at the estimate.
summary(model1)
## Call:
## ergm(formula = greys ~ edges)
## 
## Maximum Likelihood Results:
## 
##       Estimate Std. Error MCMC % z value Pr(>|z|)    
## edges  -3.6546     0.1008      0  -36.26   <1e-04 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##      Null Deviance: 5552.1  on 4005  degrees of freedom
##  Residual Deviance:  942.8  on 4004  degrees of freedom
##  
## AIC: 944.8  BIC: 951.1  (Smaller is better. MC Std. Err. = 0)

The coefficient is difficult to make sense of, so we can use the plogis function to exponentiate the coefficient and provide a probability of there being a tie. We can see that this is exactly the same as network density using network.density in statnet.

plogis(coef(model1))
##      edges 
## 0.02521848
network.density(greys)
## [1] 0.02521848

There is a .025 probability of their being a tie between two nodes in this graph which is the same as saying that the graph has a density of .025.

Fitting an ERGM: Exogenous Factors

ERGM fits a logistic regression in this case as dyads are independent of one another.

But we know this network is nonrandom and other factors likely contribute to network structure.

First, let’s add node attributes.

In statnet we use nodefactor() for categorical and nodecov() for ratio/interval variables.

We may wonder if women are more likely to form romantic ties in the Grey’s Anatomy network.

model2 <- ergm(greys ~ edges + nodefactor("gender"))
## Starting maximum pseudolikelihood estimation (MPLE):
## Evaluating the predictor and response matrix.
## Maximizing the pseudolikelihood.
## Finished MPLE.
## Stopping at the initial estimate.
## Evaluating log-likelihood at the estimate.
summary(model2)
## Call:
## ergm(formula = greys ~ edges + nodefactor("gender"))
## 
## Maximum Likelihood Results:
## 
##                      Estimate Std. Error MCMC % z value Pr(>|z|)    
## edges                 -3.4963     0.1667      0 -20.979   <1e-04 ***
## nodefactor.gender.M   -0.1527     0.1442      0  -1.058     0.29    
## nodefactor.gender.NB  -0.9084     1.0130      0  -0.897     0.37    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##      Null Deviance: 5552.1  on 4005  degrees of freedom
##  Residual Deviance:  940.8  on 4002  degrees of freedom
##  
## AIC: 946.8  BIC: 965.7  (Smaller is better. MC Std. Err. = 0)

The variable is not significant and therefore gender doesn’t seem to matter for this network by itself.

We can evaluate the maximum likelihood (higher is better).

model1$mle.lik[1]
## [1] -471.4136
model2$mle.lik[1]
## [1] -470.3988

And we can see that the performance is not much better for model 2.

Now, let’s look first at potential homophilous effects using nodematch().

First, let’s look at gender homophily.

model3 <- ergm(greys ~ edges + nodematch("gender"))
## Starting maximum pseudolikelihood estimation (MPLE):
## Evaluating the predictor and response matrix.
## Maximizing the pseudolikelihood.
## Finished MPLE.
## Stopping at the initial estimate.
## Evaluating log-likelihood at the estimate.
summary(model3)
## Call:
## ergm(formula = greys ~ edges + nodematch("gender"))
## 
## Maximum Likelihood Results:
## 
##                  Estimate Std. Error MCMC % z value Pr(>|z|)    
## edges             -3.0790     0.1072      0 -28.718   <1e-04 ***
## nodematch.gender  -2.1816     0.3347      0  -6.519   <1e-04 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##      Null Deviance: 5552.1  on 4005  degrees of freedom
##  Residual Deviance:  871.8  on 4003  degrees of freedom
##  
## AIC: 875.8  BIC: 888.4  (Smaller is better. MC Std. Err. = 0)
exp(coef(model3))
##            edges nodematch.gender 
##       0.04600607       0.11285703

This is a pretty heteronormative show, so we would expect that same-sex relationships would be rare.

This is indeed the case as the matching on sex reduces the likelihood of having a tie by about 89%.

model1$mle.lik[1]
## [1] -471.4136
model3$mle.lik[1]
## [1] -435.8813

We are still modeling logistic regressions as network configuration terms are not included (e.g. our model of tie formation does not depend on other ties).

Let’s look at race nodematch(“race”) as well to see two main competing exogenous factors.

model3a <- ergm(greys ~ edges + nodematch("race"))

summary(model3a)
## Call:
## ergm(formula = greys ~ edges + nodematch("race"))
## 
## Maximum Likelihood Results:
## 
##                Estimate Std. Error MCMC % z value Pr(>|z|)    
## edges           -4.0910     0.1841      0 -22.222  < 1e-04 ***
## nodematch.race   0.6992     0.2201      0   3.177  0.00149 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##      Null Deviance: 5552  on 4005  degrees of freedom
##  Residual Deviance:  932  on 4003  degrees of freedom
##  
## AIC: 936  BIC: 948.5  (Smaller is better. MC Std. Err. = 0)
exp(coef(model3a))
##          edges nodematch.race 
##     0.01672241     2.01222749

Homophily by race doubles the likelihood of edge formation in this graph.

Let’s look at both.

model3b <- ergm(greys ~ edges + nodematch("gender") + nodematch("race"))

summary(model3b)
## Call:
## ergm(formula = greys ~ edges + nodematch("gender") + nodematch("race"))
## 
## Maximum Likelihood Results:
## 
##                  Estimate Std. Error MCMC % z value Pr(>|z|)    
## edges             -3.5121     0.1882      0 -18.663  < 1e-04 ***
## nodematch.gender  -2.1790     0.3348      0  -6.508  < 1e-04 ***
## nodematch.race     0.6928     0.2216      0   3.127  0.00177 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##      Null Deviance: 5552.1  on 4005  degrees of freedom
##  Residual Deviance:  861.2  on 4002  degrees of freedom
##  
## AIC: 867.2  BIC: 886.1  (Smaller is better. MC Std. Err. = 0)

Note that when including nodematch terms one also would often include nodefactor as well (“main effects” vs.”interaction effects”).

Fitting an ERGM: Endogenous Terms

Let’s add network configurations to this model. A list of ergm terms is located here.

The first thing that we might be interested in is a monogamy effect. We can find that by evaluating the effect of having a degree of 1 or the network term degree(1) and this could look like this: model4 <- ergm(greys ~ edges + degree(1)) if we included no other terms.

Unfortunately, this model fails to converge. So we try again. As our model includes the match terms, let’s keep them. Nested models here aren’t as interpretable as in canonical regression due to the interdependence of model terms. Endogenous terms like degree(1) can work, but doesn’t with these data. Prior research show that four-cycles are rare in most relationship networks and we can use cycle(4) to evaluate. 4-cycles are square subgraphs where each node has two alters and indicates local clustering.

Unfortunately, the model including the 4-cycle doesn’t work. So we can try some other network terms. The triangle and the other clustering variables don’t work very well. This is likely again due to the lack of closure between triads due to heteronormative norms around dating friends exes.

We can evaluate the other side of this equation by turning to 2-stars (a node sends to two other nodes with no connection between them (e.g. an intransitive triad)).

It turns out that the endogenous terms aren’t playing nicely in this case. This is common and they have developed some parameters to help. First, we will add the gwesp (geometrically weighted edgewise shared partners) a measure of local clustering. This includes a decay function discounting parters beyond the first. The closer to 0, the more dramatic the discounting.

model4a <- ergm(greys ~ edges + gwesp(decay=.5, fixed=TRUE) + nodematch("gender")                 + nodematch("race"), control=control.ergm(seed=50))
## Starting maximum pseudolikelihood estimation (MPLE):
## Evaluating the predictor and response matrix.
## Warning in mple.existence(pl): The MPLE does not exist!
## Maximizing the pseudolikelihood.
## Warning in ergm.mple(nw, fd, m, MPLEtype = MPLEtype, init = init, control =
## control, : glm.fit: fitted probabilities numerically 0 or 1 occurred
## Finished MPLE.
## Starting Monte Carlo maximum likelihood estimation (MCMLE):
## Iteration 1 of at most 60:
## Model statistics 'gwesp.fixed.0.5' are not varying. This may indicate that the observed data occupies an extreme point in the sample space or that the estimation has reached a dead-end configuration.
## Warning: 'glpk' selected as the solver, but package 'Rglpk' is not available;
## falling back to 'lpSolveAPI'. This should be fine unless the sample size and/or
## the number of parameters is very big.
## Optimizing with step length 1.0000.
## The log-likelihood improved by 0.1058.
## Estimating equations are not within tolerance region.
## Iteration 2 of at most 60:
## Model statistics 'gwesp.fixed.0.5' are not varying. This may indicate that the observed data occupies an extreme point in the sample space or that the estimation has reached a dead-end configuration.
## Optimizing with step length 1.0000.
## The log-likelihood improved by 0.0158.
## Convergence test p-value: 0.0001. Converged with 99% confidence.
## Finished MCMLE.
## Evaluating log-likelihood at the estimate. Fitting the dyad-independent submodel...
## Bridging between the dyad-independent submodel and the full model...
## Setting up bridge sampling...
## Using 16 bridges: 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 .
## Bridging finished.
## This model was fit using MCMC.  To examine model diagnostics and check
## for degeneracy, use the mcmc.diagnostics() function.
summary(model4a)
## Call:
## ergm(formula = greys ~ edges + gwesp(decay = 0.5, fixed = TRUE) + 
##     nodematch("gender") + nodematch("race"), control = control.ergm(seed = 50))
## 
## Monte Carlo Maximum Likelihood Results:
## 
##                  Estimate Std. Error MCMC % z value Pr(>|z|)    
## edges             -3.4798     0.1823      0 -19.090  < 1e-04 ***
## gwesp.fixed.0.5   -4.7845         NA     NA      NA       NA    
## nodematch.gender  -2.1285     0.3251      0  -6.547  < 1e-04 ***
## nodematch.race     0.6806     0.2154      0   3.160  0.00158 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##      Null Deviance: 5552.1  on 4005  degrees of freedom
##  Residual Deviance:  859.1  on 4001  degrees of freedom
##  
## AIC: 867.1  BIC: 892.3  (Smaller is better. MC Std. Err. = 0.2803)

And now we can use geometrically weighted degree as a popularity measure. gwdegree estimates the change in the likelihood of a tie as a function of the degree of the involved nodes with a decay penalizing increases in degree. Negative effects indicate popularity is an important factor, while positive estimates indicates that the edges are more spread out through the graph.

model4b <- ergm(greys ~ edges + gwdegree(decay = 1, fixed = TRUE) + nodematch("gender") + nodematch("race"), control=control.ergm(seed=50))
## Starting maximum pseudolikelihood estimation (MPLE):
## Evaluating the predictor and response matrix.
## Maximizing the pseudolikelihood.
## Finished MPLE.
## Starting Monte Carlo maximum likelihood estimation (MCMLE):
## Iteration 1 of at most 60:
## Optimizing with step length 0.8179.
## The log-likelihood improved by 3.1556.
## Estimating equations are not within tolerance region.
## Iteration 2 of at most 60:
## Optimizing with step length 1.0000.
## The log-likelihood improved by 1.7014.
## Estimating equations are not within tolerance region.
## Iteration 3 of at most 60:
## Optimizing with step length 1.0000.
## The log-likelihood improved by 0.2533.
## Estimating equations are not within tolerance region.
## Iteration 4 of at most 60:
## Optimizing with step length 1.0000.
## The log-likelihood improved by 0.0564.
## Convergence test p-value: 0.7493. Not converged with 99% confidence; increasing sample size.
## Iteration 5 of at most 60:
## Optimizing with step length 1.0000.
## The log-likelihood improved by 0.0163.
## Convergence test p-value: 0.4582. Not converged with 99% confidence; increasing sample size.
## Iteration 6 of at most 60:
## Optimizing with step length 1.0000.
## The log-likelihood improved by 0.0226.
## Convergence test p-value: 0.0057. Converged with 99% confidence.
## Finished MCMLE.
## Evaluating log-likelihood at the estimate. Fitting the dyad-independent submodel...
## Bridging between the dyad-independent submodel and the full model...
## Setting up bridge sampling...
## Using 16 bridges: 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 .
## Bridging finished.
## This model was fit using MCMC.  To examine model diagnostics and check
## for degeneracy, use the mcmc.diagnostics() function.
summary(model4b)
## Call:
## ergm(formula = greys ~ edges + gwdegree(decay = 1, fixed = TRUE) + 
##     nodematch("gender") + nodematch("race"), control = control.ergm(seed = 50))
## 
## Monte Carlo Maximum Likelihood Results:
## 
##                  Estimate Std. Error MCMC % z value Pr(>|z|)    
## edges             -2.4959     0.3277      0  -7.617  < 1e-04 ***
## gwdeg.fixed.1     -1.0687     0.3057      0  -3.496 0.000473 ***
## nodematch.gender  -2.1607     0.3348      0  -6.453  < 1e-04 ***
## nodematch.race     0.5351     0.1978      0   2.705 0.006827 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##      Null Deviance: 5552.1  on 4005  degrees of freedom
##  Residual Deviance:  851.2  on 4001  degrees of freedom
##  
## AIC: 859.2  BIC: 884.3  (Smaller is better. MC Std. Err. = 0.3963)

And we can quickly check model improvement

lls <- data.frame(model1=model1$mle.lik[1], model3b=model3b$mle.lik[1], model4a=model4a$mle.lik[1], model4b=model4b$mle.lik[1])

knitr::kable(lls, caption="Comparing Models")
Comparing Models
model1 model3b model4a model4b
-471.4136 -430.6248 -429.5477 -425.5815

Fitting an ERGM: Model Fit

We have been concerned with model improvement, but should actually have been also considering model fit. We specifically want to know how well properties of the simulated model fit with the observed network. This gives us confidence in the model itself. The technique for doing this takes advantage of a unique property of our models. We modeled local properties of networks with the theory that local properties aggregate into the network as a whole. So, no whole network statistics were used to model our graph. We can use network-wide statistics to compare the simulated and observed networks.

We would expect that the observed network statistics would be in the confidence intervals of the simulated networks.

The gof function gets us the goodness-of-fit statistics and plot those statistics in statnet.

gof4b <- gof(model4b)

plot(gof4b)

And we can see that the observed statistics are generally well within the intervals. The model fits the data.

We should always also check for degeneracy. We want to make sure that the simulation process worked. We have some indication that it did as the models actually spit something out, but we should check and make sure.

The estimates should be normally distributed and centered on 0. The mcm.diagnostics function will give us the information we need especially a series of graphs of the covariates. We will see the normal distribution (or not) on the left side of the graphs. On the right side we will see a plot of the “noise” around 0 and it should look like almost like a solid blue horizontal line with jagged edges. If the peak of the distribution (right side) is very different than 0 or if the plot of a ratio variable has multiple peaks (right side) or if the left side has a slope or if the blue ban is very narrow, then you likely have degeneracy problems.

mcmc.diagnostics(model4b)

## Sample statistics summary:
## 
## Iterations = 18944:368640
## Thinning interval = 512 
## Number of chains = 1 
## Sample size per chain = 684 
## 
## 1. Empirical mean and standard deviation for each variable,
##    plus standard error of the mean:
## 
##                     Mean     SD Naive SE Time-series SE
## edges            -3.2573 16.857   0.6445         1.2506
## gwdeg.fixed.1    -2.8080 14.548   0.5563         1.0151
## nodematch.gender -0.5205  3.476   0.1329         0.1764
## nodematch.race   -2.5439 12.808   0.4897         0.9615
## 
## 2. Quantiles for each variable:
## 
##                    2.5%    25%   50%   75% 97.5%
## edges            -32.92 -15.00 -4.00 8.000 32.00
## gwdeg.fixed.1    -30.13 -12.49 -2.92 7.294 25.67
## nodematch.gender  -6.00  -3.00 -1.00 2.000  7.00
## nodematch.race   -26.00 -11.00 -3.00 6.250 23.00
## 
## 
## Are sample statistics significantly different from observed?
##                   edges gwdeg.fixed.1 nodematch.gender nodematch.race
## diff.      -3.257309942  -2.807996019     -0.520467836   -2.543859649
## test stat. -2.604581405  -2.766161763     -2.950100290   -2.645661684
## P-val.      0.009198657   0.005672039      0.003176708    0.008153131
##                (Omni)
## diff.              NA
## test stat. 11.8851570
## P-val.      0.0202743
## 
## Sample statistics cross-correlations:
##                      edges gwdeg.fixed.1 nodematch.gender nodematch.race
## edges            1.0000000     0.9725273        0.5064192      0.9127881
## gwdeg.fixed.1    0.9725273     1.0000000        0.4818725      0.8640597
## nodematch.gender 0.5064192     0.4818725        1.0000000      0.4760833
## nodematch.race   0.9127881     0.8640597        0.4760833      1.0000000
## 
## Sample statistics auto-correlation:
## Chain 1 
##              edges gwdeg.fixed.1 nodematch.gender nodematch.race
## Lag 0    1.0000000    1.00000000      1.000000000     1.00000000
## Lag 512  0.5797856    0.53762859      0.134112576     0.58754937
## Lag 1024 0.3549752    0.31039471      0.074991464     0.37069075
## Lag 1536 0.2282492    0.19518648      0.103251773     0.24775752
## Lag 2048 0.1375828    0.11737799      0.007290324     0.16039749
## Lag 2560 0.0511338    0.04773588      0.004745570     0.06025093
## 
## Sample statistics burn-in diagnostic (Geweke):
## Chain 1 
## 
## Fraction in 1st window = 0.1
## Fraction in 2nd window = 0.5 
## 
##            edges    gwdeg.fixed.1 nodematch.gender   nodematch.race 
##        0.4045548        0.1517573        0.8581428        0.5937843 
## 
## Individual P-values (lower = worse):
##            edges    gwdeg.fixed.1 nodematch.gender   nodematch.race 
##        0.6858048        0.8793784        0.3908136        0.5526564 
## Joint P-value (lower = worse):  0.4674514 
## 
## Note: MCMC diagnostics shown here are from the last round of
##   simulation, prior to computation of final parameter estimates.
##   Because the final estimates are refinements of those used for this
##   simulation run, these diagnostics may understate model performance.
##   To directly assess the performance of the final model on in-model
##   statistics, please use the GOF command: gof(ergmFitObject,
##   GOF=~model).

These are distributed around 0 and so we can be confident that our models aren’t degenerate. Note that the nodematch.sex is bumpy because of it consists of a few categories.

Let’s return to the covariates.

It is common to talk about results of an ergm by focusing on the statistical significance and magnitude of the estimates (b). For example, we can see that there is a large negative effect of gender homophily suggesting that different-sex partnerships are common in this network. Workplace dramas, like Grey’s Anatomy, feed of the overlapping romantic interests of main characters. However, it makes sense that drama on the margins is less useful. Keeping track of tertiary conflicts becomes cumbersome and many romantic relationships center on a few characters. We see that the narrative follows this logic as the geometrically weighted degree parameter is negative indicating that popularity is playing a role. Racial homophily is not as strong, but is strong consistent with racialized presentations of heternormative sexuality on tv.

We can also discuss the probablities for tie formation of our significant variables of interest (here, popularity, racial homophily, and gender homophily).

Remember that we can use plogis() to exponentiate.

plogis(coef(model4b))
##            edges    gwdeg.fixed.1 nodematch.gender   nodematch.race 
##       0.07614955       0.25565529       0.10333865       0.63067189

Gwdegree is difficult to interpret in an exponentiated form, but indicates that popularity plays a role in the network. The graph is heternormative as the gender factor is very small. The probability of gender homophily is .10. And racial homophily is an important positive factor with a probability of .63.