0

I am making a coefficient plot with multiple models using coefplot::multiplot in R; the below image is what I currently have. enter image description here

And this is the code I used to create it:

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

I would like to have the coefficients colour coded by variable (eg, Temperature) rather than by model, but don't know how to. Any advice on how to do this appreciated.

agorapotatoes
  • 311
  • 2
  • 16

1 Answers1

0

In case this is useful to anyone, I have created the figure I wanted using ggplot using the below code; it is probably not the most efficient way of doing it but it works. I have included another term submodel that allows you to separate into two types of model (circles and triangles on the below picture).

First create a function that pulls the desired information from a model and stores it in a dataframe:

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

Then send all the models you wish to use to this function and bind them together:

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

At this point you should have a dataframe that contains the model coefficients and everything else you need, including the lb and ub which are upper and lower bounds that you use to make errorbars. Now make the plot using 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()

This produces the following plot:

enter image description here

agorapotatoes
  • 311
  • 2
  • 16