Bayesian Example on Dosage-Group Rates of Liver Cancer in Lab Animals
===================================================================== 11/8/17

data(aflatoxin)

Format

A data frame with 6 observations on the following 3 variables.

dose   = dose (of aflatoxin) in ppb
total  = number of test animals
tumor  = number with liver cancer

## Consider Bayesian analysis with 6 groupwise probability parameters p1 .. p6
## with prior  Beta(0.5,0.5) for each independently. Problem is to understand
## the distribution of the fitted intercept and slope for the line 
## log(1+aflatoxin$dose) ~ log((0.5+aflatoxin$tumor)/(aflatoxin$total+1)).

> summary(lm(log(1+aflatoxin$dose) ~ log((0.5+aflatoxin$tumor)/(aflatoxin$total+1))))$coef[,1:2]
                                                   Estimate Std. Error
(Intercept)                                        4.400240  0.4456295
log((0.5 + aflatoxin$tumor)/(aflatoxin$total + 1)) 1.207899  0.2056075

## and corresponding Rsq = 0.896.
## Of course this "SE" is based on normal theory, so not helpful here.

## The posterior densities of p1 .. p6 are Beta with parameters
> AB = rbind(alpha=0.5+aflatoxin$tumor, beta=0.5+aflatoxin$total-aflatoxin$tumor)
      [,1] [,2] [,3] [,4] [,5] [,6]
alpha  0.5  2.5  1.5  4.5 20.5 28.5
beta  18.5 20.5 21.5 17.5  5.5  0.5

> phat = aflatoxin$tumor/aflatoxin$total             #### MLE
> round(phat,4)
[1] 0.0000 0.0909 0.0455 0.1905 0.8000 1.0000

## To obtain Bayesian and PB distributions (based on data) for the desired quantities:

> ParmArr = array(0, c(5000,2,2), dimnames=list(NULL, c("Bayes","ML"), 
                 c("Int","Slope")))
  dosfcn = log(1+aflatoxin$dose)
  dosSSQ = 5*var(dosfcn)
  dosav = mean(dosfcn)
  dosdif = dosfcn-dosav
  tot = aflatoxin$total
  totp = tot+1
  DoseRegr = function(mvec) {
       yv = log((0.5+mvec)/totp)
       Slop = sum(yv*dosdif)/dosSSQ
       Int = mean(yv)-Slop*dosav
       c(Int=Int,Slop=Slop) }

  set.seed(57575)
  MatBayes = array(rbinom(5000*6, tot, rbeta(5000*6, AB[1,], AB[2,])), c(6,5000))
  MatML = array(rbinom(5000*6, tot, phat), c(6,5000))

  for(i in 1:5000) {
       ParmArr[i,1,] = DoseRegr(MatBayes[,i])
       ParmArr[i,2,] = DoseRegr(MatML[,i]) }

   hist(ParmArr[,1,2], nclass=40, prob=T, xlab="Slope",
      main="Slopes from Bayesian Predictive Datasets")
   lines(density(ParmArr[,1,2]), col="red",lty=3)
   abline(v=quantile(ParmArr[,1,2],prob=.05), col="blue")
   abline(v=quantile(ParmArr[,1,2],prob=.95), col="blue")

### Recall that original dataset gave:
> DoseRegr(aflatoxin$tumor)
       Int       Slop 
-3.4450578  0.7418989 

   hist(ParmArr[,2,2], nclass=40, prob=T, xlab="Slope",
      main="Slopes from Parametric-Bootstrap Datasets")
   lines(density(ParmArr[,2,2]), col="blue",lty=3)
   abline(v=quantile(ParmArr[,2,2],prob=.05), col="red")
   abline(v=quantile(ParmArr[,2,2],prob=.95), col="red")