6
$\begingroup$

I'm working with a dataset for which I only have means, standard deviations, and sample sizes for different levels of a continuous predictor.

E.G.

Y  X SD_Y N_Y
5  1   3   4
10 2   6   2
15 3   2   8

I would like to determine the regression line that fits this data. I'm wracking my brains to remember how data points should be weighted in a linear regression (I'm also interested in using a generalized linear model as well)- by sample size, variance, SD?

Any pointers?

$\endgroup$
  • $\begingroup$ Does the Y column give the means? $\endgroup$ – Rob Hyndman Nov 2 '10 at 0:12
  • $\begingroup$ Yes. It does give the means. $\endgroup$ – jebyrnes Nov 23 '10 at 0:15
5
$\begingroup$

Weight each mean by the number of points that went into computing it. You can later use the estimated standard deviations to test the hypothesis of homoscedasticity that this approach assumes. If the Ns are as small as in your example then this test probably wouldn't have much power unless the SDs vary greatly.

$\endgroup$
  • $\begingroup$ I'm marking this as correct as it is short and sweet. However, see whuber's lovely longer explication as to WHY it is correct. And, indeed, in fooling around with the problem, weighted by n produces the correct parameter estimates and relatively similar SE and p values to working with simulated data. $\endgroup$ – jebyrnes Nov 3 '10 at 4:48
8
$\begingroup$

This is Analysis of Variance.

After all, consider one of the $y$'s, with standard deviation $s$, and let its predicted value (which depends on the corresponding $x$) be $f$. The original aim is to vary $f$ (within constraints depending on the model; often $f$ is required to be a linear function of $x$) to minimize the sum of squared residuals. Suppose we had the original dataset available. Let the values summarized by a particular $y$ be $y_1, y_2, \ldots, y_k$, so that $y$ is their mean and $s$ is their standard deviation. Their contribution to the sum of squares of residuals equals

$$\sum_{i=1}^k{\left( y_i - f \right)^2} = k \left(y - f \right)^2 + k' s^2 \text{.}$$

(I have written $k'$ because its value depends on how you compute your standard deviations: it is $k$ for one convention and $k-1$ for another.) Because the last term does not depend on $f$, it does not affect the minimization: we can neglect it.

The other term on the right hand side shows that you want to perform a weighted least squares calculation with weights equal to the counts $k$ (the "N_Y" column of the data). Equivalently, you can create a synthetic dataset by making $N_Y$ copies of each datum $(X,Y)$ and performing ordinary least squares regression.

Note that this analysis presumes nothing about the form of the prediction function: it can include any explanatory variables you like and have any form, even nonlinear ones.

Note also that the weighting does not depend on the standard deviations. This is because implicitly we have assumed that the variance of the y's is constant, so that all of the differences among the observed standard deviations are attributed to random fluctuations. This hypothesis can be tested in the usual ways (e.g., with F-tests). For the example data it holds up: those standard deviations do not vary significantly.

Edit I see, in retrospect, that this answer merely reiterates @onestop's pithy response. I'm leaving it up because it demonstrates why @onestop is correct.

$\endgroup$
  • 1
    $\begingroup$ I didn't have time for more than a pithy answer this morning. I was planning to expand on it this evening but you've done the the job very nicely, whuber. $\endgroup$ – onestop Nov 2 '10 at 18:49
  • $\begingroup$ Am I correct in thinking that if the variances did vary significantly, you should weight inversely by the variances? In the case of equal variance this reduces to weighting by the number of observations since the variance of n IID observations is (variance of one observation / n). $\endgroup$ – Chris Taylor Feb 7 '11 at 12:20
  • $\begingroup$ @Chris (This replaces an earlier reply, which mistakenly referenced a different question.) Yes, that's correct. One must always weight in direct proportion to the number of observations, though (assuming the observational errors in each group are independent). I don't think there's an easy or general rule concerning whether or not one should use variance-weighted regression when the variances don't differ by much. (The range of SD's of 6:2 in this case, with only 2 - 8 degrees of freedom, is not wide enough to provide convincing evidence of heteroscedasticity.) $\endgroup$ – whuber Feb 7 '11 at 15:09
  • $\begingroup$ @whuber Dear Whuber, I was hoping maybe you could help me out with this question , you always give very nice answers to methodology-type questions. Hope to hear from you $\endgroup$ – user929304 Apr 15 '16 at 13:26
3
$\begingroup$

Let the disaggregrate model be:

$Y_{ia} = X_a \beta + \epsilon_i$

where

$\epsilon_i \sim N(0,\sigma^2)$

Your aggregate model is given by:

$Y_a = \frac{\sum_i(Y_{ia})}{n_a}$

where,

$n_a$ is the number of observations you have corresponding to the $a$ index.

Therefore, it follows that:

$Y_a = X_a \beta + \epsilon_a$

where

$\epsilon_a \sim N(0, \frac{\sigma^2}{n_a} )$ and

$a=1, 2, ... A$

Therefore, the OLS estimate would be given by minimizing:

$\sum_a(Y_a - X_a \beta)^2$

Which yields the usual solution. So, I do not think there is any difference as far as the estimate for the slope parameters are concerned.

Edit 1

Here is a small simulation in R which illustrates the above idea (apologies for the flaky code as I am using questions such as the above to learn R).

set.seed(1);
n <- c(4,2,8);
x <- c(1,2,3);
data <- matrix(0,14,2)
mean_data <- matrix(0,3,2)
index <- 1;
for (i in 1 : 3)
{
    for(obs in 1:n[i])
    {
        data[index,1] <- x[i];
        data[index,2] <- x[i]*8 + 1.5*rnorm(1);

        mean_data[i,1] = x[i];
        mean_data[i,2] = mean_data[i,2] + data[index,2];

        index = index + 1;  
    }
    mean_data[i,2] = mean_data[i,2] / n[i];
}
beta <- lm(mean_data[,2] ~ mean_data[,1]);

The above code yields the output when you type beta:

Call:
lm(formula = mean_data[, 2] ~ mean_data[, 1])

Coefficients:
   (Intercept)  mean_data[, 1]  
      -0.03455         7.99326  

Edit 2

However, OLS is not efficient as error variances are not equal. Thus, using MLE ideas, we need to minimize:

$\sum_a{n_a (Y_a - X_a \beta)^2}$

In other words, we want to minimize:

$\sum_a{(\sqrt{n_a} Y_a - \sqrt{n_a} X_a \beta)^2}$

Thus, the MLE can be written as follows:

Let $W$ be a diagonal matrix with the $\sqrt{n_a}$ along the diagonal. Thus, the MLE estimate can be written as:

$(X' X)^{-1} X' Y$

where,

$Y = W [Y_1,Y_2,...Y_A]'$ and

$X = W [X_1,X_2,...X_A]'$

Another way to think about this is:

Consider the variance of $Y$. The transformation given above for $Y$ ensures that the variance of the individual values of $Y$ are identical thus satisfying the conditions of the Gauss-Markov theorem that OLS is BLUE.

$\endgroup$
0
$\begingroup$

I think I would calculate normalized variables (z=(x-mean(x))/(sd(x)), and run the regression. Or you could work out a way to generate samples in a bootstrap. I'm not shure if this would be the textbook solution, but intuitively it should work.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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