Failing to get correct results for GLS regression

Question

I'm trying to run the most basic GLS regression but I'm not getting the expected results. My GLS uses the identity matrix as the error covariance matrix.

In R, the GLS regression is:

set.seed(0)
n = 50
x = as.matrix(seq(1, n))
X = as.matrix(cbind(rep(1, n), x))
y = 2 + 3 * x + rnorm(n)

V = diag(n)

m = glm(y~1+x, weights=diag(V))
s = summary(m)

So far so good. Now, I do the above by hand instead of using the glm package. From wikipedia:

$$ \hat{\beta} = (X'V^{-1}X)^{-1} X' V^{-1} y \\ Cov(\hat{\beta}) = (X'V^{-1}X)^{-1} $$

Then define $\beta$'s covariance matrix as $C = Cov(\hat{\beta}) = (X'V^{-1}X)^{-1}$. Then the diagonal elements $Var(\hat{\beta}_i) = c_{ii}$ are the variances. This allows me to calculate the standard errors:

$$ se(\hat{\beta}_i) = \sqrt{Var(\hat{\beta}_i)} = \sqrt{c_{ii}} \\ \text{tvalue} = \frac{\hat{\beta}_i}{se(\hat{\beta}_i)} $$

b = solve(t(X) %*% solve(V) %*% X) %*% t(X) %*% solve(V) %*% y

var_b = solve(t(X) %*% solve(V) %*% X)
st_err = sqrt(diag(var_b))
tvalues = b/st_err

I get the coefficients correctly, but my standard errors and tvalues are different. I can also successfully calculate $\hat{\sigma}^2$ so I'm clearly overlooking something about $Cov(\hat{\beta})$.

Answer 1

By default, glm for a gaussian family will use $$\widehat{\mathrm{cov}}[\hat\beta]=\hat\sigma^2\left(X^TV^{-1}X\right)^{-1}$$

To use a known dispersion parameter instead of estimating, specify it as an option to summary

> summary(m, dispersion=1)

Call:
glm(formula = y ~ 1 + x, weights = diag(V))

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)   2.1409     0.2871   7.456 8.91e-14 ***
x             2.9954     0.0098 305.657  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for gaussian family taken to be 1)

which agrees with your calculations