如何在 R 中使用 nlme 设置每组的 phi？

Question

当使用 指定值作为 phi 时fixed = TRUE，如何为每个主题设置固定值（例如，主题 1 的值 = 0.7，主题 2 的值 = 0.5 等）？

library(nlme)
mod <- gls(rate ~ pressure,
           data = Dialyzer,
           corr = corAR1(form = ~ 1 | Subject, value = 0.7, fixed = TRUE))

Answer 1

通过一些工作，我可以完成这件事glmmTMB。我不知道如何在 中做到这一点nlme，并且我怀疑这是可能的（但不确定）。

首先拟合原始模型（ 的单个固定值phi）：

library(nlme)
mod <- gls(rate ~ pressure,
           data = Dialyzer,
           corr = corAR1(form = ~ 1 | Subject, value = 0.7, fixed = TRUE))

通过这样做和比较来热身glmmTMB：

library(glmmTMB)
trans <- function(rho) rho/sqrt(1-rho^2)
g1 <- glmmTMB(rate ~ pressure + ar1(0+factor(index)|Subject), 
        dispformula = ~ 0, 
        map = list(theta = factor(c(1, NA))), 
        start = list(theta = c(0, trans(0.7))),
        data = Dialyzer)

您需要了解什么：

• 公式中的术语ar1(0+factor(index)|Subject)是特定于主题的 AR1 术语（有关其工作原理的详细说明，请参阅相关的 glmmTMB 插图[“结构化协方差矩阵的构造”部分]，并且通常了解此答案背后的许多技术细节）
• trans()是从相关尺度到相关系数内部参数化的变换glmmTMB
• dispformula = ~ 0指定我们没有额外的残差方差项（如果我们添加它，这将是一个金块效应）
• map指定哪些参数在其起始值处保持恒定，或设置为彼此相等：theta是随机效应参数的向量（作为随机效应glmmTMB处理ar1()），第一个参数是 log-SD，第二个参数（保持恒定）是自相关参数
• start指定起始值（特别是其map元素设置为 的任何参数的固定值NA）

glmmTMB和模型的结果gls看起来足够接近（固定效应参数相同，协方差矩阵在2%以内，残差SD在1%以内）

现在我们如何破解它以获得特定于主题的固定phi值？我们将生成 20 个不同的 ar1()术语，每个术语仅适用于单个主题，方法是将每个术语乘以 0/1 虚拟（指示符）变量：

首先用特定于主题的术语拟合模型，但估计 SD 和 phi，并且可以在主题之间自由变化：

dummy <- lme4::dummy
ar1_terms <- sprintf("ar1(0+dummy(Subject, %d):factor(index)|Subject)",
                     seq_len(length(levels(Dialyzer$Subject))))
form <- reformulate(c("pressure", ar1_terms), response = "rate")
g2 <- glmmTMB(form, dispformula = ~ 0, data = Dialyzer)

部分结果：

Subject    dummy(Subject, 1):factor(index)1  12.940   0.72 (ar1)
Subject.1  dummy(Subject, 2):factor(index)1  11.088   0.64 (ar1)
Subject.2  dummy(Subject, 3):factor(index)1   8.193   0.56 (ar1)
Subject.3  dummy(Subject, 4):factor(index)1   8.782   0.55 (ar1)
Subject.4  dummy(Subject, 5):factor(index)1   8.258   0.54 (ar1)

这些似乎是合理的（SD 在 11.3 左右变化，phi 在 0.6 左右变化）。

现在需要进行一些设置来设置map和start向量，它们是与每个主题的参数相对应的 20 对 (log-SD, trans(phi))。每对中的第一个元素map应为 1（以便将所有受试者的 log-SD 参数设置为相等），每对中的第二个元素应为NA（以固定值）；每对中的第一个元素start设置为 0（log-SD 的默认起始值​​），每start对中的第二个元素为trans(phi_i)。为了便于说明，我将固定phi向量设置为(0.4, 0.5, 0.6, ... 0.9, 0.4, ...)。

Subject    dummy(Subject, 1):factor(index)1  12.940   0.72 (ar1)
Subject.1  dummy(Subject, 2):factor(index)1  11.088   0.64 (ar1)
Subject.2  dummy(Subject, 3):factor(index)1   8.193   0.56 (ar1)
Subject.3  dummy(Subject, 4):factor(index)1   8.782   0.55 (ar1)
Subject.4  dummy(Subject, 5):factor(index)1   8.258   0.54 (ar1)

现在用这个规格重新安装模型：

phivec <- rep((4:9)/10, length.out = 20)

## trick to get alternating (0, phi_i) pairs
startvec <- c(rbind(0, trans(phivec)))
## set all log-sd values equal, all phi values fixed to start values
mapvec <- factor(rbind(rep(1, 20), NA_integer_))

结果似乎符合我们的预期。

Groups     Name                              Std.Dev. Corr      
Subject    dummy(Subject, 1):factor(index)1  11.66    0.40 (ar1)
Subject.1  dummy(Subject, 2):factor(index)1  11.66    0.50 (ar1)
Subject.2  dummy(Subject, 3):factor(index)1  11.66    0.60 (ar1)
Subject.3  dummy(Subject, 4):factor(index)1  11.66    0.70 (ar1)
Subject.4  dummy(Subject, 5):factor(index)1  11.66    0.80 (ar1)
Subject.5  dummy(Subject, 6):factor(index)1  11.66    0.90 (ar1)
...