[R] Analyzing an unbalanced AB/BA cross-over design

[Peter Dalgaard BSA]

martin wolfsegger <wolfseggerm at gmx.at> writes:

> I am looking for help to analyze an unbalanced AB/BA cross-over design by
> requesting the type III SS !
> 
> # Example 3.1 from S. Senn (1993). Cross-over Trials in Clinical
> Research
> outcome<-c(310,310,370,410,250,380,330,270,260,300,390,210,350,365,370,310,380,290,260,90,385,400,410,320,340,220)
> subject<-as.factor(c(1,4,6,7,10,11,14,1,4,6,7,10,11,14,2,3,5,9,12,13,2,3,5,9,12,13))
> period<-as.factor(c(1,1,1,1,1,1,1,2,2,2,2,2,2,2,1,1,1,1,1,1,2,2,2,2,2,2))
> treatment<-as.factor(c('F','F','F','F','F','F','F','S','S','S','S','S','S','S','S','S','S','S','S','S','F','F','F','F','F','F'))
> sequence<-as.factor(c(1,1,1,1,1,1,1,1,1,1,1,1,1,1,2,2,2,2,2,2,2,2,2,2,2,2))
> example<-data.frame(outcome,subject,period,treatment,sequence)
> 
> The recommended SAS code equals
> 
> PROC GLM DATA=example;
> CLASS subject period treatment sequence;
> MODEL outcome = treatment sequence period subject(sequence);
> RANDOM subject(sequence);
> RUN;
> 
> For PROC GLM, the random effects are treated in a post hoc fashion after the
> complete fixed effect model is fit. This distinction affects other features
> in the GLM procedure, such as the results of the LSMEANS and ESTIMATE
> statements. Looking only on treatment, period, sequence and subject effects, the
> random statement can be omitted.
> 
> The R code for type I SS equals
> 
> example.lm<-lm(outcome~treatment+period+sequence+subject%in%sequence,
> data=example)
> anova(example.lm)
> 
> Response: outcome
>                  Df Sum Sq Mean Sq F value    Pr(>F)    
> treatment         1  13388   13388 17.8416  0.001427 ** 
> period            1   1632    1632  2.1749  0.168314    
> sequence          1    335     335  0.4467  0.517697    
> sequence:subject 11 114878   10443 13.9171 6.495e-05 ***
> Residuals        11   8254     750                     
> 
> According to the unbalanced design, I requested the type III SS which
> resulted in an error statement
> 
> library(car)
> Anova(example.lm, type="III")
> 
> Error in linear.hypothesis.lm(mod, hyp.matrix, summary.model = sumry,  : 
>         One or more terms aliased in model.
>  
> 
> by using glm I got results with 0 df for the sequence effect !!!!
> 
> example.glm<-glm(outcome~treatment+period+sequence+subject%in%sequence,
> data=example, family=gaussian)
> library(car)
> Anova(example.glm,type="III",test.statistic="F")
> 
> Anova Table (Type III tests)
> 
> Response: outcome
>                      SS Df       F    Pr(>F)    
> treatment         14036  1 18.7044  0.001205 ** 
> period             1632  1  2.1749  0.168314    
> sequence              0  0                      
> sequence:subject 114878 11 13.9171 6.495e-05 ***
> Residuals          8254 11            
> 
> 
> The questions based on this output are
> 
> 1) Why was there an error statement requesting type III SS based on lm ?
> 2) Why I got a result by using glm with 0 df for the period effect ?
> 3) How can I get the estimate, the StdError for the constrast (-1,1) of the
> treatment effect ?

The type I/III issues my be better left to people who actually
understand them, but note that sequence is really a embedded in
subject (7 subjects had one sequence and 6 subjects the other) so
sequence:subject is equivalent to just subject, and this is likely
also the reason that sequence is aliased within subject. Basically,
this is a fixed-effects model, and if each subject has his own
parameter, you can't estimate the sequence effect. 

I'd do this kind of design with an explicit mixed-effects model:

> summary(aov(outcome~treatment+period+sequence+Error(subject)))

Error: subject
          Df Sum Sq Mean Sq F value Pr(>F)
sequence   1    335     335  0.0321  0.861
Residuals 11 114878   10443

Error: Within
          Df  Sum Sq Mean Sq F value   Pr(>F)
treatment  1 13388.5 13388.5 17.8416 0.001427 **
period     1  1632.1  1632.1  2.1749 0.168314
Residuals 11  8254.5   750.4

possibly substituting treatment:period for sequence (it's the same
thing). 

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