4
$\begingroup$

I have a question very similar to the question asked here: is it possible to calculate standard errors (specifically, the standard error of the intercept) for generalized least squares regression coefficients if we only have sample means and sample covariance?

For example, suppose we have a sample covariance matrix (not necessarily 2 by 2):

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

and we also know sample means of each variable, as well as sample size. From this, the slopes, the intercept, and standard errors of the slope can be calculated as done in the R code example below (adapted from the answer given on Standard error of regression coefficient without raw data). It works exactly as it does for OLS regression, EXCEPT for the intercept standard error. This is because for GLS, standard errors are calculated as $\sigma^2(X'\Sigma X)^{-1}$ rather than $\sigma^2(X'X)^{-1}$ as in OLS.

So, my questions are:

  1. If we do not have access to the raw data $X$ and $Y$ and we can't access $\Sigma$, but instead only have access to means and sample covariance (and sample size), is it possible to correctly calculate intercept standard error in the context of GLS?
  2. Along similar lines, is it possible to fit a GLS regression through the origin (i.e., no intercept), and calculate standard errors for the slopes using only the sample covariance matrix, means, and sample size?
  3. If the above situations are not possible, would it be possible if we also have access to $\Sigma$?

Fully reproducible example in R:

require(phytools)
require(nlme)

n <- 50 # number of individuals
p <- 2 # number of predictors

# generate data and covariance structure
tree <- pbtree(n=n)
data <- sim.corrs(tree,vcv = matrix(c(1,.8,.8,.8,1,.8,.8,.8,1),3,3))
colnames(data) <- c("V1","V2","V3")
Sigma <- vcv(tree)

### sample means
m <- apply(data,2,function(X) ace(X,tree,method="pic")$ace[1])
m_mat <- matrix(1,nrow(data),1) %*% m

### sample covariance matrix
v <- t(data-m_mat) %*% solve(Sigma) %*% (data-m_mat) / (n-1)

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

results <- summary(gls(V3~V1+V2,data=as.data.frame(data),correlation=corBrownian(1,tree)))[[18]]
se <- sqrt(diag(solve(xx) * c(s2)))

### Notice that the se_calc intercept is wrong
print(data.frame(gls_beta=results[,1],beta_calc=b,gls_se=results[,2],se_calc=se))

### This, on the other hand, gives the correct intercept standard error
### But what if we only have means and sample covariance?
sqrt(diag(solve(t(cbind(1,data[,1:2])) %*% solve(Sigma) %*% as.matrix(cbind(1,data[,1:2])))*c(s2)))
$\endgroup$
  • $\begingroup$ I'm afraid the answers will be in the negative--certainly so when $p\ge 2$. What would work would be to have the coefficients of $X^\prime \Sigma X$ (or its inverse) instead of the covariances of the data, because this plays the analogous role in GLS. $\endgroup$ – whuber Feb 6 '15 at 21:06

Your Answer

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

Browse other questions tagged or ask your own question.