[R-sig-ME] How to analyze a 2x2 crossover with baseline in R?

Ben Bolker bbolker at gmail.com
Mon Jun 10 20:09:16 CEST 2013


Emmanuel Curis <emmanuel.curis at ...> writes:

> 
> Hello,
> 
> I have to test for the effect of a treatment in a crossover design
> with the two sequences TP/PT (T = treatment, P = placebo), but with
> also a single baseline measurement, before period 1. Sample data are
> given at the end of this message.
> 
> I can do the analysis without the baseline measurement, but when I try
> to include in it in the model, lme fails with the message
> 

> Erreur dans MEEM(object, conLin, control$niterEM) : 
>   Singularité rencontrée en résolution inverse au niveau 0, bloc 1
> 
> which should be in English something like « Error in MEEM(...).
> Singularity encountred when backsolving at level 0, bloc 1 ».
> 
> lmer also fails complaing for a non-definite matrix.
> 
> I guess this is related to the fact that period J0 is confounded with
> treatment = Baseline, as seen in 
> 
> > with( donnees, table( Periode, Traitement ) )
>        Traitement
> Periode Baseline Treatment Placebo
>      J0       16         0       0
>      M1        0         8       8
>      M2        0         8       7
> 
> as may also suggest the fact that not using both Periode and
> Traitement makes lme succeeds.
> 
> But how to write the model to fit, trying to reproduce the whole
> time-evolution of the patient outcome knowing the time and the
> condition? Does it make sense anyway, or fit periode 1 vs baseline and
> period 2 vs periode 1 the only way to to something?
> 
> When I use lm, using the patient as a fixed effect, it seems to work
> (some coefficients are NA, in agreement with the above problem I
> guess), and results suggests there is no period or sequence effect,
> especially if considering
> but I am not very satisfied with it...
> 

  The basic recipe here is to collapse your analysis into
a one-way analysis and then (hopefully) to reconstruct the
desired the contrasts.  I started to do this (below), but didn't
have quite enough time to finish it ...

donnees <- read.table(header=TRUE,
text="X7_SRF Periode Sequence Traitement Numero
548      J0 Groupe.A   Baseline      1
351      M1 Groupe.A  Treatment      1
418      M2 Groupe.A    Placebo      1
300      J0 Groupe.A   Baseline     10
180      M1 Groupe.A  Treatment     10
153      M2 Groupe.A    Placebo     10
652      J0 Groupe.A   Baseline     11
638      M1 Groupe.A  Treatment     11
491      J0 Groupe.B   Baseline     12
488      M1 Groupe.B    Placebo     12
443      M2 Groupe.B  Treatment     12
376      J0 Groupe.B   Baseline     14
466      M1 Groupe.B    Placebo     14
447      M2 Groupe.B  Treatment     14
638      J0 Groupe.A   Baseline     15
628      M1 Groupe.A  Treatment     15
872      M2 Groupe.A    Placebo     15
459      J0 Groupe.A   Baseline     16
191      M1 Groupe.A  Treatment     16
206      M2 Groupe.A    Placebo     16
731      J0 Groupe.B   Baseline     17
536      M1 Groupe.B    Placebo     17
584      M2 Groupe.B  Treatment     17
462      J0 Groupe.B   Baseline      2
373      M1 Groupe.B    Placebo      2
400      M2 Groupe.B  Treatment      2
197      J0 Groupe.B   Baseline      3
291      M1 Groupe.B    Placebo      3
160      M2 Groupe.B  Treatment      3
768      J0 Groupe.B   Baseline      4
260      M1 Groupe.B    Placebo      4
238      M2 Groupe.B  Treatment      4
606      J0 Groupe.A   Baseline      5
629      M1 Groupe.A  Treatment      5
350      M2 Groupe.A    Placebo      5
424      J0 Groupe.A   Baseline      6
304      M1 Groupe.A  Treatment      6
439      M2 Groupe.A    Placebo      6
538      J0 Groupe.A   Baseline      7
511      M1 Groupe.A  Treatment      7
501      M2 Groupe.A    Placebo      7
524      J0 Groupe.B   Baseline      8
601      M1 Groupe.B    Placebo      8
298      M2 Groupe.B  Treatment      8
511      J0 Groupe.B   Baseline      9
513      M1 Groupe.B    Placebo      9
444      M2 Groupe.B  Treatment      9")

library("nlme")
try(m1 <- lme( X7_SRF ~ Traitement + Periode, random = ~ 1|Numero,
       data = donnees,
       control = list( msVerbose = TRUE ) ))

donnees <- transform(donnees,TP=interaction(Traitement,Periode,drop=TRUE))

m2 <- lme( X7_SRF ~ TP, random = ~ 1|Numero,
       data = donnees,
       control = list( msVerbose = TRUE ) )
summary(m2)$tTable

I didn't succeed in the rest of this, but see
http://www.math.mcmaster.ca/~bolker/classes/s4c03/hw/hw2_s4c03.pdf
OR:
http://stackoverflow.com/questions/9335708/
   contrasts-for-lm-using-contrast-package
(broken URL)

for more information ...

contr1 <- matrix(
                 c(1,0,0,0,0,             # baseline
                    -1,1/2,1/2,0,0,       # baseline vs M1
                    0,-1/2,-1/2,1/2,1/2, # M1 vs M2
                   -1,1/2,0,1/2,0,       # baseline vs placebo
                    0,-1/2,1/2,-1/2,1/2), # placebo vs treatment
                    ## 0,1/2,-1/2,-1/2,1/2), # interaction:
                                          # (T2-P2)-(T1-P1) = T2+P1-T1-P2
              byrow=TRUE,ncol=5,
      dimnames=list(c("b","b_vs_M1","M1_vs_M2","b_vs_P","P_vs_T"),
                      c("B","PM1","TM1","PM2","TM2")))
                  
solve(contr1)



More information about the R-sig-mixed-models mailing list