---
title: "SISG Bayes: Session 4 INLA example"
author: "Ken Rice"
date: "`r Sys.Date()`"
output: pdf_document
---

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

Implementing the LHON example via INLA

```{r}
lhon <- data.frame( y=c(1,1,1,0,0,0), x=c(0,1,2,0,1,2), n=c(6,8,75,10,66,163) )
dmat <- rbind(lhon$n[4:6], lhon$n[1:3])
names(dmat) <- c("CC","CT","TT")

barplot(dmat, names.arg=c("CC","CT","TT"), legend.text=c("case","control")[2:1],
args.legend=list(x="topleft", bty="n"))
```

Setting this up for a non-Bayesian analysis, fitting a GLM:

```{r}
# setup data
cc.dat <- data.frame(x=c(0,1,2), success=c(6,8,75), fail=c(10,66,163))

# estimates and confidence interval from GLM analysis
logitmod <- glm(cbind(success,fail)~x,family="binomial", data=cc.dat)
coef(summary(logitmod))
confint.default(logitmod)
```

Now using INLA -- note the similar syntax, and for now using its default priors:

```{r}
library("INLA")
cc.mod <- inla(success~x,family="binomial",data=cc.dat,Ntrials=success+fail)
summary(cc.mod)
```

Now implementing prior $N(0,10)$ for the intercept (i.e. log odds for $X=0$) and -- informatively -- $N(0,\sigma^2)$ for the log odds ratio, with $\sigma$ such that the upper 95\% point of the distribution is at $OR=5$

```{r}
Upper95 <- log(5) # note log(5) = 1.61
sigma <- Upper95/qnorm(0.95)
cc.inla <- inla(success~x,family="binomial",data=cc.dat,Ntrials=success+fail,
   control.fixed=list(mean.intercept=c(0),prec.intercept=c(1/10),
                      mean=c(0),prec=c(1/sigma^2)))
summary(cc.inla)
```

For a more informative prior, repeat but set the upper 97.5\% point of the prior to $OR=1.5$:

```{r}
Upper975 <- 1.5 # note log(1.5)=0.41
sigma <- log(Upper975)/qnorm(0.975)
cc.inf.inla <- inla(success~x,family="binomial",data=cc.dat,Ntrials=success+fail,
   control.fixed=list(mean.intercept=c(0),prec.intercept=c(1/10),
                      mean=c(0),prec=c(1/sigma^2)))
summary(cc.inf.inla)
```

Some graphs to compare them, on the log odds ratio scale:

```{r}
#pdf("inlalhonpost.pdf", w=8, h=4.5)
par(mar=c(4,4,0,0)+0.2)
curve(inla.dmarginal(x, cc.inf.inla$marginals.fixed$x), -1,1.5, 
	col="forestgreen", ylab="density", lty=2, xlab=expression(beta[1]))
curve(inla.dmarginal(x, cc.inla$marginals.fixed$x), -1,1.5, 
	col="forestgreen", ylab="density", lty=1, add=TRUE)
curve(dnorm(x,0,log(Upper975)/qnorm(0.975)), -1, 1.5, col=4, lty=2, add=TRUE)
curve(dnorm(x,0,Upper95/qnorm(0.95)), -1, 1.5, col=4, lty=1, add=TRUE)
lines(y=c(0,0), x=confint.default(logitmod,2), lwd=2, col="goldenrod")
points(y=0, x=coef(logitmod)[2], lwd=2, col="goldenrod")
legend("topleft", lty=1:2, c("diffuse","informative"), bty="n")
legend("topright", lty=1, bty="n", col=c(4,"forestgreen","goldenrod"), c("prior","posterior","non-Bayes CI"))
#dev.off()
```

Add on the odds ratio scale:

```{r}
#pdf("inlalhonpostexp.pdf", w=8, h=4.5)
par(mar=c(4,4,0,0)+0.2)
plot( inla.tmarginal(exp, cc.inla$marginals.fixed$x) , ylab="density",
type="l", xlab=expression("Odds ratio, "*e^beta[1]), col="forestgreen", 
ylim=c(0,2.2), xlim=c(0.5, 3.5))
lines( inla.tmarginal(exp, cc.inf.inla$marginals.fixed$x) , 
col="forestgreen", type="l", lty=2)
curve(dlnorm(x,0,log(Upper975)/qnorm(0.975)), exp(-1), exp(1.5), col=4, lty=2, add=TRUE)
curve(dlnorm(x,0,Upper95/qnorm(0.95)), exp(-1), exp(1.5), col=4, lty=1, add=TRUE)
lines(y=c(0,0), x=exp(confint.default(logitmod,2)), lwd=2, col="goldenrod")
points(y=0, x=exp(coef(logitmod)[2]), lwd=2, col="goldenrod")
legend("topright", lty=1:2, c("diffuse","informative"), bty="n")
legend("right", lty=1, bty="n", col=c(4,"forestgreen","goldenrod"), c("prior","posterior","non-Bayes CI"))
#dev.off()
```