1
$\begingroup$

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)
$\endgroup$
2
  • $\begingroup$ Is there a particular reason you're using nlme rather than its successor, lme4? $\endgroup$
    – Macro
    Jun 10, 2012 at 16:16
  • $\begingroup$ I was hoping to get some p-values and pvals.fnc() will only work with the most simple random effects structure. $\endgroup$
    – jokel
    Jun 11, 2012 at 19:46

1 Answer 1

2
$\begingroup$

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.

$\endgroup$
2
  • $\begingroup$ Thanks for your helpful answer! I was hoping to find a way to use lme() in order to report p-values. Unfortunately it seems that pvals.fnc() will only work with a lmer() model 'when there is more than one random intercept in the model'(R Documentation). So the only random effects specification I would be left with is the first one you mention: lmer(memory ~ sess*treat + blood + (1|id), data = data). So it seems I have to give up on p-values and use likelihood ratio tests instead? $\endgroup$
    – jokel
    Jun 11, 2012 at 19:44
  • $\begingroup$ You probably mean will NOT work with more than one random intercept. But why not use the most simple model? $\endgroup$
    – Henrik
    Jun 12, 2012 at 7:32

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.