Correlation between in- and out degree

[MikeGacho]

I have a directed, binary network:

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

Vertex attribute names:
background_colony IF_bee_neighbour IF_neighbour_bee vertex.names

Edge attribute names not shown

And I'm fitting the following model :

model1 <- ergm(net~ edges + gwidegree(3.26, fixed = T) + mutual + gwesp(0.1, fixed=T) +
nodematch("background_colony", diff=TRUE) +
nodecov("IF_bee_neighbour") +
nodecov("IF_neighbour_bee") +
edgecov(net, "codetection"),
reference = ~Bernoulli,
control = control.ergm(#MCMLE.steplength = 0.5,
MCMC.burnin = 20000,
MCMC.samplesize = 5000,
parallel = 4, parallel.type = "PSOCK"),
verbose = TRUE)

Model1 converges successfully when I add either gwodegree or gwidegree. When I try to fit a model2 adding both terms I get:

Error in ergm.MCMLE(init, s, s.obs, control = control, verbose = verbose, :
Number of edges in a simulated network exceeds that in the observed by a factor of more than 20. This is a strong indicator of model degeneracy or a very poor starting parameter configuration. If you are reasonably certain that neither of these is the case, increase the MCMLE.density.guard control.ergm() parameter.

I increased MCMLE.density.guard.min to 300000 but model2 remains stuck at iteration1

I suppose the issue is due to the strong positive correlation between in- and out-degree in my network (individuals with more incoming ties also have more outgoing ties). How should I deal with this?
Basically I want to understand if the effects obtained for mutuality, triangles, nodal attributes, homophily in model1 are reliable when only gwodegree or gwidegree are included separately. Does model1 takes into account the correlation?

Any help appreciated!

[krivit]

Can you provide the formula for model2?

[MikeGacho]

model2 <- ergm(net~ edges + gwidegree(3.26, fixed = T) + gwodegree(3.26, fixed = T) + mutual + gwesp(0.1, fixed=T) +
nodematch("background_colony", diff=TRUE) +
nodecov("IF_bee_neighbour") +
nodecov("IF_neighbour_bee") +
edgecov(net, "codetection"),
reference = ~Bernoulli,
control = control.ergm(
MCMC.burnin = 20000,
MCMC.samplesize = 5000,
parallel = 4, parallel.type = "PSOCK"),
verbose = TRUE)

The way I estimated gw*degree decays was by fitting first two simple models:

model_ind <- ergm(net~ edges + gwidegree(3, fixed = FALSE, cutoff=200),
control = control.ergm(
MCMC.burnin = 20000,
MCMC.samplesize = 5000,
parallel = 4, parallel.type = "PSOCK"),
verbose = TRUE)

model_outd <- ergm(net~ edges + gwodegree(3, fixed = FALSE, cutoff=200),
control = control.ergm(
MCMC.burnin = 20000,
MCMC.samplesize = 5000,
parallel = 4, parallel.type = "PSOCK"),
verbose = TRUE)

[krivit]

Can you try estimating the degrees using curved modelling (i.e., just gwodegree()/gwidegree())? It works surprisingly well in some situations.

[MikeGacho]

It works well and fast with gwodegree but it hangs with gwidegree. Possible reasons?

gwodegree only

model_outd <-ergm(net~ edges + gwodegree(cutoff=120))
Starting maximum pseudolikelihood estimation (MPLE):
Obtaining the responsible dyads.
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.3382.
The log-likelihood improved by 3.3514.
Estimating equations are not within tolerance region.
Iteration 2 of at most 60:
Optimizing with step length 0.6287.
The log-likelihood improved by 2.9700.
Estimating equations are not within tolerance region.
Iteration 3 of at most 60:
Optimizing with step length 1.0000.
The log-likelihood improved by 1.1984.
Estimating equations are not within tolerance region.
Iteration 4 of at most 60:
Optimizing with step length 1.0000.
The log-likelihood improved by 0.2450.
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.1972.
Estimating equations are not within tolerance region.
Iteration 6 of at most 60:
Optimizing with step length 1.0000.
The log-likelihood improved by 0.0257.
Convergence test p-value: 0.5963. 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.0582.
Convergence test p-value: 0.4161. 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.0216.
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.

summary(model_outd)
ergm(formula = net ~ edges + gwodegree(cutoff = 120))

Monte Carlo Maximum Likelihood Results:

Estimate Std. Error MCMC % z value Pr(>|z|)
edges -1.96959 0.03960 0 -49.73 <1e-04 ***
gwodegree -3.29971 0.02610 0 -126.42 <1e-04 ***
gwodegree.decay 3.26087 0.03223 0 101.17 <1e-04 ***

Signif. codes: 0 ‘’ 0.001 ‘’ 0.01 ‘’ 0.05 ‘.’ 0.1 ‘ ’ 1

Null Deviance: 537197 on 387506 degrees of freedom
Residual Deviance: 70415 on 387503 degrees of freedom

AIC: 70421 BIC: 70453 (Smaller is better. MC Std. Err. = 19.84)

gwidegree only. This gets stuck at iteration 2 for hours so I stopped

model_ind <-ergm(net~ edges + gwidegree(cutoff=250))
Starting maximum pseudolikelihood estimation (MPLE):
Obtaining the responsible dyads.
Evaluating the predictor and response matrix.
Maximizing the pseudolikelihood.
Finished MPLE.
Starting Monte Carlo maximum likelihood estimation (MCMLE):
Iteration 1 of at most 60:
Model statistics ‘gwidegree’ are not varying. This may indicate that the observed data occupies an extreme point in the sample space or that the estimation has reached a dead-end configuration.
Warning: Model statistics ‘gwidegree’ are linear combinations of some set of preceding statistics at the current stage of the estimation. This may indicate that the model is nonidentifiable.
Optimizing with step length 0.0066.
The log-likelihood improved by < 0.0001.
Estimating equations are not within tolerance region.
Iteration 2 of at most 60:
Warning message:
In ergm_MCMC_sample(s, control, theta = mcmc.init, verbose = max(verbose - :
Unable to reach target effective size in iterations alotted.

I tried the same with a different network also having correlation between in- and out-degree an the behavior is the same. Able to estimate gwodegree but not gwidegree.