nlme: ajustement du modèle mixte en utilisant le modèle de covariance CSH

Question

J'essaie d'ajuster un modèle mixte avec des mesures répétées (MMRM) dans R en utilisant le paquetage nlme. La structure des données est la suivante: Chaque patient appartient à l'un des trois groupes (grp) et est affecté à un groupe de traitement (trt). Les résultats des patients (y) sont mesurés au cours de 6 visites (visite). Je veux utiliser un modèle de symétrie composé avec variance hétérogène sur les différentes visites (comme dans le type CSH pour PROC MIXED de SAS, https://support.sas.com/documentation/cdl/en/statug/63347/HTML/default/viewer.htm#statug_mixed_sect020.htm). Pour ce faire, j'ai utilisé le paramètre de corrélation dans lme pour définir la structure de corrélation à CS (corCompSymm) et le paramètre de poids de sorte que la variance est une fonction de visite.

J'ai aussi essayé d'ajouter la visite au paramètre form de corCompSymm lui-même. Le problème que j'ai: il semble que j'obtiens le même résultat que je mette ou non le paramètre poids dans l'appel à lme (en d'autres termes, il semble que je reçois le modèle CS plutôt que le CSH modèle).

En exécutant le code ci-dessous, vous remarquerez que la diagonale de la matrice de covariance des estimations des paramètres du modèle est identique quel que soit le modèle utilisé, ce qui suggère que le paramètre poids est ignoré.

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)) 

La question: Comment puis-je nlme pour s'adapter à différents paramètres de la variance pour les différentes visites?

Answer 1

Après quelques impasses, il semble que le problème est de s'assurer que le bon paramètre est défini dans l'appel de varIdent.

La bonne façon de le faire semble être:

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) 

Il semble même, mais il faut noter que le paramètre passé à varIdent a été explicitement identifié comme « forme ». Je m'attendais à un crash si c'était interprété autrement, mais j'avais tort.