Calculate sample variances from linear regression model for meta analysis?

I want to perform a meta analysis in two levels using linear regression in R (with lm() ). I understand that in the second level I should perform a weighted regression with the weights being the inverse of the sample variances in the first level.

How do you calculate these inverse sample variances from the output of the first level regression model (lm() model). Can you use the variances from the variance-covariance matrix as shown by vcov()? Why (not)? If not, how to calculate the sample variances?

To illustrate my problem I provide a simplified toy problem below.

In this problem I want to measure the effect of a study method used in different schools by different teachers. Teachers can teach in multiple schools. In the first level I summarize the effect for each school (eliminate teacher effect) and in the second level I want to measure the effect of the method.

I'm aware you can solve this with mixed models but I want to know how to solve it through two level regression using lm() using inverse variance weights.

EDIT Note that the actual data I'm working with is rather technical and very large, hence the toy example. The existing method of analysis on my data involves a mixed effect model but this is very slow given the amount of data. There is a solution proposed to speed this up by splitting the mixed effect model in a two level fixed effect model (using 2lm()) This approach is less performant but still satisfactory and much faster. I'm looking to increase the performance of the two level analysis by including the sample variances from level 1 lm() as weights in level 2 lm().

## generate some data
set.seed(101)
dat = data.frame(
method = factor(c('a','a','a','a','a' ,'b','b','b','b','b'))
, school = factor(c('one', 'one', 'two', 'two', 'two', 'three', 'three', 'three',  'four','four'))
, teacher = factor(c('Joe', 'Jane', 'Joe', 'Jane', 'Jef', 'Joe', 'Jane', 'Jef', 'Joe','Jef'))
)
dat$$outcome = rnorm(nrow(dat), mean = as.numeric(dat$$method) + as.numeric(dat$teacher)) dat #> method school teacher outcome #> 1 a one Joe 3.673964 #> 2 a one Jane 2.552462 #> 3 a two Joe 3.325056 #> 4 a two Jane 2.214359 #> 5 a two Jef 3.310769 #> 6 b three Joe 6.173966 #> 7 b three Jane 3.618790 #> 8 b three Jef 3.887266 #> 9 b four Joe 5.917028 #> 10 b four Jef 3.776741 level1 = lm(outcome ~ 0 + school + teacher, dat, contrasts = list(teacher = 'contr.sum')) summary(level1) #> #> Call: #> lm(formula = outcome ~ 0 + school + teacher, data = dat, contrasts = list(teacher = "contr.sum")) #> #> Residuals: #> 1 2 3 4 5 6 7 8 #> -0.30048 0.30048 -0.65925 -0.04748 0.70673 0.57972 -0.25300 -0.32672 #> 9 10 #> 0.38001 -0.38001 #> #> Coefficients: #> Estimate Std. Error t value Pr(>|t|) #> schoolfour 4.5028 0.5183 8.687 0.000967 *** #> schoolone 2.9402 0.5183 5.672 0.004765 ** #> schoolthree 4.5600 0.3992 11.422 0.000335 *** #> schooltwo 2.9501 0.3992 7.390 0.001788 ** #> teacher1 -0.6882 0.3441 -2.000 0.116099 #> teacher2 -0.3460 0.3441 -1.006 0.371493 #> --- #> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 #> #> Residual standard error: 0.6915 on 4 degrees of freedom #> Multiple R-squared: 0.9882, Adjusted R-squared: 0.9706 #> F-statistic: 55.99 on 6 and 4 DF, p-value: 0.0008178 ## data is aggregated per school, removing the teacher effect dat2 = subset(dat, teacher == 'Joe') dat2$outcome = predict(level1,dat2)
level2 =lm(outcome ~ method, dat2)
summary(level2)
#>
#> Call:
#> lm(formula = outcome ~ method, data = dat2)
#>
#> Residuals:
#>        1        3        6        9
#> -0.00493  0.00493  0.02862 -0.02862
#>
#> Coefficients:
#>             Estimate Std. Error t value Pr(>|t|)
#> (Intercept)  3.97937    0.02053  193.81 2.66e-05 ***
#> methodb      1.58626    0.02904   54.63 0.000335 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> Residual standard error: 0.02904 on 2 degrees of freedom
#> Multiple R-squared:  0.9993, Adjusted R-squared:  0.999
#> F-statistic:  2984 on 1 and 2 DF,  p-value: 0.0003349

## use 1 / estimated variances from parameters in level1 model as sample variance
dat2$$w = 1/diag(vcov(level1))[paste0('school',dat2$$school)]
dat2
#>   method school teacher  outcome        w
#> 1      a    one     Joe 3.974444 3.722073
#> 3      a    two     Joe 3.984304 6.274351
#> 6      b  three     Joe 5.594250 6.274351
#> 9      b   four     Joe 5.537017 3.722073
level2_w =lm(outcome ~ method, dat2, weights = w)
summary(level2_w)
#>
#> Call:
#> lm(formula = outcome ~ method, data = dat2, weights = w)
#>
#> Weighted Residuals:
#>         1         3         6         9
#> -0.011940  0.009196  0.053379 -0.069304
#>
#> Coefficients:
#>             Estimate Std. Error t value Pr(>|t|)
#> (Intercept)  3.98063    0.01985  200.51 2.49e-05 ***
#> methodb      1.59231    0.02808   56.72 0.000311 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> Residual standard error: 0.06277 on 2 degrees of freedom
#> Multiple R-squared:  0.9994, Adjusted R-squared:  0.9991
#> F-statistic:  3217 on 1 and 2 DF,  p-value: 0.0003107

• Perhaps metafor-project.org/doku.php/tips:rma_vs_lm_lme_lmer will help you. – mdewey Dec 11 '18 at 9:25
• @mdewey thanks for helping. In the link you provided. they calculate sample variance with escalc from the metafor package. I'm not sure how it's don inside the function. I really want to know, if possible, how to calculate it from a lm() model. – statastic Dec 11 '18 at 12:51