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

Question

使用

fixed = TRUE

将值指定为 phi 时，如何为每个主题设置固定值（例如，主题 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 vignette [“结构化协方差矩阵的构造”部分]，并且通常对于许多技术这个答案的详细信息）

• 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, ...)

。

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

现在用此规格改装模型：

g3 <- update(g2, start = list(theta = startvec), map = list(theta = mapvec))

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

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

Answer 2

我不是 R 方面的专家，但我相信您可以使用

corFixed

参数和

corClasses

函数来实现您正在寻找的东西，所以让我用一个例子来解释。

首先，假设您要定义一个名为

custom_corAR1

的自定义相关结构函数，该函数采用固定值向量作为参数，在函数内部，您可以使用 AR1 相关结构创建

corStruct

并设置使用

value

参数固定值。

要使用

corClasses

指定相关结构，您可以将自定义

corAR1

函数与固定值向量一起传递给

corAR1

参数，此外，您可以通过将

corFixed

设置为来指示相关结构是固定的“corAR1”！

看看这个：

library(nlme)

fixed_values <- c(0.7, 0.5) #replace this with the desired fixed values for each subject

custom_corAR1 <- function(value) {
  corStruct(AR1(form = ~ 1 | Subject),
            value = value,
            fixed = TRUE)
}

mod <- gls(rate ~ pressure,
            data = Dialyzer,
            corr = corClasses(corAR1 = custom_corAR1(fixed_values),
                              corFixed = "corAR1"))