Some information about this (fake) dataset:
study_id = unique identifier per study
link_id = unique identifier per effect size
g = hedge’s g for age differences in racial bias (more positive = older kids exhibited more racial bias)
v_g = variance of hedge’s g
v_se = standard error of hedge’s g
N = sample size
study_quality = the coded study quality score (range from 1-10; higher scores = higher study quality)
published_ec = the effect coded indicator for whether the study was published (coded 0.5) or unpublished (coded -0.5)



Setting up the nesting variables as factors
metafor also likes the effect size id (or link id) to not restart at 1 per study, thus the following:

dat %>%
  mutate(esid = as.factor(row_number()),
         samid = as.factor(study_id)) -> dat

Analyses begin here:

summary(mod<-rma.mv (yi = g, V= v_g,  random = list (~1|esid, ~1|samid), tdist=TRUE, data = dat, slab = link_id))
## 
## Multivariate Meta-Analysis Model (k = 112; method: REML)
## 
##   logLik  Deviance       AIC       BIC      AICc   
## -71.5553  143.1106  149.1106  157.2392  149.3349   
## 
## Variance Components:
## 
##             estim    sqrt  nlvls  fixed  factor 
## sigma^2.1  0.0935  0.3058    112     no    esid 
## sigma^2.2  0.0971  0.3116     75     no   samid 
## 
## Test for Heterogeneity:
## Q(df = 111) = 741.2919, p-val < .0001
## 
## Model Results:
## 
## estimate      se     tval   df    pval   ci.lb   ci.ub      
##   0.7972  0.0504  15.8312  111  <.0001  0.6974  0.8970  *** 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
confint(mod)
## 
##           estimate  ci.lb  ci.ub 
## sigma^2.1   0.0935 0.0479 0.1757 
## sigma.1     0.3058 0.2188 0.4192 
## 
##           estimate  ci.lb  ci.ub 
## sigma^2.2   0.0971 0.0239 0.1831 
## sigma.2     0.3116 0.1545 0.4279
n <- length(dat$v_g)
list.inverse.variances <- 1 / (dat$v_g)
sum.inverse.variances <- sum(list.inverse.variances)
squared.sum.inverse.variances <- (sum.inverse.variances) ^ 2
list.inverse.variances.square <- 1 / (dat$v_g^2)
sum.inverse.variances.square <- sum(list.inverse.variances.square)
numerator <- (n - 1) * sum.inverse.variances
denominator <- squared.sum.inverse.variances -  sum.inverse.variances.square
estimated.sampling.variance <- numerator / denominator

100*(estimated.sampling.variance) / (mod$sigma2[1] + mod$sigma2[2] + estimated.sampling.variance)
## [1] 15.17067
100*(mod$sigma2[1]) / (mod$sigma2[1] + mod$sigma2[2] + estimated.sampling.variance)
## [1] 41.6055
100*(mod$sigma2[2]) / (mod$sigma2[1] + mod$sigma2[2] + estimated.sampling.variance)
## [1] 43.22383
funnel(mod)

rstudent(mod, progbar=TRUE) -> rstud.es
rstud = data.frame (resid_es = rstud.es$resid,
                    se_es = rstud.es$se,
                    z_es = rstud.es$z,
                    link_id = rstud.es$slab)

# re-run model without influential cases
rstud %>% filter (abs(z_es) >= 1.96) -> influ.df
nrow(influ.df)
## [1] 6
length(influ.df$link_id)/nrow(rstud)
## [1] 0.05357143
dat %>% filter (., !link_id %in% influ.df$link_id) %>%
  mutate(esid = as.factor(row_number()),
         samid = as.factor(study_id)) -> dat_trim

summary(mod<-rma.mv (yi = g, V= v_g,  random = list (~1|esid, ~1|samid), tdist=TRUE, data = dat_trim, slab = link_id))
## 
## Multivariate Meta-Analysis Model (k = 106; method: REML)
## 
##   logLik  Deviance       AIC       BIC      AICc   
## -53.0400  106.0800  112.0800  120.0419  112.3176   
## 
## Variance Components:
## 
##             estim    sqrt  nlvls  fixed  factor 
## sigma^2.1  0.0804  0.2836    106     no    esid 
## sigma^2.2  0.0482  0.2195     72     no   samid 
## 
## Test for Heterogeneity:
## Q(df = 105) = 519.7661, p-val < .0001
## 
## Model Results:
## 
## estimate      se     tval   df    pval   ci.lb   ci.ub      
##   0.7820  0.0428  18.2760  105  <.0001  0.6971  0.8668  *** 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(mod<-rma.mv (yi = g, V= v_g, mods = ~1 + scale(study_quality, scale = F), random = list (~1|esid, ~1|samid), 
                     tdist=TRUE, data = dat, slab = link_id))
## 
## Multivariate Meta-Analysis Model (k = 112; method: REML)
## 
##   logLik  Deviance       AIC       BIC      AICc   
## -59.3302  118.6605  126.6605  137.4624  127.0414   
## 
## Variance Components:
## 
##             estim    sqrt  nlvls  fixed  factor 
## sigma^2.1  0.0890  0.2983    112     no    esid 
## sigma^2.2  0.0488  0.2210     75     no   samid 
## 
## Test for Residual Heterogeneity:
## QE(df = 110) = 534.7440, p-val < .0001
## 
## Test of Moderators (coefficient 2):
## F(df1 = 1, df2 = 110) = 28.4593, p-val < .0001
## 
## Model Results:
## 
##                                  estimate      se     tval   df    pval   ci.lb 
## intrcpt                            0.7937  0.0427  18.5663  110  <.0001  0.7090 
## scale(study_quality, scale = F)    0.0854  0.0160   5.3347  110  <.0001  0.0537 
##                                   ci.ub      
## intrcpt                          0.8784  *** 
## scale(study_quality, scale = F)  0.1171  *** 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
psych::describe(dat$study_quality)
##    vars   n mean   sd median trimmed  mad  min  max range skew kurtosis   se
## X1    1 112 5.58 2.69   5.26    5.59 3.41 1.01 9.91   8.9 0.01    -1.23 0.25
predict(mod, 2.69)
## 
##    pred     se  ci.lb  ci.ub  pi.lb  pi.ub 
##  1.0234 0.0602 0.9040 1.1428 0.2781 1.7687
predict(mod, -2.69)
## 
##    pred     se  ci.lb  ci.ub   pi.lb  pi.ub 
##  0.5640 0.0611 0.4430 0.6851 -0.1816 1.3096
summary(mod<-rma.mv (yi = g, V= v_g, mods = ~1 + se_g, random = list (~1|esid, ~1|samid), 
                     tdist=TRUE, data = dat, slab = link_id))
## 
## Multivariate Meta-Analysis Model (k = 112; method: REML)
## 
##   logLik  Deviance       AIC       BIC      AICc   
## -71.3197  142.6393  150.6393  161.4413  151.0203   
## 
## Variance Components:
## 
##             estim    sqrt  nlvls  fixed  factor 
## sigma^2.1  0.0963  0.3103    112     no    esid 
## sigma^2.2  0.0954  0.3088     75     no   samid 
## 
## Test for Residual Heterogeneity:
## QE(df = 110) = 737.7902, p-val < .0001
## 
## Test of Moderators (coefficient 2):
## F(df1 = 1, df2 = 110) = 0.0662, p-val = 0.7974
## 
## Model Results:
## 
##          estimate      se     tval   df    pval    ci.lb   ci.ub      
## intrcpt    0.8591  0.2461   3.4912  110  0.0007   0.3715  1.3468  *** 
## se_g      -0.3244  1.2609  -0.2573  110  0.7974  -2.8233  2.1744      
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(mod<-rma.mv (yi = g, V= v_g, mods = ~1 + v_g, random = list (~1|esid, ~1|samid), 
                     tdist=TRUE, data = dat, slab = link_id))
## 
## Multivariate Meta-Analysis Model (k = 112; method: REML)
## 
##   logLik  Deviance       AIC       BIC      AICc   
## -71.3108  142.6217  150.6217  161.4236  151.0026   
## 
## Variance Components:
## 
##             estim    sqrt  nlvls  fixed  factor 
## sigma^2.1  0.0963  0.3103    112     no    esid 
## sigma^2.2  0.0954  0.3088     75     no   samid 
## 
## Test for Residual Heterogeneity:
## QE(df = 110) = 737.3860, p-val < .0001
## 
## Test of Moderators (coefficient 2):
## F(df1 = 1, df2 = 110) = 0.0718, p-val = 0.7892
## 
## Model Results:
## 
##          estimate      se     tval   df    pval    ci.lb   ci.ub      
## intrcpt    0.8285  0.1273   6.5075  110  <.0001   0.5762  1.0808  *** 
## v_g       -0.8300  3.0978  -0.2679  110  0.7892  -6.9691  5.3091      
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(mod<-rma.mv (yi = g, V= v_g, mods = ~1 + published_ec, random = list (~1|esid, ~1|samid), 
                     tdist=TRUE, data = dat, slab = link_id))
## 
## Multivariate Meta-Analysis Model (k = 112; method: REML)
## 
##   logLik  Deviance       AIC       BIC      AICc   
## -67.9724  135.9447  143.9447  154.7467  144.3257   
## 
## Variance Components:
## 
##             estim    sqrt  nlvls  fixed  factor 
## sigma^2.1  0.0921  0.3035    112     no    esid 
## sigma^2.2  0.0860  0.2932     75     no   samid 
## 
## Test for Residual Heterogeneity:
## QE(df = 110) = 700.6512, p-val < .0001
## 
## Test of Moderators (coefficient 2):
## F(df1 = 1, df2 = 110) = 6.8170, p-val = 0.0103
## 
## Model Results:
## 
##               estimate      se     tval   df    pval   ci.lb   ci.ub      
## intrcpt         0.7628  0.0504  15.1247  110  <.0001  0.6628  0.8627  *** 
## published_ec    0.2634  0.1009   2.6109  110  0.0103  0.0635  0.4632    * 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
predict(mod, -0.5)
## 
##    pred     se  ci.lb  ci.ub   pi.lb  pi.ub 
##  0.6311 0.0801 0.4724 0.7898 -0.2201 1.4823
predict(mod, 0.5)
## 
##    pred     se  ci.lb  ci.ub  pi.lb  pi.ub 
##  0.8944 0.0613 0.7729 1.0160 0.0494 1.7395