center day variable at final observation in R (for linear mixed model)

Question

I'm trying to analyze some data with repeated measurements on subjects in different treatment groups. Here's a subset of my data with observations taken on days 1, 3, and 21 (the complete dataset has additional observations on days between 3 and 21).

mydata <- data.frame(
  Subject  = c(13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 29, 30, 31, 32, 33, 
           34, 35, 36, 37, 38, 39, 40, 62, 63, 64, 65, 13, 14, 15, 16, 17, 18, 
           19, 20, 21, 22, 23, 24, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 
           40, 62, 63, 64, 65, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 
           29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 62, 63, 64, 65), 
  Day       = as.numeric(c(rep(c("1", "3", "21"), each=28))), 
  Treatment = c(rep(c("B", "A", "C", "B", "C", "A", "A", "B", "A", "C", "B", "C", 
                  "A", "A", "B", "A", "C", "B", "C", "A", "A"), each = 4)), 
  Obs       = c(6.472687, 7.017110, 6.200715, 6.613928, 6.829968, 7.387583, 7.367293, 
            8.018853, 7.527408, 6.746739, 7.296910, 6.983360, 6.816621, 6.571689, 
            5.911261, 6.954988, 7.624122, 7.669865, 7.676225, 7.263593, 7.704737, 
            7.328716, 7.295610, 5.964180, 6.880814, 6.926342, 6.926342, 7.562293, 
            6.677607, 7.023526, 6.441864, 7.020875, 7.478931, 7.495336, 7.427709, 
            7.633020, 7.382091, 7.359731, 7.285889, 7.496863, 6.632403, 6.171196, 
            6.306012, 7.253833, 7.594852, 6.915225, 7.220147, 7.298227, 7.573612, 
            7.366550, 7.560513, 7.289078, 7.287802, 7.155336, 7.394452, 7.465383, 
            6.976048, 7.222966, 6.584153, 7.013223, 7.569905, 7.459185, 7.504068, 
            7.801867, 7.598728, 7.475841, 7.511873, 7.518384, 6.618589, 5.854754, 
            6.125749, 6.962720, 7.540600, 7.379861, 7.344189, 7.362815, 7.805802, 
            7.764172, 7.789844, 7.616437, NA, NA, NA, NA))

I want to analyze my data using a linear mixed model (using nlme):

mymodel <- lme(Obs ~ Treatment * Day, random = ~1 | Subject, correlation = corAR1(form = ~1 | Subject), data=mydata, na.action=na.omit)

Since I'm most interested in differences on the final day of the experiment, I'd really like the center at Day 21, but R appears to use the the lowest value as the reference level for a numerical variable. I could just set Day 21 to 0, but that messes up the time sequence, which is important for the autocorrelation structure, and changing Day to a factor just gives me an error message:

Error in MEEM(object, conLin, control$niterEM) : Singularity in backsolve at level 0, block 1

So, how do I test the effect of Treatment with Day 21 as the reference level of Day?

Answer 1

try

mydata$Day = relevel(mydata$Day, ref="Day21")

before you run your model. (I think in your model Subject should be written instead of Sample), at least for your example, because there is no Sample in mydata.

EDIT

the original example has been edited, but the same approach should still work:

mydata$Day = factor(mydata$Day, levels=c(21, 3, 1))


mymodel <- lme(Obs ~ Treatment * Day, random = ~1 | Subject, correlation = 
                corAR1(form = ~1 | Subject), data=mydata, na.action=na.omit)
summary(mymodel)

which gives:

Fixed effects: Obs ~ Treatment * Day 
Value Std.Error DF  t-value p-value
(Intercept)      7.635586 0.1159703 46 65.84086  0.0000
TreatmentB      -0.965811 0.1682325 25 -5.74093  0.0000
TreatmentC      -0.169050 0.1682325 25 -1.00486  0.3246
Day3            -0.208276 0.1133072 46 -1.83815  0.0725
Day1            -0.452858 0.1306224 46 -3.46693  0.0012
TreatmentB:Day3  0.229415 0.1636327 46  1.40201  0.1676
TreatmentC:Day3  0.060868 0.1636327 46  0.37198  0.7116
TreatmentB:Day1  0.352958 0.1932224 46  1.82669  0.0742
TreatmentC:Day1  0.334850 0.1932224 46  1.73298  0.0898