@@ -237,7 +237,7 @@ ppc_bars_grouped <-
237237# '
238238ppc_rootogram <- function (y ,
239239 yrep ,
240- style = c(" standing" , " hanging" , " suspended" ),
240+ style = c(" standing" , " hanging" , " suspended" , " discrete " ),
241241 ... ,
242242 prob = 0.9 ,
243243 size = 1 ) {
@@ -266,6 +266,44 @@ ppc_rootogram <- function(y,
266266 }
267267 tyrep <- do.call(rbind , tyrep )
268268 tyrep [is.na(tyrep )] <- 0
269+
270+ # Discrete style
271+ pred_mean <- colMeans(tyrep )
272+ pred_quantile <- t(apply(tyrep , 2 , quantile , probs = probs ))
273+ colnames(pred_quantile ) <- c(" lower" , " upper" )
274+
275+ # prepare a table for y
276+ ty <- table(y )
277+ y_count <- as.numeric(ty [match(xpos , rownames(ty ))])
278+ y_count [is.na(y_count )] <- 0
279+
280+ if (style == " discrete" ) {
281+ obs_color <- ifelse(y_count > = pred_quantile [, " lower" ] & y_count < = pred_quantile [, " upper" ], " blue" , " red" )
282+
283+ data <- data.frame (
284+ xpos = xpos ,
285+ obs = y_count ,
286+ pred_mean = pred_mean ,
287+ lower = pred_quantile [, " lower" ],
288+ upper = pred_quantile [, " upper" ],
289+ obs_color = obs_color
290+ )
291+
292+ graph <- ggplot(data , aes(x = xpos )) +
293+ geom_point(aes(y = obs , fill = " Observed" ), size = size * 3.5 , color = obs_color , shape = 18 ) +
294+ geom_pointrange(aes(y = pred_mean , ymin = lower + (pred_mean - lower )* 0.5 , ymax = upper - (upper - pred_mean )* 0.5 , color = " Expected" ), linewidth = size , size = size , fatten = 2 , alpha = 0.6 ) +
295+ geom_linerange(aes(y = pred_mean , ymin = lower , ymax = upper , color = " Expected" ), linewidth = size , size = size , alpha = 0.4 ) +
296+ scale_y_sqrt() +
297+ scale_fill_manual(" " , values = get_color(" l" )) +
298+ scale_color_manual(" " , values = get_color(" dh" )) +
299+ labs(x = expression(italic(y )), y = " Count" ) +
300+ bayesplot_theme_get() +
301+ reduce_legend_spacing(0.25 )
302+ return (graph )
303+ }
304+
305+
306+ # Standing, hanging, and suspended styles
269307 tyexp <- sqrt(colMeans(tyrep ))
270308 tyquantile <- sqrt(t(apply(tyrep , 2 , quantile , probs = probs )))
271309 colnames(tyquantile ) <- c(" tylower" , " tyupper" )
@@ -395,7 +433,7 @@ ppc_bars_data <-
395433 data <-
396434 reshape2 :: melt(tmp_data , id.vars = " group" ) %> %
397435 count(.data $ group , .data $ value , .data $ variable ) %> %
398- tidyr :: complete(.data $ group , .data $ value , .data $ variable , fill = list (n = 0 )) %> %
436+ tidyr :: complete(.data $ group , .data $ value , .data $ variable , fill = list (n = 0 )) %> %
399437 group_by(.data $ variable , .data $ group ) %> %
400438 mutate(proportion = .data $ n / sum(.data $ n )) %> %
401439 ungroup() %> %
0 commit comments