Binomial Regression
Chapter 8

Example 8.2.1:  Two independent samples
Domestic vs international graduate students

model{
 for(i in 1:2){ y[i] ~ dbin(theta[i],N[i]) } # Likelihood
 theta[1] ~ dbeta(11.26,11.26) # The prior
 theta[2] ~ dbeta(13.32,6.28)
 beta[1] <- logit(theta[1]) # Induced posterior
 beta[2] <- logit(theta[2]) - logit(theta[1])
 lambda <- theta[1]/theta[2] # Induced posterior for ratio
 prob <- step(lambda -1)  # Posterior prob. ratio > 1.
}
list(N=c(10,10), y=c(5,8))
list(theta=c(0.5,0.5))

	 node	 mean	 sd	 MC error	2.5%	median	97.5%	start	sample
	beta[1]	-2.845E-4	0.3578	0.001833	-0.7103	8.446E-4	0.6979	10001	50000
	beta[2]	0.9834	0.5509	0.002412	-0.08084	0.979	2.082	10001	50000
	lambda	0.7038	0.1505	7.036E-4	0.4398	0.6926	1.032	10001	50000
	prob	0.03614	0.1866	8.884E-4	0.0	0.0	1.0	10001	50000
	theta[1]	0.4999	0.08678	4.45E-4	0.3295	0.5002	0.6677	10001	50000
	theta[2]	0.7202	0.081	3.468E-4	0.5485	0.7255	0.8625	10001	50000

Exercise 8.2:  Oring data

model{
 for(i in 1:n){ 
   y[i] ~ dbin(theta[i],1) # Likelihood
   logit(theta[i]) <- beta[1] + beta[2]*temp[i]
 }
 for(i in 1:2){ tildetheta[i] ~ dbeta(a[i],b[i]) } # Prior
 # Induced prior on the regression coefficients
 beta[1] <- (75/20)*logit(tildetheta[1]) - (55/20)*logit(tildetheta[2])
 beta[2] <-  (-1/20)*logit(tildetheta[1]) +  (1/20)*logit(tildetheta[2])
 for(i in 1:5){ LD[i] <- (logit(alpha[i]) - beta[1]) / beta[2] }
}
list(tildetheta=c(0.5,0.5))
 # Sample size, hyperparameters of prior, and the data
list(n=23, a = c(1.6,1), b = c(1, 1.6), alpha=c(0.90,0.75,0.50,0.25,0.10))
y[ ] temp[ ]
1	53
1	57
1	58
1	63
0	66
0	67
0	67
0	67
0	68
0	69
0	70
0	70
1	70
1	70
0	72
0	73
0	75
1	75
0	76
0	76
0	78
0	79
0	81
END

Example 8.2.3:  Trauma data
Exercise 8.4

model{ 
 for(i in 1:n){ 
   death[i] ~ dbern(theta[i])
   logit(theta[i]) <- beta[1] + beta[2]*ISS[i] + beta[3]*RTS[i] +beta[4]*AGE[i] + beta[5]*TI[i] +beta[6]*AGE[i]*TI[i] 
 }
 for(i in 1:6){ beta[i] ~ dflat() }
 junk <- ID[1]
 }
list(beta= c(0,0,0,0,0,0))
DATA
Get trauma data from website

Exercise 8.6:  Two independent samples
Domestic vs international graduate students

model{
 theta[1] ~ dbeta(11.26,11.26)
 theta[2] ~ dbeta(13.32,6.28)
 beta[1] <- logit(theta[1])
 beta[2] <- logit(theta[2]) - logit(theta[1])
}

Exercise 8.9:  Bayes factors for the Oring data

model{ 
 for (i in 1:n){ 
   y[i] ~ dbin(theta[i],1)
   logit(theta[i]) <- beta[1] + beta[2]*(temp[i]-mt)/sdtemp
 }
 for(i in 1:2){ tildetheta[i] ~ dbeta(a[i],b[i]) }
 beta[1] <- G[1,1]*logit(tildetheta[1]) + G[1,2]*logit(tildetheta[2])
 beta[2] <- G[2,1]*logit(tildetheta[1]) + G[2,2]*logit(tildetheta[2])
 #Have specified logistic model and induced prior on beta
 #Now specify probs under cloglog model
 for(i in 1:n){
   cloglog(thetastar[i]) <- beta[1] + beta[2]*(temp[i]-mt)/sdtemp
   #Now give terms that go into log lik ratio
   #comparing cloglog in num to logistic in den
   v[i] <- y[i]*(log(thetastar[i]) - log(theta[i])) +(1-y[i])*(log(1-thetastar[i])-log(1-theta[i])) 
 }
 #Finally, give corresponding log of prior ratio
 for(i in 1:2){
   cloglog(tildethetastar[i]) <- beta[1] +beta[2]*(ttemp[i]-mt)/sdtemp
   logit(tttildetheta[i]) <- beta[1] + beta[2]*(ttemp[i]-mt)/sdtemp
   u[i] <- (a[i]-1)*log(tildethetastar[i]) - a[i]*log(tttildetheta[i])
               + b[i]*(log(1-tildethetastar[i])- log(1-tttildetheta[i])) 
               + beta[1]+ beta[2]*(ttemp[i]-mt)/sdtemp 
 }
 wstar1 <- sum(v[ ]) 
 wstar2 <- sum(u[ ])
 w <- exp(wstar1+wstar2) 
}
list(tildetheta=c(0.5,0.5))
list(ttemp = c(55,75),mt = 69.56522, sdtemp = 7.05708,
      n=23, a = c(1.6,1.0), b = c(1.0,1.6))
G[,1]                  G[,2]
0.2717391 0.7282609
-0.3528540 0.3528540
END

y[ ] temp[ ]
1	53
1	57
1	58
1	63
0	66
0	67
0	67
0	67
0	68
0	69
0	70
0	70
1	70
1	70
0	72
0	73
0	75
1	75
0	76
0	76
0	78
0	79
0	81
END

Exercise 8.10:  Oring data with standardized covariates

model{ 
 for(i in 1:n){ 
   y[i] ~ dbin(theta[i],1)
   logit(theta[i]) <- gamma[1] + gamma[2]*(temp[i]-mt)/sdtemp
 }
 for(i in 1:2){ tildetheta[i] ~ dbeta(a[i],b[i]) }
 gamma[1] <- Ztinv[1,1]*logit(tildetheta[1]) + Ztinv[1,2]*logit(tildetheta[2])
 gamma[2] <- Ztinv[2,1]*logit(tildetheta[1]) + Ztinv[2,2]*logit(tildetheta[2])
}
list(tildetheta=c(0.5, 0.5))
list(n=23, mt=69.565, sdtemp=7.057, a = c(1.6,1), b = c(1,1.6))
Ztinv[,1]     Ztinv[,2]
0.2717           0.7283
-0.3528         0.3528
END
y[ ] temp[ ]
1	53
1	57
1	58
1	63
0	66
0	67
0	67
0	67
0	68
0	69
0	70
0	70
1	70
1	70
0	72
0	73
0	75
1	75
0	76
0	76
0	78
0	79
0	81
END

	 node	 mean	 sd	 MC error	2.5%	median	97.5%	start	sample
	gamma[1]	-0.9997	0.5026	0.002184	-2.056	-0.9762	-0.07456	10001	50000
	gamma[2]	-1.198	0.4972	0.003094	-2.275	-1.164	-0.316	10001	50000
	tildetheta[1]	0.7743	0.1483	9.255E-4	0.4231	0.8034	0.9741	10001	50000
	tildetheta[2]	0.1471	0.08176	4.563E-4	0.03055	0.1332	0.3394	10001	50000
Exercise 8.11:  Trauma data with standardized covariates

model{ 
 for(i in 1:n){ 
   death[i] ~ dbern(theta[i])
   logit(theta[i]) <- gamma[1] + gamma[2]*(ISS[i]-14.3)/11 + gamma[3]*(RTS[i]-7.29)/1.25 + gamma[4]*(AGE[i]-31.4)/17 + gamma[5]*TI[i] + gamma[6]*(AGE[i]-31.4)*TI[i]/17 
 }
 for(i in 1:6){ 
   tildetheta[i] ~ dbeta(a[i],b[i])
   v[i] <- log(tildetheta[i]/(1-tildetheta[i])) 
   gamma[i] <- inprod(Ztinv[ i,1:6], v[1:6]) 
 }
 junk <- ID[1]
}
list(tildetheta= c(0.5,0.5,0.5,0.5,0.5,0.5))
list(n=300,a = c(1.1,3,5.9,1.3,1.1,1.5), b= c(8.5,11,1.7,12,4.9,5.5))
Ztinv[,1]   Ztinv[,2]   Ztinv[,3]   Ztinv[,4]   Ztinv[,5]   Ztinv[,6]
-2.603          -2.460            -3.702           -3.702            6.063             7.403
-0.343          -0.343             0.343            0.343            0.686             -0.686
-2.091          -2.091            -2.091           -2.091           4.182               4.182
2.901            2.218              2.559            2.559          -5.119              -5.119
1.149            1.005              1.005            1.149          -3.154              -1.154
-5.460          -4.778            -4.778            -5.460          10.238            10.238
END

DATA
Get trauma data from website

Exercise 8.13:  Two-sample problem looking at induced priors

model{
 beta[1] ~ dnorm(0,0.001)
 beta[2] ~ dnorm(0,0.001)
 tildetheta[1]<- exp(beta[1])/(1+exp(beta[1]))
 tildetheta[2]<- exp(beta[1] + beta[2])/(1+exp(beta[1] + beta[2])) 
}

Part d
model{ 
 for(i in 1:2){ tildetheta[i] ~ dunif(0,1) }
 beta[1] <- logit(tildetheta[1])
 beta[2] <- logit(tildetheta[2]) - logit(tildetheta[1])
}

Exercise 8.14:  O-ring data with reference priors

model{ 
 for(j in 1:2){ beta[ j ] ~ dnorm(0,b) }
 for(i in 1:23){
   y[i] ~ dbin(0.5,1)
   logit(theta[i]) <- beta[1] + beta[2]*(temp[i] - mt)/sdtemp
   }
}
list(beta = c(0,0))
list(mt=69.565, sdtemp=7.057, b=1)
y[ ] temp[ ]
1	53
1	57
1	58
1	63
0	66
0	67
0	67
0	67
0	68
0	69
0	70
0	70
1	70
1	70
0	72
0	73
0	75
1	75
0	76
0	76
0	78
0	79
0	81
END

Exercise 8.16:  Two sample problem:  Partial priors

model{ 
 tildetheta[1] ~ dbeta(11.26,11.26)
 beta[1] <- logit(tildetheta[1])
 beta[2] ~ dnorm(0,b)
 tildetheta[2] <- exp(beta[1] +beta[2])/(1+ exp(beta[1]+beta[2]))
}
list(b=1)

Exercise 8.18:  Trauma data:  Partial priors

model{
 for(i in 1:n){ 
   death[i] ~ dbern(0.5)
   logit(theta[i]) <- gamma[1] + gamma[2]*(ISS[i]-14.3)/11 + gamma[3]*(RTS[i] - 7.29)/1.25 +gamma[4]*(AGE[i] -31.4)/17 + gamma[5]*TI[i] + gamma[6]*(AGE[i]-31.4)*TI[i]/17 
 }
 for(i in 1:3){ tildetheta[i] ~ dbeta(a[i],bb[i]) }
 gamma[1] <- logit(tildetheta[1])
 gamma[2] <- logit(tildetheta[2]) - gamma[1]
 gamma[3] <- (logit(tildetheta[2]) - logit(tildetheta[3]))/2
 for(j in 4:6){ gamma[j] ~ dnorm(0,b) }
 junk <- ID[1]
}
list(tildetheta = c(0.1,0.1,0.1))
list(n=300, a=c(2.06,2.7,1.9 ), bb=c(21.2,10.66,3.7), b=1)
DATA
Get trauma data from website

Exercise 8.20:  Toenail fungus data:  Mixed model

model{ 
 for(i in 1:1908){ 
   y[i] ~ dbern(p[i])
   logit(p[i]) <- gamma[ID[i]] + beta[1]*Trt[i] + beta[2]*Time[i] + beta[3]*Trt[i]*Time[i] 
 }
 # Note: ID[ ] identifies individuals,
 # 294 people listed as numbers between 1 and 383.
 for(i in 1:383){ gamma[i] ~ dnorm(mu,tau) }
 beta[1] ~ dnorm(0,0.01)
 beta[2] ~ dnorm(0,0.01)
 beta[3] ~ dnorm(0,0.01)
 mu ~ dnorm(0,0.01)
 tau ~ dgamma(2,0.5)
 junk <- Visit[1] 
}
list(beta=c(0,0,0),mu=0,tau=1)
DATA
Get toenail data from website

Exercise 8.21: Cow Abortion Data
Response is y[i] = 1 if cow aborted and 0 if not
Covariates are GR = gravidity (number of previous successful pregnancies) 
and DO = days open (the number of days between the last birth and the current pregnancy). 

model{
 meangr <- mean(GR[ ])       
 meando <- mean(DO[ ])
 stdgr <- sd(GR[ ])
 stddo <- sd(DO[ ])
 for(k in 1:13145){ 
   grstd[k] <- (GR[k] - meangr)/stdgr   
   dostd[k] <- (DO[k] - meando)/stddo   
 }
 for(k in 1:9) { gamma[k] ~ dnorm(mu, tau) }
 tau <- 1/pow(sigma,2)   
 sigma ~ dunif(0,10)
 for(k in 1:13145){
   y[k] ~ dbern(pr[k])
   logit(pr[k]) <-  beta[1]*dostd[k] + beta[2]*grstd[k] + beta[3]*dostd[k]*grstd[k] +gamma[herd[k]]   
 }
 mu ~ dnorm (-2, 0.001)
 for(i in 1:3){ beta[i] ~ dflat() }
}
list(beta=c( -0.000867, -0.2093, -0.000215), sigma = 1, mu =-2)

DATA
Get cow data from website