How do I dodge error bars and points? Issues using the position argument

Question

I am having trouble getting my error bars and points to "dodge".

I have looked at: ggplot2 - Dodge horizontal error bars with points

But this doesn't seem to be working.

My data:

structure(list(CircleID = c(312L, 355L, 287L, 235L, 411L, 379L
), clusterID = c(3L, 1L, 3L, 2L, 4L, 5L), Centroid_Y = c(7434929.783, 
7465878.542, 7475850.883, 7494449.891, 7453518.8, 7443857.934
), Centroid_X = c(569904.2122, 427795.6068, 476436.8265, 407557.5587, 
479871.428, 541163.2806), Lat = c(-23.19319, -22.91354, -22.824849, 
-22.65447, -23.02663, -23.11351), Lon = c(123.68308, 122.2959, 
122.77037, 122.10026, 122.80355, 123.402), Area = c(19.495698, 
19.49945, 19.495405, 19.495243, 19.497044, 19.500582), AlphaDiv = c(14L, 
16L, 19L, 22L, 18L, 13L), Richness = c(50L, 49L, 56L, 68L, 50L, 
59L), BetaDiv = c(2.50877193, 2.161290323, 2.027027027, 2.090909091, 
1.777777778, 3.411214953), Div = c(3.035329541, 2.831802152, 
3.012866305, 3.011564561, 3.097181771, 3.057795467), EdibleRichness = c(23L, 
26L, 27L, 31L, 35L, 17L), EdibleDiv = c(1.840690806, 1.650296103, 
1.882220841, 1.632854736, 1.960030193, 2.024585764), EdibleShrubRichness = c(6L, 
10L, 2L, 8L, 11L, 5L), EdibleShrubDiv = c(0.501333895, 0.490388177, 
0, 0.762407485, 0.772461793, 0.131197816), EdibleTreeRichness = c(10L, 
8L, 16L, 11L, 17L, 11L), EdibleTreeDiv = c(1.277456129, 1.125851256, 
1.437726452, 0.904393668, 1.606223117, 1.784590699), EphemShrubRichness = c(27L, 
23L, 24L, 24L, 23L, 28L), EphemShrubDiv = c(2.627285386, 2.371839483, 
2.753388655, 2.42111203, 2.493962379, 2.546223571), FFDiv = c(1.8583798, 
1.1289252, 1.5282626, 1.440973, 1.1089241, 1.3398642), TSFDiv = c(0.3762979, 
1.5408334, 2.0779354, 1.7571675, 0.5433299, 2.0140312), SuccDiv = c(0.156349121, 
1.077707981, 0.956184167, 0.560543868, 0.80547797, 1.020026211
), time = c(0.657399946, 0.772208489, 0.866721107, 0.886513771, 
0.901100342, 0.955033821), MaxArea = c(17.56358868, 7.571722324, 
8.460637152, 16.03766892, 14.100645, 6.761454355), MeanArea = c(3.335023315, 
1.870935628, 1.981612187, 3.216795584, 3.187732935, 1.573368285
), MaxProp = c(90.0895607, 38.66809073, 33.7053729, 82.26452431, 
72.32196327, 34.67309004), MeanProp = c(5.399550816, 5.399550816, 
2.772138805, 4.351757239, 5.399550816, 5.399550816), OneYearRain = c(50.70000005, 
64.39999938, 57.30000091, 59.60000074, 57.09999967, 51.30000043
), OneRainAvg = c(22.19444453, 20.54722206, 22.46388875, 22.38888903, 
21.7055556, 20.09166669), TwoYearRain = c(358.5999901, 333.2999961, 
376.5000069, 337.8000021, 371.9000043, 354.6999978), TwoRainAvg = c(15.10000043, 
14.35416673, 15.16249991, 14.91250032, 15.45416655, 14.07500009
), ThreeYearRain = c(813.5999915, 706.8000026, 782.5000027, 723.3000007, 
778.8000005, 804.1999913), ThreeRainAvg = c(4.441666732, 4.983333329, 
4.308333317, 4.241666694, 4.858333309, 4.966666728), AverageRain = c(235.874, 
243.247, 221.319, 232.387, 223.683, 235.739), spinsandplain = c(0.93, 
0.86, 0.82, 0.88, 0.82, 0.81), claylake = c(0, 0, 0, 0, 0, 0), 
    NyurnmaProp = c(0, 0.32, 0, 0, 0.12, 0.07), NyukuraProp = c(0.95, 
    0.48, 0.4, 0.78, 0.72, 0.6), ManguuProp = c(0.01, 0.12, 0.09, 
    0.02, 0.14, 0.16), KunarkaProp = c(0.02, 0.03, 0.3, 0.17, 
    0.01, 0.1), No.Fires = c(4L, 24L, 21L, 15L, 20L, 19L)), row.names = c(NA, 
6L), class = "data.frame")

The models:

m1 = glm(Richness~TSFDiv+scale(NyukuraProp)+scale(spinsandplain)+scale(AverageRain),na.action=na.fail,family=poisson,data=data) 
m2 = glm(Richness~TSFDiv+scale(NyukuraProp)+scale(AverageRain),na.action=na.fail,family=poisson,data=data) 
m3 = glm(Richness~TSFDiv+scale(NyukuraProp)+scale(spinsandplain),na.action=na.fail,family=poisson,data=data) 

Create data frame of model coefficients and standard errors Function to extract what we need

ce = function(model.obj) {
  extract = summary(get(model.obj))$coefficients[ ,1:2]
  return(data.frame(extract, vars=row.names(extract), model=model.obj))
}

Run function on the models and bind into single data frame

coefs = do.call(rbind, sapply(paste0("m",1:3), ce, simplify=FALSE))
names(coefs)[2] = "se" 

Remove the intercept (you may need to change this depending on the coeffs in the model)

Coeffs = coefs[!(row.names(coefs) %in% c('m1.(Intercept)')), ]
Coeffs = Coeffs[!(row.names(Coeffs) %in% c('m2.(Intercept)')), ]
Coeffs = Coeffs[!(row.names(Coeffs) %in% c('m3.(Intercept)')), ]

GGPLOT

TickLabels <- c("Average rainfall","Nyukura","Spinifex","Time-since-fire div")

ggplot(Coeffs, aes(vars, Estimate))+
    geom_hline(yintercept=0, lty=2, lwd=1, colour="black") +
    geom_errorbar(aes(ymin=Estimate-se, ymax=Estimate+se, colour=vars),lwd=1, width=0,position = position_dodge(width=1))+
    geom_point(size=3, aes(colour=vars))+
    #facet_grid(. ~ model) +
    coord_flip()+
    theme_classic()+
    scale_color_viridis(discrete=TRUE,option="D") +           #direction=-1
    scale_x_discrete(labels= TickLabels)+
    theme(axis.title = element_text(size=18)) +    
    theme(axis.line = element_line(colour="black")) +
    theme(axis.ticks = element_line()) +
    theme(axis.text = element_text(size = 18, colour="black")) +
    theme(panel.border = element_blank()) +
    theme(axis.line = element_line(colour="black"))+
    theme(legend.position = "none")+
    labs(x="", y="Estimate",title="Species richness")+
    theme(title = element_text(size=22))

I have tried:

# position=position_dodge(width = 0.90)   
# position=position_dodge(1))
# position=position_nudge(x=-1) 
# pd <- position_dodge(width = 0.4) 

The plot

What am I doing wrong?

Answer 1

I assume that you want to add a little bit of space between the outputs from different models? One option might be to specify the vertical axis as numeric, and add a small offset to for each model manually.

For example:

library(tidyverse); library(ggplot2); library(viridis)

Coeffs %>% arrange(vars, model) %>% mutate(v_offset = case_when(model == 'm1' ~ -.1,
                                                              model == 'm2' ~  0, 
                                                              model == 'm3' ~ .1)) %>%
        mutate(my_yaxis = as.numeric(as.factor(vars)) + v_offset) -> Coeffs
        

TickLabels <- c("Average rainfall","Nyukura","Spinifex","Time-since-fire div")
TickPositions =  unique(as.numeric(as.factor(Coeffs$vars)))

Modified ggplot command

   ggplot(Coeffs, aes(my_yaxis, Estimate))+
        geom_hline(yintercept=0, lty=2, lwd=1, colour="black") +
        geom_errorbar(aes(ymin=Estimate-se, ymax=Estimate+se, colour=vars),lwd=1,    width=0,position = position_dodge(width=1))+
        geom_point(size=3, aes(colour=vars))+
        #facet_grid(. ~ model) +
        coord_flip()+
        theme_classic()+
        scale_color_viridis(discrete=TRUE,option="D") +           #direction=-1
        scale_x_continuous(breaks = TickPositions, labels= TickLabels) +
        theme(axis.title = element_text(size=18)) +    
        theme(axis.line = element_line(colour="black")) +
        theme(axis.ticks = element_line()) +
        theme(axis.text = element_text(size = 18, colour="black")) +
        theme(panel.border = element_blank()) +
        theme(axis.line = element_line(colour="black"))+
        theme(legend.position = "none")+
        labs(x="", y="Estimate",title="Species richness")+
        theme(title = element_text(size=22)) 
     

Note, that the aes now refers to the newly created numerical variable instead of vars. Also scale_x_discrete has been swapped to scale_x_numeric as the axis is no longer discrete.

These questions are related and give alternative approaches: ggplot2 - Dodge horizontal error bars with points, Vertical equivalent of position_dodge for geom_point on categorical scale