################################################################### # # Example 6.2.1 from Press (2003) # ################################################################### # # Occurrence or nonoccurrence of infection following birth by # Ceasarian section. # # Response variable # # y=1 if caesarian birth resulted in an infection # y=0 if caesarian birth DID NOT result in an infection # # Explanatory variables # # x1=1 if caesarian was nonplanned # x2=1 if risk factors were present at the time of birth # x3=1 if antibiotics were given as a prophylaxis # ################################################################### # Auxiliary functions and densities # --------------------------------- like=function(beta){ prob = pnorm(x%*%beta) return(prod((prob^y)*((1-prob)^(1-y)))) } like1=function(x,beta){ prob = pnorm(x%*%beta) return(prod((prob^y)*((1-prob)^(1-y)))) } prior = function(beta){prod(dnorm(beta,0,2.236068))} rprior = function(n){matrix(rnorm(4*n,0,2.236068),n,4)} post = function(beta){like(beta)*prior(beta)} post1 = function(x,beta){like1(x,beta)*prior(beta)} q025 = function(x){quantile(x,0.025)} q975 = function(x){quantile(x,0.975)} # Dataset and exploratory analysis # -------------------------------- n = c(98,18,2,26,58,9,40) count = c(11,1,0,23,28,0,8) y = NULL for (i in 1:7) y = c(y,rep(0,n[i]-count[i]),rep(1,count[i])) x = matrix(c(rep(c(1,1,1),n[1]),rep(c(0,1,1),n[2]),rep(c(0,0,1),n[3]), rep(c(1,1,0),n[4]),rep(c(0,1,0),n[5]),rep(c(1,0,0),n[6]), rep(c(0,0,0),n[7])),ncol=3,byrow=T) x = cbind(1,x) par(mfrow=c(3,1)) plot(y,xlab="birth",ylab="",main="x1=1 if caesarian was nonplanned") lines(x[,2],col=2) plot(y,xlab="birth",ylab="",main="x2=1 if risk factors were present at the time of birth") lines(x[,3],col=2) plot(y,xlab="birth",ylab="",main="x3=1 if antibiotics were given as a prophylaxis") lines(x[,4],col=2) # Maximum likelihood estimates # ---------------------------- bhat = c(-1.093022,0.607643,1.197543,-1.904739) Vhat = matrix(0,4,4) Vhat[1,1] = 0.040745 Vhat[2,2] = 0.073101 Vhat[3,3] = 0.062292 Vhat[4,4] = 0.080788 Vhat[1,2] = -0.007038 Vhat[1,3] = -0.039399 Vhat[1,4] = 0.004829 Vhat[2,3] = -0.006940 Vhat[2,4] = -0.050162 Vhat[3,4] = -0.016803 for (i in 2:4) for (j in 1:(i-1)) Vhat[i,j] = Vhat[j,i] cVhat = t(chol(Vhat)) # MCMC scheme (Random walk Metropolis) # ------------------------------------ set.seed(1234) beta = bhat tau = 1 betas = NULL burn = 1000 M = 1000 LAG = 20 niter = burn+LAG*M for (iter in 1:niter){ betadraw = beta + tau*cVhat%*%rnorm(4) accept = min(1,post(betadraw)/post(beta)) if (runif(1)