Why do different regression models give identical parameters?

Question

I am working with data from an experiment where we measured two correlated traits in replicate individuals of the same genotype at four different sites. The two traits are partially regulated by the same pathways (they're methylation of cytosines in different sequence contexts if that helps anyone), so it's not surprising that they are correlated, but we also strongly suspect there maybe correlation with some other unmeasured confounder.

Here is a subset of the data:

# Import data into R.
df <- data.frame(
  genotype = c("6191", "6046", "8241", "8227", "9452", "6191", "9452", "6237", "6104", "6231", "6231", "6180", "6231", "6191", "6046", "6231", "6104", "6046", "8241", "6231", "6126", "6237", "6191", "6126", "6231", "8241", "5856", "9452", "6104", "9470"),
  site = c("B", "C", "C", "A", "C", "C", "D", "B", "B", "C", "B", "D", "B", "B", "D", "A", "B", "A", "D", "D", "A", "A", "B", "D", "C", "B", "A", "C", "C", "D"),
  trait1 = c(-0.4056050, -0.5276847, -0.4070822, -0.5523137, -0.3886427, -0.4417745, -0.4890365, -0.6902847, -0.6573103, -0.5832851, -0.7258661, -0.3142865, -0.7943594, -0.5554425, -0.7154426, -0.9325366, -0.5365588, -0.6726796,-0.5333700, -1.1490379, -0.7341835, -0.5674334, -0.5937559, -0.6678709, -0.9884700, -0.6867066, -0.7534590, -0.5746081, -0.6980742, -0.4437108),
  trait2 = c(-1.441314, -1.542336, -1.521273, -1.647513, -1.699523, -1.730442, -1.717863, -2.996804, -1.754690, -1.712154, -2.243960, -1.562110, -2.960188, -2.245098, -2.236504, -3.224924, -2.369385, -2.371202, -1.931709, -2.548280, -2.200219, -2.934834, -1.776365, -1.895126, -3.252790, -2.379290, -3.261246, -1.924550, -2.308284, -1.475468)
)

Trait1 and Trait2 have correlation coefficient of ~0.8:

cor(df$trait1, df$trait2, method = 'spearman')

We are especially interested in how trait 1 varies across genotypes and sites. A basic ANOVA using terms for genotype, environment, and their interaction gives the following:

> (model1 <- aov(trait1 ~ genotype * site, data=df))

Call:
   aov(formula = trait1 ~ genotype * site, data = df)

Terms:
                 genotype      site genotype:site Residuals
Sum of Squares  0.6445498 0.0515288     0.1477629 0.1287884
Deg. of Freedom        11         3             9         6

Residual standard error: 0.1465085
24 out of 48 effects not estimable
Estimated effects may be unbalanced

To account for the correlation we have tried including trait 2 as a fixed effect:

> (model2 <- aov(trait1 ~ genotype * site + trait2, data=df))

Call:
   aov(formula = trait1 ~ genotype * site + trait2, data = df)

Terms:
                 genotype      site    trait2 genotype:site Residuals
Sum of Squares  0.6445498 0.0515288 0.1208514     0.0980477 0.0576523
Deg. of Freedom        11         3         1             9         5

Residual standard error: 0.10738
24 out of 49 effects not estimable
Estimated effects may be unbalanced

When the cofactor is added the sum of squares for the genotype:site interaction and the residuals change. However, the sum of squares for the main effects are identical to seven decimal places. I find that very suspicious. The only reason I can think of for them to be identical is if the effect of trait2 is perfectly orthogonal to the fixed effects of genotype and site, but not the interaction term, which seems very unlikely.

I can't find any errors in how I have formatted the data in previous steps, and I also find it odd that the SS remain identical between two models even if I take arbitrary subsamples of the whole data. That suggests to me that it has something to do with some subtle aspect of the geometry of linear regression.

Can anyone tell me what could explain these identical effects?

PS I realise that model2 might not be the best way to analyse the data if there is a latent factor, but we are interested in understanding these models to help us work out what a possible latent factor might be.

Answer 1

This is because R computes type I sums of squares (SS) by default. With type I SS, the SS of the effects are computed in order of entry in the model formula, with no regard to effects that come after. So, in your case, the SS for genotype is computed first, without any regard to the other predictors in the model. Because it appears first in both models, its SS will be identical between them. Next comes site. The way R processes interactions is that they are always placed at the end. So in the first model, the interaction comes third, and in the second model, the interaction comes fourth, so it will have different SS in these two models.

If you want the type III SS, commonly reported by other statistical software, you can use car::Anova around your aov call. See this post and its linked posts for more information on what these SS mean.

Answer 2

Interesting question. I think you need to do a bit more detective work for your data (as already suggested in the comments to your post) and then report back here with your findings.

The first thing you need to look at for your entire dataset df is this command:

table(df$genotype, df$site)

This will tell you whether you have all possible combinations of genotype and site represented in your data. If that is the case, you would want all entries of the resulting table to be non-zero.

If the table command suggests that there are some possible combinations of gentotype and site NOT represented in your data, then that is a sign that you will NOT be able to fit a model which contains the interaction between genotype and site. If you do, you will see that some of the effects included in the model are NOT estimable:

m <- lm(trait1 ~ genotype*site, data = df)
summary(m)

(Try these last 2 commands for the data subset you shared here and you will see exactly which effects are NOT estimable.)