[R-sig-ME] Reformat: Assistance with specification of crossover design model

[Adam Smith]

Thanks very much for your thoughts, Ben. I've added some additional thoughts below
that I hope you'll correct me on.  

Best,
Adam

----------------------------------------

>
>> I'm analyzing data (variable names in brackets) generated by a
>> 4-period [period], 3-treatment [trt] crossover design. The design
>> was strongly balanced (all treatments preceded all others, including
>> itself) and uniform within period (all treatments occurred the same
>> # of times in each of the 4 periods).
>
>> This produced 18 distinct sequences of treatments (e.g., ABCC, CAAB,
>> etc.), to which a single individual [id] was randomly assigned. Two
>> covariates [cov1, cov2] were also measured on each individual.
>
>> The response [y] was continuous, and each individual was associated
> with a pre-study baseline measure [y0].
>
>> Because a single individual was assigned to each sequence, I believe
>> that a random effect for each individual will capture individual,
>> baseline, and sequence effects (they're perfectly confounded).
>
>> The question is simple: does the response differ among treatments?
>> I'm hoping the correct specification is as simple as I think it is,
>> illustrated below with mock data. My specific questions:
>
>
> Thanks for a very clear question. I'm going to take a stab at
> these, but haven't thought *very* deeply about them; hopefully
> this will inspire corrections/alternative answers.
>
>> 1 - Does this model correctly capture treatment, period, and
>> potential carryover effects? As I understand it, in a strongly
>> balanced design such as this, carryover and treatment effects are
>> not confounded, so my assumption is that I don't have to specify any
>> addition variable to capture this effect (e.g., a variable
>> indicating the treatment in the preceding period). I'm happy to be
>> corrected.
>
> I think you're right that the individual-level random effect
> controls for pre-treatment and sequence/carryover effects. However,
> I'm not sure it will capture period or treatment-by-period effects
> (i.e. residuals within the same period might be correlated, and
> residuals from individuals who received the treatment in the same
> period might be correlated).

On further thought, I think the carryover effects should be modeled
explicitly. In fact, we can distinguish between two types of carryover
effects - so-called 'self' (treatment follows itself) and 'mixed' 
(treatment follows a different treatment) carryover effects.  A common 
model to analyze in this case is:

Ydks = mu + Ak + Td(k,s) + I*Sd(k-1, s) + (1-I)*Md(k-1, s) + Us + eks

where:

Ydks = response to treatment d in period k for subject s
Ak = effect of period
Td = effect of treatment
I = indicator if treatment in period k-1 was the same as in period k
Sd = Self carryover effect of treatment d
Md = Mixed carryover effect of treatment d
   (NOTE: no carryover effect in period 1; i.e., Sd and Md = 0 in period 1)
   (I'm not sure how to distinguish this from treatment "0"; currently in 
    period 1, I've set Sd = Md = NA, but this is incorrect b/c it drops from
    consideration the observations from period 1)
Us = subject effect (random)
eks = error term

The corresponding LMM is specified, I think, as follows:

datURL <- "https://www.dropbox.com/s/qyer7qrl1ay2h22/xover_test.csv"

# Download data
dat <- repmis::source_data(datURL, sep = ",", header = TRUE)
dat <- within(dat, {
              id = factor(id)
              period = factor(period)
              cov1 = factor(cov1)
              cov2 = factor(cov2)
              selfco = factor(self * prevtrt)
              mixedco = factor((1-self) * prevtrt)
              trt = factor(trt)
              })
str(dat)


# Proposed linear mixed model analysis
require(lme4)
mod1 <- lmer(y ~ trt + period + selfco + mixedco + cov1 + cov2 + 
                (1|id), data = dat)
summary(mod1)

I tried including the baseline measurement (y0) as you suggested, but it 
reduced the variance estimate of (1|id) to zero.  I'm inclined, then, to 
drop the baseline measurement and leave the random intercept for individuals.  
In fact, the variance is very near zero even when excluding y0. 

>
> One way to see whether you might have missed something would
> be to draw (e.g.) boxplots of residuals by whatever groupings
> you might be concerned about.
>
>
>> 2 - Is an interaction between treatment and period prescribed?
>> Because the design is uniform within periods, I believe that period
>> and treatment effects are likewise not confounded.
>
> I think so.
>
>> 3 - Does the random effect for each individual captures variation
>> due to individual, sequence, and baseline measurement?
>
> I think so.
>
>> Thanks very much for your assistance,
>>
>> Adam Smith
>> Department of Natural Resources Science
>> University of Rhode Island
>>
>> # Download data
>> datURL <- "https://dl.dropboxusercontent.com/u/23278690/xover_test.csv"
>> dat <- repmis::source_data(datURL, sep = ",", header = TRUE)
>> dat <- within(dat, {
>> ? ? ?id = factor(id)
>> ? ? ?trt = factor(trt)
>> ? ? ?period = factor(period)
>> ? ? ?cov1 = factor(cov1)
>> ? ? ?cov2 = factor(cov2)
>> })
>>
>> require(lme4)
>> # Proposed linear mixed model analysis
>> mod1 <- lmer(y ~ trt + period + cov1 + cov2 + (1|id), data = dat)
>
> I think you do want trt*period (possibly as a random effect?)
> In a setting where I had access to some regularization (e.g.
> blme) I would be tempted to treat period as random -- but I wouldn't
> do it in this vanilla model, because there aren't enough periods to
> estimate it reliably.

I'm still unsure about the trt*period interaction.  It will eat several degrees of
freedom as a fixed effect.  Including it as a random effect (1|trt:period) produces
an estimate of zero variance.  What is my interpretation if treating it as a random
effect?

>
> While you don't _need_ to incorporate pre-treatment covariate,
> carryover effects, etc., it might be interesting -- I think these
> would be partitioned out of the among-individual variation ...
>
>> # Test global treatment effect, for example...
>> mod2 <- update(mod1, . ~ . - trt)
>> anova(mod1, mod2)
>>
>> sessionInfo()
>>
>> R version 3.0.3 (2014-03-06)
>>
>> Platform: x86_64-w64-mingw32/x64 (64-bit)
>
>> attached base packages:[1] stats graphics grDevices utils datasets
> methods base
>
>> other attached packages:[1] car_2.0-19 AICcmodavg_1.35 lme4_1.1-5
>> Rcpp_0.11.1 Matrix_1.1-3 plyr_1.8.1
>