Take the 2-minute tour ×
Cross Validated is a question and answer site for people interested in statistics, machine learning, data analysis, data mining, and data visualization. It's 100% free, no registration required.

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).

share|improve this question
    
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

2 Answers 2

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 
share|improve this answer
    
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

share|improve this answer

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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