5  Poisson modelleren

# Zoals altijd: de pakket moeten wel geïnstalleerd zijn.
library(tidyverse)
Warning: package 'tidyverse' was built under R version 4.1.3
-- Attaching packages --------------------------------------- tidyverse 1.3.1 --
v ggplot2 3.3.6     v purrr   0.3.4
v tibble  3.1.7     v dplyr   1.0.9
v tidyr   1.2.0     v stringr 1.4.0
v readr   2.1.2     v forcats 0.5.1
Warning: package 'ggplot2' was built under R version 4.1.3
Warning: package 'tibble' was built under R version 4.1.3
Warning: package 'tidyr' was built under R version 4.1.3
Warning: package 'readr' was built under R version 4.1.3
Warning: package 'dplyr' was built under R version 4.1.3
-- Conflicts ------------------------------------------ tidyverse_conflicts() --
x dplyr::filter() masks stats::filter()
x dplyr::lag()    masks stats::lag()
library(brms)
Warning: package 'brms' was built under R version 4.1.3
Loading required package: Rcpp
Warning: package 'Rcpp' was built under R version 4.1.3
Loading 'brms' package (version 2.17.0). Useful instructions
can be found by typing help('brms'). A more detailed introduction
to the package is available through vignette('brms_overview').

Attaching package: 'brms'
The following object is masked from 'package:stats':

    ar
library(bayesplot)
Warning: package 'bayesplot' was built under R version 4.1.3
This is bayesplot version 1.9.0
- Online documentation and vignettes at mc-stan.org/bayesplot
- bayesplot theme set to bayesplot::theme_default()
   * Does _not_ affect other ggplot2 plots
   * See ?bayesplot_theme_set for details on theme setting
library(ProbBayes)
Warning: package 'ProbBayes' was built under R version 4.1.3
Loading required package: LearnBayes

Attaching package: 'LearnBayes'
The following object is masked from 'package:brms':

    rdirichlet
Loading required package: gridExtra

Attaching package: 'gridExtra'
The following object is masked from 'package:dplyr':

    combine
Loading required package: shiny
Warning: package 'shiny' was built under R version 4.1.3

5.1 Poisson Modelleren

Hier observeren we tellingen \(y_1,...,y_n\) verdeeld volgens een Poisson distributie met gemiddelde \(\lambda\).

Schrijf het model in termen van het logaritme van het gemiddelde:

\[\theta=log\lambda\] is het model met \(N(\mu,\sigma)\) prior op de log gemiddelde parameter \(\theta\).

5.2 Leren over website tellingen

In het ProbBayes pakket, omvat de variabele Count in de dataset web_visits tellingen van dagelijkse bezoeken aan een blog_website. We zijn geïnteresseerd in leren van het gemiddelde aantal bezoeken \(\lambda\).

We plaatsen een \(N(0, 10)\) prior op \(\theta=log\lambda\) die de beperkte informatie van de prior over de lokatie van deze parameter uitdrukt.

5.2.1 Bayesiaans fitten

In deze run van de brm() functie gaan we uit van Poisson sampling en een normaal prior met gemiddelde 0 en standaard deviatie 10 geplaatst op het loggemiddelde \(\theta=log\lambda\).

fit <- brm(Count ~ 0 + Intercept, 
            data = web_visits,
            family = poisson,
             refresh = 0,
           prior = prior(normal(0, 10), 
                         class = b,
                         coef = "Intercept"))
Compiling Stan program...
Start sampling

We bevestigen de prior met de prior_summary() functie.

prior_summary(fit)
         prior class      coef group resp dpar nlpar lb ub  source
        (flat)     b                                       default
 normal(0, 10)     b Intercept                                user

De summary() functie geeft de samenvattingen van de posterior van \(\theta\).

summary(fit)
 Family: poisson 
  Links: mu = log 
Formula: Count ~ 0 + Intercept 
   Data: web_visits (Number of observations: 28) 
  Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
         total post-warmup draws = 4000

Population-Level Effects: 
          Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept     4.58      0.02     4.54     4.62 1.00     1431     1936

Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

De posterior_samples() functie geeft de posterior simulaties van \(\theta\).

post <- posterior_samples(fit)
Warning: Method 'posterior_samples' is deprecated. Please see ?as_draws for
recommended alternatives.
head(post)
  b_Intercept    lprior      lp__
1    4.615509 -3.328038 -136.9109
2    4.600793 -3.327360 -135.9627
3    4.608702 -3.327724 -136.3967
4    4.601794 -3.327406 -136.0080
5    4.567200 -3.325820 -136.0463
6    4.562010 -3.325583 -136.3331

5.2.2 Posterior predictieve model checks

Feitelijk is dit een slecht model voor deze gegevens. Men kan dat zien door verschillende posterior predictive checks.

De pp_check() toont dichtheidsplots van 10 gerepliceerde datasets van de posterior voorspellende verdeling. Merk op dat deze gerepliceerde datasets er anders uitzien (kleinere variatie) dan de geobserveerde data.

pp_check(fit)
Using 10 posterior draws for ppc type 'dens_overlay' by default.

De pp_check() functie zal een posterior voorspellende controle implementeren met behulp van verschillende controle functies. Hier gebruiken we \((y^_,s)\) als een bivariate controlefunctie. De scatterplot geeft de waarden van \((y^_,s)\) weer. Hieruit blijkt dat de waargenomen gegevens meer variatie vertonen dan het Poisson model voorspelt.

pp_check(fit, type = "stat_2d")
Using all posterior draws for ppc type 'stat_2d' by default.