Multilevel-Wachstumsmodell mit nlme/ggplot2 vs lme4/ggplot2 visualisieren

Question

Ich versuche, die Ergebnisse eines nlme Objekts ohne Erfolg zu visualisieren. Wenn ich dies mit einem lmer Objekt mache, wird das richtige Diagramm erstellt. Mein Ziel ist es, nlme zu verwenden und eine angepasste Wachstumskurve für jedes Individuum mit ggplot2 zu visualisieren. Die predict() Funktion scheint mit nlme und lmer Objekten anders zu funktionieren.

Modell:

#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 Visualisierung Versuch:

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) 

wenn ich benutze:

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

es ret Urnen die gleichen angepassten Werte für jedes Individuum und so erzeugt ggplot für jedes die gleiche Wachstumskurve. Running test <-data.frame(ranef(autoregressive, level=1)) gibt verschiedene Abschnitte und Steigungen nach Patient ID zurück. Interessanterweise, wenn ich das Modell mit lmer anpasse und den untenstehenden Code ausführe, wird die korrekte Darstellung zurückgegeben. Warum funktioniert predict() anders mit nlme und lmer Objekten?

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

Answer 1

in reproduzierbares Beispiel erstellen, fand ich, dass der Fehler nicht in predict() vorkommende wurde noch in ggplot() sondern im lme Modell.

Daten:

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

Modell 1 ohne NA-Wert paßt erfolgreich.

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

Einführung NA Ursachen invertierbare Koeffizientenmatrix Fehler im Modell 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)) 

aber nicht in Modell 2, die eine stark vereinfach within-Gruppe AR (1) Korrelationsstruktur aufweist.

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

Jedoch ändert opt="optim" zu opt="nlminb" passt Modell 1 erfolgreich.

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

Der folgende Code visualisiert erfolgreich Modell 3 (ehemals Modell 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) 

Bitte beachte, dass ich bin nicht ganz sicher, was das Optimierungs "optim"-"nlminb" Ändern tut und warum es funktioniert.