I have two factors j (with 3 levels A, B, C) and k (with 3 levels M, N, O), with k nested within j. Level A of j is the reference level and it has only one k level, M, in it.
What I want to test with the model is: (1) are the means for each non-control j level significantly different than the mean for j="A"? (2) is the mean of each j:k interaction significantly different from the j mean?
Here I generate a synthetic data set where every j:k level-pair is significant but every level of j is not (according to the above criteria), then use the data to fit the model "x ~ j/k" (k is nested in j), and look at the summary data:
> set.seed(0)
> j = gl(3, 30, labels=LETTERS[1:3])
> k = gl(3, 10, length=90, labels=LETTERS[13:15])
> k[1:30] = "M" # Reference j="A" has K="M" only
> x = 100*(2^as.numeric(k)) # 200,400,800 for k=M,N,O
> x[1:30] = mean(c(200,400,800)) # Reference j="A" has x=same mean as "B" and "C"
> x = x + rnorm(90, 0, 4) # A little randomness
> x = x + rep(c(0, 1, 1.2), each=30) # A slight non-significant deviation by j.
> df = data.frame(x=x, j=j, k=k)
> f = lm(x ~ j/k, data=df)
> summary(f)
Call:
lm(formula = x ~ j/k, data = df)
Residuals:
Min 1Q Median 3Q Max
-8.3711 -2.7546 -0.0538 2.2673 10.2899
Coefficients: (2 not defined because of singularities)
Estimate Std. Error t value Pr(>|t|)
(Intercept) 100.0878 0.6635 150.84 <2e-16 ***
jB 101.5278 1.3270 76.51 <2e-16 ***
jC 100.6014 1.3270 75.81 <2e-16 ***
jA:kN NA NA NA NA
jB:kN 198.9841 1.6253 122.43 <2e-16 ***
jC:kN 200.3907 1.6253 123.30 <2e-16 ***
jA:kO NA NA NA NA
jB:kO 598.8600 1.6253 368.46 <2e-16 ***
jC:kO 601.9267 1.6253 370.35 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 3.634 on 83 degrees of freedom
Multiple R-squared: 0.9998, Adjusted R-squared: 0.9998
F-statistic: 8.132e+04 on 6 and 83 DF, p-value: < 2.2e-16
At first glance it appears that part of what I wanted to test is present in the table. I should see six significant entries besides the intercept, and I do. But closer inspection shows that none of the summary table entries, showing results of testing the 7 coefficients generated by the model, is actually testing what I want to test!
Rows jB and jC are significant. I expected these would be testing the jB and jC means against the jA mean and thus would not be significant. But that is not what they are testing! Row jB:kN and the other three interaction rows are significant. I expected these would be testing the jB:kN (etc.) means against the jB (etc.) mean, but they are not!
Notice that jB and jC coefficients in the "Estimate" column are significantly different than 0. This indicates the summary table jB and jC entries are not showing the test of a coefficient representing the jB and jC means minus the jA mean, since these were constructed to be insignificantly different. Actually, these coefficients are comparing the jB:kM and jc:kM nested group means to the control group jA mean, as verified by looking at the difference of the "jB:kM" and "jC:kM" means and the intercept:
> mean(df$x[df$j == "B" & df$k == "M"]) - mean(df$x[df$j == "A"])
[1] 101.5278
> mean(df$x[df$j == "C" & df$k == "M"]) - mean(df$x[df$j == "A"])
[1] 100.6014
These are exactly the jB and jC coefficient values. Thus coefficient jB is significant if the jB:kM group mean differs significantly from the jA control group mean. That is not the coefficient or test I expected and is not one in which I am interested. Instead, I want to see a coefficient that is the entire j level's mean minus the control mean, so the test tells if the j level differs from the control. Also, I want two more coefficients that are the jB:kM and jC:kM nested level's means minus the j level mean, so the test tells if the nested level differs from the mean of the level it is nested within. How can I respecify the model (or contrast matrix?) to make this happen?
Are the jB:kN coefficient and other three interaction coefficients equal to the means within that nested level minus the containing level's mean? They are not:
> mean(df$x[df$j == "B" & df$k == "N"]) - mean(df$x[df$j == "B" & df$k == "M"])
[1] 198.9841
This value is the jB:kN coefficent, so it is the difference between the jB:kN group mean and the jB:kM group mean! Again that is not the coefficient or test I expected and is not one in which I am interested. Instead, I want to see a coefficient that is the jB nested level's mean minus the j level mean, so the test tells if the nested level differs from the mean of the level it is nested within.
In summary, I have three problems:
I want the reference level mean compared to the mean of each outer j level, NOT to the mean of the first nested k level within the j level.
I want my nested levels compared to the mean of the j level that contains them, NOT to the reference level mean or the "pseudo-reference" level that is the first nested level within each j level.
I want additional tests to appear in the summary table, comparing the last nested k level means to the mean of the j level that contains them.
How can I do this?
