Introduction

In this lab, we will be testing a number of hypotheses about a network’s structure using exponential random graph modeling (ERGM) techniques using the statnet package in R.11 Mark S. Handcock, David R. Hunter, Carter T. Butts, Steven M. Goodreau, and Martina Morris (2003). statnet: Software tools for the Statistical Modeling of Network Data. statnetproject.org; see also ?? statnet. For more information about ERGMs, see generally D. Lusher, J. Koskinen, & G. Robins (2012) Exponential Random Graph Models for Social Networks. statnet provides a comprehensive framework for ERGM-based network modeling, including tools for model estimation, model evaluation, model-based network simulation, and network visualization. This functionality is powered by a central Markov chain Monte Carlo Maximum Likelihood Estimation (MCMCMLE) algorithm.22 For a great introduction to MCMC grounded in graph theory, see Jeremy Kun, Markov Chain Monte Carlo Without All the Bullshit.

statnet resources:

● Developer website

● User guide

● Tutorial

Data

We will analyze the communication behaviors within a team of seventeen members who were involved in designing military installations.

• CRIeq.txt: each team member’s communication to retrieve information from other team members on the topic of environmental quality (eq). This is a directed, binary relation.

• CAIeq.txt: each team member’s communication to allocate information to other team members on the topic of environmental quality (eq). This is a directed, binary relation.

• EXeq_cons.txt: each team member’s expertise on the topic of environmental quality (eq) as perceived on average by all team members. This is a continuous attribute.

Hypotheses

We will test various hypotheses based on the Theory of Transactive Memory.33 See Monge & Contractor (2003) Theories of Communication Networks, 198—203.

Hypothesis 1: Individuals are less likely to retrieve information from those who retrieve information from them.

Hypothesis 2a: Information retrieval tends to be transitive. That is, if individual i retrieves information from individual k, and individual k retrieves information from individual j, then individual i is more likely to retrieve information from individual j.

Hypothesis 2b: Transitivity increases at a sub-linear rate as a function of the number of ties.

Hypothesis 3a: Individuals tend to retrieve information from other members with high expertise.

Hypothesis 3b: Individuals with low expertise tend to retrieve information from many others.

Hypothesis 4: Individuals tend to retrieve information from members to whom they allocate information.

Part I. Building & Visualizing the Networks (30 pts)

The analysis will use three files: the CRIeq.txt as the network file, EXeq_cons.txt as the attribute file, and CAIeq.txt as the co-variate network file. To begin, we must convert the data files into matrices, transform those matrices into networks, and attach the attribute file to our base network.

Understanding the Base Network

Let’s begin by looking at the summary of our base network.

Visualization

Before we conduct further analysis, let’s visualize our base network. Similar to our approach in Lab 2, we will begin by establishing set coordinates for our nodes in order to simplify visual comparisons.

Next we can visualize our base network.

Next, we will size the nodes by the their expertise value.

Let’s compare this visualization to sizing by in-degree centrality.

Understanding the Covariate Network

Let’s explore the summary statistics of our co-variate network.

## Network attributes:
##   vertices = 17
##   directed = TRUE
##   hyper = FALSE
##   loops = FALSE
##   multiple = FALSE
##   bipartite = FALSE
##  total edges = 23 
##    missing edges = 0 
##    non-missing edges = 23 
##  density = 0.08455882 
## 
## Vertex attributes:
##   vertex.names:
##    character valued attribute
##    17 valid vertex names
## 
## No edge attributes
## 
## Network edgelist matrix:
##       [,1] [,2]
##  [1,]    7    2
##  [2,]    6    4
##  [3,]    7    4
##  [4,]    9    4
##  [5,]   14    4
##  [6,]    6    5
##  [7,]   17    5
##  [8,]    2    6
##  [9,]   14    6
## [10,]    7    8
## [11,]    2    9
## [12,]    4    9
## [13,]    5    9
## [14,]    6    9
## [15,]   14    9
## [16,]   17    9
## [17,]    4   10
## [18,]    7   11
## [19,]    4   13
## [20,]    5   16
## [21,]    6   16
## [22,]    4   17
## [23,]    5   17

Visualization

We will repeat the visualization process for our co-variate network. Observe the location and distribution of edges in the following visualization.

Next, we will size the nodes by their expertise scores.

Part II: Constructing & Analyzing the ERGM Model (70 pts)

Next, we’re going to construct an Exponential Random Graph Model. Note that model construction is integral to this process; ERGM is not a single method of analysis, but a type of modelling that requires a theoretical grounding specific to the network and the hypotheses posed by the researcher. While the base ERGM simulation is set to take up to twenty iterations of simulations to fit the model by estimating parameters, it will stop at fewer if the generated networks converge on estimates of our coefficient values or paramters. ERGMs are used to predict ties as a function of individual covariates (i.e. attribute data, like “EX” in our example) or network structure.

We will start with a very basic model, looking only at the probability of edge formation, otherwise known as density, to demonstrate how co-efficients can be translated into odds ratios.44 The term for tie density (edges) is often used similarly to an intercept term in a linear regression or other linear model such as R’s glm. Keep in mind that the base network is about information retrieval. Our model will primarily allow us to ask what the probability is that an information retrieval relationship will form between two nodes.

model0 <- ergm(CRIeq ~ edges)

## Evaluating log-likelihood at the estimate.

summary(model0)

## 
## ==========================
## Summary of model fit
## ==========================
## 
## Formula:   CRIeq ~ edges
## 
## Iterations:  5 out of 20 
## 
## Monte Carlo MLE Results:
##       Estimate Std. Error MCMC % p-value    
## edges  -1.7288     0.1695      0  <1e-04 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##      Null Deviance: 377.1  on 272  degrees of freedom
##  Residual Deviance: 230.6  on 271  degrees of freedom
##  
## AIC: 232.6    BIC: 236.3    (Smaller is better.)

What do coefficients mean? Coefficients are the change in the (log-odds) likelihood of a tie for a unit change in a predictor. In our basic model above, our only predictor is information on the number of edges in the network. We can see that the coefficient estimate is negative, suggesting that a tie is more likely not to form than form (i.e. density is less than .5.) To get a better sense of how less likely it is for a tie to form, we can translate our log-odds into a probability.

Here coefficients mean that for a unit increase in the density of network, the log-odds of a tie in network are 1.72 times lower than the log-odds of not having a tie in the network. To translate to percentage, it’s P = odds/(1+odds) = 0.151, the probability of having a tie in network is 15.1%

To translate our estimated coefficient for the edges parameter into a probability, we can take the inverse-logit. We are thereby finding the probability that a tie will form in the network, looking only at the number of ties in the base network to build our model. We can see below that this probability is equal to our network density.

plogis(coef(model0)[[1]])

## [1] 0.1507353

network.density(CRIeq)

## [1] 0.1507353

Conceptually, this should be fairly easy to follow: if all we know about a network is the number of ties we have and we’re attempting to predict the probability that an edge will exist solely based on that information, then the probability of an edge existing at any point in the network equals the density of the network as a whole.55 As we move on from this toy model, keep in mind that the estimate for our edges parameter will change as we add additional predictors or network statistics as additional terms will partially explain tie formation.

Base Network ERGM

As we build a network, we can evaluate whether individual network statistics or node attributes prove our hypotheses and whether they do so in a way that is significantly different from random chance. Later, we will evaluate whether the model does a good job of explaining our observed network.

model1 <- ergm(CRIeq ~ edges      # Set the base term based on density/edge formation. 
               + mutual           # H1
               + transitive       # H2a: Transitive triads ( type 120D, 030T, 120U, or 300) # What if the tie is part of a transitive triad?
               + nodeicov("EX")   # H3a
               + nodeocov("EX")   # H3b
               + edgecov(CAIeq)   # H4
) 

## Starting maximum likelihood estimation via MCMLE:

## Iteration 1 of at most 20:

## Optimizing with step length 1.

## The log-likelihood improved by 0.1837.

## Step length converged once. Increasing MCMC sample size.

## Iteration 2 of at most 20:

## Optimizing with step length 1.

## The log-likelihood improved by 0.004038.

## Step length converged twice. Stopping.

## Evaluating log-likelihood at the estimate. Using 20 bridges: 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 .
## This model was fit using MCMC.  To examine model diagnostics and check for degeneracy, use the mcmc.diagnostics() function.

summary(model1)

## 
## ==========================
## Summary of model fit
## ==========================
## 
## Formula:   CRIeq ~ edges + mutual + transitive + nodeicov("EX") + nodeocov("EX") + 
##     edgecov(CAIeq)
## 
## Iterations:  2 out of 20 
## 
## Monte Carlo MLE Results:
##               Estimate Std. Error MCMC % p-value    
## edges          -6.9987     1.1704      0 < 1e-04 ***
## mutual         -1.4200     1.0017      0 0.15746    
## transitive      0.2795     0.1224      0 0.02323 *  
## nodeicov.EX     9.6094     2.0908      0 < 1e-04 ***
## nodeocov.EX     1.2373     1.6358      0 0.45008    
## edgecov.CAIeq   2.1239     0.6769      0 0.00189 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##      Null Deviance: 377.1  on 272  degrees of freedom
##  Residual Deviance: 123.8  on 266  degrees of freedom
##  
## AIC: 135.8    BIC: 157.4    (Smaller is better.)

kable(plogis(coef(model1)))

edges	0.0009122
mutual	0.1946590
transitive	0.5694190
nodeicov.EX	0.9999329
nodeocov.EX	0.7751006
edgecov.CAIeq	0.8932078

ERGM Model 2

Now we will change our model slightly to avoid convergence problems that lead to degeneracy. The terms transitive and dgwesp both rely on triangle formations so including both of them in the model leads to a situation similar to colinearity in a generalized linear model. While they measure slightly different things (take a look at the ergm-terms documentation to understand more about what’s happening under the hood), they’re typically used interchangeably. As a result, we need to use two separate models to test the relevance of these parameters.

model2 <- ergm(CRIeq ~ edges 
               + mutual                           # H1
               + dgwesp(0.5, fixed=T, type="OTP") # H2b: OTP "transitive shared partner" ordered pair (i,j) iff i->k->j.
               + nodeicov("EX")                   # H3a
               + nodeocov("EX")                   # H3b
               + edgecov(CAIeq)                   # H4
) 

## Starting maximum likelihood estimation via MCMLE:

## Iteration 1 of at most 20:

## Optimizing with step length 1.

## The log-likelihood improved by 0.1527.

## Step length converged once. Increasing MCMC sample size.

## Iteration 2 of at most 20:

## Optimizing with step length 1.

## The log-likelihood improved by 0.007917.

## Step length converged twice. Stopping.

## Evaluating log-likelihood at the estimate. Using 20 bridges: 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 .
## This model was fit using MCMC.  To examine model diagnostics and check for degeneracy, use the mcmc.diagnostics() function.

summary(model2) 

## 
## ==========================
## Summary of model fit
## ==========================
## 
## Formula:   CRIeq ~ edges + mutual + dgwesp(0.5, fixed = T, type = "OTP") + 
##     nodeicov("EX") + nodeocov("EX") + edgecov(CAIeq)
## 
## Iterations:  2 out of 20 
## 
## Monte Carlo MLE Results:
##                     Estimate Std. Error MCMC %  p-value    
## edges                -6.9226     1.2205      0  < 1e-04 ***
## mutual               -1.5643     1.0477      0 0.136584    
## gwesp.OTP.fixed.0.5   0.6182     0.3600      0 0.087128 .  
## nodeicov.EX           9.0243     2.3208      0 0.000127 ***
## nodeocov.EX           1.3914     1.7717      0 0.432963    
## edgecov.CAIeq         2.1138     0.6780      0 0.002022 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##      Null Deviance: 377.1  on 272  degrees of freedom
##  Residual Deviance: 124.3  on 266  degrees of freedom
##  
## AIC: 136.3    BIC: 157.9    (Smaller is better.)

kable(plogis(coef(model2)))

edges	0.0009843
mutual	0.1730248
gwesp.OTP.fixed.0.5	0.6498149
nodeicov.EX	0.9998796
nodeocov.EX	0.8008127
edgecov.CAIeq	0.8922413

Model Diagnostics

Next, judge convergence of the MCMC processes of the first model, using the mcmc.diagnostics() function. The function will plot the change of model statistics during the last iteration of the MCMC estimation procedure.77 Note that although the edge graphs appear to be periodic, the dips between whole numbers are due to the fact that edges are always whole numbers. For each model statistic, the left hand side plot gives the change of the statistic with iterations, and the right hand side plot is a histogram of the statistic values. Both are normalized, so the observed data locate at 0.

mcmc.diagnostics(model1)      # Performs the markov chain monte carlo diagnostics

## Sample statistics summary:
## 
## Iterations = 16384:4209664
## Thinning interval = 1024 
## Number of chains = 1 
## Sample size per chain = 4096 
## 
## 1. Empirical mean and standard deviation for each variable,
##    plus standard error of the mean:
## 
##                   Mean     SD Naive SE Time-series SE
## edges          0.05029  7.075  0.11054        0.15078
## mutual        -0.02881  2.373  0.03707        0.04868
## transitive     0.36719 29.434  0.45991        0.66936
## nodeicov.EX    0.03732  3.796  0.05931        0.08193
## nodeocov.EX   -0.03853  2.375  0.03712        0.05091
## edgecov.CAIeq  0.02905  1.797  0.02808        0.03401
## 
## 2. Quantiles for each variable:
## 
##                  2.5%     25%     50%    75%  97.5%
## edges         -12.000  -5.000  0.0000  5.000 15.000
## mutual         -4.000  -2.000  0.0000  1.000  6.000
## transitive    -40.000 -21.000 -5.0000 15.000 76.625
## nodeicov.EX    -6.588  -2.588 -0.2353  2.471  8.176
## nodeocov.EX    -4.353  -1.706 -0.1765  1.412  5.118
## edgecov.CAIeq  -4.000  -1.000  0.0000  1.000  3.000
## 
## 
## Sample statistics cross-correlations:
##                   edges    mutual transitive nodeicov.EX nodeocov.EX
## edges         1.0000000 0.7291288  0.9108447   0.9905801   0.9172898
## mutual        0.7291288 1.0000000  0.8584708   0.7329587   0.8646695
## transitive    0.9108447 0.8584708  1.0000000   0.9175585   0.9336619
## nodeicov.EX   0.9905801 0.7329587  0.9175585   1.0000000   0.9054735
## nodeocov.EX   0.9172898 0.8646695  0.9336619   0.9054735   1.0000000
## edgecov.CAIeq 0.5319630 0.4756379  0.5163785   0.5182302   0.5649074
##               edgecov.CAIeq
## edges             0.5319630
## mutual            0.4756379
## transitive        0.5163785
## nodeicov.EX       0.5182302
## nodeocov.EX       0.5649074
## edgecov.CAIeq     1.0000000
## 
## Sample statistics auto-correlation:
## Chain 1 
##                 edges      mutual  transitive  nodeicov.EX nodeocov.EX
## Lag 0     1.000000000 1.000000000 1.000000000  1.000000000  1.00000000
## Lag 1024  0.247867613 0.223412510 0.330347377  0.264160565  0.26215333
## Lag 2048  0.114867170 0.092421595 0.137483230  0.118401772  0.11281823
## Lag 3072  0.023164791 0.045190866 0.044179047  0.025828238  0.02508956
## Lag 4096  0.017483954 0.015038755 0.031692227  0.015022567  0.01327444
## Lag 5120 -0.006790644 0.007322916 0.002891676 -0.002736252 -0.01474736
##          edgecov.CAIeq
## Lag 0      1.000000000
## Lag 1024   0.150202745
## Lag 2048   0.061661892
## Lag 3072   0.030062910
## Lag 4096  -0.009328426
## Lag 5120  -0.009108246
## 
## Sample statistics burn-in diagnostic (Geweke):
## Chain 1 
## 
## Fraction in 1st window = 0.1
## Fraction in 2nd window = 0.5 
## 
##         edges        mutual    transitive   nodeicov.EX   nodeocov.EX 
##       0.35313      -0.05042       0.05859       0.45225      -0.04138 
## edgecov.CAIeq 
##      -0.11956 
## 
## Individual P-values (lower = worse):
##         edges        mutual    transitive   nodeicov.EX   nodeocov.EX 
##     0.7239890     0.9597893     0.9532810     0.6510881     0.9669915 
## edgecov.CAIeq 
##     0.9048278 
## Joint P-value (lower = worse):  0.8083045 .

## Warning in formals(fun): argument is not a function

## 
## 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).

Repeat the process for the second model.

mcmc.diagnostics(model2)      # Performs the markov chain monte carlo diagnostics

## Sample statistics summary:
## 
## Iterations = 16384:4209664
## Thinning interval = 1024 
## Number of chains = 1 
## Sample size per chain = 4096 
## 
## 1. Empirical mean and standard deviation for each variable,
##    plus standard error of the mean:
## 
##                          Mean     SD Naive SE Time-series SE
## edges                0.095459  5.464  0.08537        0.09239
## mutual               0.003418  1.716  0.02681        0.02746
## gwesp.OTP.fixed.0.5  0.037285 10.293  0.16083        0.17605
## nodeicov.EX          0.063290  2.938  0.04590        0.04995
## nodeocov.EX         -0.056913  1.697  0.02652        0.02856
## edgecov.CAIeq       -0.025635  1.662  0.02597        0.02825
## 
## 2. Quantiles for each variable:
## 
##                        2.5%    25%        50%   75%  97.5%
## edges               -10.000 -4.000  0.000e+00 4.000 11.000
## mutual               -3.000 -1.000  0.000e+00 1.000  3.000
## gwesp.OTP.fixed.0.5 -19.453 -7.257 -3.800e-02 7.034 20.415
## nodeicov.EX          -5.647 -1.941  1.176e-01 2.059  5.824
## nodeocov.EX          -3.412 -1.235 -5.000e-09 1.118  3.235
## edgecov.CAIeq        -3.000 -1.000  0.000e+00 1.000  3.000
## 
## 
## Sample statistics cross-correlations:
##                         edges    mutual gwesp.OTP.fixed.0.5 nodeicov.EX
## edges               1.0000000 0.5349603           0.9041329   0.9833636
## mutual              0.5349603 1.0000000           0.7193093   0.5417068
## gwesp.OTP.fixed.0.5 0.9041329 0.7193093           1.0000000   0.9207986
## nodeicov.EX         0.9833636 0.5417068           0.9207986   1.0000000
## nodeocov.EX         0.8686584 0.7423456           0.8979661   0.8418218
## edgecov.CAIeq       0.4015638 0.3621196           0.4195749   0.3807982
##                     nodeocov.EX edgecov.CAIeq
## edges                 0.8686584     0.4015638
## mutual                0.7423456     0.3621196
## gwesp.OTP.fixed.0.5   0.8979661     0.4195749
## nodeicov.EX           0.8418218     0.3807982
## nodeocov.EX           1.0000000     0.4511698
## edgecov.CAIeq         0.4511698     1.0000000
## 
## Sample statistics auto-correlation:
## Chain 1 
##                 edges        mutual gwesp.OTP.fixed.0.5  nodeicov.EX
## Lag 0     1.000000000  1.0000000000        1.0000000000  1.000000000
## Lag 1024  0.078775604  0.0237120185        0.0900419710  0.084280262
## Lag 2048  0.005724991 -0.0058077377        0.0013584399  0.003751389
## Lag 3072  0.012390431  0.0008266370       -0.0004036169  0.014634936
## Lag 4096 -0.017039156  0.0329981011        0.0053092201 -0.025310686
## Lag 5120  0.007574065 -0.0008314628        0.0120795806  0.006554031
##           nodeocov.EX edgecov.CAIeq
## Lag 0    1.0000000000   1.000000000
## Lag 1024 0.0737313539   0.083802358
## Lag 2048 0.0058394419   0.017793639
## Lag 3072 0.0114812248   0.004184903
## Lag 4096 0.0212565017  -0.005882316
## Lag 5120 0.0008786673  -0.010749413
## 
## Sample statistics burn-in diagnostic (Geweke):
## Chain 1 
## 
## Fraction in 1st window = 0.1
## Fraction in 2nd window = 0.5 
## 
##               edges              mutual gwesp.OTP.fixed.0.5 
##            -0.44914            -0.19654            -0.01433 
##         nodeicov.EX         nodeocov.EX       edgecov.CAIeq 
##            -0.58727             0.30896             0.30764 
## 
## Individual P-values (lower = worse):
##               edges              mutual gwesp.OTP.fixed.0.5 
##           0.6533333           0.8441867           0.9885681 
##         nodeicov.EX         nodeocov.EX       edgecov.CAIeq 
##           0.5570236           0.7573544           0.7583558 
## Joint P-value (lower = worse):  0.3837084 .

## Warning in formals(fun): argument is not a function

## 
## 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).

Model Evaluation

To evaluate the goodness-of-fit for our model, we need to simulate many variations of the model.99 See ?? simulate for more information.

Let’s visually inspect two of our random networks based on our first model.

ggnet2(sim[[1]], mode = baseLayout, label = TRUE, arrow.size = 8, edge.color="blue", arrow.gap = 0.03) + theme(panel.background = element_rect(fill = "#fffff8")) + guides(color = FALSE, size = FALSE)

ggnet2(sim[[10]], mode = baseLayout, label = TRUE, arrow.size = 8, edge.color="orange", arrow.gap = 0.03) + theme(panel.background = element_rect(fill = "#fffff8")) + guides(color = FALSE, size = FALSE)

Next we’re going to extract the number of triangles from each of the 100 samples, create a histogram of that model, and place a red arrow at the value of the observed network.

model.tridist <- sapply(1:100, function(x) summary(sim[[x]] ~triangle)) # Extracts the tiangle data from the simulated networks
hist(model.tridist, xlim=c(0,140), breaks = 20)                                     # Plots that triangle distribution as a histogram
CRIeq.tri <- summary(CRIeq ~ triangle)                                  # Saves the CRIeq triangle data from the summary to the CRI.eq variable
arrows(CRIeq.tri,20, CRIeq.tri, 5, col="red", lwd=3)                    # Adds an arrow to the plotted histogram

Goodness of Fit

Next, we will calculate the Goodness of Fit for several of the parameters in our model. A p-value closer to one is better; this represents the difference between the observed networks and simulations.

gof

## 
## Goodness-of-fit for in-degree 
## 
##    obs min mean max MC p-value
## 0    9   7 9.16  11       1.00
## 1    3   0 1.78   4       0.44
## 2    0   0 0.65   3       0.94
## 3    0   0 0.49   3       1.00
## 4    0   0 0.48   3       1.00
## 5    2   0 0.62   3       0.24
## 6    0   0 0.65   3       0.96
## 7    0   0 0.66   3       0.98
## 8    1   0 0.99   4       1.00
## 9    1   0 0.69   2       1.00
## 10   0   0 0.46   3       1.00
## 11   1   0 0.23   2       0.44
## 12   0   0 0.11   1       1.00
## 13   0   0 0.03   1       1.00
## 
## Goodness-of-fit for out-degree 
## 
##   obs min mean max MC p-value
## 0   4   0 2.02   7       0.30
## 1   1   0 3.23   8       0.30
## 2   3   0 3.80   8       0.92
## 3   4   0 3.52   8       0.96
## 4   3   0 2.79   7       1.00
## 5   2   0 1.34   4       0.80
## 6   0   0 0.29   2       1.00
## 7   0   0 0.01   1       1.00
## 
## Goodness-of-fit for edgewise shared partner 
## 
##      obs min  mean max MC p-value
## esp0  10   1 10.87  20       0.98
## esp1  13   5 12.48  20       1.00
## esp2  11   2 10.20  22       0.92
## esp3   5   0  5.55  20       1.00
## esp4   2   0  1.72  11       0.82
## esp5   0   0  0.24   3       1.00
## 
## Goodness-of-fit for minimum geodesic distance 
## 
##     obs min   mean max MC p-value
## 1    41  22  41.06  54       1.00
## 2    19  10  35.80  60       0.08
## 3     4   0  12.18  35       0.22
## 4     0   0   2.91  23       0.60
## 5     0   0   0.48  10       1.00
## 6     0   0   0.04   3       1.00
## Inf 208 112 179.53 240       0.14

# -------------------------------------------------------------------------------------------------
# Test the goodness of fit of the model
# Compiles statistics for these simulations as well as the observed network, and calculates p-values 
# -------------------------------------------------------------------------------------------------

par(mfrow=c(2,2))   # Separate the plot window into a 2 by 2 orientation
plot(gof)           # Plot the goodness of fit

Submitting the Lab (5 pts)

After knitting your file to RPubs, copy the URL and paste it into the comment field of the Lab 2 Assignment on Canvas. Save this .Rmd file and submit it in the file portion of your Canvas assignment. Make sure to review your file and its formatting. Run spell check (built into RStudio) and proofread your answers before submitting. If you can’t publish to RPubs, save your HTML file as a PDF and submit that instead.1111 There are many different ways to do this with different browsers. Google it.