ergm model running problems after update of the R version

[tmorguet]

Dear authors of the ergm package,

we have enthusiastically used your ergm package in R version 3.6.2 of 19/10/20 and RStudio version 1.2.5033 for our empirical studies and gained valuable insights.

In our most comprehensive ergm model including structural effects (~edges+mutual+gwidegree+gwodegree+gwesp+gwdsp), actor relation effects (+absdiff("")+nodeocov("")+nodeicov("")+nodematch("")), and control=control.ergm(MCMC.burnin = 40000, MCMC.samplesize = 4000, MCMC.interval = 4000, MCMLE.maxit = 21)), the step length converged twice after 16 iterations and the run time is 1.2 hours.
Unfortunately, after an update of R to version 4.1.2 it happened that the model did not start the first iteration despite it started the estimation. The same was the case when testing the codes on other computers with R versions 4.1.0 and 4.1.2.

Do you perhaps have an explanation why the estimation of our comprehensive ergm model does not run anymore or a hint what to consider when updating R regarding the ergm package? Does the ergm code have to be adapted in any way to run on newer R versions without problems?

The reason why we have to use a current R version is the desire to calculate further analyses based on our ergm models, specifically multilevel ergms (package mlergm), process analyses for ergms (package ergMargins), and frailty ergms (package fergm). Any additional advice on what to look for technically when calculating these models?

Any feedback would be greatly appreciated!

@handcock @drh20drh20 @CarterButts @sgoodreau @krivit @martinamorris

[mbojan]

Thanks @tmorguet . We will need more details about the analysis you are running in order to reproduce the problem. Please provide:

1. Exact call to ergm() you are running.
2. Output of sessionInfo() (we need to know the package versions)
3. Can you provide the network object you are trying to model?

[tmorguet]

Thank you @mbojan for your fast response and your willingness to help! I'm happy to provide you with all requested detailed information regarding our model:

Model_03 <- ergm(Edgelist_Network~edges # Arc
		+mutual #Reciprocity 
                 +gwidegree(decay=0.2, fixed=T, attr=NULL,cutoff=30,levels=NULL) #Popularity spread (indegree)
                 +gwodegree(decay=0.2, fixed=T, attr=NULL,cutoff=30,levels=NULL) #Acitivty spread (outdegree)
                 +gwesp(decay=0.2, fixed=T, cutoff=30) #Transitive Closure
                 +gwdsp(decay=0.2, fixed=T, cutoff=30) # Mutiple 2-path (multiple connectivity)
                 
                 #Control variable
                 +nodematch("var1")
                 
                 #Actor relation effects
                 +absdiff("var2")+nodeocov("var2")+nodeicov("var2")
                 +absdiff("var3")+nodeocov("var3")+nodeicov("var3")
                 +absdiff("var4")+nodeocov("var4")+nodeicov("var4")
                 +absdiff("var5")+nodeocov("var5")+nodeicov("var5")
                 , control=control.ergm(MCMC.burnin = 40000, MCMC.samplesize = 4000, MCMC.interval = 4000, MCMLE.maxit = 21))

R version 4.1.2 (2021-11-01)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 10 x64 (build 19043)

Matrix products: default

locale:
[1] LC_COLLATE=German_Germany.1252  LC_CTYPE=German_Germany.1252    LC_MONETARY=German_Germany.1252 LC_NUMERIC=C                   
[5] LC_TIME=German_Germany.1252    

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] fergm_1.1.4           rstan_2.21.2          ggplot2_3.3.5         StanHeaders_2.21.0-7  ergMargins_0.1.3      mlergm_0.8           
 [7] statnet_2019.6        tsna_0.3.5            sna_2.6               statnet.common_4.5.0  ergm.count_4.0.2      tergm_4.0.2          
[13] networkDynamic_0.11.0 ergm_4.1.2            network_1.17.1        readxl_1.3.1          reshape2_1.4.4        haven_2.4.3          

loaded via a namespace (and not attached):
 [1] jsonlite_1.7.2          RcppParallel_5.1.4      stats4_4.1.2            cellranger_1.1.0        robustbase_0.93-9       Rttf2pt1_1.3.9         
 [7] pillar_1.6.4            lattice_0.20-45         glue_1.5.1              extrafontdb_1.0         RColorBrewer_1.1-2      colorspace_2.0-2       
[13] cowplot_1.1.1           Matrix_1.3-4            plyr_1.8.6              lpSolve_5.6.15          pkgconfig_2.0.3         purrr_0.3.4            
[19] scales_1.1.1            processx_3.5.2          tzdb_0.2.0              tibble_3.1.6            generics_0.1.1          ellipsis_0.3.2         
[25] cachem_1.0.6            withr_2.4.3             cli_3.1.0               magrittr_2.0.1          crayon_1.4.2            memoise_2.0.1          
[31] ps_1.6.0                GGally_2.1.2            fansi_0.5.0             nlme_3.1-153            MASS_7.3-54             forcats_0.5.1          
[37] pkgbuild_1.3.0          tools_4.1.2             loo_2.4.1               prettyunits_1.1.1       hms_1.1.1               lifecycle_1.0.1        
[43] matrixStats_0.61.0      stringr_1.4.0           V8_3.6.0                trust_0.1-8             munsell_0.5.0           callr_3.7.0            
[49] compiler_4.1.2          rlang_0.4.12            grid_4.1.2              gtable_0.3.0            codetools_0.2-18        inline_0.3.19          
[55] reshape_0.8.8           curl_4.3.2              R6_2.5.1                lpSolveAPI_5.5.2.0-17.7 gridExtra_2.3           rle_0.9.2              
[61] dplyr_1.0.7             extrafont_0.17          fastmap_1.1.0           utf8_1.2.2              readr_2.1.1             stringi_1.7.6          
[67] parallel_4.1.2          Rcpp_1.0.7              vctrs_0.3.8             DEoptimR_1.0-9          tidyselect_1.1.1        coda_0.19-4  

function (x, vertex.attr = NULL, vertex.attrnames = NULL, directed = TRUE, 
    hyper = FALSE, loops = FALSE, multiple = FALSE, bipartite = FALSE, 
    ...) 
{
    g <- as.network(x, directed = directed, hyper = hyper, loops = loops, 
        multiple = multiple, bipartite = bipartite, ...)
    if (!is.null(vertex.attr)) {
        if (is.null(vertex.attrnames)) {
            if (!is.null(names(vertex.attr))) 
                vertex.attrnames <- names(vertex.attr)
            else {
                vertex.attrnames <- 1:length(vertex.attr)
                warning("Vertex attribute names not given; making some up.")
            }
        }
        for (i in 1:length(vertex.attr)) g <- set.vertex.attribute(g, 
            vertex.attrnames[[i]], vertex.attr[[i]])
    }
    g
}

Please let me know if I can provide you with more details. Otherwise I would be grateful for any hints to improve the model specification of our model_03 so that it runs.

Thank you again and kind regards

[mbojan]

Thanks @tmorguet . For item 3 I meant can you please upload the network data object you are trying to model (Edgelist_network in your code).

Other than that, that's quite a complex model you are fitting so there are a lot of potential issues with respect to model specification versus data that can lead to the behavior you are describing. These are not necessarily technically a bug in the software. Did you try fitting simpler models? Did they converge?

[mbojan]

Converting to discussion. Will create issue(s) as soon as one is identified.

[CarterButts]

Just as an observation (based on ongoing experience with the new 4.x platform), but it does appear that the new code is (1) faster (and perhaps higher quality?) when it works, but (2) overzealous about both hanging (when it doesn't think that the system has converged) and giving up (when it concludes that mixing is too slow and/or there is too much multicollinearity). Changing the convergence detection back to Hotelling can sometimes help, but that still doesn't reproduce 3.x behavior. My own experience is that switching to stochastic approximation (see the control.ergm) often produces reasonable results, although the default SA settings may require tuning. The settings associated with the old Robins-Monro option used to work extremely well, but unfortunately that seems to be broken right now. Another option to consider, for folks who really want Geyer-Thompson-Hummel solutions rather than SA, is to (1) fit using SA with reasonable settings (be sure to turn off the bridge sampler, which wastes time here), and then (2) use the SA solution with the init argument to control.ergm to now fit the model using the usual G-T-H MCMCMLE scheme. SA can be very fast, and G-T should converge quickly once one is near the MLE. This may not work for every model, but if your problem is that the MCMC is getting bogged down during the stepping process (where it is moving you from the MPLE/CD estimate towards the MLE), using SA as an additional approximation step to get near the MLE may fix it. In other cases, SA may simply be good enough as-is. There's nothing inherently wrong with it (indeed, this is what p-net has always used), we just haven't put as much time and effort into refining the ergm implementation as has gone into the implementation of G-T-H.