Why won’t this simple model on a small network converge?

Part 1: fitting with the latest ergm version

# install.packages("ergm")
library(ergm)
# install.packages("ergm")
packageVersion('ergm')
## [1] '4.7.5'

First let’s confirm this network exists and is not extreme.

nw <- network.initialize(118, directed=FALSE)
nw <- san(nw~edges + degree(0), target.stats = c(173, 43))
summary(nw~edges + degree(0))
##   edges degree0 
##     173      43
plot(nw)

Yup, it exists. And one can clearly see that it could have many more or fewer isolates at this density, or much higher or lower density with this number of isolates.

Fit the desired model with all default settings

set.seed(0)
nw <- network.initialize(n = 118, directed=F)
fit <- ergm(nw~edges + degree(0), target.stats = c(173, 43)) 
## Error in ergm.MCMLE(init, s, s.obs, control = control, verbose = verbose, : Unconstrained MCMC sampling did not mix at all. Optimization cannot continue.

We get an error. The full output is suppressed, but it is filled with lines like: The log-likelihood improved by 7655885619538640582404466006606.0000.

What happens if we add in some of the fitting parameters that have been discussed with other users’ ergm 4. convergence issues?

mycontrol <- control.ergm(
              MCMLE.maxit = 500,  
              main.method = c("Stochastic-Approximation"),
              MCMC.interval = 1e5,
              MCMC.samplesize = 1e6,
              MCMLE.termination = "Hotelling",
              MCMC.effectiveSize=NULL,
              SAN = control.san(
                SAN.maxit = 500, 
                SAN.nsteps = 1e8
              )
            )

fit <- ergm(nw~edges + degree(0), target.stats = c(173, 43), control = mycontrol)
## Starting maximum pseudolikelihood estimation (MPLE):
## Obtaining the responsible dyads.
## Evaluating the predictor and response matrix.
## Maximizing the pseudolikelihood.
## Finished MPLE.
## Stochastic approximation algorithm with theta_0 equal to:
##     edges   degree0 
## -2.694787  4.711083 
## Starting burnin of 1600000 steps
## Phase 1: 200 steps (interval = 100000)
## Stochastic Approximation estimate:
##     edges   degree0 
## -493.6500 -344.2563
## Phase 3:  1000 iterations (interval=1e+05)
## Error in snearPD(V): Matrix 'x' has negative elements on the diagonal.

The Phase 3 coefficients make it clear that the estimation has gone off track, even before the error appears. The error itself is hard to interpret (at least for me).

How do we do with fewer isolates?

The examples below consider the same model, but with 10, 20, 30, and 35 isolates, instead of the observed 43. I explored these cases to get a sense of how this particular model behaves as it moves further away from the expected level of isolates. Summary: 10 and 20 do good, 30 maybe starts drifting away a bit, 35 returns the same negative matrix error as 43.

fit <- ergm(nw~edges + degree(0), target.stats = c(173, 10), control = mycontrol)
## Starting maximum pseudolikelihood estimation (MPLE):
## Obtaining the responsible dyads.
## Evaluating the predictor and response matrix.
## Maximizing the pseudolikelihood.
## Finished MPLE.
## Stochastic approximation algorithm with theta_0 equal to:
##      edges    degree0 
## -3.6268625  0.2226657 
## Starting burnin of 1600000 steps
## Phase 1: 200 steps (interval = 100000)
## Stochastic Approximation estimate:
##     edges   degree0 
## -3.580928  0.637440
## Phase 3:  1000 iterations (interval=1e+05)
## 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.
fit
## 
## Call:
## ergm(formula = nw ~ edges + degree(0), target.stats = c(173, 
##     10), control = mycontrol)
## 
## Stochastic Approximation Maximum Likelihood Coefficients:
##   edges  degree0  
## -3.5755   0.6749
goffit <- gof(fit, control=control.gof.ergm(nsim=1e2))
goffit$summary.model
##         obs min   mean max MC p-value
## edges   173 132 174.61 207       0.92
## degree0  10   2   9.73  21       0.96
plot(goffit)

Fits quite well.

fit <- ergm(nw~edges + degree(0), target.stats = c(173, 20), control = mycontrol)
## Starting maximum pseudolikelihood estimation (MPLE):
## Obtaining the responsible dyads.
## Evaluating the predictor and response matrix.
## Maximizing the pseudolikelihood.
## Finished MPLE.
## Stochastic approximation algorithm with theta_0 equal to:
##     edges   degree0 
## -3.372197  1.566150 
## Starting burnin of 1600000 steps
## Phase 1: 200 steps (interval = 100000)
## Stochastic Approximation estimate:
##     edges   degree0 
## -3.357320  1.824547
## Phase 3:  1000 iterations (interval=1e+05)
## 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.
fit
## 
## Call:
## ergm(formula = nw ~ edges + degree(0), target.stats = c(173, 
##     20), control = mycontrol)
## 
## Stochastic Approximation Maximum Likelihood Coefficients:
##   edges  degree0  
##  -3.346    1.845
goffit <- gof(fit, control=control.gof.ergm(nsim=1e2))
goffit$summary.model
##         obs min   mean max MC p-value
## edges   173 115 177.18 217       0.76
## degree0  20   9  19.47  37       0.90
plot(goffit)

Still does pretty good, although that edge count might be worth keeping an eye on?

fit <- ergm(nw~edges + degree(0), target.stats = c(173, 30), control = mycontrol)
## Starting maximum pseudolikelihood estimation (MPLE):
## Obtaining the responsible dyads.
## Evaluating the predictor and response matrix.
## Maximizing the pseudolikelihood.
## Finished MPLE.
## Stochastic approximation algorithm with theta_0 equal to:
##     edges   degree0 
## -3.119518  2.306999 
## Starting burnin of 1600000 steps
## Phase 1: 200 steps (interval = 100000)
## Stochastic Approximation estimate:
##     edges   degree0 
## -3.193341  2.394727
## Phase 3:  1000 iterations (interval=1e+05)
## 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.
fit
## 
## Call:
## ergm(formula = nw ~ edges + degree(0), target.stats = c(173, 
##     30), control = mycontrol)
## 
## Stochastic Approximation Maximum Likelihood Coefficients:
##   edges  degree0  
##  -3.169    2.480
goffit <- gof(fit, control=control.gof.ergm(nsim=1e2))
goffit$summary.model
##         obs min   mean max MC p-value
## edges   173   0 178.67 253       0.66
## degree0  30   9  27.34 118       0.48
plot(goffit)

Getting off a bit more?

fit <- ergm(nw~edges + degree(0), target.stats = c(173, 35), control = mycontrol)
## Unable to match target stats. Using MCMLE estimation.
## Starting maximum pseudolikelihood estimation (MPLE):
## Obtaining the responsible dyads.
## Evaluating the predictor and response matrix.
## Warning in mple.existence(pl): The MPLE does not exist!
## Maximizing the pseudolikelihood.
## Finished MPLE.
## Stochastic approximation algorithm with theta_0 equal to:
##     edges   degree0 
## -2.974803 18.060191 
## Starting burnin of 1600000 steps
## Phase 1: 200 steps (interval = 100000)
## Stochastic Approximation estimate:
##     edges   degree0 
## -2.269612 17.721862
## Phase 3:  1000 iterations (interval=1e+05)
## Error in snearPD(V): Matrix 'x' has negative elements on the diagonal.

That pesky negative matrix error again.