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?
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:
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