# # Chapter 4 -- Bayesian Computation With R # Bioassay Experiment # # Remove all objects just to be safe # rm(list=ls(all=TRUE)) # library(LearnBayes) # # Deaths of Animals in Drug Experiment # # Dose Deaths Sample Size # -0.86 0 5 # -0.30 1 5 # -0.05 3 5 # 0.73 5 5 # # # Deaths ~ Binomial(n,p) # p = exp(beta*X)/(1 + exp(beta*X) = 1/(1 + exp(-beta*X)) # x <- c(-0.86, -0.3, -0.05, 0.73) n <- c(5, 5, 5, 5) y <- c(0, 1, 3, 5) data <- cbind(x, n, y) response <- cbind(y, n - y) # results <-glm(response ~ x,family=binomial(link=logit)) summary(results) #&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& # Call: # glm(formula = response ~ x, family = binomial(link = logit)) # # Deviance Residuals: # 1 2 3 4 # -0.17236 0.08133 -0.05869 0.12237 # # Coefficients: # Estimate Std. Error z value Pr(>|z|) # (Intercept) 0.8466 1.0191 0.831 0.406 # x 7.7488 4.8728 1.590 0.112 # # (Dispersion parameter for binomial family taken to be 1) # # Null deviance: 15.791412 on 3 degrees of freedom # Residual deviance: 0.054742 on 2 degrees of freedom # AIC: 7.9648 # # Number of Fisher Scoring iterations: 7 #&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& # # Prior belief about regression parameters -- if dose level = -0.7, # then the median and 90th percentile of the probability of death # are 0.2 and 0.5. These numbers are fed to the beta.select function # in the LearnBayes Library to find the proper BETA distribution # priorlowdoselevel <- beta.select(list(p=.5,x=.2),list(p=.9,x=.5)) # [1] 1.12 3.56 # # Prior belief continued -- if does level is high = 0.6, then the belief # is that the median and 90th percentile of the probability of death # are 0.8 and 0.98. Again, the beta.select function is used: # priorhighdoselevel <- beta.select(list(p=.5,x=.8),list(p=.9,x=.98)) # [1] 2.10 0.74 # # Deaths of Animals FROM PRIOR Drug Experiment # # Dose Deaths Sample Size # -0.70 1.12 1.12+3.56=4.68 (using the fact that the two parameters are # 0.60 2.10 2.10+0.74=2.84 the number of deaths and survivors! Add # to get the total) # # When these two prior "data points" are combined with the orginal sample, # then we have a POSTERIOR distribution. # # Bind the new rows with rbind(...) -- NOTE THE TYPO ON p.72!!!!! He has this: # prior <- rbind(c(-0.7, 4.68, 1.12),c(0.6, 2.10, 0.74)) This is clearly incorrect # because with a dose of 0.6 almost everyone should DIE -- he has numbers backwards. # prior <- rbind(c(-0.7, 4.68, 1.12),c(0.6, 2.84, 2.10)) data.new <- rbind(data, prior) # # Plot the Posterior Distribution -- the function logisticpost is in the LearnBayes Library, # it computes the log posterior density. The data argument is a matrix organized as # above -- with dose, # successes, sample size, as the columns. # -- the function mycontour is also in the LearnBayes Library, # it facilitates the use of the R contour command # windows() mycontour(logisticpost,c(-3,3,-1,11),data.new,xlab="",ylab="",main="",col="red",font=2) # Main title mtext("Contour Plot of the Posterior Distribution of Beta0 and Beta1\nFor the Bioassay Example",side=3,line=1.00,cex=1.2,font=2) # x-axis title mtext("Beta0",side=1,line=2.75,font=2,cex=1.2) # y-axis title mtext("Beta1",side=2,line=2.75,font=2,cex=1.2) # # Simcontour is a function in the LearnBayes Library that simulates pairs of (Beta0, Beta1) from # Above Posterior Density # s <- simcontour(logisticpost,c(-2,3,-1,11),data.new,1000) points(s, pch=21, col = "blue", cex = 1.1,font=2) # windows() plot(density(s$y),xlab="",ylab="",main="",col="blue",font=2) # Main title mtext("Density of Posterior Distribution of Beta1\nFor the Bioassay Example",side=3,line=1.00,cex=1.2,font=2) # x-axis title mtext("Beta1",side=1,line=2.75,font=2,cex=1.2) # y-axis title mtext("Density",side=2,line=2.75,font=2,cex=1.2) # # Dose such that Death = 0.5 is -Beta0/Beta1 (see Poole, 2006 for math demonstration) # theta <- -s$x/s$y windows() hist(theta,xlab="",ylab="",main="",col="red",font=2) # Main title mtext("Histogram of Simulated Values of the LD-50 Parameter\nProbablity of Death = 0.5",side=3,line=1.00,cex=1.2,font=2) # x-axis title mtext("LD-50 -- Death Probability of 1/2",side=1,line=2.75,font=2,cex=1.2) # y-axis title mtext("Frequency",side=2,line=2.75,font=2,cex=1.2) #