Intro

We work with 2 panels of the Sampson data to:

Verify the values of the nw.stats for all of the operators
Fit and simulate from the separable and joint models, showing that both produce an equilibrium target edgecount equal to summary(nw.series ~ edges)

Setup & data

# loadedNamespaces()
# unloadNamespace("ergm")
# unloadNamespace("statnet.common")

master_lib <- "~/R/GHmaster-library"

# #Make lib 1st position in libPath
defaultPaths <- .libPaths()
.libPaths(c(master_lib, defaultPaths))
.libPaths()

[1] "C:/Users/Martina Morris/Documents/R/GHmaster-library"
[2] "C:/Users/Martina Morris/Documents/R/win-library/4.0" 
[3] "C:/Program Files/R/R-4.0.5/library"

#devtools::dev_mode(on=T, path = .libPaths()[1])
library(network)
library(ergm)
library(tergm)
library(tsna)

sessioninfo::package_info(pkgs = c("statnet.common", "network", "rle",
                                   "ergm", "tergm"),
                          dependencies = FALSE)

 package        * version    date       lib
 ergm           * 4.0-6512   2021-06-22 [1]
 network        * 1.17.1-685 2021-06-16 [1]
 rle              0.9.2-234  2021-06-16 [1]
 statnet.common   4.5.0-362  2021-06-16 [1]
 tergm          * 4.0-2301   2021-06-22 [1]
 source                                 
 Github (statnet/ergm@763ddda)          
 Github (statnet/network@253413a)       
 Github (statnet/rle@d08b185)           
 Github (statnet/statnet.common@e3c23c1)
 Github (statnet/tergm@eac9f3b)         

[1] C:/Users/Martina Morris/Documents/R/GHmaster-library
[2] C:/Users/Martina Morris/Documents/R/win-library/4.0
[3] C:/Program Files/R/R-4.0.5/library

Types of networks

first network

data(samplk)
samplk1

 Network attributes:
  vertices = 18 
  directed = TRUE 
  hyper = FALSE 
  loops = FALSE 
  multiple = FALSE 
  bipartite = FALSE 
  total edges= 55 
    missing edges= 0 
    non-missing edges= 55 

 Vertex attribute names: 
    cloisterville group vertex.names 

No edge attributes

list

# summary shows the nodes & edges at each slice
samp.list <- list(samplk1,samplk2)
samp.list

[[1]]
 Network attributes:
  vertices = 18 
  directed = TRUE 
  hyper = FALSE 
  loops = FALSE 
  multiple = FALSE 
  bipartite = FALSE 
  total edges= 55 
    missing edges= 0 
    non-missing edges= 55 

 Vertex attribute names: 
    cloisterville group vertex.names 

No edge attributes

[[2]]
 Network attributes:
  vertices = 18 
  directed = TRUE 
  hyper = FALSE 
  loops = FALSE 
  multiple = FALSE 
  bipartite = FALSE 
  total edges= 57 
    missing edges= 0 
    non-missing edges= 57 

 Vertex attribute names: 
    cloisterville group vertex.names 

No edge attributes

NetSeries

# summary shows sum of cross-sectional node/edge counts after net1
samp.series <- NetSeries(samp.list)
samp.series

 Combined 1 networks on '.NetworkID':
  1: n = 18, directed = TRUE, bipartite = FALSE, loops = FALSE

 Network attributes:
  vertices = 18 
  directed = TRUE 
  hyper = FALSE 
  loops = FALSE 
  multiple = FALSE 
  bipartite = FALSE 
  ergm:
    constraints: ~blockdiag(".NetworkID") + discord(.PrevNet)
  total edges= 57 
    missing edges= 0 
    non-missing edges= 57 

 Vertex attribute names: 
    .NetworkID .NetworkName cloisterville group vertex.names 

No edge attributes

samp.series  %ergmlhs% "constraints"

~blockdiag(".NetworkID") + discord(.PrevNet)
<environment: 0x000000001ca18c08>

networkDynamic

Collapsed version shows the unique total node and edgecount.

# summary shows the unique # nodes and edges in both cases
# we start at 1 to make the series index consistent with the other
# network storage mechanisms

samp.dyn <- networkDynamic(network.list = samp.list, start=1)

Onsets and termini not specified, assuming each network in network.list should have a discrete spell of length 1
Argument base.net not specified, using first element of network.list instead
Created net.obs.period to describe network
 Network observation period info:
  Number of observation spells: 1 
  Maximal time range observed: 1 until 3 
  Temporal mode: discrete 
  Time unit: step 
  Suggested time increment: 1

samp.dyn

NetworkDynamic properties:
  distinct change times: 3 
  maximal time range: 1 until  3 

Includes optional net.obs.period attribute:
 Network observation period info:
  Number of observation spells: 1 
  Maximal time range observed: 1 until 3 
  Temporal mode: discrete 
  Time unit: step 
  Suggested time increment: 1 

 Network attributes:
  vertices = 18 
  directed = TRUE 
  hyper = FALSE 
  loops = FALSE 
  multiple = FALSE 
  bipartite = FALSE 
  net.obs.period: (not shown)
  total edges= 77 
    missing edges= 0 
    non-missing edges= 77 

 Vertex attribute names: 
    active cloisterville group vertex.names 

 Edge attribute names: 
    active

samp.collapse <- network.collapse(samp.dyn)
samp.collapse

 Network attributes:
  vertices = 18 
  directed = TRUE 
  hyper = FALSE 
  loops = FALSE 
  multiple = FALSE 
  bipartite = FALSE 
  total edges= 77 
    missing edges= 0 
    non-missing edges= 77 

 Vertex attribute names: 
    cloisterville group vertex.names 

No edge attributes

Descriptives

Traditional

Cross sectional stats and form/diss changes from tsna

samp.stats <- cbind(edges = tErgmStats(samp.dyn, "~edges"),
                    form = tEdgeFormation(samp.dyn),
                    diss = tEdgeDissolution(samp.dyn))

samp.stats

Time Series:
Start = 1 
End = 3 
Frequency = 1 
  edges form diss
1    55   55    0
2    57   22   20
3     0    0   57

Term counts

The ones with operator names replicate the edge counts performed by the operators.

start = network.edgecount(samplk1)
cross = network.edgecount(samplk2)
change = network.edgecount(xor(samplk1, samplk2))
form = network.edgecount(samplk1 | samplk2)
persist = network.edgecount(samplk1 & samplk2)
diss = -persist
dissolved = network.edgecount(samplk1) - persist

cbind(start=start, cross=cross, change=change, 
      form=form, persist=persist, diss=diss, dissolved=dissolved
      )

     start cross change form persist diss dissolved
[1,]    55    57     42   77      35  -35        20

Edges only models

Will work with the NetSeries object here.

Summary

summary(samp.series ~
          Cross(~edges) +
          Change(~edges) +
          Form(~edges) +
          Diss(~edges))

       edges Change~edges   Form~edges   Diss~edges 
          57           42           77          -35

Fit

Separable

sep.fit <- tergm(
  samp.series ~
    Form(~edges) +
    Diss(~edges),
  estimate = "CMLE"
)

summary(sep.fit)

Call:
tergm(formula = samp.series ~ Form(~edges) + Diss(~edges), estimate = "CMLE")

Conditional Maximum Likelihood Results:

           Estimate Std. Error MCMC % z value Pr(>|z|)    
Form~edges  -2.3427     0.2232      0 -10.497   <1e-04 ***
Diss~edges  -0.5596     0.2803      0  -1.996   0.0459 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

     Null Deviance: 424.2  on 306  degrees of freedom
 Residual Deviance: 221.2  on 304  degrees of freedom
 
AIC: 225.2  BIC: 232.7  (Smaller is better. MC Std. Err. = 0)

Joint

jnt.fit <- tergm(
  samp.series ~
    Cross(~edges) +
    Change(~edges),
  estimate = "CMLE"
)

summary(jnt.fit)

Call:
tergm(formula = samp.series ~ Cross(~edges) + Change(~edges), 
    estimate = "CMLE")

Conditional Maximum Likelihood Results:

             Estimate Std. Error MCMC % z value Pr(>|z|)    
edges         -0.8915     0.1792      0  -4.976   <1e-04 ***
Change~edges  -1.4511     0.1792      0  -8.100   <1e-04 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

     Null Deviance: 424.2  on 306  degrees of freedom
 Residual Deviance: 221.2  on 304  degrees of freedom
 
AIC: 225.2  BIC: 232.7  (Smaller is better. MC Std. Err. = 0)

Simulation

Separable

sep.sim <- simulate(sep.fit, nsim = 1, 
                    time.slices = 1000, 
                    nw.start = "first",
                    monitor = ~edges,
                    stats = T,
                    output = "stats")

Monitor stats: edges

cbind(obs = summary(samp.series ~ edges),
      sim = colMeans(sep.sim$stats))

      obs    sim
edges  57 58.779

target = summary(samp.series ~ edges)
plot(as.numeric(sep.sim$stats), type="l",
     main = "Cross-sectional edgecount in simulated network timeseries",
     sub = "Model: Form(~edges) + Diss(~edges)",
     ylab="edges") 
abline(h=target, col="red")

Model stats

cbind(obs = sep.fit$nw.stats,
      sim = colMeans(sep.sim$stats.gen))

           obs     sim
Form~edges  77  80.188
Diss~edges -35 -37.363

plot(sep.sim$stats.gen)

Joint

jnt.sim <- simulate(jnt.fit, nsim = 1, 
                    time.slices = 1000, 
                    nw.start = "first",
                    monitor = ~edges,
                    stats = T,
                    output = "stats")

Monitor stats: edges

cbind(obs = summary(samp.series ~ edges),
      sim = colMeans(jnt.sim$stats))

      obs    sim
edges  57 59.379

target = summary(samp.series ~ edges)
plot(as.numeric(jnt.sim$stats), type="l",
     main = "Cross-sectional edgecount in simulated network timeseries",
     sub = "Model: Cross(~edges) + Change(~edges)",
     ylab="edges") 
abline(h=target, col="red")

Model stats

cbind(obs = jnt.fit$nw.stats,
      sim = colMeans(jnt.sim$stats.gen))

             obs    sim
edges         57 59.379
Change~edges  42 43.348

plot(jnt.sim$stats.gen)

Add mutuals

Summary

summary(samplk1 ~edges + mutual)

 edges mutual 
    55     14

summary(samp.series ~
          Cross(~edges + mutual) +
          Change(~edges + mutual) +
          Form(~edges + mutual) +
          Diss(~edges + mutual))

        edges        mutual  Change~edges Change~mutual    Form~edges 
           57            15            42             5            77 
  Form~mutual    Diss~edges   Diss~mutual 
           24           -35            -8

Fit

Separable

sep.fit <- tergm(
  samp.series ~
    Form(~edges + mutual) +
    Diss(~edges + mutual),
  estimate = "CMLE",
  snctrl = list(CMLE.MCMC.interval = 1e5)
)

summary(sep.fit)

Call:
tergm(formula = samp.series ~ Form(~edges + mutual) + Diss(~edges + 
    mutual), estimate = "CMLE", snctrl = list(CMLE.MCMC.interval = 1e+05))

Monte Carlo Conditional Maximum Likelihood Results:

            Estimate Std. Error MCMC % z value Pr(>|z|)    
Form~edges   -3.0735     0.3429      0  -8.962   <1e-04 ***
Form~mutual   2.2136     0.5318      0   4.162   <1e-04 ***
Diss~edges   -0.2063     0.3484      0  -0.592    0.554    
Diss~mutual  -1.1674     0.7995      0  -1.460    0.144    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

     Null Deviance: 424.2  on 306  degrees of freedom
 Residual Deviance: 201.6  on 302  degrees of freedom
 
AIC: 209.6  BIC: 224.5  (Smaller is better. MC Std. Err. = 1.333)

mcmc.diagnostics(sep.fit, which="plots")


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

Joint

jnt.fit <- tergm(
  samp.series ~
    Cross(~edges + mutual) +
    Change(~edges + mutual),
  estimate = "CMLE",
  snctrl = list(CMLE.MCMC.interval = 1e5)
)

summary(jnt.fit)

Call:
tergm(formula = samp.series ~ Cross(~edges + mutual) + Change(~edges + 
    mutual), estimate = "CMLE", snctrl = list(CMLE.MCMC.interval = 1e+05))

Monte Carlo Conditional Maximum Likelihood Results:

              Estimate Std. Error MCMC % z value Pr(>|z|)    
edges         -1.46342    0.24938      0  -5.868  < 1e-04 ***
mutual         1.81748    0.51683      0   3.517 0.000437 ***
Change~edges  -1.28808    0.24528      0  -5.251  < 1e-04 ***
Change~mutual -0.05336    0.64386      0  -0.083 0.933953    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

     Null Deviance: 424.2  on 306  degrees of freedom
 Residual Deviance: 211.9  on 302  degrees of freedom
 
AIC: 219.9  BIC: 234.8  (Smaller is better. MC Std. Err. = 1.028)

Simulation

Separable

sep.sim <- simulate(sep.fit, nsim = 1, 
                    time.slices = 1000, 
                    nw.start = "first",
                    monitor = ~edges,
                    stats = T,
                    output = "stats")

Monitor stats: edges

cbind(obs = summary(samp.series ~ edges),
      sim = colMeans(sep.sim$stats))

      obs    sim
edges  57 59.939

target = summary(samp.series ~ edges)
plot(as.numeric(sep.sim$stats), type="l",
     main = "Cross-sectional edgecount in simulated network timeseries",
     sub = "Model: Form(~edges + mutual) + Diss(~edges + mutual)",
     ylab="edges") 
abline(h=target, col="red")

Model stats

cbind(obs = sep.fit$nw.stats,
      sim = colMeans(sep.sim$stats.gen),
      ratio = with(data.frame(obs = sep.fit$nw.stats,
                              sim = colMeans(sep.sim$stats.gen)),
                     obs/sim))

            obs     sim     ratio
Form~edges   77  81.209 0.9481708
Form~mutual  24  26.011 0.9226866
Diss~edges  -35 -38.669 0.9051178
Diss~mutual  -8  -9.201 0.8694707

plot(sep.sim$stats.gen)

Joint

jnt.sim <- simulate(jnt.fit, nsim = 1, 
                    time.slices = 1000, 
                    nw.start = "first",
                    monitor = ~edges,
                    stats = T,
                    output = "stats")

Monitor stats: edges

cbind(obs = summary(samp.series ~ edges),
      sim = colMeans(jnt.sim$stats))

      obs    sim
edges  57 59.776

target = summary(samp.series ~ edges)
plot(as.numeric(jnt.sim$stats), type="l",
     main = "Cross-sectional edgecount in simulated network timeseries",
     sub = "Model: Cross(~edges + mutual) + Change(~edges + mutual)",
     ylab="edges") 
abline(h=target, col="red")

Model stats

cbind(obs = jnt.fit$nw.stats,
      sim = colMeans(jnt.sim$stats.gen),
      ratio = with(data.frame(obs = jnt.fit$nw.stats,
                              sim = colMeans(jnt.sim$stats.gen)),
                     obs/sim))

              obs    sim     ratio
edges          57 59.776 0.9535600
mutual         15 16.399 0.9146899
Change~edges   42 42.754 0.9823642
Change~mutual   5  5.301 0.9432183

plot(jnt.sim$stats.gen)

Temporal net data structures and comparisons

Intro

Setup & data

Types of networks

first network

list

NetSeries

networkDynamic

Descriptives

Traditional

Term counts

Edges only models

Summary

Fit

Separable

Joint

Simulation

Separable

Monitor stats: edges

Model stats

Joint

Monitor stats: edges

Model stats

Add mutuals

Summary

Fit

Separable

Joint

Simulation

Separable

Monitor stats: edges

Model stats

Joint

Monitor stats: edges

Model stats