R help - May 2006 - post hoc comparison in repeated measure

Hi, I have a simple dataset with repeated measures.
one factor is treatment with 3 levels (treatment1,
treatment2 and control), the other factor is time (15
time points). Each treatment group has 10 subjects
with each followed up at each time points, the
response variable is numeric, serum protein amount. So
the between subject factor is treatment, and the
within subject factor is time. I ran a 2-way ANOVA
with repeated measures considering time as the within
subject factor:
aov(response~treat*time+Error(subject/time),dat)

The results told me that the treatment is marginally
significant (p=0.04). I would like to know where that
significance came from, so I did ALL pairwise t tests
(treat1 vs. control, treat2 vs. control, treat1 vs.
treat2) at each of the time point. There are 2 ways I
can do these t tests, using the MSE from the ANOVA
(the MSE used for the treatement effect in the ANOVA,
i.e. the treatment by time interaction) as the t test
error, or simply ran ordinary t tests using only the
data of the treatment levels in comparison. What I
found is that using the first approach, I couldn't
find any pairwise comparison statistically significant
which I thought I should find at least one
significant, because the ANOVA treatment effect is
mariginally significant (p = 0.04). Using the second
approach, I did find some pariwise comparisons
significant. Can anyone explain to me why?

BTW, is there a R function that can do post hoc
comparison on repeated measure ANOVA (from avo() with
Error term)?

I haven't heard from anyone with my previous post. I
guess I should post my dataset and my code here:

The dataset has one factor as treatment with 4 levels
(treatment1, treatment2, treatment3 and control), and
another factor as time (36 time points). On average
Each treatment group has 10 subjects with each
followed up at each time points. The response variable
is numeric, serum protein amount. So the between
subject factor is treatment, and the within subject
factor is time. I ran a 2-way ANOVA without
interaction with repeated measures considering time as
the within subject factor:

dat<-read.table("dat.txt", sep='\t', header=T,
row.names=1)

attach(dat)
## run aov:
fit<-aov(x~as.factor(treatment)+as.factor(Time)+Error(as.factor(animal)))
summary(fit)

## it should produce a ANOVA table as follow:

Error: as.factor(animal)
                     Df  Sum Sq Mean Sq F value 
Pr(>F)  
as.factor(treatment)  3 1405900  468633  3.0584
0.04051 *
Residuals            36 5516274  153230               
  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05
'.'
0.1 ' ' 1 

Error: Within
                  Df   Sum Sq  Mean Sq F value   
Pr(>F)    
as.factor(Time)   35  7298866   208539  22.796 <
2.2e-16 ***
Residuals       1365 12487181     9148                
     
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05
'.'
0.1 ' ' 1 

As you can see, the treatnment effect is marginally
significant. So I decided to do all pairwise
comparisons on treatments at each time point using t
tests. I tried 2 ways to do the t tests: one is the
regular t tests, i.e. only using the data in
comparison at each time point for t test computation,
then yes, I do see a lot of significant comparisons;
the 2nd approach is to use the MSE from ANOVA (i.e.
153230) as the error term for t test calculation; then
I don't see any significant comparisons which I think
is inconsistent with the marginally significant p
value (0.04) for treatment main effect. Below is my
code, did I do something wrong here, or my assumption
above is not right: significant main effect doesnot
necessarily mean at least one pairwise comparisons
must be significant.

## to do pairwise regular t test:
pval<-list()
for (j in seq(5,180,5) ) {
pval[[as.character(j)]]<-pairwise.t.test(dat2$x[dat2$Time==j],dat2$treatment[dat2$Time==j],pool.sd=F,var.equal=T,p.adjust.method='none')$p.value
}
sort(unlist(pval)[!is.na(unlist(pval))])

## as you can see, there are indeed many significant
pairwise comparisons

## to do pairwise t tests using MSE from ANOVA
pval <- NULL
for (i in 1:3) {
    for (k in (i+1):4) {
    for (j in seq(5,180,5) ) {
       
        ## ANOVA post test
        pval <- c(pval,
2*(1-pt(abs(mean(x[Time==j&treatment==k])-mean(x[Time==j&treatment==i]))/sqrt(153230*(1/sum(treatment[Time==j]==i)+1/sum(treatment[Time==j]==k))),
df=sum(treatment[Time==j]==i)+sum(treatment[Time==j]==k)-2)))

    }
    }
}   
sort(pval)
## As you can see, none of the p values is significant

Can anyone explain?

Thanks

BTW, is there a R function that can do post hoc
comparison on repeated measure ANOVA (from avo()
with Error term)?



__________________________________________________


-------------- next part --------------
An embedded and charset-unspecified text was scrubbed...
Name: dat.txt
Url:
https://stat.ethz.ch/pipermail/r-help/attachments/20060510/a4bcaa69/attachment-0003.txt

This is how I would do it.


arraychip <- read.table("arraychip.dat", sep='\t',
                        header=T, row.names=1)
for (i in 2:4) arraychip[[i]] <- factor(arraychip[[i]])


## run aov:
fit <- aov(x ~ treatment + Time + Error(animal), data=arraychip)
summary(fit)

## single stratum for the same ANOVA
fit2 <- aov(x ~ treatment + animal + Time, data=arraychip)
summary(fit2)

trt.means <- model.tables(fit2, cterms="treatment",
type="means")
treat.means <- as.vector(trt.means$tables$treatment)
treat.n <- as.vector(trt.means$n$treatment)
names(treat.means) <- levels(arraychip$treatment)


## this expression works in R and S-Plus
mi.mj <- outer(treat.means, treat.means, "-")
s2 <- summary(fit2)$"Mean Sq"[2]             ## S-Plus
s2 <- summary(fit2)[[1]]["animal","Mean Sq"] ## R
s2.n <- s2 / treat.n
si.sj <- sqrt(outer(s2.n, s2.n, "+"))
q.tukey <- qtukey(.95, 4 ,36) /sqrt(2)

mi.mj
lower <- mi.mj - q.tukey*si.sj
upper <- mi.mj + q.tukey*si.sj
lower
upper


## in S-Plus you can also use
multicomp.default(treat.means,
                  df.residual=summary(fit2)$Df[2],
                  vmat=diag(
                    ## rep(
                    summary(fit2)$"Mean Sq"[2]
                    ##   ,  length(trt.means))
                    / treat.n
                     ),
                  method="tukey")