Intro

We work with 2 panels of the Sampson data to:

Fit and simulate from the single Operator(~edges) models, to see what they produce in equilibrium
Verify that the estimated coefs = logit(p(component change))

Setup

library(tergm)
#library(tsna)
library(dplyr)

sessioninfo::package_info(pkgs = c("network",
                                   "ergm", "tergm"),
                          dependencies = FALSE)

 package * version date       lib source        
 ergm    * 4.0.1   2021-06-21 [1] CRAN (R 4.0.5)
 network * 1.17.1  2021-06-14 [1] CRAN (R 4.0.5)
 tergm   * 4.0.1   2021-06-24 [1] CRAN (R 4.0.5)

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

Data

data(samplk)
samp.series <- NetSeries(list(samplk1,samplk2))
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.dyn <- networkDynamic(network.list = list(samplk1,samplk2), 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

Descriptives

Cross sectional stats and form/diss changes from tsna

samp.stats <- cbind(edges = tsna::tErgmStats(samp.dyn, "~edges"),
                    form = tsna::tEdgeFormation(samp.dyn),
                    diss = tsna::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

dyads = network.dyadcount(samplk1)
type <- c("form", "diss", "cross", "change")

prob.form = samp.stats[2,"form"]/(dyads-samp.stats[1,"edges"])
prob.diss = samp.stats[2,"diss"]/network.edgecount(samplk1)
prob.cross = network.edgecount(samplk2)/dyads
prob.change = network.edgecount(xor(samplk1, samplk2))/dyads

def.form = "Fraction of empty dyads at t1 that formed a tie by t2"
def.diss = "Fraction of ties at t1 that dissolved by t2"
def.cross = "Density at t2 (probability of a tie)"
def.change = "Fraction of all dyads at t1 that changed state by t2"

data.frame(type = type,
           defn = c(def.form, def.diss, def.cross, def.change),
           prob = c(prob.form, prob.diss, prob.cross, prob.change)) %>%
  kableExtra::kable(caption = "Probabilities of tie change components", 
        digits=3) %>%
  kableExtra::kable_styling(full_width=F,
                            bootstrap_options = c("striped"))

Probabilities of tie change components
type	defn	prob
form	Fraction of empty dyads at t1 that formed a tie by t2	0.088
diss	Fraction of ties at t1 that dissolved by t2	0.364
cross	Density at t2 (probability of a tie)	0.186
change	Fraction of all dyads at t1 that changed state by t2	0.137

# target edgecount
target.edges <- summary(samp.series ~ edges)

Models

Form

This model sets the tie formation rates based on the observed probability of tie formation from t1 to t2. But the other free parameter is not set, so the effective tie dissolution coef is 0, which means ties dissolve with probability 0.5 (compared to the observed probability of 0.36).

Fit

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

summary(form.fit)

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

Conditional Maximum Likelihood Results:

           Estimate Std. Error MCMC % z value Pr(>|z|)    
Form~edges  -2.3427     0.2232      0   -10.5   <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: 225.4  on 305  degrees of freedom
 
AIC: 227.4  BIC: 231.1  (Smaller is better. MC Std. Err. = 0)

Simulate

The net result is a lower predicted equilibrium tie count, compared to t2 (57)

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

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

print(form.mon <- cbind(obs = target.edges,
                        sim = colMeans(form.sim$stats)))

      obs    sim
edges  57 46.245

Diss

This model sets the tie dissolution rates based on the observed probability of tie formation from t1 to t2. But the other free parameter is not set, so the effective tie formation coef is 0, which means ties form with probability 0.5. (compared to the observed probability of 0.09).

Fit

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

summary(diss.fit)

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

Conditional Maximum Likelihood Results:

           Estimate Std. Error MCMC % z value Pr(>|z|)  
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: 420.1  on 305  degrees of freedom
 
AIC: 422.1  BIC: 425.8  (Smaller is better. MC Std. Err. = 0)

Simulate

The net result is a much higher predicted equilibrium tie count, compared to t2 (57).

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


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

print(diss.mon <- cbind(obs = target.edges,
                        sim = colMeans(diss.sim$stats)))

      obs     sim
edges  57 177.268

Cross

This model sets the target cross-sectional tie count to the t2 observed value. Both the formation and dissolution coefficients are effectively set to 0 in this case, so the rate of tie formation, and dissolution, is 0.5.

Fit

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

summary(cross.fit)

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

Conditional Maximum Likelihood Results:

      Estimate Std. Error MCMC % z value Pr(>|z|)    
edges  -1.4744     0.1468      0  -10.04   <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: 294.2  on 305  degrees of freedom
 
AIC: 296.2  BIC: 300  (Smaller is better. MC Std. Err. = 0)

Simulate

Here, the equilibrium tie count hovers around the t2 value.

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

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

print(cross.mon <- cbind(obs = target.edges,
                         sim = colMeans(cross.sim$stats)))

      obs    sim
edges  57 57.206

Change

This model sets the dyad change rates based on the observed probability of of a dyad changing state from t1 to t2. There is no bias towards formation or dissolution, so if the starting network has density less than 0.5, you have more potential formations, and if greater, you have more potential dissolutions.

Fit

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

summary(change.fit)

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

Conditional Maximum Likelihood Results:

             Estimate Std. Error MCMC % z value Pr(>|z|)    
Change~edges  -1.8383     0.1661      0  -11.07   <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: 244.8  on 305  degrees of freedom
 
AIC: 246.8  BIC: 250.5  (Smaller is better. MC Std. Err. = 0)

Simulate

The equilibrium density for this model is 0.5, well above the observed density at t2 (r round(`prob.cross, 2)).

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


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

print(change.mon <- cbind(obs = target.edges,
                          sim = colMeans(change.sim$stats)))

      obs     sim
edges  57 153.256

Compare coefficients

Here, we verify that the estimated exp(coefficient) values are equal to the odds of each tie component as calculated by hand.

They are equal, which means we can interpret the coefficients in terms of the probabilities of each change component. The coefficients are on the log-odds scale, so to obtain the corresponding probability, we transform using p(tie change component) = 1/(1+exp(coef))).

As with ergm coefficients, this can be interpreted in terms of the conditional probability at the dyad level.

In the dynamic context, this is the probability of a dyad changing state, if that change satisfies the definition of this component.

odds.form = prob.form / (1-prob.form)
odds.diss = prob.diss / (1-prob.diss)
odds.cross = prob.cross / (1-prob.cross)
odds.change = prob.change / (1-prob.change)

ecoef.form = exp(coef(form.fit))
ecoef.diss = exp(coef(diss.fit))
ecoef.cross = exp(coef(cross.fit))
ecoef.change = exp(coef(change.fit))

data.frame(
  exp.coef = c(ecoef.form, ecoef.diss, ecoef.cross, ecoef.change),
  handcalc.odds = c(odds.form, odds.diss, odds.cross, odds.change),
  handcalc.prob = c(prob.form, prob.diss, prob.cross, prob.change)) %>%
  kableExtra::kable(
    caption = "Estimated exp(coef) vs. hand calculated odds and probabilities", 
    digits=3) %>%
  kableExtra::kable_styling(bootstrap_options = c("striped"))

Estimated exp(coef) vs. hand calculated odds and probabilities
	exp.coef	handcalc.odds	handcalc.prob
Form~edges	0.096	0.096	0.088
Diss~edges	0.571	0.571	0.364
edges	0.229	0.229	0.186
Change~edges	0.159	0.159	0.137

coef <- c(NA_real_,
          coef(form.fit),
          coef(diss.fit),
          coef(cross.fit),
          coef(change.fit))
mon.edges <- c(form.mon[1],
               form.mon[2],
               diss.mon[2],
               cross.mon[2],
               change.mon[2])
names(coef)[1] <- "obs"

data.frame(coef=round(coef,1), edges=round(mon.edges)) %>% 
  kableExtra::kable(caption = "Compare coefs and obs to equilibrium edgecounts", 
        digits=3) %>%
  kableExtra::kable_styling(bootstrap_options = c("striped"))

Compare coefs and obs to equilibrium edgecounts
	coef	edges
obs	NA	57
Form~edges	-2.3	46
Diss~edges	-0.6	177
edges	-1.5	57
Change~edges	-1.8	153

Single Operator(~edges) CMLE

Intro

Setup

Data

Descriptives

Models

Form

Fit

Simulate

Diss

Fit

Simulate

Cross

Fit

Simulate

Change

Fit

Simulate

Compare coefficients