2017-04-13 166 views
0

我正嘗試使用nlme軟件包在R中使用重複測量(MMRM)模型擬合混合模型。nlme:使用CSH協方差模型擬合混合模型

數據的結構如下: 每個患者屬於三個組(grp)之一,並被分配到一個治療組(trt)。 患者結果(y)在6次訪問(訪問)期間測量。

我想在不同訪問中使用具有異構差異的複合對稱模型(如SAS的PROC MIXED的CSH類型,https://support.sas.com/documentation/cdl/en/statug/63347/HTML/default/viewer.htm#statug_mixed_sect020.htm)。

爲此,我使用lme中的相關性參數將相關結構設置爲CS(corCompSymm)和權重參數,因此方差是訪問的函數。

我也嘗試添加訪問corCompSymm本身的窗體參數。

我有這個問題:看起來我得到了相同的結果,不管我是否在lme調用中設置了權重參數(換句話說,看起來我得到CS模型而不是CSH模型)。

執行下面的代碼,您會注意到模型參數估計的協方差矩陣的對角線無論使用什麼模型都是相同的,這表明權重參數被忽略。

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以適應不同的訪問不同的變化的參數?

回答

0

經過幾次死路,看來問題在於確保在調用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參數被明確認定爲「形式」。如果這種解釋有任何其他的方式,我預計會發生崩潰,但我錯了。