# Partial, adjusted or Type-III SS of factorial experiment in R

The partial SS for two factor (N and P factors) factorial experiment with interaction can be calculated as:

\begin{eqnarray*} \textrm{SS}\left(N_{i}|\mu,P_{j},\left(NP\right)_{ij}\right) & = & \textrm{SS}\left(\mu,N_{i},P_{j},\left(NP\right)_{ij}\right)-\textrm{SS}\left(\mu,P_{j},\left(NP\right)_{ij}\right)\\ \textrm{SS}\left(P_{j}|\mu,N_{i},\left(NP\right)_{ij}\right) & = & \textrm{SS}\left(\mu,N_{i},P_{j},\left(NP\right)_{ij}\right)-\textrm{SS}\left(\mu,N_{i},\left(NP\right)_{ij}\right)\\ \textrm{SS}\left(\left(NP\right)_{ij}|\mu,N_{i},P_{j}\right) & = & \textrm{SS}\left(\mu,N_{i},P_{j},\left(NP\right)_{ij}\right)-\textrm{SS}\left(\mu,N_{i},P_{j}\right) \end{eqnarray*}

I coded these partial SS for two factor factorial experiment with interaction in R. I'm not getting the right answer and would appreciate if you point out my mistake. Thanks

# Partial SS
y <- c(55, 56, 57, 53, 54, 55, 51, 52, 53, 61, 62, 63)
Trt <- gl(  n = 4, k = 3, length = 4 * 3
, labels = c("N0P0", "N0P1", "N1P0", "N1P1")
, ordered = FALSE)
N <- gl(n = 2, k = 6, length = 2 * 6
, labels = c("Low", "High")
, ordered = FALSE)
P <- gl(n = 2, k = 3, length = 2 * 6
, labels = c("Low", "High")
, ordered = FALSE)
Data <- data.frame(y, Trt, N, P)
Fit1 <- aov(formula = y ~ N * P, data = Data)
ANOVASummary1 <- summary(Fit1)
print(ANOVASummary1)

library(MASS)

X <-
cbind(
model.matrix(object = y ~ 1,        data = Data)
, model.matrix(object = y ~ -1 + N,   data = Data)
, model.matrix(object = y ~ -1 + P,   data = Data)
, model.matrix(object = y ~ -1 + N:P, data = Data)
)

H  <- X  %*% ginv(t(X) %*% X) %*% t(X)

X1 <-
cbind(
model.matrix(object = y ~ -1 + N ,  data = Data)
, model.matrix(object = y ~ -1 + P,   data = Data)
, model.matrix(object = y ~ -1 + N:P, data = Data)
)
H1  <- X1  %*% ginv(t(X1) %*% X1)  %*% t(X1)
SSM1  <- t(y) %*% (H - H1) %*% y
SSM1

X2 <-
cbind(
model.matrix(object = y ~  1,       data = Data)
, model.matrix(object = y ~ -1 + P,   data = Data)
, model.matrix(object = y ~ -1 + N:P, data = Data)
)
H2  <- X2  %*% ginv(t(X2) %*% X2)  %*% t(X2)
SSM2  <- t(y) %*% (H - H2) %*% y
SSM2

X3 <-
cbind(
model.matrix(object = y ~  1,       data = Data)
, model.matrix(object = y ~ -1 + N,   data = Data)
, model.matrix(object = y ~ -1 + N:P, data = Data)
)

H3  <- X3  %*% ginv(t(X3) %*% X3)  %*% t(X3)
SSM3  <- t(y) %*% (H - H3)  %*% y
SSM3

X4 <-
cbind(
model.matrix(object = y ~  1,     data = Data)
, model.matrix(object = y ~ -1 + N, data = Data)
, model.matrix(object = y ~ -1 + P, data = Data)
)

H4  <- X4  %*% ginv(t(X4) %*% X4)  %*% t(X4)
SSM4  <- t(y) %*% (H - H4)  %*% y
SSM4

print(ANOVASummary1)


Edit

The print(ANOVASummary1) gives the ANOVA. I want to get the SS in ANOVA using Linear Model approach. Using the Linear Model approach I got the correct Sequential SS but could not figure out how to get the Partial, Adjusted or Type-III SS ( See here and here).

-
What answer do you believe is correct and why? Could you please reduce the length of this code to something as small and simple as possible that still exhibits the problem? –  whuber Dec 24 '12 at 16:00
Thanks @whuber for showing interest in my problem. I'm not able to get the right partial SS for N and P. Any idea! –  MYaseen208 Dec 24 '12 at 16:07
I don't quite follow your question. (Presumably, by "under" in the 1st sentence, you mean understand, & by "partial SS", you mean type I (sequential) Sums of Squares.) In addition to @whuber's good suggestions, could you elaborate on your question? Is this a real situation or for your own learning? What is it about type I SS you are trying to understand? Can you comment your code? What exactly are you trying to do & getting wrong? –  gung Dec 25 '12 at 17:26

Try this:

> library(car)
> Anova(Fit1, type="III")
Anova Table (Type III tests)

Response: y
Sum Sq Df F value    Pr(>F)
(Intercept)   9408  1    9408 1.425e-13 ***
N               24  1      24  0.001195 **
P                6  1       6  0.039969 *
N:P            108  1     108 6.364e-06 ***
Residuals        8  8
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

-
Thanks @Stat for your answer. How to get these SS of squares using Linear Model (matrix approach in this case) approach?. Thanks –  MYaseen208 Dec 26 '12 at 2:02

Type III SS are available from lm() via drop1()

> lmNP<-lm(y~N*P, data=Data)
> drop1(lmNP, y ~ N * P)
Single term deletions
Model:
y ~ N * P
Df Sum of Sq RSS    AIC
<none>                8  3.134
N       1        24  32 17.770
P       1         6  14  7.850
N:P     1       108 116 33.224


Compare with output from car::Anova() in the previous example. I got this from W Venables article Exegesis on linear models http://www.stats.ox.ac.uk/pub/MASS3/Exegeses.pdf

-