Fishing

rm(list = ls())
library(brms)
library(dplyr)
library(ggplot2)
library(grid)
library(gridExtra)
library(bayesplot)
zinb <- read.csv("http://stats.idre.ucla.edu/stat/data/fish.csv")
zinb$camper <- factor(zinb$camper, labels = c("no", "yes"))
head(zinb)
skimr::skim(zinb)

ggplot(zinb, 
       aes(x=xb, y=zg, color=as.factor(livebait)))+
  geom_point()+
  facet_wrap(~camper)+
  labs(color="LiveBait")+
  theme_bw()

ggplot(zinb, 
       aes(x=xb, y=zg, 
           size=count,
           color=as.factor(livebait)))+
  geom_point()+
  facet_wrap(~camper)+
  labs(color="LiveBait")+
  theme_bw()

ggplot(zinb, 
       aes(x=count, 
           fill=camper))+
  geom_bar(position="dodge")+
  theme_bw()

ggplot(zinb, 
       aes(x=count, 
           fill=camper))+
  geom_bar(position="dodge")+
  theme_bw() + facet_wrap(~persons, 
                          scales="free")

ggplot(zinb, 
       aes(y=count, 
           fill=as.factor(child)))+
  geom_boxplot(alpha=0.5)+
  ylim(0,50)+
  theme_bw()+facet_wrap(~nofish)

get_prior(count ~ persons + child + camper, 
          data = zinb,
          family = zero_inflated_poisson("log"))
summary(zinb$count)
fit_zinb1 <- brm(count ~ persons + 
                   child + 
                   camper, 
                 data = zinb,
                 prior=c(
                   prior(normal(34, 5),  class = Intercept),
                   prior(normal(40, 3),  class = b, coef=camperyes),
                   prior(normal(20, 4),  class = b, coef=child),
                   prior(normal(20, 8),  class = b, coef=persons),
                   prior(gamma(5, 5),   class = zi)),
                 family = zero_inflated_poisson("log"), 
                 chains=4, cores=4)
fit_zinb1$prior
summary(fit_zinb1)
plot(fit_zinb1)
pp_check(fit_zinb1, ndraws=100)
pp_check(fit_zinb1, type = "rootogram")
pp_check(fit_zinb1, type = "error_hist", ndraws = 11)
pp_check(fit_zinb1, type = "scatter_avg", ndraws = 100)
pp_check(fit_zinb1, type = "stat_2d")
pp_check(fit_zinb1, type = "loo_pit")

y<-zinb$count
yrep<-posterior_predict(fit_zinb1 )
ppc_boxplot(y, yrep[1:8, ])
ppc_dens(y, yrep[200:202, ])
ppc_freqpoly(y, yrep[1:3,], alpha = 0.1, size = 1, binwidth = 5)

group<-zinb$camper
ppc_dens_overlay(y, yrep[1:25,])
ppc_ecdf_overlay(y, yrep[sample(nrow(yrep), 25), ])
ppc_hist(y, yrep[1:8, ])
ppc_boxplot(y, yrep[1:8, ])
ppc_dens(y, yrep[200:202, ])
ppc_freqpoly(y, yrep[1:3,], alpha = 0.1, size = 1, binwidth = 5)
ppc_freqpoly_grouped(y, yrep[1:3,], group) + yaxis_text()
ppc_dens_overlay_grouped(y, yrep[1:25, ], group = group)
ppc_ecdf_overlay_grouped(y, yrep[1:25, ], group = group)
ppc_violin_grouped(y, yrep, group, size = 1.5)
ppc_violin_grouped(y, yrep, group, alpha = 0, y_draw = "points", y_size = 1.5)
ppc_violin_grouped(y, yrep, group, alpha = 0, y_draw = "both",
                   y_size = 1.5, y_alpha = 0.5, y_jitter = 0.33)

ppc_stat(y, yrep)
q25 <- function(y) quantile(y, 0.25)
ppc_stat(y, yrep, stat = "q25")
ppc_stat(y, yrep, stat = function(y) quantile(y, 0.25))
ppc_stat_grouped(y, yrep, group)
ppc_stat_grouped(y, yrep, group) + yaxis_text()

bayesplot_theme_set(ggplot2::theme_linedraw())
color_scheme_set("viridisE")
ppc_stat_2d(y, yrep, stat = c("mean", "sd"))

bayesplot_theme_set(ggplot2::theme_grey())
color_scheme_set("brewer-Paired")
ppc_stat_2d(y, yrep, stat = c("median", "mad"))


theme_set(theme_bw())
color_scheme_set("brightblue")
ppc_intervals(y, yrep)
ppc_intervals(y, yrep, size = 1.5, fatten = 0)
ppc_ribbon(y, yrep)
ppc_ribbon(y, yrep, y_draw = "points")
ppc_ribbon(y, yrep, y_draw = "both")

ppc_error_hist(y, yrep[1:3, ])
ppc_error_hist_grouped(y, yrep[1:3, ], group)
ppc_error_scatter(y, yrep[10:14, ])
ppc_error_scatter_avg(y, yrep)
ppc_scatter_avg_grouped(y, yrep, group, facet_args = list(scales = "free_x"))

ppc_bars(y, yrep)
ppc_bars_grouped(y, yrep, group, prob = 0.5, freq = FALSE)
ppc_rootogram(y, yrep)
ppc_rootogram(y, yrep, prob = 0)

loo1 <- loo(fit_zinb1, save_psis = TRUE, cores = 4)
psis1 <- loo1$psis_object
lw <- weights(psis1)
ppc_loo_pit_overlay(y, yrep, lw = lw)
ppc_loo_pit_qq(y, yrep, lw = lw)
ppc_loo_pit_qq(y, yrep, lw = lw, compare = "normal")
keep_obs <- 1:50
ppc_loo_intervals(y, yrep, psis_object = psis1, subset = keep_obs)
ppc_loo_intervals(y, yrep, psis_object = psis1, subset = keep_obs,
                  order = "median")

ppc_ribbon_grouped(y, yrep, x = zinb$persons, group, y_draw = "both") +
  ggplot2::scale_x_continuous(breaks = pretty)

ce<-conditional_effects(fit_zinb1)
p1<-plot(ce, rug=TRUE, theme=theme_bw())
p1$persons
grid.arrange(p1$persons,
             p1$child, 
             p1$camper,
             ncol=3)

get_prior(count ~ persons + child + camper,
          zi ~ child, 
          data = zinb, family = zero_inflated_poisson())




fit_zinb2 <- brm(bf(count ~ persons + 
                      child + 
                      camper, 
                      zi ~ child), 
                 data = zinb, chains=4, cores=4)
fit_zinb2 <- brm(bf(count ~ persons + 
                      child + 
                      camper, 
                    zi ~ child),
                 prior=c(
                   prior(normal(34, 5),  class = Intercept),
                   prior(normal(40, 3),  class = b, coef=camperyes),
                   prior(normal(20, 4),  class = b, coef=child),
                   prior(normal(20, 8),  class = b, coef=persons),
                   prior(gamma(5, 5),    class = zi )),
                 family = zero_inflated_poisson("log"), 
                 data = zinb, chains=4, cores=4)
fit_zinb2$prior
summary(fit_zinb2)
ce2<-conditional_effects(fit_zinb1)
p2<-plot(ce2, rug=TRUE, theme=theme_bw())
grid.arrange(p1$persons,
             p1$child, 
             p1$camper,
             p2$persons,
             p2$child, 
             p2$camper,
             ncol=3)
LOO(fit_zinb1, fit_zinb2)