Same SE for all coefficients of a linear model

Question

I am trying to find shellfish densities for a field-based aquaculture program and wanted to understand the impact of reef structure and site on the response. This experiment only has a sample size of 36. I looked at the distribution and log-transformed it for normality before running a linear model:

model<-lm(data$logged ~ data$reef + data$site)

My output looks like this:

             Estimate Std. Error t value Pr(>|t|)    
(Intercept)  10.2593     0.1599  64.179  < 2e-16 ***
data$reef2   -1.2147     0.1958  -6.204 1.06e-06 ***
data$reef3   -1.4936     0.1958  -7.629 2.61e-08 ***
data$reef4   -2.6405     0.1958 -13.487 9.03e-14 ***
data$reef5   -1.0867     0.1958  -5.550 6.18e-06 ***
data$reef6   -0.8922     0.1958  -4.557 9.30e-05 ***
data$siteB   -0.6314     0.1384  -4.561 9.21e-05 ***
data$siteC   0.2567     0.1384   1.854   0.0743 .  

Residual standard error: 0.3391 on 28 degrees of freedom
Multiple R-squared:  0.8943,    Adjusted R-squared:  0.8679 
F-statistic: 33.85 on 7 and 28 DF,  p-value: 4.937e-12

What I don't understand is why my SE value are the same for all of my coefficients. Can someone help me understand how to fix this issue?

Thanks.

Answer 1

Mathematically, if writing the linear model as $y = X\beta + \epsilon$, where $X = \begin{bmatrix}e & x_1 & x_2 & \cdots & x_p\end{bmatrix}$, $\beta = \begin{bmatrix}\beta_0 & \beta_1 & \beta_2 & \cdots & \beta_p\end{bmatrix}'$. The standard error of $\hat{\beta}_j, 1 \leq j \leq p$, denoted by $\hat{\sigma}_{\hat{\beta}_j}$, is then given by \begin{align} \hat{\sigma}\sqrt{e_j'(X'X)^{-1}e_j}, \tag{1} \end{align} where $e_j$ is a $(p+1)$-long column vector whose $(j + 1)$-st entry is $1$ and all the other entries $0$, $\hat{\sigma} = (n - p - 1)^{-1}y'(I - H)y$, $H = X(X'X)^{-1}X'$. In other words, $\hat{\sigma}_{\hat{\beta}_j}$ is the square root of the $(j + 1, j + 1)$ diagonal entry of the matrix $\hat{\sigma}^2(X'X)^{-1}$. Therefore, if all the diagonal entries of $(X'X)^{-1}$ are the same (or block-wise same), then you would see identical (or block-wise identical) standard errors of OLS estimates.

Although you did not include details of the input data in your post, the output implies that your data probably came from a designed experiment, where the columns of $X$ are mutually orthogonal. In this case, the diagonal entries of $(X'X)^{-1}$ could be identical. For example, suppose \begin{align} X = \begin{bmatrix} 1_9 & 1_9 & 0 & 0 \\ 1_9 & 0 & 1_9 & 0 \\ 1_9 & 0 & 0 & 1_9 \\ 1_9 & 0 & 0 & 0 \end{bmatrix}, \end{align} where $1_9$ is a $9$-long column vector consisting of all ones. It is then easy to verify that \begin{align} X'X = \begin{bmatrix} 36 & 9 & 9 & 9 \\ 9 & 9 & 0 & 0 \\ 9 & 0 & 9 & 0 \\ 9 & 0 & 0 & 9 \end{bmatrix}, \quad (X'X)^{-1} = \frac{1}{9}\begin{bmatrix} 1 & -1 & -1 & -1 \\ -1 & 2 & 1 & 1 \\ -1 & 1 & 2 & 1 \\ -1 & 1 & 1 & 2 \end{bmatrix}. \end{align} It can be seen that all the $(2, 2), (3, 3), (4, 4)$ diagonal entries of $(X'X)^{-1}$ equal to $\frac{2}{9}$, yielding the same standard errors of $\hat{\beta}_1, \hat{\beta}_2$ and $\hat{\beta}_3$.

To analyze your particular case, you can use model.matrix command to print out your design matrix $X$ and compute $(X'X)^{-1}$ to verify the reason.

Answer 2

To accompany Doctor Milt's answer (+1) [and Zhanxiong's answer, also +1, which was written while I was posting mine], here is a simulation showing the same thing happening with random data of the same size and structure as yours. You can play around and see whether the std errors are equal as you change the sample size or add/remove individual observations:

df <- expand.grid(reef=factor(rep(1:6)), site=factor(c("A", "B", "C")))
df <- rbind(df, df)
df$outcome <- rnorm(n=nrow(df))
summary(lm(outcome ~ reef + site, data=df))

And the output is:

> dim(df)
[1] 36  3
> summary(lm(outcome ~ reef + site, data=df))

Call:
lm(formula = outcome ~ reef + site, data = df)

Residuals:
     Min       1Q   Median       3Q      Max 
-1.92710 -0.61335  0.04716  0.48348  1.44549 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)  
(Intercept)  0.01300    0.41824   0.031   0.9754  
reef2       -0.83613    0.51223  -1.632   0.1138  
reef3       -0.74835    0.51223  -1.461   0.1552  
reef4       -0.30631    0.51223  -0.598   0.5547  
reef5       -0.19428    0.51223  -0.379   0.7073  
reef6       -0.08205    0.51223  -0.160   0.8739  
siteB        0.66652    0.36220   1.840   0.0764 .
siteC        0.22147    0.36220   0.611   0.5458  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.8872 on 28 degrees of freedom
Multiple R-squared:  0.2265,    Adjusted R-squared:  0.03312 
F-statistic: 1.171 on 7 and 28 DF,  p-value: 0.3502

You can try things like summary(lm(outcome ~ reef + site, data=df[-5, ])) to see the impact of removing single rows of data.

Answer 3

It doesn't look like there's an issue.

Your standard errors are equal because you presumably have a balanced data set, i.e. two observations for each of the 18 possible combinations of reef and site.

Answer 4

The (usual) equation for the variance of a regression parameter (not the intercept) $\hat\beta_j$ is as follows.

$$ \widehat{\text{var}}\left( \hat\beta_j \right) = \dfrac{ s^2 }{ (n-1)s^2_{X_j} }\times\dfrac{1}{ 1-R^2_j } $$

$s^2$ is the residual variance.

$n$ is the sample size

$s^2_{X_j}$ is the variance of the feature $j$

$R^2_j$ is the variance-inflation factor (VIF) for feature $j$, which is calculated by running a regression of feature $j$ on the other features and calculating the $R^2$ of that regression.

To get your observed result, one of the following must have occurred.

1. The change in $s^2_{X_j}$ perfectly offset the change in the VIF for every variable

2. All $s^2_{X_j}$ and all $R^2_j$ are equal.

(Perhaps there could be a combination of the two where some changes in $R^2_j$ are offset by changes in $s^2_{X_j}$ and the rest of the features have the same $s^2_{X_j}$ and $R^2_j$.)

Possibility #1 requires a lot of coincidences. On the other hand, for a designed experiment where all of the features have equal variance and are independent of each other, possibility #2 will occur.

And since not all of your slope coefficients are equal, apply the above logic to the groups of coefficients that do have equal coefficients.

The equation for $ \widehat{\text{var}}\left( \hat\beta_j \right) $ is available on the VIF Wikipedia article, which has references to other material.

Answer 5

Thanks to the answer by Zhanxiong, despite orthogonal variables, there is another possible situation that makes the $\hat\sigma_{\hat\beta_j}$ the same. For example, variables j=1 and 2, when $\sum_{i=1}^{n}x_{ij}x_{ij}' = \sum_{i=1}^{n}x_{ij}^2$ (for elements on the diagonal) are the same for j=1 and 2.

This will happen if we merge the data in some ways. Taking the example in migration matrix, for each ordered pair (p1,p2), we have 2 variables, number of migrants from p1 to p2, and number of migrants from p2 to p1. For each row (obs as one place), these two variables are different; and these two variables are not orthogonal. But for $\sum_{i=1}^{n}x_{ij}^2$, they are the same, since each "number of migrants from p1 to p2" can be find as "number of migrants from p2 to p1" for the pair (p2,p1).