library('igraph')
library('GGally')
library('ggplot2')
library('tergm')
library('ergm')
library('network')
library('networkDynamic')
library('ergm.count')
library('sna')
library('statnet.common')
library('tsna')
#### Setup ####
# Read in the edge lists of the three films
edgL4 <- read.csv('SW4_links.csv')
edgL5 <- read.csv('SW5_links.csv')
edgL6 <- read.csv('SW6_links.csv')

# Turn them into networks
sw4Net <- as.network(edgL4, matrix.type = 'edgelist',  
                    directed = FALSE, vertex.names = NULL)
sw5Net <- as.network(edgL5, matrix.type = 'edgelist', 
                     directed = FALSE, vertex.names = NULL)
sw6Net <- as.network(edgL6, matrix.type = 'edgelist', 
                     directed = FALSE, vertex.names = NULL)

# Read in the character attributes
charAtr4 <- as.matrix(read.table('SW4.txt', header = TRUE))
charAtr5 <- as.matrix(read.table('SW5.txt', header = TRUE))
charAtr6 <- as.matrix(read.table('SW6.txt', header = TRUE))

# Attach the attributes to the networks
sw4Net %v% 'Name' <- charAtr4[,2]
sw5Net %v% 'Name' <- charAtr5[,2]
sw6Net %v% 'Name' <- charAtr6[,2]

sw4Net %v% 'Faction' <- charAtr4[,3]
sw5Net %v% 'Faction' <- charAtr5[,3]
sw6Net %v% 'Faction' <- charAtr6[,3]

sw4Net %v% 'FactionId' <- as.numeric(charAtr4[,4])
sw5Net %v% 'FactionId' <- as.numeric(charAtr5[,4])
sw6Net %v% 'FactionId' <- as.numeric(charAtr6[,4])

sw4Net %v% 'Force' <- charAtr4[,5]
sw5Net %v% 'Force' <- charAtr5[,5]
sw6Net %v% 'Force' <- charAtr6[,5]

sw4Net %v% 'ForceId' <- as.numeric(charAtr4[,6])
sw5Net %v% 'ForceId' <- as.numeric(charAtr5[,6])
sw6Net %v% 'ForceId' <- as.numeric(charAtr6[,6])

sw4Net %v% 'ForceCol' <- charAtr4[,7]
sw5Net %v% 'ForceCol' <- charAtr5[,7]
sw6Net %v% 'ForceCol' <- charAtr6[,7]

# Make an igraph and matrix of all this for potential later use
sw4Ig <- asIgraph(sw4Net)
sw5Ig <- asIgraph(sw5Net)
sw6Ig <- asIgraph(sw6Net)

sw4Mat <- as.matrix(sw4Net)
sw5Mat <- as.matrix(sw5Net)
sw6Mat <- as.matrix(sw6Net)
#### Plotting ####
# Summarise and plot the networks with force and faction colours
summary(sw4Net)
## Network attributes:
##   vertices = 20
##   directed = FALSE
##   hyper = FALSE
##   loops = FALSE
##   multiple = FALSE
##   bipartite = FALSE
##  total edges = 58 
##    missing edges = 0 
##    non-missing edges = 58 
##  density = 0.3052632 
## 
## Vertex attributes:
## 
##  Faction:
##    character valued attribute
##    attribute summary:
##   Bounty_Hunter Galactic_Empire            Jedi         Neutral  Rebel_Alliance 
##               4               2               2               2               9 
##            Sith 
##               1 
## 
##  FactionId:
##    numeric valued attribute
##    attribute summary:
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##     1.0     1.0     2.0     2.8     5.0     6.0 
## 
##  Force:
##    character valued attribute
##    attribute summary:
##      None Sensitive 
##        16         4 
## 
##  ForceCol:
##    character valued attribute
##    attribute summary:
##   black hotpink 
##      16       4 
## 
##  ForceId:
##    numeric valued attribute
##    attribute summary:
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##     1.0     1.0     1.0     1.2     1.0     2.0 
## 
##  Name:
##    character valued attribute
##    attribute summary:
##    the 10 most common values are:
##        BERU       BIGGS       C_3PO   CHEWBACCA DARTH_VADER     DODONNA 
##           1           1           1           1           1           1 
##   GOLD_FIVE GOLD_LEADER      GREEDO         HAN 
##           1           1           1           1 
##   vertex.names:
##    character valued attribute
##    20 valid vertex names
## 
## Edge attributes:
## 
##  weight:
##    integer valued attribute
##    58values
## 
## Network edgelist matrix:
##       [,1] [,2]
##  [1,]    1   16
##  [2,]    2   16
##  [3,]    3   16
##  [4,]    4   16
##  [5,]    5   16
##  [6,]    6   16
##  [7,]    7   16
##  [8,]    8   16
##  [9,]    9   16
## [10,]    1    6
## [11,]    2    1
## [12,]    1    4
## [13,]    1    9
## [14,]    1    7
## [15,]    8    4
## [16,]   10    7
## [17,]    3    4
## [18,]    3    5
## [19,]    3    2
## [20,]    4    5
## [21,]    2    4
## [22,]    2    5
## [23,]    2    7
## [24,]    7    4
## [25,]    3    7
## [26,]    4    6
## [27,]    2    6
## [28,]    7    6
## [29,]   11   17
## [30,]   10   11
## [31,]   10   17
## [32,]    9    6
## [33,]    9    4
## [34,]   12    9
## [35,]    9   18
## [36,]    2    9
## [37,]    7   11
## [38,]    7   17
## [39,]    9    7
## [40,]   10    6
## [41,]   13   14
## [42,]   13   19
## [43,]   13    4
## [44,]   14   19
## [45,]   14    4
## [46,]    4   19
## [47,]    8    7
## [48,]    7   15
## [49,]    4   15
## [50,]    8   15
## [51,]    8    2
## [52,]    2   15
## [53,]   15   19
## [54,]   14   15
## [55,]    8   19
## [56,]   15   20
## [57,]    8   14
## [58,]    4   20
ggnet2(sw4Net, mode = 'circle', color.palette = 'Pastel1', edge.color = 'grey',
       node.color = 'Faction', node.size = 20, legend.position = 'none',
       label = 'Name', label.color = 'ForceCol', label.size = 7)

summary(sw5Net)
## Network attributes:
##   vertices = 20
##   directed = FALSE
##   hyper = FALSE
##   loops = FALSE
##   multiple = FALSE
##   bipartite = FALSE
##  total edges = 54 
##    missing edges = 0 
##    non-missing edges = 54 
##  density = 0.2842105 
## 
## Vertex attributes:
## 
##  Faction:
##    character valued attribute
##    attribute summary:
##   Bounty_Hunter Galactic_Empire            Jedi  Rebel_Alliance            Sith 
##               2               3               3              10               2 
## 
##  FactionId:
##    numeric valued attribute
##    attribute summary:
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    1.00    1.00    1.50    2.20    3.25    5.00 
## 
##  Force:
##    character valued attribute
##    attribute summary:
##      None Sensitive 
##        14         6 
## 
##  ForceCol:
##    character valued attribute
##    attribute summary:
##   black hotpink 
##      14       6 
## 
##  ForceId:
##    numeric valued attribute
##    attribute summary:
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##     1.0     1.0     1.0     1.3     2.0     2.0 
## 
##  Name:
##    character valued attribute
##    attribute summary:
##    the 10 most common values are:
##   BOBA_FETT       C-3PO   CHEWBACCA        DACK DARTH_VADER      DERLIN 
##           1           1           1           1           1           1 
##     EMPEROR         HAN       LANDO        LEIA 
##           1           1           1           1 
##   vertex.names:
##    character valued attribute
##    20 valid vertex names
## 
## No edge attributes
## 
## Network edgelist matrix:
##       [,1] [,2]
##  [1,]    1   14
##  [2,]    2   14
##  [3,]    3   14
##  [4,]    4   14
##  [5,]    5   14
##  [6,]    6   14
##  [7,]    2    1
##  [8,]    7    1
##  [9,]    1    3
## [10,]    4    3
## [11,]    2    3
## [12,]    5    3
## [13,]    8    3
## [14,]    9    3
## [15,]    6    3
## [16,]    1    4
## [17,]    2    4
## [18,]    9    4
## [19,]    7    4
## [20,]    1   15
## [21,]    4   15
## [22,]    5   15
## [23,]    2   15
## [24,]    1    5
## [25,]    4    5
## [26,]    2    5
## [27,]   10    5
## [28,]    6    5
## [29,]    7    5
## [30,]    9    5
## [31,]    7    2
## [32,]    2   10
## [33,]    3   13
## [34,]    4   13
## [35,]    4   16
## [36,]    3   16
## [37,]   11   17
## [38,]    9   17
## [39,]   12   17
## [40,]    7   17
## [41,]    9   11
## [42,]    1    9
## [43,]    7    9
## [44,]    2    9
## [45,]    3   18
## [46,]    9   12
## [47,]    9   19
## [48,]    3   20
## [49,]   13   20
## [50,]    1    6
## [51,]    4    6
## [52,]    2    6
## [53,]    9    6
## [54,]    7    6
ggnet2(sw5Net, mode = 'circle', color.palette = 'Pastel1', edge.color = 'grey',
       node.color = 'Faction', node.size = 20, legend.position = 'none',
       label = 'Name', label.color = 'ForceCol', label.size = 7)

summary(sw6Net)
## Network attributes:
##   vertices = 19
##   directed = FALSE
##   hyper = FALSE
##   loops = FALSE
##   multiple = FALSE
##   bipartite = FALSE
##  total edges = 53 
##    missing edges = 0 
##    non-missing edges = 53 
##  density = 0.3099415 
## 
## Vertex attributes:
## 
##  Faction:
##    character valued attribute
##    attribute summary:
##   Bounty_Hunter Galactic_Empire            Jedi  Rebel_Alliance            Sith 
##               3               2               3               9               2 
## 
##  FactionId:
##    numeric valued attribute
##    attribute summary:
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   1.000   1.000   2.000   2.316   3.500   5.000 
## 
##  Force:
##    character valued attribute
##    attribute summary:
##      None Sensitive 
##        13         6 
## 
##  ForceCol:
##    character valued attribute
##    attribute summary:
##   black hotpink 
##      13       6 
## 
##  ForceId:
##    numeric valued attribute
##    attribute summary:
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   1.000   1.000   1.000   1.316   2.000   2.000 
## 
##  Name:
##    character valued attribute
##    attribute summary:
##    the 10 most common values are:
## ADMIRAL_ACKBAR    BIB_FORTUNA         BOUSHH          C-3PO      CHEWBACCA 
##              1              1              1              1              1 
##    DARTH_VADER        EMPEROR            HAN          JABBA      JERJERROD 
##              1              1              1              1              1 
##   vertex.names:
##    character valued attribute
##    19 valid vertex names
## 
## No edge attributes
## 
## Network edgelist matrix:
##       [,1] [,2]
##  [1,]    1   14
##  [2,]    2   14
##  [3,]    3   14
##  [4,]    4   14
##  [5,]    5   14
##  [6,]    6   14
##  [7,]    7   14
##  [8,]    8   14
##  [9,]    2    1
## [10,]    1    3
## [11,]    1    6
## [12,]    1    5
## [13,]    1    7
## [14,]    1    8
## [15,]    9   15
## [16,]   10    2
## [17,]    2    4
## [18,]    2    5
## [19,]    4    5
## [20,]   11    4
## [21,]   11    2
## [22,]   11    6
## [23,]   11    7
## [24,]    6    7
## [25,]    2    6
## [26,]    2    7
## [27,]   10    5
## [28,]   10    4
## [29,]    6    5
## [30,]    7    5
## [31,]    3    5
## [32,]    6    3
## [33,]    3    7
## [34,]    9   13
## [35,]    5   16
## [36,]    5   17
## [37,]    6    8
## [38,]   12    6
## [39,]    3    8
## [40,]   12    3
## [41,]    2    3
## [42,]   12    8
## [43,]    2    8
## [44,]    7    8
## [45,]    5    8
## [46,]   12    2
## [47,]   12    7
## [48,]   12    5
## [49,]    9   18
## [50,]    9    5
## [51,]   13    5
## [52,]    3   19
## [53,]   12   19
ggnet2(sw6Net, mode = 'circle', color.palette = 'Pastel1', edge.color = 'grey',
       node.color = 'Faction', node.size = 20, legend.position = 'none',
       label = 'Name', label.color = 'ForceCol', label.size = 7)

#### Hypothesis Test ####
# ERGMs Let's Go, add diff=TRUE to nodematch for differential homophily
sw4ERGM <- ergm(sw4Net ~ edges 
                + nodematch('Faction')
                + nodefactor('Force')
                + gwesp(decay = 0.5, fixed = TRUE))
##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.4429.
##The log-likelihood improved by 1.7119.
##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.0735.
##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.0529.
##Convergence test p-value: 0.5832. Not converged with 99% confidence; increasing sample size.
## Iteration 4 of at most 60:
## Optimizing with step length 1.0000.
## The log-likelihood improved by 0.0252.
## Convergence test p-value: 0.0168. 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.0657. 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.0135.
## Convergence test p-value: 0.2334. Not converged with 99% confidence; increasing sample size.
## Iteration 7 of at most 60:
## Optimizing with step length 1.0000.
## The log-likelihood improved by 0.0048.
## 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.
sw5ERGM <- ergm(sw5Net ~ edges 
                + nodematch('Faction')
                + nodefactor('Force')
                + gwesp(decay = 0.5, fixed = TRUE))
## 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.9045.
## The log-likelihood improved by 1.8434.
## Estimating equations are not within tolerance region.
## Iteration 2 of at most 60:
## Optimizing with step length 1.0000.
## The log-likelihood improved by 0.2599.
## 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.0223.
## Convergence test p-value: 0.0228. Not converged with 99% confidence; increasing sample size.
## Iteration 4 of at most 60:
## Optimizing with step length 1.0000.
## The log-likelihood improved by 0.0423.
## Convergence test p-value: 0.4334. 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.0134.
## Convergence test p-value: 0.0457. 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.0055.
## Convergence test p-value: 0.0062. 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.
sw6ERGM <- ergm(sw6Net ~ edges 
                + nodematch('FactionId')
                + nodefactor('Force')
                + gwesp(decay = 0.5, fixed = TRUE))
## 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.5407.
## The log-likelihood improved by 2.0323.
## Estimating equations are not within tolerance region.
## Iteration 2 of at most 60:
## Optimizing with step length 1.0000.
## The log-likelihood improved by 0.1183.
## 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.0972.
## Convergence test p-value: 0.5714. Not converged with 99% confidence; increasing sample size.
## Iteration 4 of at most 60:
## Optimizing with step length 1.0000.
## The log-likelihood improved by 0.2279.
## Estimating equations are not within tolerance region.
## Iteration 5 of at most 60:
## Optimizing with step length 1.0000.
## The log-likelihood improved by 0.2414.
## Estimating equations are not within tolerance region.
## Estimating equations did not move closer to tolerance region more than 1 time(s) in 4 steps; increasing sample size.
## Iteration 6 of at most 60:
## Optimizing with step length 1.0000.
## The log-likelihood improved by 0.0283.
## Convergence test p-value: 0.0192. Not converged with 99% confidence; increasing sample size.
## Iteration 7 of at most 60:
## Optimizing with step length 1.0000.
## The log-likelihood improved by 0.0064.
## Convergence test p-value: 0.0136. Not converged with 99% confidence; increasing sample size.
## Iteration 8 of at most 60:
## Optimizing with step length 1.0000.
## The log-likelihood improved by 0.0427.
## Convergence test p-value: 0.0116. Not converged with 99% confidence; increasing sample size.
## Iteration 9 of at most 60:
## Optimizing with step length 1.0000.
## The log-likelihood improved by 0.0308.
## Convergence test p-value: 0.0155. Not converged with 99% confidence; increasing sample size.
## Iteration 10 of at most 60:
## Optimizing with step length 1.0000.
## The log-likelihood improved by 0.0062.
## Convergence test p-value: 0.0006. 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.
#### Goodness of fit ####
# Do GoF for each and then plot them, simples
mcmc.diagnostics(sw4ERGM)
## Sample statistics summary:
## 
## Iterations = 200704:2123776
## Thinning interval = 2048 
## Number of chains = 1 
## Sample size per chain = 940 
## 
## 1. Empirical mean and standard deviation for each variable,
##    plus standard error of the mean:
## 
##                               Mean     SD Naive SE Time-series SE
## edges                      -1.1532 14.479   0.4722         1.0322
## nodematch.Faction          -0.3426  4.397   0.1434         0.2458
## nodefactor.Force.Sensitive -0.3000  7.640   0.2492         0.4214
## gwesp.fixed.0.5            -1.8003 23.767   0.7752         1.6280
## 
## 2. Quantiles for each variable:
## 
##                              2.5%    25%   50%   75% 97.5%
## edges                      -35.57  -9.25 1.000  9.00 21.00
## nodematch.Faction           -9.00  -3.00 0.000  3.00  8.00
## nodefactor.Force.Sensitive -18.00  -4.00 0.000  5.00 12.00
## gwesp.fixed.0.5            -55.61 -16.19 1.856 14.76 35.37
## 
## 
## Are sample statistics significantly different from observed?
##                 edges nodematch.Faction nodefactor.Force.Sensitive
## diff.      -1.1531915        -0.3425532                 -0.3000000
## test stat. -1.1172704        -1.3936010                 -0.7119719
## P-val.      0.2638787         0.1634381                  0.4764822
##            gwesp.fixed.0.5 Overall (Chi^2)
## diff.           -1.8003322              NA
## test stat.      -1.1058820       4.1746335
## P-val.           0.2687776       0.3862184
## 
## Sample statistics cross-correlations:
##                                edges nodematch.Faction
## edges                      1.0000000         0.7323134
## nodematch.Faction          0.7323134         1.0000000
## nodefactor.Force.Sensitive 0.7913563         0.4930686
## gwesp.fixed.0.5            0.9959739         0.7281236
##                            nodefactor.Force.Sensitive gwesp.fixed.0.5
## edges                                       0.7913563       0.9959739
## nodematch.Faction                           0.4930686       0.7281236
## nodefactor.Force.Sensitive                  1.0000000       0.7946953
## gwesp.fixed.0.5                             0.7946953       1.0000000
## 
## Sample statistics auto-correlation:
## Chain 1 
##               edges nodematch.Faction nodefactor.Force.Sensitive
## Lag 0     1.0000000        1.00000000                 1.00000000
## Lag 2048  0.5007570        0.29525394                 0.39153187
## Lag 4096  0.3250061        0.17895958                 0.24655491
## Lag 6144  0.2384037        0.13262799                 0.14850980
## Lag 8192  0.2125321        0.13056602                 0.10016967
## Lag 10240 0.1412955        0.08256301                 0.09667606
##           gwesp.fixed.0.5
## Lag 0           1.0000000
## Lag 2048        0.4707481
## Lag 4096        0.3045870
## Lag 6144        0.2206387
## Lag 8192        0.1939215
## Lag 10240       0.1290108
## 
## Sample statistics burn-in diagnostic (Geweke):
## Chain 1 
## 
## Fraction in 1st window = 0.1
## Fraction in 2nd window = 0.5 
## 
##                      edges          nodematch.Faction 
##                    -1.0700                    -1.3366 
## nodefactor.Force.Sensitive            gwesp.fixed.0.5 
##                    -0.8296                    -1.1125 
## 
## Individual P-values (lower = worse):
##                      edges          nodematch.Faction 
##                  0.2846340                  0.1813583 
## nodefactor.Force.Sensitive            gwesp.fixed.0.5 
##                  0.4067641                  0.2659398 
## Joint P-value (lower = worse):  0.2940664 .
## 
## 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).

gof4 <- gof(sw4ERGM ~ degree + distance)
## Warning in gof.formula(object = object$formula, coef = coef, GOF = GOF, : No
## parameter values given, using 0.
par(mfrow = c(2, 2))
plot(gof4, main = 'Episode 4 ERGM Goodness of Fit')

mcmc.diagnostics(sw5ERGM)

## Sample statistics summary:
## 
## Iterations = 9984:132096
## Thinning interval = 256 
## Number of chains = 1 
## Sample size per chain = 478 
## 
## 1. Empirical mean and standard deviation for each variable,
##    plus standard error of the mean:
## 
##                                Mean     SD Naive SE Time-series SE
## edges                      -0.44351  8.290   0.3792         0.6608
## nodematch.Faction          -0.06485  4.211   0.1926         0.3351
## nodefactor.Force.Sensitive -0.72176  7.497   0.3429         0.5723
## gwesp.fixed.0.5            -0.76745 14.622   0.6688         1.0933
## 
## 2. Quantiles for each variable:
## 
##                              2.5%    25%     50%   75% 97.5%
## edges                      -17.00  -6.00  0.0000 5.000 14.07
## nodematch.Faction           -9.00  -3.00  0.0000 3.000  7.00
## nodefactor.Force.Sensitive -16.00  -6.00  0.0000 5.000 13.00
## gwesp.fixed.0.5            -29.91 -10.57 -0.4927 9.966 25.75
## 
## 
## Are sample statistics significantly different from observed?
##                 edges nodematch.Faction nodefactor.Force.Sensitive
## diff.      -0.4435146       -0.06485356                 -0.7217573
## test stat. -0.6712055       -0.19354611                 -1.2612474
## P-val.      0.5020896        0.84653130                  0.2072197
##            gwesp.fixed.0.5 Overall (Chi^2)
## diff.           -0.7674467              NA
## test stat.      -0.7019591       1.9478353
## P-val.           0.4827047       0.7500375
## 
## Sample statistics cross-correlations:
##                                edges nodematch.Faction
## edges                      1.0000000         0.6408628
## nodematch.Faction          0.6408628         1.0000000
## nodefactor.Force.Sensitive 0.7649779         0.2274658
## gwesp.fixed.0.5            0.9802194         0.6529219
##                            nodefactor.Force.Sensitive gwesp.fixed.0.5
## edges                                       0.7649779       0.9802194
## nodematch.Faction                           0.2274658       0.6529219
## nodefactor.Force.Sensitive                  1.0000000       0.7616518
## gwesp.fixed.0.5                             0.7616518       1.0000000
## 
## Sample statistics auto-correlation:
## Chain 1 
##                edges nodematch.Faction nodefactor.Force.Sensitive
## Lag 0     1.00000000        1.00000000                 1.00000000
## Lag 256   0.50382094        0.50249456                 0.47075861
## Lag 512   0.22708482        0.25667402                 0.25048343
## Lag 768   0.11644912        0.12208842                 0.11227958
## Lag 1024  0.01682002        0.05019640                 0.03008083
## Lag 1280 -0.04838246        0.05230081                -0.06401783
##          gwesp.fixed.0.5
## Lag 0         1.00000000
## Lag 256       0.45456400
## Lag 512       0.20240182
## Lag 768       0.09590914
## Lag 1024      0.03178156
## Lag 1280     -0.02034737
## 
## Sample statistics burn-in diagnostic (Geweke):
## Chain 1 
## 
## Fraction in 1st window = 0.1
## Fraction in 2nd window = 0.5 
## 
##                      edges          nodematch.Faction 
##                      2.545                      1.517 
## nodefactor.Force.Sensitive            gwesp.fixed.0.5 
##                      1.924                      2.734 
## 
## Individual P-values (lower = worse):
##                      edges          nodematch.Faction 
##                0.010922353                0.129282433 
## nodefactor.Force.Sensitive            gwesp.fixed.0.5 
##                0.054297684                0.006264907 
## Joint P-value (lower = worse):  0.479471 .

## 
## 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).
gof5 <- gof(sw5ERGM ~ degree + distance)
## Warning in gof.formula(object = object$formula, coef = coef, GOF = GOF, : No
## parameter values given, using 0.
par(mfrow = c(2,2))
plot(gof5, main = 'Episode 5 ERGM Goodness of Fit')

mcmc.diagnostics(sw6ERGM)

## Sample statistics summary:
## 
## Iterations = 49408:258048
## Thinning interval = 256 
## Number of chains = 1 
## Sample size per chain = 816 
## 
## 1. Empirical mean and standard deviation for each variable,
##    plus standard error of the mean:
## 
##                                 Mean     SD Naive SE Time-series SE
## edges                       0.007353  9.326   0.3265         0.6291
## nodematch.FactionId        -0.290441  3.648   0.1277         0.2235
## nodefactor.Force.Sensitive  0.183824  7.914   0.2770         0.5062
## gwesp.fixed.0.5             0.143392 16.411   0.5745         1.1543
## 
## 2. Quantiles for each variable:
## 
##                              2.5%    25%    50%   75% 97.5%
## edges                      -19.00  -6.00 0.0000  6.00 17.00
## nodematch.FactionId         -7.00  -3.00 0.0000  2.00  7.00
## nodefactor.Force.Sensitive -15.00  -5.00 0.0000  6.00 16.00
## gwesp.fixed.0.5            -32.36 -11.18 0.4667 11.36 30.22
## 
## 
## Are sample statistics significantly different from observed?
##                  edges nodematch.FactionId nodefactor.Force.Sensitive
## diff.      0.007352941          -0.2904412                  0.1838235
## test stat. 0.011687169          -1.2997479                  0.3631131
## P-val.     0.990675201           0.1936874                  0.7165204
##            gwesp.fixed.0.5 Overall (Chi^2)
## diff.            0.1433924              NA
## test stat.       0.1242213        5.386586
## P-val.           0.9011401        0.256506
## 
## Sample statistics cross-correlations:
##                                edges nodematch.FactionId
## edges                      1.0000000           0.5915482
## nodematch.FactionId        0.5915482           1.0000000
## nodefactor.Force.Sensitive 0.7321176           0.2676883
## gwesp.fixed.0.5            0.9877398           0.5837440
##                            nodefactor.Force.Sensitive gwesp.fixed.0.5
## edges                                       0.7321176       0.9877398
## nodematch.FactionId                         0.2676883       0.5837440
## nodefactor.Force.Sensitive                  1.0000000       0.7168744
## gwesp.fixed.0.5                             0.7168744       1.0000000
## 
## Sample statistics auto-correlation:
## Chain 1 
##              edges nodematch.FactionId nodefactor.Force.Sensitive
## Lag 0    1.0000000          1.00000000                  1.0000000
## Lag 256  0.5753105          0.44286357                  0.5386298
## Lag 512  0.3503821          0.20657481                  0.3020993
## Lag 768  0.2439780          0.15131336                  0.2137235
## Lag 1024 0.1807586          0.11628507                  0.1627960
## Lag 1280 0.1106107          0.08958448                  0.1482910
##          gwesp.fixed.0.5
## Lag 0         1.00000000
## Lag 256       0.53079601
## Lag 512       0.31595205
## Lag 768       0.22581277
## Lag 1024      0.16016001
## Lag 1280      0.09303864
## 
## Sample statistics burn-in diagnostic (Geweke):
## Chain 1 
## 
## Fraction in 1st window = 0.1
## Fraction in 2nd window = 0.5 
## 
##                      edges        nodematch.FactionId 
##                      1.349                      1.075 
## nodefactor.Force.Sensitive            gwesp.fixed.0.5 
##                      1.274                      1.344 
## 
## Individual P-values (lower = worse):
##                      edges        nodematch.FactionId 
##                  0.1772140                  0.2824550 
## nodefactor.Force.Sensitive            gwesp.fixed.0.5 
##                  0.2025443                  0.1790002 
## Joint P-value (lower = worse):  0.5265426 .

## 
## 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).
gof6 <- gof(sw6ERGM ~ degree + distance)
## Warning in gof.formula(object = object$formula, coef = coef, GOF = GOF, : No
## parameter values given, using 0.
par(mfrow = c(2,2))
plot(gof6, main = 'Episode 6 ERGM Goodness of Fit')

#### Get the results, we done ####
summary(sw4ERGM)
## Call:
## ergm(formula = sw4Net ~ edges + nodematch("Faction") + nodefactor("Force") + 
##     gwesp(decay = 0.5, fixed = TRUE))
## 
## Monte Carlo Maximum Likelihood Results:
## 
##                            Estimate Std. Error MCMC % z value Pr(>|z|)    
## edges                      -4.89999    0.77679      0  -6.308   <1e-04 ***
## nodematch.Faction          -0.00506    0.34180      0  -0.015   0.9882    
## nodefactor.Force.Sensitive  0.40367    0.22065      0   1.829   0.0673 .  
## gwesp.fixed.0.5             2.10081    0.47326      0   4.439   <1e-04 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##      Null Deviance: 263.4  on 190  degrees of freedom
##  Residual Deviance: 193.1  on 186  degrees of freedom
##  
## AIC: 201.1  BIC: 214.1  (Smaller is better. MC Std. Err. = 0.2443)
summary(sw5ERGM)
## Call:
## ergm(formula = sw5Net ~ edges + nodematch("Faction") + nodefactor("Force") + 
##     gwesp(decay = 0.5, fixed = TRUE))
## 
## Monte Carlo Maximum Likelihood Results:
## 
##                            Estimate Std. Error MCMC % z value Pr(>|z|)    
## edges                       -3.9554     0.6209      0  -6.370  < 1e-04 ***
## nodematch.Faction            1.5203     0.3786      0   4.015  < 1e-04 ***
## nodefactor.Force.Sensitive   0.6039     0.2512      0   2.404  0.01623 *  
## gwesp.fixed.0.5              1.1348     0.3593      0   3.158  0.00159 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##      Null Deviance: 263.4  on 190  degrees of freedom
##  Residual Deviance: 180.4  on 186  degrees of freedom
##  
## AIC: 188.4  BIC: 201.4  (Smaller is better. MC Std. Err. = 0.3271)
summary(sw6ERGM)
## Call:
## ergm(formula = sw6Net ~ edges + nodematch("FactionId") + nodefactor("Force") + 
##     gwesp(decay = 0.5, fixed = TRUE))
## 
## Monte Carlo Maximum Likelihood Results:
## 
##                            Estimate Std. Error MCMC % z value Pr(>|z|)    
## edges                       -3.3548     0.7182      0  -4.671  < 1e-04 ***
## nodematch.FactionId         -0.1960     0.3573      0  -0.549 0.583267    
## nodefactor.Force.Sensitive  -0.2175     0.1954      0  -1.113 0.265642    
## gwesp.fixed.0.5              1.4047     0.3917      0   3.587 0.000335 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##      Null Deviance: 237.1  on 171  degrees of freedom
##  Residual Deviance: 189.3  on 167  degrees of freedom
##  
## AIC: 197.3  BIC: 209.9  (Smaller is better. MC Std. Err. = 0.3083)