simulated_zipoisson.Rmd

---
title: "ZeroInflatedPoisson"
author: "Anders Sundelin"
date: "2022-10-18"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
library(rethinking)
```

## Simulating issue prediction with Bayesian Models
```{r init_sim}
N <- 350
set.seed(700716)

# probability of a event not happening
prob_untouched <- 0.2

# on average, we add 10 rows per week
change_size <- rnorm(N, 10, 25)
# tells how "scouty" our average developer are. High value -> more likely to remove issues, less likely to add new ones
scoutiness <- rbinom(N, 10, .66)

SZ <- (change_size - mean(change_size))/sd(change_size)
SC <- (scoutiness - mean(scoutiness))/sd(scoutiness)

# The log of the rate is a linear model - rates are always positive
# Thus, the rate is the exponentiation of the linear model
rate_introd <- exp(2.4 + .33 * SZ - SC)

untouched <- rbinom(N, 1, prob_untouched)

# simulate issues introduced
introd <- (1-untouched) * rpois(N, rate_introd)
```

```{r}
simplehist(introd, xlab="issues introduced", lwd=4)
```


### brms model (zi NegBin)

```{r}
# given that var(y) != mean(y) a Poisson might not be the right thing to use.
library(brms) # rethinking does not have support for zi-negbin

introd <- as.numeric(introd)
size <- scale(change_size)
scout <- scale(scoutiness)

d <- data.frame(y = introd, size = size, scout = scout)

p <- get_prior(bf(
  y ~ 1 + size + scout,
  zi ~ 1),
  family = zero_inflated_negbinomial(),
  data = d)

p[1,1] <- "normal(0, 0.5)"
p[4,1] <- "normal(0, 1)"

m1 <- brm(bf(
  y ~ 1 + size + scout,
  zi ~ 1),
  family = zero_inflated_negbinomial(),
  data = d,
  prior = p
)

summary(m1)
```

### Recovering the model parameters (zi-poisson)

```{r}
issue_model <- ulam(
  alist(
    y ~ dzipois(p, lambda),
    logit(p) <- ap, # assume simplest possible model
    log(lambda) <- al + bsize * size + bscout * scout,
    ap ~ dnorm(0, 1.5),
    al ~ dnorm(3, .5),
    bsize ~ dnorm(0, 0.4),
    bscout ~ dnorm(0, 0.4)
  ), data=list(y=introd, size=SZ, scout=SC), chains=4)
```

```{r}
precis(issue_model)
```

With the corrected lambda (making log(lambda) into a linear model), we can easily recover the parameters of the model.

```{r}
post <- extract.samples(issue_model)
mean(inv_logit(post$ap))
```
```{r}
ns <- 100
SZ_seq <- seq(from=-3, to=3, length.out=ns)
SC_seq <- seq(from=-3, to=3, length.out=ns)
lambda <- link(issue_model, data=data.frame(size=SZ_seq, scout=SC_seq))
lmu <- apply(lambda$lambda, 2, mean)
lci <- apply(lambda$lambda, 2, PI)
plot(SZ, introd, xlim=range(SZ_seq), ylim=c(0,150), xlab="normalized size", ylab="issues introduced")
lines(SZ_seq, lmu, lty=1, lwd=1.5)
shade(lci, SZ_seq, xpd=TRUE)
```
```{r}
plot(SC, introd, xlim=range(SC_seq), ylim=c(0,150), xlab="normalized scoutiness", ylab="issues introduced")
lines(SC_seq, lmu, lty=1, lwd=1.5)
shade(lci, SC_seq, xpd=TRUE)
```


Question for my Bayesian buddies, including Prof. Torkar:

* In my analysis, I standardized the data, using the complete data set, then used that data to recover the parameters.

A: That is a first good step to see if your model can recover the basic (ground truth) parameters.

* My hunch is that this introduces a risk of overfitting - making the model trust the data too much.

A: Aah, yes, indeed. The second step could be to do what you propose next below.

* I guess we should have sampled some data, standardized using that dataset, and then run the model on the remainder - akin to ML train/test/validation sets?

A: But remember also that we rely on an epistemological reasoning here, i.e., information theory (LOO). What LOO does is approximating leave-one-out (LOO) CV. It does it remarkably well and the results are very reliable (including the fact that we also get diagnostics w/ LOO). I have over time fallen back on relying more and more on LOO and not focus on scaling the sample and use that before I use the test/validation set.


### Prior selection, incl. plots

Choosing a "flat prior" in log space can be challenging (and usually requires standardizing variables).
If log(lambda) has a normal distribution, then lambda has a log-normal distribution.
Plotting the priors
```{r}
curve(dlnorm(x, 30, 30), from=0, to=100, n=200)
```
That is a prior that is extremely biased towards 0. So should not be used in log space.

Seeking a more sensible prior (we expect about 20-40 issues to be able to be introduced):
```{r}
curve(dlnorm(x, 3, .5), from=0, to=100, n=200)
```

Prior for beta:

Choosing blindly, we chose a norm(20,40) prior for beta
```{r}
Nbeta <- 100
a <- rnorm(Nbeta, 3, .5)
b <- rnorm(Nbeta, 20, 40)
plot(NULL, xlim=c(-2,2), ylim=c(0,100))
for (i in 1:Nbeta) curve(exp(a[i] + b[i]*x), add=TRUE, col=grau())
```

Wild linear relationships - we need to tame the prior.


```{r }
Nbeta <- 100
a <- rnorm(Nbeta, 3, 0.5)
b <- rnorm(Nbeta, 0, 0.4)
plot(NULL, xlim=c(-2,2), ylim=c(0,100))
for (i in 1:Nbeta) curve(exp(a[i] + b[i]*x), add=TRUE, col=grau())
```

The Normal(0, 0.4) prior allows some strong relationships, while keeping most of the lines relatively straight (skeptical of the data).

### Viewing the priors on the natural outcome scale

```{r}
x_seq <- seq(from=-3, to=3, length.out=100)
lambda <- sapply(x_seq, function(x) exp(a+b*x))
plot(NULL, xlim=range(x_seq), ylim=c(0,500), xlab="size", ylab="issues introduced")
for (i in 1:Nbeta) lines(x_seq, lambda[i,], col=grau(), lwd=1.5)

```

Most of the priors are relatively flat, but a few show strong coupling. This would allow such trends to appear, if they were present in the data. But in general, the priors assume "no relation" - this is because we want the data to show these relationships - we do not want to bias our analysis via the prior.