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
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
|
################################
# Section 9.2.6 An Example
################################
library(LearnBayes)
data(birdextinct)
attach(birdextinct)
logtime=log(time)
plot(nesting,logtime)
out = (logtime > 3)
text(nesting[out], logtime[out], label=species[out], pos = 2)
S=readline(prompt="Type <Return> to continue : ")
windows()
plot(jitter(size),logtime,xaxp=c(0,1,1))
S=readline(prompt="Type <Return> to continue : ")
windows()
plot(jitter(status),logtime,xaxp=c(0,1,1))
##### Least-squares fit
fit=lm(logtime~nesting+size+status,data=birdextinct,x=TRUE,y=TRUE)
summary(fit)
##### Sampling from posterior
theta.sample=blinreg(fit$y,fit$x,5000)
S=readline(prompt="Type <Return> to continue : ")
windows()
par(mfrow=c(2,2))
hist(theta.sample$beta[,2],main="NESTING",
xlab=expression(beta[1]))
hist(theta.sample$beta[,3],main="SIZE",
xlab=expression(beta[2]))
hist(theta.sample$beta[,4],main="STATUS",
xlab=expression(beta[3]))
hist(theta.sample$sigma,main="ERROR SD",
xlab=expression(sigma))
apply(theta.sample$beta,2,quantile,c(.05,.5,.95))
quantile(theta.sample$sigma,c(.05,.5,.95))
S=readline(prompt="Type <Return> to continue : ")
###### Estimating mean extinction times
cov1=c(1,4,0,0)
cov2=c(1,4,1,0)
cov3=c(1,4,0,1)
cov4=c(1,4,1,1)
X1=rbind(cov1,cov2,cov3,cov4)
mean.draws=blinregexpected(X1,theta.sample)
c.labels=c("A","B","C","D")
windows()
par(mfrow=c(2,2))
for (j in 1:4)
hist(mean.draws[,j],
main=paste("Covariate set",c.labels[j]),xlab="log TIME")
S=readline(prompt="Type <Return> to continue : ")
######## Predicting extinction times
cov1=c(1,4,0,0)
cov2=c(1,4,1,0)
cov3=c(1,4,0,1)
cov4=c(1,4,1,1)
X1=rbind(cov1,cov2,cov3,cov4)
pred.draws=blinregpred(X1,theta.sample)
c.labels=c("A","B","C","D")
windows()
par(mfrow=c(2,2))
for (j in 1:4)
hist(pred.draws[,j],
main=paste("Covariate set",c.labels[j]),xlab="log TIME")
S=readline(prompt="Type <Return> to continue : ")
######### Model checking via posterior predictive distribution
pred.draws=blinregpred(fit$x,theta.sample)
pred.sum=apply(pred.draws,2,quantile,c(.05,.95))
par(mfrow=c(1,1))
ind=1:length(logtime)
windows()
matplot(rbind(ind,ind),pred.sum,type="l",lty=1,col=1,
xlab="INDEX",ylab="log TIME")
points(ind,logtime,pch=19)
out=(logtime>pred.sum[2,])
text(ind[out], logtime[out], label=species[out], pos = 4)
S=readline(prompt="Type <Return> to continue : ")
######### Model checking via bayes residuals
prob.out=bayesresiduals(fit,theta.sample,2)
windows()
par(mfrow=c(1,1))
plot(nesting,prob.out)
out = (prob.out > 0.35)
text(nesting[out], prob.out[out], label=species[out], pos = 4)
|
