Specifying mixed-linear model in lme() for pre-/post-treatment design

Question

Dear crossvalidated community!

I am analyzing a data of a study in which every participant's performance in a memory task was measured twice (once under placebo, once under drug treatment in random order). Moreover every subject was tested in a blood test measuring the level of stable blood marker. I would like to analyze the data with respect to drug-effects and blood marker effects using lme() from the nlme() package.

'treatment' codes placebo/drug 'session' codes the first or second testing in order to account for learning/re-test effects 'blood' codes for the blood marker level (continous variable), which should be entered as a covariate in the analysis 'id' codes for the subjects

Would the following equation correct for the described data:

# some data
data <- structure(list(id = structure(c(1L, 12L, 14L, 15L, 16L, 17L, 
18L, 19L, 20L, 2L, 3L, 4L, 5L, 6L, 7L, 8L, 9L, 10L, 11L, 13L, 
1L, 12L, 14L, 15L, 16L, 17L, 18L, 19L, 20L, 2L, 3L, 4L, 5L, 6L, 
7L, 8L, 9L, 10L, 11L, 13L), .Label = c("1", "10", "11", "12", 
"13", "14", "15", "16", "17", "18", "19", "2", "20", "3", "4", 
"5", "6", "7", "8", "9"), class = "factor"), sess = structure(c(1L, 
2L, 2L, 1L, 1L, 2L, 2L, 1L, 1L, 2L, 2L, 1L, 1L, 2L, 2L, 1L, 1L, 
2L, 2L, 1L, 1L, 2L, 2L, 1L, 1L, 2L, 2L, 1L, 1L, 2L, 2L, 1L, 1L, 
2L, 2L, 1L, 1L, 2L, 2L, 1L), .Label = c("a", "b"), class = "factor"), 
    treat = structure(c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
    1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 
    2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L
    ), .Label = c("drug", "placebo"), class = "factor"), memory = structure(c(2L, 
    12L, 17L, 18L, 1L, 20L, 19L, 15L, 14L, 16L, 8L, 7L, 5L, 9L, 
    4L, 6L, 10L, 3L, 11L, 13L, 2L, 12L, 17L, 18L, 1L, 20L, 19L, 
    15L, 14L, 16L, 8L, 7L, 5L, 9L, 4L, 6L, 10L, 3L, 11L, 13L), .Label = c("1", 
    "10", "11", "12", "13", "14", "15", "16", "17", "18", "19", 
    "2", "20", "3", "4", "5", "6", "7", "8", "9"), class = "factor"), 
    blood = structure(c(3L, 11L, 17L, 18L, 20L, 19L, 9L, 10L, 
    6L, 16L, 5L, 2L, 12L, 15L, 4L, 13L, 14L, 8L, 7L, 1L, 3L, 
    11L, 17L, 18L, 20L, 19L, 9L, 10L, 6L, 16L, 5L, 2L, 12L, 15L, 
    4L, 13L, 14L, 8L, 7L, 1L), .Label = c("1", "100", "17", "25", 
    "38", "39", "40", "44", "55", "6", "62", "66", "70", "77", 
    "81", "83", "85", "86", "89", "94"), class = "factor")), .Names = c("id", 
"sess", "treat", "memory", "blood"), row.names = c(NA, -40L), class = "data.frame")

# the model
m1 <- lme(memory~sess*treat*blood, random=~1|id/treat, data=data)
summary(m1)

Answer 1

A few notes:

1. If you look at your data you obtain:

   > str(data)
   'data.frame':   40 obs. of  5 variables:
   $ id    : Factor w/ 20 levels "1","10","11",..: 1 12 14 15 16 17 18 19 20 2 ...
   $ sess  : Factor w/ 2 levels "a","b": 1 2 2 1 1 2 2 1 1 2 ...
   $ treat : Factor w/ 2 levels "drug","placebo": 1 1 1 1 1 1 1 1 1 1 ...
   $ memory: Factor w/ 20 levels "1","10","11",..: 2 12 17 18 1 20 19 15 14 16 ...
   $ blood : Factor w/ 20 levels "1","100","17",..: 3 11 17 18 20 19 9 10 6 16 ...
   

   The dependent variable memory and the covariate blood shouldn't be factors, but numeric:

   > data$memory <- as.numeric(data$memory)
   > data$blood <- as.numeric(data$blood)
   

2. I am no expert for regression designs, but why do you add the covariate and all it's interactions? Wouldn't it be sufficient to have something like memory ~ sess*treat + blood? I think everything on top of this needs justification.

3. The random effect needs to be stratified for both within-subjects factor sess and treat something that is cumbersome in lme and lme4 is preferable (for the cost of loosing straight forward p-values, but you can still get them via bootstrap). Something between the following extremes should work:

   lmer(memory ~ sess*treat + blood + (1|id), data = data)
   lmer(memory ~ sess*treat + blood + (1+sess*treat|id), data = data)
   

   Have a look at the paper from Baayen and Milin (2010) for different random effects structures. And to obtain p-values from a lmer model you can use the pvals.fnc function from the languageR package.

Baayen, R. H., & Milin, P. (2010). Analyzing Reaction Times. International Journal of Psychological Research, 3(2), 12–28.