bcwr/Chapter.10.3.R at master · rghan/bcwr · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
#############################################################
# Section 10.3  Binary Response Regression with a Probit Link
#############################################################

#################################################
# Section 10.3.1. Missing data and Gibbs sampling
#################################################

library(LearnBayes)

 data(donner)
 attach(donner)
 X=cbind(1,age,male)

 fit=glm(survival~X-1,family=binomial(link=probit))
 summary(fit)

 m=10000
 fit=bayes.probit(survival,X,m)

 apply(fit$beta,2,mean)

 apply(fit$beta,2,sd)

 a=seq(15,65)
 X1=cbind(1,a,1)
 p.male=bprobit.probs(X1,fit$beta)

 plot(a,apply(p.male,2,quantile,.5),type="l",ylim=c(0,1),
   xlab="age",ylab="Probability of Survival")
 lines(a,apply(p.male,2,quantile,.05),lty=2)
 lines(a,apply(p.male,2,quantile,.95),lty=2)

###################################################
# Section 10.3.2  Proper priors and model selection
###################################################

library(LearnBayes)
data(donner)
y=donner$survival
X=cbind(1,donner$age,donner$male)

beta0=c(0,0,0); c0=100
P0=t(X)%*%X/c0

bayes.probit(y,X,1000,list(beta=beta0,P=P0))$log.marg

bayes.probit(y,X[,-2],1000,
   list(beta=beta0[-2],P=P0[-2,-2]))$log.marg

bayes.probit(y,X[,-3],1000,
   list(beta=beta0[-3],P=P0[-3,-3]))$log.marg

bayes.probit(y,X[,-c(2,3)],1000,
   list(beta=beta0[-c(2,3)],P=P0[-c(2,3),-c(2,3)]))$log.marg