使用nlme/ggplot2 vs lme4/ggplot2可视化多级增长模型

问题描述：

我试图将没有成功的nlme对象的结果可视化。当我使用lmer对象时，会创建正确的绘图。我的目标是使用nlme并通过ggplot2可视化每个人的拟合增长曲线。 predict()函数似乎与nlme和lmer对象的工作方式不同。使用nlme/ggplot2 vs lme4/ggplot2可视化多级增长模型

型号：

#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))通过患者ID返回不同的截距和斜率。有趣的是，当我将模型与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)))

通过“可视化模型估计的随机效应”是否意味着绘制每个人的拟合增长曲线？我想你可以改变'geom_line（aes（y = fitted（autoregressive），group = id））'或者开始加入'data $ fit Niek

谢谢你回复@Niek。我尝试使用'fitted（）'，但是它为每个人返回相同的拟合值。我更新了我的问题。谢谢！ –

你有一个可重复的例子吗？没有你的数据我不能看一看。尝试用随机生成的数据或公共数据集重现问题。 –

答

在创建重现的实例中，我发现错误未在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)

无NA值模型1成功地适合。在模型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

###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的原因可逆系数矩阵的错误，它具有更简单的组内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"做，为什么它的工作原理。

使用nlme/ggplot2 vs lme4/ggplot2可视化多级增长模型

相关推荐