Same SE for all coefficients of a linear model 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.

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}
e_9 & e_9   & 0   & 0 \\
e_9 & 0     & e_9 & 0 \\
e_9 & 0     & 0   & e_9 \\
e_9 & 0     & 0   & 0 
\end{bmatrix},
\end{align}
where $e_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 equal to $2$, 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.
