问题

使用fixed = TRUE来指定phi的值时，如何为每个受试者设置固定值（例如，Subject 1的值为0.7，Subject 2的值为0.5等）？

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

希望这对你有帮助。

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

答案1

得分: 7

我可以在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() 处理为随机效应），第一个参数是对数标准差，第二个（将保持恒定）是自相关参数
• start 指定起始值（特别是对于其中 map 元素设置为 NA 的任何参数的固定值）

glmmTMB 和 gls 模型的结果看起来足够接近（固定效应参数相同，协方差矩阵在2%内，残差标准差在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 向量，这是对应于每个Subject参数的20对（log-SD，trans(phi)）。
每个map 对中的第一个元素应该是1（以便所有主体的对数标准差参数都相等），每对中的第二个元素应该是 NA（以固定值）；
每个start 对中的第一个元素设置为0（对数标准差的默认起始值），每对中的第二个元素是 trans(phi_i)。
我将固定的 phi 向量设为 (0.4, 0.5, 0.6, ... 0.9, 0.4, ...) 以作示范。

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

## 用以获取交替的 (0, phi_i) 对的技巧
startvec <- c(rbind(0, trans(phivec)))
## 将所有对数标准差值设为相等，将所有 phi 值固定为起始值
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

<details>
<summary>英文:</summary>

With some work I can get this done in `glmmTMB`. I don&#39;t know how to do it in `nlme`, and I doubt it&#39;s possible (but don&#39;t know for sure).

Start by fitting the original model (a single fixed value of `phi`):

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

Warming up by doing this in glmmTMB and comparing:

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)

What you need to know:

• the ar1(0+factor(index)|Subject) term in the formula is a subject-specific AR1 term (see the relevant glmmTMB vignette for a detailed explanation of how this works ["Construction of structured covariance matrices" section], and generally for many of the technical details underlying this answer)
• trans() is the transformation from the correlation scale to the internal parameterization of correlation coefficients in glmmTMB
• dispformula = ~ 0 specifies that we have no additional residual variance term (this would be a nugget effect if we added it)
• map specifies which parameters to hold constant at their starting values, or set equal to each other: theta is the vector of random effect parameters (glmmTMB handles ar1() as a random effect), the first parameter is a log-SD, the second (to be held constant) is the autocorrelation parameter
• start specifies starting values (and, in particular, the fixed values for any parameters with their map element set to NA)

The results of the glmmTMB and gls models seem close enough (fixed effect parameters the same, covariance matrix within 2%, residual SD within 1%)

Now how do we hack this to get subject-specific fixed phi values? We'll generate 20 different ar1() terms, each of which applies only to a single subject, by multiplying each term by a 0/1 dummy (indicator) variable:

First fit the model with subject-specific terms, but with both SD and phi estimated and free to vary among subjects:

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)

Partial results:

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)

These seem plausible (SDs varying around 11.3, phi varying around 0.6).

Now little bit of hacking to set up the map and start vectors, which are 20 pairs of (log-SD, trans(phi)) corresponding to each Subject's parameters. The first element in each map pair should be 1 (so that the log-SD parameters for all subjects are set equal), and the second in each pair should be NA (to fix the value); the first element in each start pair is set to 0 (the default starting value for log-SD) and the second element in each start pair is trans(phi_i). I set the fixed phi vector to (0.4, 0.5, 0.6, ... 0.9, 0.4, ...) for illustration.

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

Now refit the model with this specification:

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

The results seem to be what we expected.

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

答案2

得分: 0

我不是R语言的专家，但我相信您可以使用corFixed参数和corClasses函数来实现您想要的效果，让我用一个例子来解释一下。

首先，假设您想要定义一个自定义的相关性结构函数，名为custom_corAR1，这个函数以一个固定值的向量作为参数，在函数内部，您可以使用AR1相关性结构创建corStruct，并使用value参数设置固定值。

要使用corClasses指定相关性结构，您可以将自定义的corAR1函数与固定值向量一起传递给corAR1参数，此外，您可以通过将corFixed设置为"corAR1"来指示相关性结构是固定的！

请查看以下代码：

library(nlme)

fixed_values <- c(0.7, 0.5) # 用所需的每个主体的固定值替换这些值

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

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