[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