R:在 nlme 函数中使用因子变量

Question

library(nlme)
model <- nlme(height ~ (R0) + 1,
              data = Loblolly,
              fixed = list(R0 ~ 1),
              random = list(Seed = pdDiag(list(R0 ~ 1))),
              start = list(fixed = c(R0 = -8.5)))

这是一个只有 1 个固定效应参数的简单模型。这个模型很合适，但是当我想引入一个因子水平协变量(即年龄)时，我遇到了以下错误。

Loblolly$age2 <- as.factor(ifelse(Loblolly$age < 12.5, 0, 1))
model2 <- nlme(height ~ (R0 + age2) + 1,
              data = Loblolly,
              fixed = list(R0 ~ 1 + (age2)),
              random = list(Seed = pdDiag(list(R0 ~ 1))),
              start = list(fixed = c(R0 = -8.5, age2 = 1)))

Error in chol.default((value + t(value))/2) : 
  the leading minor of order 1 is not positive definite
In addition: Warning messages:
1: In Ops.factor(R0, age2) : ‘+’ not meaningful for factors
2: In Ops.factor(R0, age2) : ‘+’ not meaningful for factors
3: In Ops.factor(R0, age2) : ‘+’ not meaningful for factors

这似乎是一个语法错误，但我不确定如何修复它。

Answer 1

首先，您的模型规范不正确:因为您将固定效应定义为 RO在 fixed = list(R0 ~ 1 + (age2)) ，您必须在模型定义中使用此定义。

模型拟合指令则变为:

model2 <- nlme(height ~ (R0) + 1,
          data = Loblolly,
          fixed = list(R0 ~ 1 + (age2)),
          random = list(Seed = pdDiag(list(R0 ~ 1))),
          start = list(fixed = c(R0 = -8.5, age2 = 1)))

现在这会导致一条新的错误消息:

Error in nlme.formula(height ~ (R0) + 1, data = Loblolly, fixed = list(R0 ~  : 
  step halving factor reduced below minimum in PNLS step

请注意 nlme有一个 verbose论点(在我们的案例中没有那么丰富)。

但似乎在没有收敛时会发生此错误。
在这种情况下，这是由于您的起始值不再适合此模型规范。

我只是尝试了一组不同的值，例如:

model2 <- nlme(height ~ (R0) + 1,
               data = Loblolly,
               fixed = list(R0 ~ 1 + (age2)),
               random = list(Seed = pdDiag(list(R0 ~ 1))),
               start = list(fixed = c(R0 = 0, age2 = 30)), verbose=TRUE)

一个收敛并提供一个模型

> model2
Nonlinear mixed-effects model fit by maximum likelihood
  Model: height ~ (R0) + 1 
  Data: Loblolly 
  Log-likelihood: -305.1093
  Fixed: list(R0 ~ 1 + (age2)) 
R0.(Intercept)       R0.age21 
      12.96167       36.80548 

Random effects:
 Formula: R0 ~ 1 | Seed
        R0.(Intercept) Residual
StdDev:   0.0002761926 9.145988

Number of Observations: 84
Number of Groups: 14