rsnastatnetergm

Obtain average homophily coefficients from ERGM for meta-regression


I have run the following ERGM model on several classes:

mFriend<-ergm(networks$`network_id`~edges+mutual+ 
                  nodeofactor('sex') + nodeofactor('minor') +
                  nodemix ('sex', base = c(2,3))+nodemix ('minor', base = c(2,3)) + absdiff('gpa'),
                control = control.ergm(seed=1, MCMC.burnin=50000, MCMC.interval=5000), verbose=TRUE)

And have the following summary for one class:

Monte Carlo Maximum Likelihood Results:

                    Estimate Std. Error MCMC % z value Pr(>|z|)    
edges                -3.3709     0.7549      0  -4.466  < 1e-04 ***
mutual                2.7828     0.4667      0   5.963  < 1e-04 ***
nodeofactor.sex.1    -0.8513     1.0916      0  -0.780 0.435492    
nodeofactor.minor.1   0.2849     0.9734      0   0.293 0.769728    
mix.sex.0.0           2.2518     0.6215      0   3.623 0.000291 ***
mix.sex.1.1           3.2358     0.8129      0   3.980  < 1e-04 ***
mix.minor.0.0         0.1571     0.5563      0   0.282 0.777615    
mix.minor.1.1           -Inf     0.0000      0    -Inf  < 1e-04 ***
absdiff.gpa          -0.5670     0.3446      0  -1.645 0.099939 .  

The next step of my work is to do meta-regression with these classes and add covariates. To do so I need to obtain average homophily coefficients and std.errors for each class and add them to meta-regression:

rma(homophily_estimates, homophily_errors, mods = ~ minor_perc+class_size+mean_income, data=dfmeta2)

I am, however, not sure on how to proceed since this is my first time doing meta-regression. Do you have any advice?


Solution

  • I presume you are aiming at an analysis similar to e.g. Lubbers & Snijders in https://doi.org/10.1016/j.socnet.2007.03.002 ? I'm not sure what precisely you mean by "how to proceed?", but here's and attempt at guessing it :) I'll edit my answer if the guess is not exactly right. The steps are essentially:

    1. Fit the model(s) to all the networks.
    2. Create a data frame with model coefficients
    3. Join with the class-level variables that you have (income, class size etc.)
    4. Fit the meta-regression

    For example something like that:

    library(ergm)
    #> Loading required package: network
    library(purrr)
    library(dplyr)
    
    # Let's pretend the three Sampson's networks are three classes
    data(samplk)
    networks <- list(samplk1, samplk2, samplk3)
    
    # Fit the models to every network
    l <- map(networks, ~ ergm(.x ~ edges + nodematch("group")))
    #> Starting maximum pseudolikelihood estimation (MPLE):
    # (...)
    #> Evaluating log-likelihood at the estimate.
    
    # Pretend the list of networks is named which allows
    # to bind_rows() with `.id` argument
    names(l) <- paste0("cl", seq(along = l))
    
    # Data frame of coefficients
    d <- map(l, broom::tidy) |> 
      bind_rows(.id = "cl")
    d
    #> # A tibble: 6 × 7
    #>   cl    term            estimate std.error mcmc.error statistic  p.value
    #>   <chr> <chr>              <dbl>     <dbl>      <dbl>     <dbl>    <dbl>
    #> 1 cl1   edges              -2.11     0.212          0     -9.98 1.79e-23
    #> 2 cl1   nodematch.group     1.73     0.318          0      5.45 5.06e- 8
    #> 3 cl2   edges              -2.42     0.239          0    -10.1  5.86e-24
    #> 4 cl2   nodematch.group     2.47     0.334          0      7.40 1.34e-13
    #> 5 cl3   edges              -2.48     0.245          0    -10.1  6.28e-24
    #> 6 cl3   nodematch.group     2.53     0.338          0      7.48 7.34e-14
    
    # Artificial class-level data to be added
    classd <- data.frame(
      cl = c("cl1", "cl2", "cl3"),
      income = rnorm(3),
      class_size = c(15, 21, 18)
    )
    
    # Join model coefs with class-level data
    dd <- left_join(d, classd, by = "cl")
    
    # Fit the meta-regression
    metafit <- metafor::rma(
      estimate ~ income, 
      std.error^2, 
      data = dd,
      subset = term == "nodematch.group"
    )
    summary(metafit)
    #> 
    #> Mixed-Effects Model (k = 3; tau^2 estimator: REML)
    #> 
    #>   logLik  deviance       AIC       BIC      AICc   
    #>  -0.5423    1.0846    7.0846    1.0846   31.0846   
    #> 
    #> tau^2 (estimated amount of residual heterogeneity):     0.0625 (SE = 0.2450)
    #> tau (square root of estimated tau^2 value):             0.2500
    #> I^2 (residual heterogeneity / unaccounted variability): 36.08%
    #> H^2 (unaccounted variability / sampling variability):   1.56
    #> R^2 (amount of heterogeneity accounted for):            35.50%
    #> 
    #> Test for Residual Heterogeneity:
    #> QE(df = 1) = 1.5646, p-val = 0.2110
    #> 
    #> Test of Moderators (coefficient 2):
    #> QM(df = 1) = 1.3664, p-val = 0.2424
    #> 
    #> Model Results:
    #> 
    #>          estimate      se    zval    pval    ci.lb   ci.ub      
    #> intrcpt    2.0185  0.3004  6.7190  <.0001   1.4297  2.6073  *** 
    #> income     0.3102  0.2654  1.1689  0.2424  -0.2099  0.8304      
    #> 
    #> ---
    #> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1