Cross Validated is a question and answer site for people interested in statistics, machine learning, data analysis, data mining, and data visualization. It only takes a minute to sign up.
Sign up to join this community
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?
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.
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.
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.
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.
