nlme: fit mixed model using CSH covariance model

Question

I am trying to fit a mixed model with repeated measurements (MMRM) model in R using the nlme package.

The structure of the data is as follows: Each patient belongs to one of three groups (grp) and is assigned to a treatment group (trt). Patients outcomes (y) are measured during 6 visits (visit).

I want to use a compound symmetry model with heterogeneous variance over the different visits (as in the CSH type for SAS's PROC MIXED, https://support.sas.com/documentation/cdl/en/statug/63347/HTML/default/viewer.htm#statug_mixed_sect020.htm).

To do that I used the correlation parameter in lme to set the correlation structure to CS (corCompSymm) and the weights parameter so the variance is a function of visit.

I have also tried adding visit to the form parameter of corCompSymm itself.

The problem I have: it appears that I get the same results whether or not I set the weights parameter in the call to lme (in other words, it appears that I am getting the CS model rather than the CSH model).

Executing the code below, you will notice that the diagonal of the covariance matrix of the model parameter estimates are identical no matter what model is used which suggests the weight parameter is being ignored.

remove(list = objects())
library(nlme)

set.seed(55)

npatients     = 200; 
nvisits       = 6;

#---
# Generate some data:
subject_table = data.frame(subject = sprintf("S%03d", 1:npatients),
                           trt     = sample(x = c("P", "D"),       replace = T, size = npatients),
                           grp     = sample(x = c("A", "B", "C"),  replace = T, size = npatients))
subject_table = merge(subject_table, 
                      data.frame(visit.number = 1:6))
subject_table = transform(subject_table, 
                          visit = sprintf("V%02d", visit.number),
                          y     = rnorm(nrow(subject_table), mean = 0, sd = visit.number^2))
subject_table = transform(subject_table, 
                          visit   = factor(visit),
                          subject = factor(subject, ordered = T, levels =     sort(unique(as.character(subject)))),
                          grp     = factor(grp),
                          trt     = factor(trt))
#---
# Fit MMRM model to data using nlme
cs_model       = lme(y ~ trt*visit*grp,                              # fixed     effects 
                     random      = ~1|subject,                       # random effects 
                     data        = subject_table,                    # data
                     correlation = corCompSymm(form=~1|subject))     # CS correlation matrix within patient

csh_model_v1   = lme(y ~ trt*visit*grp,                              # fixed effects 
                     random      = ~1|subject,                       # random effects 
                     data        = subject_table,                    # data
                     weights     = varIdent(~1|visit),               # different "weight" within each visit (I think)
                     correlation = corCompSymm(form=~1|subject))     # CS correlation matrix within patient

csh_model_v2   = lme(y ~ trt*visit*grp,                              # fixed effects 
                     random      = ~1|subject,                       # random effects 
                     data        = subject_table,                    # data
                     weights     = varIdent(~visit|subject),         # different "weight" within each visit (I think)
                     correlation = corCompSymm(form=~1|subject))     # CS correlation matrix within patient

csh_model_v3   = lme(y ~ trt*visit*grp,                              # fixed effects 
                     random      = ~1|subject,                       # random effects 
                     data        = subject_table,                    # data
                     correlation = corCompSymm(form=~visit|subject)) # CS correlation matrix within patient

diag(vcov(cs_model))
diag(vcov(csh_model_v1))
diag(vcov(csh_model_v2))
diag(vcov(csh_model_v3))

The question: How do I get nlme to fit different variance parameters for the different visits?

Answer 1

After a few dead ends, it appears that the problem is making sure that the correct parameter is being set in the call to varIdent.

The correct way to do it appears to be:

csh_model_right = lme(y ~ trt*visit*grp,                          # fixed effects 
                  random      = ~1|subject,                   # random effects 
                  data        = subject_table,                # data
                  weights     = varIdent(form=~1|visit),      # different "weight" within each visit (I know)
                  correlation = corCompSymm(),                # CS correlation matrix within subject per random statement above
                  control     = lme.control) 

It looks the same, but notice that the parameter passed to varIdent was explicitly identified as "form". I had expected a crash if this was interpreted any other way, but I was wrong.