Как раскрасить по коэффициенту, а не по модели в R coefplot :: multiplot

Question

Я делаю график коэффициентов с несколькими моделями, используя coefplot::multiplot в R; ниже изображение то, что у меня есть в настоящее время.

И вот код, который я использовал для его создания:

coefplot::multiplot(sc.mod.env.msrh, sc.mod.env.msrhmu, sc.mod.env.msrhat, sc.mod.env.msrhatmu, 
    coefficients=c("(Intercept)",'MeanSpeed', 'RH', 'MeanUpdraft', 'TKE','AirTemp'), 
    newNames=c(MeanSpeed='Horizontal Wind Speed', RH='Relative Humidity', MeanUpdraft='Vertical Wind Speed', AirTemp='Temperature'), 
    single=FALSE, 
    ncol=2, 
    names=c(sc.mod.env.msrhatmu="a) Global model w/ horizontal wind speed", sc.mod.env.tkerhatmu="b) Global model w/ TKE", sc.mod.env.msatmu="c) Global model w/ horizontal wind speed, \n RH removed", sc.mod.env.tkeatmu="d) Global model w/ TKE, \n RH removed"))+
theme_bw()+
theme(legend.position="none")+
ggtitle("")

Я бы хотел, чтобы цвета коэффициентов кодировались по переменной (например, по температуре), а не по модели, но я не знаю, как это сделать. Любые советы о том, как это сделать, приветствуются.

Answer 1

На случай, если это кому-нибудь пригодится, я создал нужную фигуру, используя ggplot, используя приведенный ниже код; это, вероятно, не самый эффективный способ сделать это, но это работает. Я включил еще один термин подмодель, который позволяет разделить модель на два типа (кружки и треугольники на рисунке ниже).

Сначала создайте функцию, которая извлекает нужную информацию из модели и сохраняет ее в кадре данных:

get_estimates_for_coefplot<- function(mod, time, modname){
test2<-data.frame(summary(mod)$coefficients)
test2['term']<-names(coef(mod))
test2['model']<-modname
test2['estimate']<-test2$Estimate
test2['submodel']<-time
test2['std.error']<-test2$Std..Error 
test2['ub']<-test2$estimate+test2$std.error
test2['lb']<-test2$estimate-test2$std.error
Newdata <- test2 %>%
filter(!grepl(".*s.*",term))
return(Newdata)
}

Затем отправьте все модели, которые вы хотите использовать, для этой функции и свяжите их вместе:

df<-get_estimates_for_coefplot(sc.mod.env.msrhatmu, 'Time of day \n not included', 'a) Horizontal wind speed, \n RH included')
df<-rbind(df, get_estimates_for_coefplot(sc.mod.env.tkerhatmu, 'Time of day \n not included', 'b) TKE, RH included'))
modcoeff<-df[c('term','estimate','model','std.error','submodel','lb','ub')]
modcoeff <- modcoeff %>% 
  relabel_predictors(c(MeanSpeed = "Horizontal \n wind speed",
                     RH = "Relative \n Humidity",
                     AirTemp = "Temperature",
                     "I(hour^2)" = "Time of Day^2",
                     hour = "Time of Day",
                     MeanUpdraft = "Vertical \n wind speed",
                     TKE = "TKE"))
modcoeff$term<- factor(modcoeff$term, levels = c("Time of Day^2","Time of Day","Vertical \n wind speed","Temperature", "Relative \n Humidity","TKE", "Horizontal \n wind speed"))
modcoeff$submodel<- factor(modcoeff$submodel, levels = c('Time of day \n not included', 'Time of day \n included'))

На данный момент у вас должен быть фрейм данных, который содержит коэффициенты модели и все остальное, что вам нужно, включая lb и ub, которые представляют собой верхнюю и нижнюю границы, которые вы используете для создания панелей ошибок. Теперь сделайте сюжет используя ggplot.

pd <- position_dodge(width=0.5)
ggplot(modcoeff, 
   aes(x=term,
       y=estimate, 
       color=term,
       group=std.error,
       ymin=lb,
       ymax=ub)) +
  geom_hline(yintercept = 0, color='darkgrey') +
  geom_point(aes(shape=submodel),size=2.7, position=pd, stat="identity") + 
  geom_errorbar(aes(width=0.2),position=pd) +
  facet_wrap(~model)+
  xlab("")+
  ylab('Standardised coefficient value')+
  theme_bw()+
  theme(legend.position = c(0.83, 0.25))+
  theme(legend.title = element_blank())+
  guides(colour=FALSE)+
  coord_flip()

Это дает следующий сюжет: