Визуализируйте многоуровневую модель роста с помощью nlme/ggplot2 и lme4/ggplot2

Я безуспешно пытаюсь визуализировать результаты объекта nlme. Когда я делаю это с объектом lmer, создается правильный график. Моя цель — использовать nlme и визуализировать подогнанную кривую роста для каждого человека с помощью ggplot2. Функция predict() работает по-разному с объектами nlme и lmer.

модель:

#AR1 with REML
autoregressive <- lme(NPI ~ time,
                  data = data,
                  random = ~time|patient,
                  method = "REML",
                  na.action = "na.omit",
                  control = list(maxlter=5000, opt="optim"),
                  correlation = corAR1())

nlme попытка визуализации:

data <- na.omit(data)

data$patient <- factor(data$patient,
                   levels = 1:23)

ggplot(data, aes(x=time, y=NPI, colour=factor(patient))) +
    geom_point(size=1) +
    #facet_wrap(~patient) +
    geom_line(aes(y = predict(autoregressive,
                              level = 1)), size = 1)

неправильная визуализация

когда я использую:

data$fit<-fitted(autoregressive, level = 1) 
geom_line(aes(y = fitted(autoregressive), group = patient))

он возвращает одинаковые подогнанные значения для каждого человека, поэтому ggplot создает одинаковую кривую роста для каждого. Выполнение test <-data.frame(ranef(autoregressive, level=1)) возвращает различные точки пересечения и наклоны в зависимости от идентификатора пациента. Интересно, что когда я сопоставляю модель с lmer и запускаю приведенный ниже код, он возвращает правильный график. Почему predict() работает по-разному с объектами nlme и lmer?

timeREML <- lmer(NPI ~ time + (time | patient), 
                 data = data,
                 REML=T, na.action=na.omit)

ggplot(data, aes(x = time, y = NPI, colour = factor(patient))) +
    geom_point(size=3) +
    #facet_wrap(~patient) +
    geom_line(aes(y = predict(timeREML)))

правильный сюжет


r ggplot2 nlme
person Santi Allende    schedule 13.04.2017    source источник
comment
Под визуализацией случайных эффектов, оцениваемых моделью, вы имеете в виду построение подогнанной кривой роста для каждого человека? Я думаю, вы могли бы изменить geom_line(aes(y = fitted(autoregressive), group = id)) или начать с добавления data$fit<-fitted(autoregressive)   -  person Niek    schedule 13.04.2017
comment
Спасибо за ответ @Niek. Я пытался использовать fitted(), но он возвращает одинаковые подходящие значения для каждого человека. Я обновил свой вопрос выше. Спасибо!   -  person Santi Allende    schedule 13.04.2017
comment
У вас есть воспроизводимый пример? Без ваших данных я не могу посмотреть. Попробуйте воспроизвести проблему со случайно сгенерированными данными или общедоступным набором данных.   -  person Bishops_Guest    schedule 13.04.2017

Ответы (1)


arrow_upward
0
arrow_downward

При создании воспроизводимого примера я обнаружил, что ошибка возникает не в модели predict() и не в ggplot(), а в модели lme.

Данные:

###libraries
library(nlme)
library(tidyr)
library(ggplot2)

###example data
df <- data.frame(replicate(78, sample(seq(from = 0, 
            to = 100, by = 2), size = 25, 
            replace = F)))

##add id
df$id <- 1:nrow(df)

##rearrange cols
df <- df[c(79, 1:78)]

##sort columns
df[,2:79] <- lapply(df[,2:79], sort)

##long format
df <- gather(df, time, value, 2:79)

##convert time to numeric
df$time <- factor(df$time)
df$time <- as.numeric(df$time)

##order by id, time, value
df <- df[order(df$id, df$time),]

##order value
df$value <- sort(df$value)

Модель 1 без значений NA подходит успешно.

###model1
model1 <- lme(value ~ time,
                  data = df,
                  random = ~time|id,
                  method = "ML",
                  na.action = "na.omit",
                  control = list(maxlter=5000, opt="optim"),
                  correlation = corAR1(0, form=~time|id,
                                       fixed=F))

Введение NA вызывает ошибку матрицы обратимых коэффициентов в модели 1.

###model 1 with one NA value
df[3,3] <- NA

model1 <- lme(value ~ time,
                  data = df,
                  random = ~time|id,
                  method = "ML",
                  na.action = "na.omit",
                  control = list(maxlter=2000, opt="optim"),
                  correlation = corAR1(0, form=~time|id,
                                       fixed=F))

Но не в модели 2, которая имеет более упрощенную структуру корреляции AR(1) внутри группы.

###but not in model2
model2 <- lme(value ~ time,
                  data = df,
                  random = ~time|id,
                  method = "ML",
                  na.action = "na.omit",
                  control = list(maxlter=2000, opt="optim"),
                  correlation = corAR1(0, form = ~1 | id))

Однако изменение opt="optim" на opt="nlminb" успешно подходит для модели 1.

###however changing the opt to "nlminb", model 1 runs 
model3 <- lme(value ~ time,
          data = df,
          random = ~time|id,
          method = "ML",
          na.action = "na.omit",
          control = list(maxlter=2000, opt="nlminb"),
          correlation = corAR1(0, form=~time|id,
                               fixed=F))

Приведенный ниже код успешно визуализирует модель 3 (ранее модель 1).

df <- na.omit(df)

ggplot(df, aes(x=time, y=value)) +
    geom_point(aes(colour = factor(id))) +
    #facet_wrap(~id) +
    geom_line(aes(y = predict(model3, level = 0)), size = 1.3, colour = "black") +
    geom_line(aes(y = predict(model3, level=1, group=id), colour = factor(id)), size = 1)

Обратите внимание, что я не совсем уверен, что делает изменение оптимизатора с "optim" на "nlminb" и почему это работает.

person Santi Allende    schedule 15.04.2017