[R-sig-ME] Crossover design, degrees of freedom for t test

[Michael Young]

Hello,
I'm fitting a model for a crossover trial and I'm concerned that the
degrees of freedom given by lme4/lmer/lmerTest and also by nlme/lme
are wrong, or else I'm fitting the model incorrectly.

There are two treatments (U or F). Each participant was treated with
the U treatment twice on two separate days and the F treatment once on
another day. Order was randomized, and only one participant was
treated per day (so there's no overall random effect for day, but
possibly a random effect for the replicated days of the U treatment
*within* participant).

The outcome was measured on each day, immediately prior to treatment
then at a few time points after treatment. The pre-treatment
measurement allows the removal of a participant's daily variation in
that measure.

I've created a toy dataset to represent this design. This dataset has
10 participants, each observed on three different days (day A,B, C).
Within each day there were four measurements (time0,time1,time2,time3)
where time0 is the baseline (pre-treatment) measurement. Note that day
A corresponds to the F treatment and day B and C correspond to the U
treatment.

x <- structure(list(time = c("time0", "time1", "time2", "time3", "time0",
"time1", "time2", "time3", "time0", "time1", "time2", "time3",
"time0", "time1", "time2", "time3", "time0", "time1", "time2",
"time3", "time0", "time1", "time2", "time3", "time0", "time1",
"time2", "time3", "time0", "time1", "time2", "time3", "time0",
"time1", "time2", "time3", "time0", "time1", "time2", "time3",
"time0", "time1", "time2", "time3", "time0", "time1", "time2",
"time3", "time0", "time1", "time2", "time3", "time0", "time1",
"time2", "time3", "time0", "time1", "time2", "time3", "time0",
"time1", "time2", "time3", "time0", "time1", "time2", "time3",
"time0", "time1", "time2", "time3", "time0", "time1", "time2",
"time3", "time0", "time1", "time2", "time3", "time0", "time1",
"time2", "time3", "time0", "time1", "time2", "time3", "time0",
"time1", "time2", "time3", "time0", "time1", "time2", "time3",
"time0", "time1", "time2", "time3", "time0", "time1", "time2",
"time3", "time0", "time1", "time2", "time3", "time0", "time1",
"time2", "time3", "time0", "time1", "time2", "time3", "time0",
"time1", "time2", "time3"), day = c("A", "A", "A", "A", "B",
"B", "B", "B", "C", "C", "C", "C", "A", "A", "A", "A", "B", "B",
"B", "B", "C", "C", "C", "C", "A", "A", "A", "A", "B", "B", "B",
"B", "C", "C", "C", "C", "A", "A", "A", "A", "B", "B", "B", "B",
"C", "C", "C", "C", "A", "A", "A", "A", "B", "B", "B", "B", "C",
"C", "C", "C", "A", "A", "A", "A", "B", "B", "B", "B", "C", "C",
"C", "C", "A", "A", "A", "A", "B", "B", "B", "B", "C", "C", "C",
"C", "A", "A", "A", "A", "B", "B", "B", "B", "C", "C", "C", "C",
"A", "A", "A", "A", "B", "B", "B", "B", "C", "C", "C", "C", "A",
"A", "A", "A", "B", "B", "B", "B", "C", "C", "C", "C"), ID = c(1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L,
3L, 3L, 3L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 5L,
5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 6L, 6L, 6L, 6L, 6L,
6L, 6L, 6L, 6L, 6L, 6L, 6L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L,
7L, 7L, 7L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 9L,
9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 10L, 10L, 10L, 10L,
10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L), treat = c("F", "F",
"F", "F", "U", "U", "U", "U", "U", "U", "U", "U", "F", "F", "F",
"F", "U", "U", "U", "U", "U", "U", "U", "U", "F", "F", "F", "F",
"U", "U", "U", "U", "U", "U", "U", "U", "F", "F", "F", "F", "U",
"U", "U", "U", "U", "U", "U", "U", "F", "F", "F", "F", "U", "U",
"U", "U", "U", "U", "U", "U", "F", "F", "F", "F", "U", "U", "U",
"U", "U", "U", "U", "U", "F", "F", "F", "F", "U", "U", "U", "U",
"U", "U", "U", "U", "F", "F", "F", "F", "U", "U", "U", "U", "U",
"U", "U", "U", "F", "F", "F", "F", "U", "U", "U", "U", "U", "U",
"U", "U", "F", "F", "F", "F", "U", "U", "U", "U", "U", "U", "U",
"U"), outcome = c(96.5, 94.7, 77.4, 92.4, 98.1, 102.2, 104.8,
100.9, 106, 104.8, 105.7, 114.5, 78.2, 81.1, 72.8, 82.3, 88.6,
82.1, 80.2, 79.3, 64.6, 69, 72, 75.4, 98.8, 92.7, 79.9, 82.5,
76.8, 75.9, 82.9, 86.7, 89.1, 89.1, 85.9, 82.7, 93.3, 83, 94.6,
85, 72.6, 83.1, 90.3, 94.5, 93.1, 76.4, 76.5, 81.5, 97.7, 88.9,
76.4, 81.8, 96.6, 93.5, 99.2, 97.7, 103.9, 100.6, 87.5, 86.2,
77.5, 84.7, 85.2, 89.9, 84.6, 74.6, 81.8, 86.8, 71, 75.3, 85.6,
85.1, 80.9, 80.8, 70.3, 72.3, 79.1, 67.2, 71.6, 75.7, 82.8, 74.2,
76.8, 76.1, 101, 104.5, 95.1, 105.7, 81.9, 83.3, 88.3, 92.5,
102.5, 92.6, 94, 92.6, 105.8, 103.4, 79.4, 84.9, 103.2, 100.4,
103.1, 107.6, 92.7, 90.8, 95.1, 104.6, 78.1, 96.5, 97.3, 99.7,
64.1, 73.2, 86.7, 83.2, 80.1, 78.4, 107.3, 100.5)), .Names = c("time",
"day", "ID", "treat", "outcome"), class = "data.frame", row.names = c(NA,
-120L))

#I've fit this model:
library(lme4)
library(lmerTest)
m1 <- lmer(outcome~treat*time + (1|ID/day),data=x)
summary(m1)

#or in nlme:
library(nlme)
m1_nlme <- lme(outcome~treat*time, random=~1|ID/day,data=x)
summary(m1_nlme)

I'm confused by the degrees of freedom indicated here for the t test
corresponding to the interaction effect of any specific timepoint.
This interaction is the causal effect of treatment relative to
baseline comparing treatment F with treatment U. It seems to me that
there are truly only 10 participants, so I'm not sure why lmerTest (as
well as nlme) is indicating there are 84 degrees of freedom for this
test.  At most, there are 20*3 observations for any specific
time-point by exposure interaction (because each of the 10
participants each received one baseline and one specific
post-treatment observation on three separate days).

There should be a rough correspondence here to a difference in
difference t test (which requires averaging across replicated days of
the U treatment), and that difference in difference t test has 9
degrees of freedom:

###subtract time0 (baseline) from a specific timepoint eg time2
seperately for treatments
time2F <- x[time=="time2"&treat=="F",][order(ID),]$outcome -
x[time=="time0"&treat=="F",][order(ID),]$outcome

#Since U is replicated, do this seperately for each replicate then average.
time2U_B <- x[time=="time2"&treat=="U"&day=="B",][order(ID),]$outcome -
  x[time=="time0"&treat=="U"&day=="B",][order(ID),]$outcome

time2U_C <- x[time=="time2"&treat=="U"&day=="C",][order(ID),]$outcome -
  x[time=="time0"&treat=="U"&day=="C",][order(ID),]$outcome

time2U <- rowMeans(cbind(time2U_B,time2U_C))

t.test(time2F,time2U,paired=TRUE)

The p values are different in part because the standard errors are
different but also because of the substantially different degrees of
freedom being used values.

Note that in my actual data I have even more timepoints--ie I actually
measure the outcome about 10 times after treatment on each day. This
leads to very large degrees of freedom in the mixed effects model
despite the fact that this is really only 10 participants. Does
measuring the outcome more times in the same people really increase
our ability to estimate the distribution of the test statistic for a
specific time by exposure interaction? I'm not convinced that really
makes sense, because if I had measured the outcome 1000 times after
treatment in just three individuals I would probably find a
significant effect even if just one of the participants responded to
treatment. I have no way of accounting for the fact that most
individuals might not have a different response between treatments but
I just "got lucky" and selected a single individual who does respond.


It seems like I'm implying that I'm concerned there might be
participant-specific effects. So to allow participant-specific effects
of treatment on change from baseline measurement, how would I do that?

The following gives me a warning message:

m2 <- lmer(outcome~treat*time + (1 + treat|ID/day),data=x)
#Warning message:
#In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv,  :
#  Model is nearly unidentifiable: large eigenvalue ratio
 #- Rescale variables?

summary(m2)
#Error in calculation of the Satterthwaite's approximation. The output
of lme4 package is returned
#summary from lme4 is returned
#some computational error has occurred in lmerTest

This fits in nlme, but still gives the same degrees of freedom (84, in
this case)

m2_nlme <- lme(outcome~treat*time, random=~treat|ID/day,data=x)
summary(m2_nlme)


In any case, I'm not convinced this random slope is even allowing a
participant specific effect of treatment on change from baseline so
much as a participant-specific effect of treatment day (including
observations on the treatment day which preceded treatment).


It seems that I actually want something like this, but this doesn't fit:
m3 <- lmer(outcome~treat*time + (1 + time:treat|ID/day),data=x)

I'm wondering if it's even possible to estimate participant-specific
time by treatment effects without replicates measures within each
day's timepoint?

That's as far as I've gotten. If anyone could shed some light it would
be much appreciated.

Thanks,
Michael