Run a Poisson regression model with varying prior precisions on the linear effect

inla_prior_poisson.R

if (!require(INLA)){
    install.packages("INLA", 
                     repos=c(getOption("repos"),
                             INLA="https://inla.r-inla-download.org/R/testing"), 
                     dep=TRUE)
}

library(INLA)
library(tidyverse)
library(broom)

## Poisson, depends on x
set.seed(1)
n <- 50
x <- rnorm(n)
E <- runif(n)
beta <- 2
y <- rpois(n, lambda = E*exp(beta*x))

df <- data.frame(x=x, y=y, E=E)

## Formula involving Y as a matrix
formula = y ~ x - 1


fit_model <- function(x, formula, prec){
    inla(formula = formula,
         family = c("poisson"),
         control.fixed=list(prec = prec),
         data = x,
         E = E)
}



tidy.inla <- function(x){
    # https://gist.github.com/njtierney/aff88d1cb182c959b436154f48faf1e6
    # x = model_inla
    term_names <- rownames(x$summary.fixed)
    
    tibble::as_tibble(x$summary.fixed) %>%
        dplyr::mutate(terms = term_names) %>%
        dplyr::select(terms,
                      dplyr::everything())
}

scientific_10 <- function(x) {
    parse(text=gsub("e", " %*% 10^", scales::scientific_format()(x)))
}

fitted_models <- as.list(10^seq(-2,3,by=0.5)) %>%
    set_names(. , .) %>%
    map(~fit_model(df, formula, .x))


fitted_models %>%
    map_df(~tidy(.x), .id = "prec") %>%
    mutate(prec = parse_number(prec)) %>%
    ggplot(data=. , aes(x = prec, y = mean)) +
    geom_pointrange(aes(ymin = `0.025quant`,
                        ymax = `0.975quant`)) +
    scale_x_log10(label = scientific_10) +
    theme_bw() +
    annotation_logticks(sides = "b") +
    theme(panel.grid.minor.x = element_blank()) +
    xlab("Prior precision of linear effect") +
    ylab("Posterior distribution of linear effect")