модель смешанного эффекта с взаимодействием в фиксированных эффектах в ggplot - PullRequest
Новая
       90

модель смешанного эффекта с взаимодействием в фиксированных эффектах в ggplot

0 голосов
/ 29 января 2020

Я строю взаимодействие фиксированных эффектов в модели смешанных эффектов на основе объекта lmer(). Для этого я прогнозирую новые значения на основе моей модели. Это прекрасно работает, за исключением того, что из-за того, как я их генерирую, прогнозы растягиваются на весь возможный диапазон оси X. Теперь я могу ограничить предсказанные линии регрессии диапазоном их соответствующей переменной группировки, определив new.dat на основе al oop (изменить значения max и min в зависимости от переменной группировки "Variety") et c., Но - Есть ли более элегантное / простое решение для построения этого? Я что-то упускаю (я относительно новичок в R)?

Данные: 

library(datasets)
data("Oats")

# manipulate data so it resembles more my actual data
Oats <- Oats %>%
  filter((Variety == "Golden Rain"  & nitro>=0.2) | (Variety == "Marvellous" & nitro <=0.4) | (Variety == "Victory" & nitro<=0.4 & nitro>=0.2))  #%>%

Модель и график: 

mod2 <- lmer(yield ~ nitro * Variety + (1| Variety), data=Oats)

new.dat <- data.frame(nitro=seq(min(Oats$nitro),max(Oats$nitro), length.out = 48), Variety= Oats$Variety)
new.dat$pred<-predict(mod2,newdata=new.dat,re.form=~0)

ggplot(data=Oats, aes(x=nitro, y=yield, col = Variety)) +
  geom_point() +
  geom_line(data=new.dat, aes(y=pred)) +
  geom_point(data=new.dat, aes(y=pred))

resulting plot

Большое спасибо за каждый намек!

1 Ответ

1 голос
/ 29 января 2020

Вы можете получить это, рассчитав мин / макс для каждой группы, а затем вычислив последовательности по группам. Ведение с Tidyverse, поскольку ваш код уже использует его: 

library(tidyverse)
library(pairwiseCI)
#> Loading required package: MCPAN
#> Loading required package: coin
#> Loading required package: survival
library(lme4)
#> Loading required package: Matrix
#> 
#> Attaching package: 'Matrix'
#> The following objects are masked from 'package:tidyr':
#> 
#>     expand, pack, unpack

data("Oats")

 ## manipulate data so it resembles more my actual data
Oats <-
  Oats %>%
  filter((Variety == "Golden Rain"  & nitro>=0.2) | (Variety == "Marvellous" & nitro <=0.4) | (Variety == "Victory" & nitro<=0.4 & nitro>=0.2))  #%>%

mod2 <- lmer(yield ~ nitro * Variety + (1| Variety), data=Oats)
#> Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
#> unable to evaluate scaled gradient
#> Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
#> Hessian is numerically singular: parameters are not uniquely determined

## Calculate min/max by group
all_vals <-
  Oats %>%
  group_by(Variety) %>%
  summarize(min_nitro = min(nitro),
            max_nitro = max(nitro))
## Calculate sequence for each group
new.dat <-
  all_vals %>%
  group_split(Variety) %>%
  map_dfr(~ data.frame(Variety = .x$Variety, nitro = seq(.x$min_nitro, .x$max_nitro, length.out = 20)))


new.dat$pred<-predict(mod2,newdata=new.dat,re.form=~0)

ggplot(data=Oats, aes(x=nitro, y=yield, col = Variety)) +
  geom_point() +
  geom_line(data=new.dat, aes(y=pred)) +
  geom_point(data=new.dat, aes(y=pred))

img

...