R: using factor variables in an nlme function

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

Here is a simple model with one fixed effect parameter. This model fits well, but when I wanted to introduce a factor level covariance (i.e. age), I ran into the following error.

 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 

This is apparently a syntax error, but I'm not sure how to fix it.

Answer 1

Firstly, your model specification is incorrect: since you define fixed effects as RO in fixed = list(R0 ~ 1 + (age2)) , you should use this definition in the model definition.

Then, the installation instructions for the model will look like this:

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

Now this results in a new error message:

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

Please note that nlme has a verbose argument (in our case, it is not so informative).

But it seems that this error occurs when there is no convergence. In this case, this is due to your initial values ​​that no longer meet this model specification.

I just tried a different set of values, for example:

 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) 

that one converges and provides a model

 > 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