Day 12: Logistic Models

Consider the following Binomial data (from the Bioassay Experiment example):

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)

Using the Frequentist approach.

A standard classical analysis fits the model by maximum likelihood. The R function glm is used to do this fitting, and the summary output presents the estimates and the associated standard errors.

response = cbind(y, n - y)
results = glm(response ~ x, family = binomial)
summary(results)
## 
## Call:
## glm(formula = response ~ x, family = binomial)
## 
## 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

To see how well the resulting logistic model fit the data, we could create the following plot

p.hats=y/n
plot(x,p.hats)
curve(exp(0.85+7.75*x)/(1+exp(0.85+7.75*x)),-1,1,add=T)

#points(x,exp(0.85+7.75*x)/(1+exp(0.85+7.75*x)),pch=19)

Using the Bayesian approach.

Suppose that the user has prior beliefs about the regression parameters that she inputs through the following conditional means prior. This prior is constructed by thinking about the probability of death at two different dose levels, \(x_L\) and \(x_H\). When the dose level is \(x_L = −0.7\), the median and 90th percentile of the probability of death \(p_L\) are respectively 0.2 and 0.5. When the dose level is \(x_H = 0.6\), the user believes that the median and 90th percentile of the probability of death \(p_H\) are given respectively by 0.8 and 0.98. One matches this information with a beta prior using the beta.select function.

library(LearnBayes)
beta.select(list(p=.5,x=.2),list(p=.9,x=.5))
## [1] 1.12 3.56
beta.select(list(p=.5,x=.8),list(p=.9,x=.98))
## [1] 2.10 0.74
prior=rbind(c(-0.7, 4.68, 1.12), c(0.6, 2.84, 0.74))
data.new=rbind(data, prior)

We note that this prior has the same functional form as the likelihood, where the beta parameters can be viewed as the numbers of deaths and survivals in a prior experiment performed at two dose levels used in constructing the priors. The function ‘logisticpost’ and ‘simcontour’ can be used to obtain simulated values from the resulting posterior distribution.

mycontour(logisticpost,c(-3,3,-1,9),data.new, xlab="beta0", ylab="beta1")

mycontour(logisticpost,c(-3,3,-1,9),data.new, xlab="beta0", ylab="beta1")
s=simcontour(logisticpost,c(-2,3,-1,11),data.new,1000)
points(s)

plot(density(s$y),xlab="beta1")

mean(s$x); mean(s$y)
## [1] -0.2068083
## [1] 2.621943
p.hats=y/n
plot(x,p.hats)
curve(exp(mean(s$x)+mean(s$y)*x)/(1+exp(mean(s$x)+mean(s$y)*x)),-1,1,add=T)

#points(x,exp(0.85+7.75*x)/(1+exp(0.85+7.75*x)),pch=19)

In this setting, one parameter of interest is the LD-50, the value of the dose \(x\) such that the probability of death is equal to one-half. Note that the LD-50 is equal to \(\theta=-\beta_0/\beta_1\).

theta=-s$x/s$y
hist(theta,xlab="LD-50",breaks=20)

The 95% probability interval for LD-50, \(\theta\)

quantile(theta,c(.025,.975))
##       2.5%      97.5% 
## -0.3061804  0.7010341