R: using factor variables in 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 just 1 fixed effect parameter. This model fits fine, but when I wanted to introduce a factor level covariate (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 appears to be a syntax error, but I'm not sure how to go about fixing it.

Answer 1

First, your model specification is not correct: as you define fixed effects as RO in fixed = list(R0 ~ 1 + (age2)), you must use this definition in model definition.

Model fitting instruction then becomes:

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 leads to 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

Note that nlme has a verbose argument (is not so informative in our case).

But it appears that this error happens when there is no convergence. In this case, this is due to your starting values, that are not adequate anymore with this model specification.

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

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