1
$\begingroup$

Dataset

The data is a typical dataset encounted in two-way classification model: enter image description here

Situation

I applied two-way anova by running R on this dataset to look at the sum of squares

fullwithinteract <- lm(torque~hole+barrel+hole:barrel,df)
anova(fullwithinteract)

And the results is as follows: enter image description here

The problem arises when I tried to compute Sum Sq for barrel using a design matrix. I used the design matrix:

List item

Then let $Y$ be the response vector, and $X$ be the matrix above, I go on to compute in R:

sum(( Y -   X %*% solve(t(X)%*%X) %*% t(X) %*% Y)^2)

which gives an answer of 14010.5, while according to the anova table output by R, the answer should be 26565.1+1910.5=28475.6.

Question

What have I misunderstood? Is this even calculatable using a design matrix approach?

REMARKS I hope to remark that using the following two matrices (left for calculating the SSE associated with the full model; right for calculating the SSE associated with the reduced model that nullify the interaction effects), the SSEs computed coincide with that computed by R. Now the main problem is the discrepancy between the SSEs computed using a design matrix and R for the factors barrel and hole. enter image description here

REMARKS 2 I would like to give some details to how $Y$ is computed. Calling the original dataset df, $Y$ is given by

df$torque

And the data and the design matrix used in the situation above can be found at: https://www.notion.so/hephaes/Data-e4ed9065085e49a8bc91c760ebba1bf3

$\endgroup$
6
  • $\begingroup$ Your design matrix needs an intercept. $\endgroup$ – Noah May 12 '20 at 15:18
  • $\begingroup$ Thx for the response. Is the first column of all 1 in the design matrix the intercept column? If it is not how exactly should I add an intercept? $\endgroup$ – hephaes May 12 '20 at 15:20
  • $\begingroup$ Yes; X <- cbind(1, X) should work. $\endgroup$ – Noah May 12 '20 at 15:20
  • $\begingroup$ @Noah Then there will be two columns of $1$s next to each other. $\endgroup$ – Dave May 12 '20 at 15:21
  • 2
    $\begingroup$ You're right, my mistake. I saw three columns and assumed there was an interaction and two dummies, where the rest of the data simply had been cut off. $\endgroup$ – Noah May 12 '20 at 15:40
4
$\begingroup$

The reason is that R uses type 1 sums of squares (SS) by default and the formula you used to calculate SS for barrel relies on using type 3 SS.

For type 1 SS, the order of the terms entered into the model matters; the SS for a predictor is the difference between the SS error (SSE) for the model with only the terms entered before the predictor and the SSE for the model with the terms entered before the predictor and the predictor itself. Using your method of calculating SSE, we get:

mm <- model.matrix(torque ~ hole+barrel+hole:barrel, data = df)

anova(fullwithinteract)
#> Analysis of Variance Table
#> 
#> Response: torque
#>             Df  Sum Sq Mean Sq F value   Pr(>F)   
#> hole         1  1035.1  1035.1  2.1672 0.214959   
#> barrel       1 26565.1 26565.1 55.6192 0.001727 **
#> hole:barrel  1    55.1    55.1  0.1154 0.751151   
#> Residuals    4  1910.5   477.6                    
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

First, we can calculate the SSE for the model that contains only intercept and hole:

X1 <- mm[,1:2] #Intercept and hole
sum((Y - X1%*% solve(t(X1)%*%X1) %*% t(X1) %*%Y)^2)
#> [1] 28530.75

Then we can calculate the SSE for the model that contains intercept, hole, and barrel:

X2 <- mm[,1:3] #Intercept, hole, and barrel
sum((Y - X2%*% solve(t(X2)%*%X2) %*% t(X2) %*%Y)^2)
#> [1] 1965.625

Below we see that the difference between these is the SS that R reports for barrel using anova().

28530.75 - 1965.625
#> [1] 26565.12

For type 3 SS, the SS for a predictor is the difference between the SSE for the model with all terms other than the predictor and the SSE for the model with all terms, which is how you calculated your SS for barrel manually. We can use the car package to get the type 3 SS and verify that your calculation is correct.

car::Anova(fullwithinteract, type = 3)
#> Anova Table (Type III tests)
#> 
#> Response: torque
#>             Sum Sq Df  F value   Pr(>F)    
#> (Intercept)  68450  1 143.3133 0.000279 ***
#> hole           306  1   0.6412 0.468141    
#> barrel       12100  1  25.3337 0.007317 ** 
#> hole:barrel     55  1   0.1154 0.751151    
#> Residuals     1910  4                      
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

First, we can calculate the SSE for the model that contains all terms except barrel:

X1 <- mm[,-3] #All except barrel
sum((Y - X1%*% solve(t(X1)%*%X1) %*% t(X1) %*%Y)^2)
#> [1] 14010.5

Then we can calculate the SSE for the model that contains all terms:

X2 <- mm #All
sum((Y - X2%*% solve(t(X2)%*%X2) %*% t(X2) %*%Y)^2)
#> [1] 1910.5

Here we see that the difference matches the type 3 SS for barrel reported by car::Anova().

14010.5 - 1910.5
#> [1] 12100
$\endgroup$
1
  • $\begingroup$ Thank you for the detailed answer! I also searched for type I, II, III SS just before your answer was posted, and indeed this is the only correct explanation. I thought the SS calculated in my lecture notes is the only way (but clear this is not now). The reference is here: mcfromnz.wordpress.com/2011/03/02/… So future readers may consider consulting this for more information. $\endgroup$ – hephaes May 12 '20 at 17:39

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.