8
$\begingroup$

After searching here:

Perform simple regression without raw data

I am still curious about this.

Is it possible to derive the standard error of a regression coefficient from summary data alone?

E.g., assume we are given the following variance-covariance matrix.

$\begin{bmatrix}Var(X) & Cov(X,Y)\\Cov(X,Y) & Var(Y)\end{bmatrix}$

We can derive the regression coefficient $\beta_{XY} = Cov(X,Y) / Var(X)$.

Given a specific n, is it possible to derive the standard error of $\beta$ as well? If so, which formula is being used? It appears that all formulas for regression standard errors that I could find assume that you know the variance of residuals of the regression, which we don't know from summary data alone.

---UPDATE--- Based on feedback by Greg and whuber, I am now at a point, where I know that

$\sigma^2(b) = \sigma^2(X'X)^{-1}=\sigma^2(Var(X)(n-1))^{-1}$,

but typically matrix formulas for the standard error of $b$ then note that $\sigma^2$, the unknown variance of the errors, is estimated with MSE.

Based on comments by whuber, sum of squares residuals is:

$Y'Y - bX'X$, which according to Greg can be spelled out as

$Var(Y)(n-1) - b(Var(x)(n-1))$

And to go from here to MSE, I think (?) all that is left is divide by n-2,

so $(Var(Y)(n-1) - b(Var(x)(n-1))) / (n-2)$

And plugging this all in yields:

$\sigma^2(b) =((Var(Y)(n-1) - b(Var(X)(n-1))) / (n-2))\times(Var(X)(n-1))^{-1}$

I will try to see if this works out in R, and will report back. If anyone sees something blatantly wrong, I appreciate the heads up.

$\endgroup$
  • $\begingroup$ You know the sample size, presumably? $\endgroup$ – Jeremy Miles Jan 27 '15 at 18:08
  • $\begingroup$ Yes, presumably the sample size n would be known. $\endgroup$ – Felix Thoemmes Jan 27 '15 at 18:35
  • 1
    $\begingroup$ In response to your last line, the sum of squares of the residuals is equal to $Y^\prime Y - b^\prime X^\prime Y$ where $b$ is the vector of estimated coefficients. These terms can be found in the $\text{Cov}(X,Y)$ and $\text{Var}(Y)$ summaries (combined with the means). $\endgroup$ – whuber Jan 27 '15 at 18:51
10
$\begingroup$

With the usual notation, arrange the independent variables by columns in an $n\times (p+1)$ array $X$, with one of them filled with the constant $1$, and arrange the dependent values $Y$ into an $n$-vector (also a column). Assuming the appropriate inverses exist, the desired formulas are

$$b = (X^\prime X)^{-1} X^\prime Y$$

for the $p+1$ parameter estimates $b$,

$$s^2 = (Y^\prime Y - b^\prime X^\prime Y)/(n-p-1)$$

for the residual standard error, and

$$V = (X^\prime X)^{-1} s^2$$

for the variance-covariance matrix of $b$.

When only the summary statistics $\newcommand{\m}{\mathrm m} \m_X$ (the means of the columns of $X$, forming a $p+1$-covector which includes a mean of $1$ for the constant), $m_Y$ (the mean of $Y$), $\text{Cov}(Y)$, $\text{Var}(Y)$, and $\text{Cov}(X,Y)$ are available, first recover the needed values via

$$(X^\prime X)_0 = (n-1)\text{Cov}(X) + n \m_X \m_X^\prime,$$

$$Y^\prime Y = (n-1)\text{Var}(Y) + n m_Y,$$

$$(X^\prime Y)_0 = (n-1)\text{Cov}(X,Y) + n m_Y \m_X.$$

Then to obtain $X^\prime X$, border $(X^\prime X)_0$ symmetrically by a vector of column sums (given by $n \m_X$) with the value $n$ on the diagonal; and to obtain $X^\prime Y$, augment the vector $(X^\prime Y)_0$ with the sum $n m_Y$. For instance, when using the first column for the constant term, these bordered matrices will look like

$$X^\prime X = \pmatrix{ n & n\m_X \\ n\m_X^\prime & (X^\prime X)_0 }$$

and

$$X^\prime Y = \left(n m_Y, (X^\prime Y)_0\right)^\prime$$

in block-matrix form.

If the means are not available--the question did not indicate they are--then replace them all with zeros. The output will estimate an "intercept" of $0$, of course, and its standard error of the intercept will likely be incorrect, but the remaining coefficient estimates and standard errors will be correct.


Code

The following R code generates data, uses the preceding formulas to compute $b$, $s^2$, and the diagonal of $V$ from only the means and covariances of the data (along with the values of $n$ and $p$ of course), and compares them to standard least-squares output derived from the data. In all examples I have run (including multiple regression with $p\gt 1$) agreement is exact to the default output precision (about seven decimal places).

For simplicity--to avoid doing essentially the same set of operations three times--this code first combines all the summary data into a single matrix v and then extracts $X^\prime X$, $X^\prime Y$, and $Y^\prime Y$ from its entries. The comments note what is happening at each step.

n <- 24
p <- 3
beta <- seq(-p, p, length.out=p)# The model
set.seed(17)
x <- matrix(rnorm(n*p), ncol=p) # Independent variables
y <- x %*% beta + rnorm(n)      # Dependent variable plus error
#
# Compute the first and second order data summaries.
#
m <- rep(0, p+1)                # Default means
m <- colMeans(cbind(x,y))       # If means are available--comment out otherwise
v <- cov(cbind(x,y))            # All variances and covariances
# 
# From this point on, only the summaries `m` and `v` are used for the calculations
# (along with `n` and `p`, of course).
#
m <- m * n                      # Compute column sums
v <- v * (n-1)                  # Recover sums of squares of residuals
v <- v + outer(m, m)/n          # Adjust to obtain the sums of squares
v <- rbind(c(n, m), cbind(m, v))# Border with the sums and the data count
xx <- v[-(p+2), -(p+2)]         # Extract X'X
xy <- v[-(p+2), p+2]            # Extract X'Y
yy <- v[p+2, p+2]               # Extract Y'Y
b <- solve(xx, xy)              # Compute the coefficient estimates
s2 <- (yy - b %*% xy) / (n-p-1) # Compute the residual variance estimate
#
# Compare to `lm`.
#
fit <- summary(lm(y ~ x))
(rbind(Correct=coef(fit)[, "Estimate"], From.summary=b))    # Coeff. estimates
(c(Correct=fit$sigma, From.summary=sqrt(s2)))               # Residual SE
#
# The SE of the intercept will be incorrect unless true means are provided.
#
se <- sqrt(diag(solve(xx) * c(s2))) # Remove `diag` to compute the full var-covar matrix
(rbind(Correct=coef(fit)[, "Std. Error"], From.summary=se)) # Coeff. SEs
| cite | improve this answer | |
$\endgroup$
  • 2
    $\begingroup$ I can confirm that this is indeed correct and working -- and impressive to say the least. What an outstanding and thorough answer. $\endgroup$ – Felix Thoemmes Jan 28 '15 at 15:59
4
$\begingroup$

One quick and dirty approach when you have the summaries but not the raw data is to generate a dataset with the specified summaries and sample size and then run the simulated data through your regular routines to compute whatever you want. The mvrnorm function in the MASS package for R will generate random normal data with a given mean vector and covariance matrix. I am sure other programs will as well (or you can create your own by generating data and then multiplying by the appropriate matrix).

From the theoretical side you can start by thinking about having all your variables centered (mean subtracted so that the current mean is exactly 0, this really only affects the intercept) then $x'x=var(x)\times (n-1)$, $x'y=cov(x,y)\times(n-1)$, and $y'y=var(y)\times(n-1)$. So just plug those values into the matrix version of the formula for the standard error that you want.

| cite | improve this answer | |
$\endgroup$
  • $\begingroup$ Hi Greg, I was aware of the quick and dirty simulation approach that you described, but was also really wanting to get an analytic solution. In my particular problem, I have about 60,000 of these covariance matrices, and rather not simulate all of them. I will double-check your proposed analytic solution, and if it pans out, will report back here, and mark your answer as correct. I appreciate the help so far. $\endgroup$ – Felix Thoemmes Jan 27 '15 at 19:43
  • 1
    $\begingroup$ @felize2000, I have been planning on writing an R function to do this (it would be part of the ObsSens package), but it has been low on the priority list. Maybe I need to revisit this and bump it up on the list. $\endgroup$ – Greg Snow Jan 27 '15 at 21:41
  • $\begingroup$ was not aware of the ObsSens package - just downloaded it. Highly useful for me since I am often doing observational studies in which unobserved confounding cannot be ruled out even after adjustment on observables. $\endgroup$ – Felix Thoemmes Jan 28 '15 at 16:01
  • $\begingroup$ I wasn't aware that mvrnorm can reproduce a "given mean vector and covariance matrix." I thought that it drew independent random values from a multivariate normal distribution, causing its actual output to deviate from the population parameters due to chance. Is there some kind of option that forces it to simulate a desired mean and covariance matrix in the data exactly? $\endgroup$ – whuber Jan 28 '15 at 16:20
  • 1
    $\begingroup$ @whuber, specify empirical=TRUE in the call to mvrnorm and the mean and variance will be as specified (within round off error anyways). $\endgroup$ – Greg Snow Jan 28 '15 at 17:14

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.