построить данные наблюдений и прогнозировать данные по двум моделям (лм и лм) на одном и том же графике

Question

Как можно отобразить данные наблюдений и результаты разных моделей (lm и lme) на одном графике? Я попробовал код ниже, но он работал только для очков. Я хотел бы добавить данные, предсказанные встроенными моделями разных цветов.

#Data
d <- runif(160,0,100)#data
y <- rnorm(16,1,0.05)*x + rnorm(16,0,0.5)#data
df = data.frame(d,y)

#Models
#linear - 1
m1 = lm(y~d, data = df)
summary(m1)
# linear - 2
m2 = lm(y~d+I(d^2), data = df)
summary(m2)

df$Class10<-with(df,ifelse(d<20,"<20",ifelse(d<30,"20-30",
 ifelse(d<40,"30-40",ifelse(d<50,"40-50",ifelse(d<60,"50-60",
 ifelse(d<70,"60-70",ifelse(d<80,"70-80",ifelse(d<90,"80-90",
 ifelse(d>=90,">90","ERROR"))))))))))
# number of classes
length(unique(df$Class10))
# classes
sort(unique(df$Class10))
# observations by class
table(df$Class10)
plot(table(df$Class10))

# b0
m10 = lme(y~d, random=~1|Class10, method="ML" ,data = df)

# b1
m10 = lme(y~d, random=~-1+d|Class10, method="ML" , data = df)

# 
m10 = lme(y~d, random=~d|Class10, method="ML" , data = df,
control = lmeControl(niterEM = 5200, msMaxIter = 5200))

#plot points - It works
plot(df$d, df$y)  
points(df$d, predict(m1), col="blue")
points(df$d, predict(m10, level=1), col="red")

#curve
plot(df$d, df$y)
curve(predict(m1,newdata=data.frame(d=x)),lwd=2, add=T)
curve(predict(m10,newdata=data.frame(d=x)),lwd=1, add=T)#error
# line
plot(df$d,df$y)  
curve(predict(m1,newdata=data.frame(d=x)),lwd=2, add=T)
lines(df$d, predict(m10, level=1),col="green")#error

Есть ли какой-нибудь способ, например, в ggplot2?

Answer 1

Вот способ! Мне нравится использовать broom и broom.mixed, чтобы получить полный тиббл с прогнозируемыми значениями для каждой модели.

library(tidyverse)
library(lme4)
library(broom)
library(broom.mixed)


df <- ChickWeight

lin <- lm(weight ~ Time,df)

mlm <- lmer(weight ~ Time + (1 | Chick),df)

df <- df %>% 
  mutate(linpred = broom::augment(lin)[,3] %>% pull(),
         mlmpred = broom.mixed::augment(mlm)[,4] %>% pull())

ggplot(df,aes(Time,weight,group = Chick)) + 
  geom_line(alpha = .2) + 
  geom_line(aes(y = linpred,color = 'Fixed Linear Effect')) + 
  geom_line(aes(y = mlmpred,color = 'Random Intercepts'), alpha = .4) +
  scale_color_manual(values = c('blue','red')) +
  labs(color = '') +
  theme_minimal()