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

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

这似乎是语法错误，但是我不确定如何解决它。

最佳答案

首先，您的模型规范不正确:在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

关于R:在nlme函数中使用因子变量，我们在Stack Overflow上找到一个类似的问题：https://stackoverflow.com/questions/46795349/