0

I have a mixed model with a non-linear term, so I would like to use the R package nlme instead of lme. However, switching to nlme, even without adding anything to the model, causes Rstudio and R to crash.

I have found that even generated data, which can easily be fitted using lme, causes this behaviour (on my computer).

Let's start by loading the libraries and setting up a data.frame with the grouping id and spatial coordinate x.

library(nlme)

nid <- 300
nx <- 10

data <- expand.grid(
  x = seq(nx),
  id = seq(nid)
)

Now, let's add correlated error and uncorrelated error as separate columns, as well as a random intercept value per id. The output of arima.sim requires a normalisation step.

data$ec <- c(
  replicate(
    nid,
    as.numeric(
      arima.sim(
        model = list(
          order = c(1, 0, 0),
          ar = 0.5
        ),
        n = nx
      )
    )
  )
)
data$ec <- data$ec / sd(data$ec)
data$eu <- rnorm(nid * nx)
data$random <- rep(rnorm(nid), each = nx)

Now, we can create 3 dependent variables, for 3 models. The first is a mixed model with uncorrelated (regular) error. The second includes an exponential (AR1) correlation structure. The third combines both. I am adding an intercept of 1, an sd of the random effect of 2 and an sd of the total residual error of 3.

data$y1 <- 1 + 2 * data$random + 3 * data$eu
data$y2 <- 1 + 2 * data$random + 3 * data$ec
data$y3 <- 1 + 2 * data$random + sqrt(8) * data$ec + sqrt(1) * data$eu

All of the following lme models fit without problem, giving the expected result.

l1 <- lme(
  fixed = y1 ~ 1,
  random = ~ 1 | id,
  data = data,
  method = "ML"
)
l2 <- lme(
  fixed = y2 ~ 1,
  random = ~ 1 | id,
  correlation = corExp(
    form = ~ x | id
  ),
  data = data,
  method = "ML"
)
l3 <- lme(
  fixed = y3 ~ 1,
  random = ~ 1 | id,
  correlation = corExp(
    form = ~ x | id,
    nugget = TRUE
  ),
  data = data,
  method = "ML"
)

As far as I know, the following nlme code specifies exactly the same models as above. The first runs without issues. But the ones with a correlation structure crash R / RStudio. No warning or error message is provided. Fiddling with the arguments and with nlmeControl does not help, though I do think nlmeControl could be the place to search for a solution.

nlme(
  model = y1 ~ b0,
  fixed = b0 ~ 1,
  random = b0 ~ 1,
  group = ~ id,
  data = data,
  start = list(
    fixed = fixed.effects(l1),
    random = setNames(random.effects(l1), "b0")
  ),
  method = "ML"
)
nlme(
  model = y2 ~ b0,
  fixed = b0 ~ 1,
  random = b0 ~ 1,
  group = ~ id,
  correlation = corExp(
    form = ~ x
  ),
  data = data,
  start = list(
    fixed = fixed.effects(l2),
    random = setNames(random.effects(l2), "b0")
  ),
  method = "ML"
)
nlme(
  model = y3 ~ b0,
  fixed = b0 ~ 1,
  random = b0 ~ 1,
  group = ~ id,
  correlation = corExp(
    form = ~ x,
    nugget = TRUE
  ),
  data = data,
  start = list(
    fixed = fixed.effects(l3),
    random = setNames(random.effects(l3), "b0")
  ),
  method = "ML"
)

Has anyone experienced this before? Does my example code give the same problem on your computer? What are good strategies to change nlmeControl to attempt to remedy this?

LBogaardt
  • 402
  • 1
  • 4
  • 25
  • I can reproduce the crash on my Win10 machine with R 4.0.5 and nlme 3.1-152. Using `data <- as.data.frame(data)` I can avoid the crash on the model with `y2~b0`. – Marco Sandri Sep 13 '21 at 11:46
  • I can also reproduce with Ubuntu and a much older R version. Such a crash (segfault?) is always a bug. You should report it. – Roland Sep 13 '21 at 12:08
  • Unfortunately, `nlme` is no longer maintained. Is there any way to see if it produces a warning message right before it crashes? – LBogaardt Sep 14 '21 at 09:15
  • 1
    Of course it is still maintained. nlme is a recommended package, which means it is maintained by R Core (the devs of R). – Roland Sep 15 '21 at 06:03

0 Answers0