bcwr/Chapter.10.4.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
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
###################################################
# Section 10.4 Estimating a Table of Means
###################################################

library(LearnBayes)

 data(iowagpa)
 rlabels = c("91-99", "81-90", "71-80", "61-70", "51-60", "41-50",
     "31-40", "21-30")
 clabels = c("16-18", "19-21", "22-24", "25-27", "28-30")
 gpa = matrix(iowagpa[, 1], nrow = 8, ncol = 5, byrow = T)
 dimnames(gpa) = list(HSR = rlabels, ACTC = clabels)
 gpa

 samplesizes = matrix(iowagpa[, 2], nrow = 8, ncol = 5, byrow = T)
 dimnames(samplesizes) = list(HSR = rlabels, ACTC = clabels)
 samplesizes

 act = seq(17, 29, by = 3)
 matplot(act, t(gpa), type = "l", lwd = 3,
  xlim = c(17, 34), col=1:8, lty=1:8)
 legend(30, 3, lty = 1:8, lwd = 3, legend = c("HSR=9", "HSR=8",
     "HSR=7", "HSR=6", "HSR=5", "HSR=4", "HSR=3", "HSR=2"), col=1:8)

 MU = ordergibbs(iowagpa, 5000)

 postmeans = apply(MU, 2, mean)
 postmeans = matrix(postmeans, nrow = 8, ncol = 5)
 postmeans=postmeans[seq(8,1,-1),]
 dimnames(postmeans)=list(HSR=rlabels,ACTC=clabels)
 round(postmeans,2)

matplot(act, t(postmeans), type = "l", lty=1:8, lwd = 3, col = 1, xlim = c(17, 34))
 legend(30, 3, lty = 1:8, lwd = 2, legend = c("HSR=9", "HSR=8",
     "HSR=7", "HSR=6", "HSR=5", "HSR=4", "HSR=3", "HSR=2"))

 postsds = apply(MU, 2, sd)
 postsds = matrix(postsds, nrow = 8, ncol = 5)
 postsds=postsds[seq(8,1,-1),]
 dimnames(postsds)=list(HSR=rlabels,ACTC=clabels)
 round(postsds,3)

 s=.65
 se=s/sqrt(samplesizes)
 round(postsds/se,2)

 FIT=hiergibbs(iowagpa,5000)

 par(mfrow=c(2,1))
 plot(density(FIT$beta[,2]),xlab=expression(beta[2]),
  main="HIGH SCHOOL RANK")
 plot(density(FIT$beta[,3]),xlab=expression(beta[3]),
  main="ACT SCORE")
 quantile(FIT$beta[,2],c(.025,.25,.5,.75,.975))

 quantile(FIT$beta[,3],c(.025,.25,.5,.75,.975))

 quantile(FIT$var,c(.025,.25,.5,.75,.975))

 posterior.means = apply(FIT$mu, 2, mean)
 posterior.means = matrix(posterior.means, nrow = 8, ncol = 5,
  byrow = T)

 par(mfrow=c(1,1))
 matplot(act, t(posterior.means), type = "l", lwd = 3, lty=1:8, col=1,
   xlim = c(17, 34))
 legend(30, 3, lty = 1:8, lwd = 2, legend = c("HSR=9", "HSR=8",
     "HSR=7", "HSR=6", "HSR=5", "HSR=4", "HSR=3", "HSR=2"))

 p=1-pnorm((2.5-FIT$mu)/.65)
 prob.success=apply(p,2,mean)

 prob.success=matrix(prob.success,nrow=8,ncol=5,byrow=T)
 dimnames(prob.success)=list(HSR=rlabels,ACTC=clabels)
 round(prob.success,3)