[R] Analyzing an unbalanced AB/BA cross-over design
Peter Dalgaard BSA
p.dalgaard at biostat.ku.dk
Wed Jan 29 12:30:40 CET 2003
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
More information about the R-help
mailing list