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
56
57
|
#############################################################
# 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)
S=readline(prompt="Type <Return> to continue : ")
###################################################
# 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
|
