R:在 nlme 函数中使用因子变量
R: using factor variables in nlme function
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
这似乎是一个语法错误,但我不确定如何修复它。
首先,您的模型规范不正确:当您在 fixed = list(R0 ~ 1 + (age2))
中将固定效应定义为 RO
时,您必须在模型定义中使用此定义。
模型拟合指令则变为:
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
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
这似乎是一个语法错误,但我不确定如何修复它。
首先,您的模型规范不正确:当您在 fixed = list(R0 ~ 1 + (age2))
中将固定效应定义为 RO
时,您必须在模型定义中使用此定义。
模型拟合指令则变为:
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