nlme：CSH共分散モデルを使用して混合モデルに適合

Question

nlmeパッケージを使用して、Rで繰り返し測定（MMRM）モデルを持つ混合モデルを適合させようとしています。

データ構造は次のとおりです。 各患者は3つのグループ（grp）のうちの1つに属し、治療グループ（trt）に割り当てられます。 患者アウトカム（y）は、6回の訪問（訪問）中に測定される。

（SASのPROC MIXEDのCSHタイプのように、異なる訪問で異種分散を持つ複合対称モデルを使用したい、https://support.sas.com/documentation/cdl/en/statug/63347/HTML/default/viewer.htm#statug_mixed_sect020.htm）。

これを行うには、相関構造をCS（corCompSymm）と重みパラメータに設定するためにlmeの相関パラメータを使用して、分散が訪問の関数であるようにしました。

また、corCompSymm自体のフォームパラメータへの訪問を追加しようとしました。

問題は、私が持っている

：私がLMEに電話で重みパラメータを設定するか否かに関係なく、同じ結果を得ることが表示されます（つまり、私がCSHのではなく、CSモデルを取得していますことが表示されますモデル）。

以下のコードを実行すると、重みパラメータが無視されていることを示唆するモデルが使用されているにもかかわらず、モデルパラメータ推定値の共分散行列の対角が同一であることがわかります。

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

質問： はどのように私は別の訪問のために異なる分散パラメータに合うようにnlmeを得るのですか？

Answer 1

いくつかのデッドエンドが終わった後、正しいパラメータがvarIdentの呼び出しで設定されていることを確認しているようです。

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) 

それは同じに見えますが、varIdentに渡されたパラメータが明示的に「形」として同定されたことに気づく：

それを行うための正しい方法があるように思われます。これが他の方法で解釈された場合、私はクラッシュが予想されましたが、私は間違っていました。