5
$\begingroup$

I'm trying to do a meta-analysis for the first time, comparing measurements of a simple experiment treatment against a control in a variety of species. I started by fitting a mixed-effects model to mean values collected from a set of published studies. This works fine, but it of course ignores the fact that the means are estimated with different levels of precision. It sounds like this is a problem that should be solved by weighting the means appropriately. Several of these studies published standard errors or their raw data, so I can use SEs for weighting.

My question is: What is special about 1/SE (or 1/ σ2)? Why not use 2/SE, or 10/SE, or any other such weighting factor?

I was stimulated to ask this because the SEs differ by almost two orders of magnitude due to a few of the studies using extremely precise instruments. Greater precision is obviously better, but it's not obvious to me that these precise estimates are worth as much as 100 less precisely estimates. (Weighting also introduced some algorithmic problems when fitting the model, but that is probably solvable and I'm more interested in the conceptual justification).

$\endgroup$
8
  • 3
    $\begingroup$ i feel like i dont understand your question, the weighting is assuming a diagonal (non constant) covariance online.stat.psu.edu/stat501/lesson/13/13.1 $\endgroup$
    – seanv507
    Commented May 13 at 16:04
  • 6
    $\begingroup$ Note that weights are appear in both the numerator and denominator of the WLS formulae, so using $2/\sigma_i^2$ is the same as using $1/\sigma_i^2$. Also, sometimes Wikipedia is your friend: en.wikipedia.org/wiki/Weighted_least_squares $\endgroup$
    – jbowman
    Commented May 13 at 16:12
  • 1
    $\begingroup$ Not what this question is about but if these are studies on a variety of species, don't you have to make rather strong assumptions for this meta analysis? Not only may different species react differently to this treatment but the differences between species themselves might not be "random". $\endgroup$
    – dipetkov
    Commented May 13 at 21:22
  • 1
    $\begingroup$ @dipetkov Yes, this is true of most meta-analyses in my field. Species differences are entirely confounded with paper/methodological differences. In this instance, though, it seems reasonable to use this simple approach based on what we understand about the phenomenon. $\endgroup$
    – mkt
    Commented May 14 at 11:26
  • 1
    $\begingroup$ You are referring to inverse variance weighting en.m.wikipedia.org/wiki/Inverse-variance_weighting (see derivation there) which is useful in sensor data fusion (for an estimate of a mean value with least variance, if we have independent samples) $\endgroup$
    – Ggjj11
    Commented May 14 at 18:04

2 Answers 2

6
$\begingroup$

You are asking about what is called "weighted least squares". The idea is that if for some $X_i$ values (values of the independent variable) your observations show more variation in the dependent variable $Y_i$ than is the case for other $X_i$ values, those with more variation would also have larger residuals around the regression line and hence influence the sum of squared residuals more than the observations with smaller variation around the line do. You clearly would in general not want this to be the case.

If you know the std. deviation $\sigma_i$ of your $Y$ data for each particular $X_i$ value, you can transform your $Y_i$ and $X_i$ data by dividing them through the known std. deviation $\sigma_i$. If you then regress $\frac{Y_i}{\sigma_i}$ on $\frac{X_i}{\sigma_i}$ the residuals in this new regression do have the same variance for each $X_i$ value.

Mostly, the true $\sigma_i$ are unknown, and instead the sample estimates $s_i$ are used. In the situation you describe you know the std. errors of the means in all the studies, and thus you can divide through these std. errors.

EDIT

There is a tricky detail I forgot to mention! Also, the intercept has to be divided by $s_i$ and finally, the residual $r_i$ term also. I will demonstrate this in the R script below too. So, the regression equation changes from

$y_i = b_0 + b_1x_i + r_i$

into

$\frac{y_i}{s_i} = b_0\frac{1}{s_i} + b_1\frac{x_i} {s_i} + \frac{r_i}{s_i}$

In the "transformed" regression equation, the sum of squares of the (transformed) residuals which is minimized is $\sum\frac{r_i^2}{\ s_i^2}$. Due to the squares, it is said that the "regression weights" are equal to $\frac{1}{s_i^2}$.

If you would use for the weights $\frac{9}{s_i^2}$ instead of simply $\frac{1}{s_i^2}$, then the sum of squares of the residuals would be nine times as high, and hence the residual standard error would be $\sqrt9=3$ times as high, as you can verify in the script below, comparing the results of model2 and model3. But for the regression coefficients and their std. error and p-values nothing changes!

# Generate some data: 100 cases, 10 for each x value 1 through 10 .
set.seed(123)
x <- rep(1:10, each=10)
# Generate random error, dependent on the x value.
error <- rnorm(100,0,1) * x
y <- 1+0.5*x + error

# Notice the heteroscedasticity over the x values.
plot(x,y)

# Calculate sd per x value.
sd_i <- tapply(y, x, sd)

# Repeat each sd value 10 times, as each x value appears 10 times.
sd_i <- rep(sd_i, each=10)

# Calculate the weights variable for the linear model to be aplied later.
theweights <- 1/sd_i^2

# Now instead of running a linear model on y and x using the regression weights, 
# we now (for illustrative purposes) transform the data as described and 
# run a linear model on these transformed data, NOT using the weights. 
y_w <- y / sd_i
x_w <- x / sd_i
newintercept <- 1 /sd_i 

# Run a linear model on these transformed data.
model1 <- lm(y_w ~ 0 + newintercept + x_w)
summary(model1)

# Now, run a linear model using regression weights: results are equal to those
# of model1.
model2 <- lm(y ~ x, weights=theweights)
summary(model2)

# Use another arbitrary numerator in the weights.
theweights <- 9/sd_i^2
model3 <- lm(y ~ x, weights=theweights)
summary(model3)

The results of this script are:

Scatterplot of y against x: larger variances for larger x values

> model1 <- lm(y_w ~ 0 + newintercept + x_w)
> summary(model1)

Call:
lm(formula = y_w ~ 0 + newintercept + x_w)

Residuals:
    Min      1Q  Median      3Q     Max 
-2.5047 -0.7491 -0.0990  0.7445  1.9961 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
newintercept   0.8120     0.3541   2.293    0.024 *  
x_w            0.7078     0.1228   5.765 9.48e-08 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.006 on 98 degrees of freedom
Multiple R-squared:  0.5763,    Adjusted R-squared:  0.5677 
F-statistic: 66.65 on 2 and 98 DF,  p-value: < 2.2e-16

> 
> # Now, run a linear model using regression weights: results are equal to those
> # of model1.
> model2 <- lm(y ~ x, weights=theweights)
> summary(model2)

Call:
lm(formula = y ~ x, weights = theweights)

Weighted Residuals:
    Min      1Q  Median      3Q     Max 
-2.5047 -0.7491 -0.0990  0.7445  1.9961 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   0.8120     0.3541   2.293    0.024 *  
x             0.7078     0.1228   5.765 9.48e-08 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.006 on 98 degrees of freedom
Multiple R-squared:  0.2533,    Adjusted R-squared:  0.2456 
F-statistic: 33.24 on 1 and 98 DF,  p-value: 9.475e-08

> 
> # Use another arbitrary numerator in the weights.
> theweights <- 9/sd_i^2
> model3 <- lm(y ~ x, weights=theweights)
> summary(model3)

Call:
lm(formula = y ~ x, weights = theweights)

Weighted Residuals:
   Min     1Q Median     3Q    Max 
-7.514 -2.247 -0.297  2.234  5.988 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   0.8120     0.3541   2.293    0.024 *  
x             0.7078     0.1228   5.765 9.48e-08 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 3.019 on 98 degrees of freedom
Multiple R-squared:  0.2533,    Adjusted R-squared:  0.2456 
F-statistic: 33.24 on 1 and 98 DF,  p-value: 9.475e-08
$\endgroup$
3
  • $\begingroup$ Thanks, BenP. I appreciate the overview, but my question is focused on why we use the specific weights that we do and not some variant. Your answer partly deals with this, though you appear to suggest weighting by standard deviation, while weighting by variance is most common. If you could clarify this point and expand on your answer to address the bolded part of the question (jbowman's comment is a partial answer), I'd be happy to accept. $\endgroup$
    – mkt
    Commented May 14 at 11:32
  • $\begingroup$ @mkt I will add a few things to my answer, it was too fast, sorry. $\endgroup$
    – BenP
    Commented May 14 at 14:11
  • $\begingroup$ Many thanks! And in response to your last comment - no need to apologise at all. I was merely requesting an extension; I appreciate both the original answer and the helpful update. $\endgroup$
    – mkt
    Commented May 14 at 17:24
1
$\begingroup$

To understand this precisely, we can't use SE and $\sigma^2$ interchangeably. One is the property of the analysis, the other is a property of the data and both involve assumptions about the data generating process; all the more so when the data are heteroscedastic. Consider the linear model:

$$ Y_i = \alpha + \beta X_i + \epsilon_i$$

In the usual regression model we say $\epsilon \sim $ iid $\mathcal{N}(0, \sigma^2)$ iid meaning independent, identically distributed.

For simplicity, let's retain that they are independent, but suppose they are heteroscedastic meaning that they are NOT identically distributed. Therefore there is no $\sigma^2$ to speak of, except perhaps as a weighted average of each of the $\sigma^2_i$). We don't need to know anything about the shape or form of heteroscedasticity to arrive at a general result as you will see.

There are several questions that might arise. A classical statistician might ask: what's the best estimator of model parameters $\alpha$ and $\beta$? By "best" I mean it's unbiased and achieves the Cramer Rao lower bound. To simplify, I might restrict to linear estimators of the form $\hat{Y} = cY$.

The Gauss Markov theorem addresses this problem. (Note Gauss and Markov were lived a century or two apart), the Gauss Markov theorem states that inverse variance weighted least squares model is BLUE that is the best, linear unbiased estimator of $\alpha$ and $\beta$. An immediate corollary affirms OLS as a best estimator for independent data.

That is why you frequently see inverse variance weighting as an essential component of analyzing non-stationary observations. The reference "Linear Regression Modeling" by Seber and Lee goes through the proof and details.

$\endgroup$
1
  • $\begingroup$ Thank you, much appreciated! $\endgroup$
    – mkt
    Commented May 15 at 11:47

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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