从 ERGM 获取平均同质系数以进行元回归

问题描述 投票:0回答:1

我在几个班级上运行了以下 ERGM 模型:

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)

并对一堂课进行以下总结:

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 .  

我的下一步工作是对这些类进行元回归并添加协变量。为此,我需要获得每个类的平均同质系数和 std.errors 并将它们添加到元回归中:

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

但是,我不确定如何继续,因为这是我第一次进行元回归。你有什么建议吗?

r social-networking statnet ergm
1个回答
0
投票

我认为您的目标是进行类似于以下的分析:吕贝尔斯和斯奈德斯在https://doi.org/10.1016/j.socnet.2007.03.002?我不确定你所说的“如何进行?”到底是什么意思,但这里是并尝试猜测它:)如果猜测不完全正确,我将编辑我的答案。步骤基本上是:

  1. 使模型适合所有网络。
  2. 创建具有模型系数的数据框
  3. 加入您拥有的班级级别变量(收入、班级人数等)
  4. 拟合元回归

例如类似的事情:

library(ergm)
#> Loading required package: network
#> 
#> 'network' 1.18.2 (2023-12-04), part of the Statnet Project
#> * 'news(package="network")' for changes since last version
#> * 'citation("network")' for citation information
#> * 'https://statnet.org' for help, support, and other information
#> 
#> 'ergm' 4.6.0 (2023-12-17), part of the Statnet Project
#> * 'news(package="ergm")' for changes since last version
#> * 'citation("ergm")' for citation information
#> * 'https://statnet.org' for help, support, and other information
#> 'ergm' 4 is a major update that introduces some backwards-incompatible
#> changes. Please type 'news(package="ergm")' for a list of major
#> changes.
library(purrr)
library(dplyr)
#> 
#> Attaching package: 'dplyr'
#> The following objects are masked from 'package:stats':
#> 
#>     filter, lag
#> The following objects are masked from 'package:base':
#> 
#>     intersect, setdiff, setequal, union

# 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):
#> Obtaining the responsible dyads.
#> Evaluating the predictor and response matrix.
#> Maximizing the pseudolikelihood.
#> Finished MPLE.
#> Evaluating log-likelihood at the estimate. 
#> Starting maximum pseudolikelihood estimation (MPLE):
#> Obtaining the responsible dyads.
#> Evaluating the predictor and response matrix.
#> Maximizing the pseudolikelihood.
#> Finished MPLE.
#> Evaluating log-likelihood at the estimate. 
#> Starting maximum pseudolikelihood estimation (MPLE):
#> Obtaining the responsible dyads.
#> Evaluating the predictor and response matrix.
#> Maximizing the pseudolikelihood.
#> Finished MPLE.
#> Evaluating log-likelihood at the estimate.
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
© www.soinside.com 2019 - 2024. All rights reserved.