Hello John, thank you for the fast reply.

Thanks to your answer I was able to reproduce the "between" residuals by
simply writing:
sum(residuals(mod.ok)^2)

But I must admit that how to obtain the "within" was beyond me, so some
advice here would be of great help.
(Also, it might be worth adding these residuals to the standard Anova
output, since I believe some researchers using the package could find it
useful when reporting there results)


By working with the Anova package, I stumbled into a conflict between the
Anova and the aov outputs.
The toy data used was reassembled from the (?aov) example (and filled in for
balance, and contr.sum contrasts where assigned).
When I compared the Anova (on SS type III) and the aov F (And p.value)
results - I found significant differences between them (though the SS where
identical).
If I understood correctly, for balanced designs both test should have
yielded the same results. any ideas on what the source of the difference
might be?

Here is the code:




## get the data:
utils::data(npk, package="MASS")
library(reshape) # for data manipulation
colnames(npk)[5] <- "value" # so that the data set will read proparly in
"reshape"

npk.wide = cast(block ~  N +P+ K, data = npk)

# now we fill in the NA's, to get a balanced design
is.na(npk.wide) #it has na's. we'll fill them with that factor's combination
mean, so to keep some significance around
combination.mean <- apply(npk.wide, 2, function(x){mean(x,na.rm = T)})
for(i in c(2:9))
{ npk.wide[is.na(npk.wide)[,i],i] = combination.mean[i-1] }
is.na(npk.wide) #no na's anymore

# recreating a balanced npk for the aov
npk.long <- melt(npk.wide)

head(npk.long,4)
#         block    value N P K
#X0_0_0       1 46.80000 0 0 0
#X0_0_0.1     2 51.43333 0 0 0
#X0_0_0.2     3 51.43333 0 0 0
#X0_0_0.3     4 51.43333 0 0 0
head(npk.wide,4)
#  block    0_0_0 0_0_1    0_1_0 0_1_1    1_0_0    1_0_1    1_1_0    1_1_1
#1     1 46.80000  52.0 54.33333  49.5 63.76667 57.00000 62.80000 54.36667
#2     2 51.43333  55.5 56.00000  50.5 59.80000 54.66667 57.93333 58.50000
#3     3 51.43333  55.0 62.80000  50.5 69.50000 54.66667 57.93333 55.80000
#4     4 51.43333  45.5 44.20000  50.5 62.00000 54.66667 57.93333 48.80000


# npk.wide$block # it is a factor - good.


# change into contr.sum
for(i in c(1,3:5)) npk.long[,i] <- C(npk.long[,i] , "contr.sum")
npk.wide[,1] <- C(npk.wide[,1] , "contr.sum")



## as a test, not particularly sensible statistically
npk.aovE <- aov(value~  N*P*K + Error(block), npk.long)
npk.aovE
summary(npk.aovE)

## output
Error: block
          Df  Sum Sq Mean Sq F value Pr(>F)
Residuals  5 153.147  30.629

Error: Within
          Df Sum Sq Mean Sq F value    Pr(>F)
N          1 378.56  378.56 39.1502 3.546e-07 ***
P          1  16.80   16.80  1.7378   0.19599
K          1 190.40  190.40 19.6911 8.657e-05 ***
N:P        1  42.56   42.56  4.4018   0.04319 *
N:K        1  66.27   66.27  6.8535   0.01298 *
P:K        1   0.96    0.96  0.0996   0.75415
N:P:K      1  74.00   74.00  7.6533   0.00899 **
Residuals 35 338.43    9.67
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## END output

# the residuals SS
# 338.43 + 153.147
# ==
sum(residuals(lm(value~  N*P*K , npk.long))^2)
#  491.58



idata  <- cast(N+P+K  ~ . , data = npk.long)[,c(1:3)]
idata
for(i in 1:3) idata[,i] <- C(idata[,i] , "contr.sum")

mod.ok <- lm(as.matrix(npk.wide[,-1]) ~  1, data = npk.wide)
av.ok <- Anova(mod.ok, idata=idata, idesign=~ N*P*K, type = "III")
summary(av.ok, multivariate=FALSE)


## output - with different F and p.value results
Univariate Type III Repeated-Measures ANOVA Assuming Sphericity

                SS num Df Error SS den Df         F    Pr(>F)
(Intercept) 144541      1      153      5 4719.0302 1.238e-08 ***
N              379      1       49      5   38.6145  0.001577 **
P               17      1       19      5    4.3915  0.090257 .
K              190      1       24      5   38.8684  0.001554 **
N:P             43      1       24      5    8.6888  0.031971 *
N:K             66      1       19      5   17.3195  0.008810 **
P:K              1      1       49      5    0.0983  0.766580
N:P:K           74      1      153      5    2.4161  0.180810
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## End output


# the SS residuals are the same thou
sum(residuals(mod.ok)^2)
#  491.58


















On Thu, Feb 19, 2009 at 12:24 AM, John Fox <jfox@mcmaster.ca> wrote:

> Dear Tal,
>
> I suppose that the "between" residuals would be obtained, for your example,
> by residuals(mod.ok). I'm not sure what the "within" residuals are. You
> could apply the transformation for each within-subject effect to the matrix
> of residuals to get residuals for that effect -- is that what you had in
> mind? A list of transformations is in the element $P of the Anova.mlm
> object.
>
> Regards,
>  John
>
> ------------------------------
> John Fox, Professor
> Department of Sociology
> McMaster University
> Hamilton, Ontario, Canada
> web: socserv.mcmaster.ca/jfox
>
>
> > -----Original Message-----
> > From: r-help-bounces@r-project.org [mailto:r-help-bounces@r-project.org]
> On
> > Behalf Of Tal Galili
> > Sent: February-18-09 4:04 PM
> > To: r-help@r-project.org
> > Subject: [R] [package-car:Anova] extracting residuals from Anova for Type
> > II/III Repeated Measures ?
> >
> > Hello dear R members.
> > I have been learning the Anova syntax in order to perform an SS type III
> > Anova with repeated measures designs (thank you Prof. John Fox!)
> > And another question came up: where/what are the (between/within)
> residuals
> > for my model?
> >
> >
> >
> > ############  Play code:
> >
> >
> > phase <- factor(rep(c("pretest", "posttest", "followup"), c(5, 5, 5)),
> >     levels=c("pretest", "posttest", "followup"))
> > hour <- ordered(rep(1:5, 3))
> > idata <- data.frame(phase, hour)
> > idata
> >
> > mod.ok <- lm(cbind(pre.1, pre.2, pre.3, pre.4, pre.5,
> >                      post.1, post.2, post.3, post.4, post.5,
> >                      fup.1, fup.2, fup.3, fup.4, fup.5) ~
> >  treatment*gender,
> >                 data=OBrienKaiser)
> > av.ok <- Anova(mod.ok, idata=idata, idesign=~phase*hour)
> >
> >
> > summary(av.ok, multivariate=FALSE)
> >
> > ## Univariate Type II Repeated-Measures ANOVA Assuming Sphericity
> > ##
> > ##                                  SS num Df Error SS den Df       F
> >  Pr(>F)
> > ## treatment                   211.286      2  228.056     10  4.6323
> >  0.037687
> > ## gender                       58.286      1  228.056     10  2.5558
> >  0.140974
> > ## treatment:gender            130.241      2  228.056     10  2.8555
> >  0.104469
> > ## phase                       167.500      2   80.278     20 20.8651
> > 1.274e-05
> > ## treatment:phase              78.668      4   80.278     20  4.8997
> >  0.006426
> > ## gender:phase                  1.668      2   80.278     20  0.2078
> >  0.814130
> > ## treatment:gender:phase       10.221      4   80.278     20  0.6366
> >  0.642369
> > ## hour                        106.292      4   62.500     40 17.0067
> > 3.191e-08
> > ## treatment:hour                1.161      8   62.500     40  0.0929
> >  0.999257
> > ## gender:hour                   2.559      4   62.500     40  0.4094
> >  0.800772
> > ## treatment:gender:hour         7.755      8   62.500     40  0.6204
> >  0.755484
> > ## phase:hour                   11.083      8   96.167     80  1.1525
> >  0.338317
> > ## treatment:phase:hour          6.262     16   96.167     80  0.3256
> >  0.992814
> > ## gender:phase:hour             6.636      8   96.167     80  0.6900
> >  0.699124
> > ## treatment:gender:phase:hour  14.155     16   96.167     80  0.7359
> >  0.749562
> >
> >
> >
> >
> >
> >
> >
> >
> >
> > --
> > ----------------------------------------------
> >
> >
> > My contact information:
> > Tal Galili
> > Phone number: 972-50-3373767
> > FaceBook: Tal Galili
> > My Blogs:
> > www.talgalili.com
> > www.biostatistics.co.il
> >
> >       [[alternative HTML version deleted]]
> >
> > ______________________________________________
> > R-help@r-project.org mailing list
> > https://stat.ethz.ch/mailman/listinfo/r-help
> > PLEASE do read the posting guide
> http://www.R-project.org/posting-guide.html
> > and provide commented, minimal, self-contained, reproducible code.
>
>
>


-- 
----------------------------------------------


My contact information:
Tal Galili
Phone number: 972-50-3373767
FaceBook: Tal Galili
My Blogs:
www.talgalili.com
www.biostatistics.co.il

	[[alternative HTML version deleted]]


