R help - May 2002 - nlme cross-over and fixed nested

I have problem getting the concept of a nested fixed variable into the nlme
scheme. I fear the question is very stupid.  In the past I had asked this
before, and never got a reply (in other cases, the response was within
hours). I also checked the S-list, where several similar enquiries of other
people are orphaned.

We have a cross-over design, where patient are treated two weeks with
placebo, two weeks with a drug ("verum"), and a washout phase in
between.
For each phase, there is ONE number (let's call it "glucose
level")
describing the treatment effect. Blind treatement sequence is Seq=1 for
placebo/verum, Seq=2 for verum/placebo.

In Pinheiro/Bates, cross-over trials are mentioned on page 107, but the
parameter "Side" in the example is not exactly what we have.
We used the scheme below to analyze the data (there is a small effect of
imperfect washout in the example below), but a reviewer complained that the
analysis is not correct, because Patient should be nested in Seq. As I can
see, the results look reasonable, but I admit that nesting is nowhere
indicated.
I know how to handle the case is Seq if assumed a random effect, but I don't
see why it should be random, AND I want an estimate of the effect.

So, if the question is stupid, please help me out at least an tell my WHY.


# cross-over design
library(nlme)
NSubj<-20
# Placebo Phase
# Seq: 1 for "Placebo first phase", 2 for "Verum first
phase"
glup<-data.frame(Patient=as.factor(1:NSubj),
                Seq=as.factor(rep(c(1,2),NSubj/2)),
                Treat=rep("Pla",NSubj),
                Glu=rnorm(NSubj,30,4)
   )
# Verum Phase
gluv<-data.frame(Patient=as.factor(1:NSubj),
                Seq=as.factor(rep(c(1,2),NSubj/2)),
                Treat=rep("Ver",NSubj),
                Glu=glup$Glu-rnorm(NSubj,6,2)
   )

# Simulate imperfect washout when verum was in first phase
# == Interaction of Seq and Treat
glu<-rbind(gluv,glup)
glu$Glu<-glu$Glu+as.integer((glu$Seq=="2")&(glu$Treat=="Pla"))*rnorm(NSubj,2
,3)

# Sequence as fixed effect
glu.lme<-lme(Glu~Treat*Seq,data=glu,random=~1|Patient)
print(summary(glu.lme))

# result looks perfectly valid for me
#                  Value Std.Error DF   t-value p-value
#(Intercept)   24.694617  1.602160 18 15.413324  <.0001
#TreatPla       5.106889  1.096139 18  4.658980  0.0002
#Seq2           0.466712  2.265797 18  0.205981  0.8391
#TreatPla:Seq2  3.525737  1.550174 18  2.274413  0.0354

# Sequence as random ???????
glu.lme<-lme(Glu~Treat,data=glu,random=~1|Seq/Patient)
print(summary(glu.lme))



Dieter Menne


---------------------------------------
Dr. Dieter Menne
Biomed Software
72074 T?bingen
Tel (49) (7071) 52176
FAX (49) (7071) 55 10 46
dieter.menne at menne-biomed.de
www.menne-biomed.de



-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !)  To: r-help-request at
stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._

"Dieter Menne" <dieter.menne at menne-biomed.de> writes:
> NSubj<-20
> # Placebo Phase
> # Seq: 1 for "Placebo first phase", 2 for "Verum first
phase"
> glup<-data.frame(Patient=as.factor(1:NSubj),
>                 Seq=as.factor(rep(c(1,2),NSubj/2)),
>                 Treat=rep("Pla",NSubj),
>                 Glu=rnorm(NSubj,30,4)
>    )
> # Verum Phase
> gluv<-data.frame(Patient=as.factor(1:NSubj),
>                 Seq=as.factor(rep(c(1,2),NSubj/2)),
>                 Treat=rep("Ver",NSubj),
>                 Glu=glup$Glu-rnorm(NSubj,6,2)
>    )
I don't think you mean that. One of them wants to have 

Seq=as.factor(rep(c(2,1),NSubj/2))

otherwise patient 1 has two 1st phases, pt.2 has two 2nd phases, etc.

For a balanced design, I'd suggest

summary(aov(Glu~Treat*Seq+Error(Patient),data=glu))

which should show the wellknown crossover-design feature that the
interaction is compared to the Patient variation.


Or, put differently, it makes little sense to use the Seq variable, as
you define it, as a main effect. Rather, you'd want an indicator of
whether the current observation is the first or the second one, so
> glu <- transform(glu,
first=as.factor((Seq==1)==(Treat=="Pla")))
> summary(aov(Glu~Treat*first+Error(Patient),data=glu))
Error: Patient
            Df Sum Sq Mean Sq F value Pr(>F)
Treat:first  1  14.96   14.96  0.2739 0.6071
Residuals   18 982.75   54.60               

Error: Within
          Df Sum Sq Mean Sq F value    Pr(>F)    
Treat      1 478.04  478.04 75.9305 7.106e-08 ***
first      1   2.66    2.66  0.4221    0.5241    
Residuals 18 113.32    6.30                      
---
Signif. codes:  0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` '
1


Notice that this reverses the roles of interaction and main effect
relative to using Seq directly:
> summary(aov(Glu~Treat*Seq+Error(Patient),data=glu))
Error: Patient
          Df Sum Sq Mean Sq F value Pr(>F)
Seq        1  14.96   14.96  0.2739 0.6071
Residuals 18 982.75   54.60               

Error: Within
          Df Sum Sq Mean Sq F value    Pr(>F)    
Treat      1 478.04  478.04 75.9305 7.106e-08 ***
Treat:Seq  1   2.66    2.66  0.4221    0.5241    
Residuals 18 113.32    6.30                      
---
Signif. codes:  0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` '
1





-- 
   O__  ---- Peter Dalgaard             Blegdamsvej 3  
  c/ /'_ --- Dept. of Biostatistics     2200 Cph. N   
 (*) \(*) -- University of Copenhagen   Denmark      Ph: (+45) 35327918
~~~~~~~~~~ - (p.dalgaard at biostat.ku.dk)             FAX: (+45) 35327907
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !)  To: r-help-request at
stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._