0
$\begingroup$

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
$\endgroup$
  • $\begingroup$ Perhaps metafor-project.org/doku.php/tips:rma_vs_lm_lme_lmer will help you. $\endgroup$ – mdewey Dec 11 '18 at 9:25
  • $\begingroup$ @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. $\endgroup$ – statastic Dec 11 '18 at 12:51
0
$\begingroup$

Under different circumstances, I would post this as a comment and not an answer...but as this is an important issue, I'm going to opt for the more "permanent" approach.

The answer is that you can't use lm() for the 2-level analysis (with or without the variance-adjusted weights)...or, perhaps it is better to say that you shouldn't want to use lm().

Briefly, if you include the 2nd-level variables in an lm model (with or without weights, but in the scenario I'm about to describe, you should actually want to use the weights), you can obtain a separate measure for the effect of each of the distinct 2-level units. However, you entirely lose the ability to generalize to the larger population. And, as this is the implied goal with any meta-analysis, you really shouldn't want to do this. Said another way, if you use a fixed effects approach (lm()), then you can state something about how these particular 2nd-level units are/aren't distinct from each other. However, you cannot generalize anything to the larger population. It is the mixed-effects approach that allows you to make this generalization.

Happy to provide clarification as needed.

$\endgroup$
  • $\begingroup$ Thanks for your answer! I'm not interested in generalisation to a population but in inference only on the experimental units. I really want to try to calculate these sample variances in a two level fixed effect model. I added a paragraph to the question explaining why I want to do that :) $\endgroup$ – statastic Dec 14 '18 at 9:49

Your Answer

By clicking "Post Your Answer", you acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.

Not the answer you're looking for? Browse other questions tagged or ask your own question.