Jim — Oct 15, 2013, 3:50 PM
# Career Trajectory Example
# Exchangeable Modeling
library(LearnBayes)
# Model Description
modelString = "
model {
for (i in 1:n){
y[i] ~ dnorm (y.hat[i], tau.y[i])
tau.y[i] <- pa[i]*pow(sigma.y, -2)
y.hat[i] <- a[player[i]] + b[player[i]]*x[i] + c[player[i]]*x[i]*x[i]
}
sigma.y ~ dunif (0, 100)
for (j in 1:J){
a[j] <- B[j,1]
b[j] <- B[j,2]
c[j] <- B[j,3]
B[j,1:3] ~ dmnorm (mu.beta[], Tau.B[,])
}
mu.beta[1:3] ~ dmnorm(mean[1:3],prec[1:3 ,1:3 ])
Tau.B[1:3 , 1:3] ~ dwish(Omega[1:3 ,1:3 ], 3)
}
"
writeLines(modelString, con="hoftraj.bug")
# Use sluggerdata data frame from LearnBayes packages
# y is batting average AVG = H / AB
# player is number of slugger (1 through 10)
b.data <- list(y=sluggerdata$H / sluggerdata$AB,
x=sluggerdata$Age,
player=as.integer(sluggerdata$Player),
pa=sluggerdata$AB,
n=199, J=10,
mean = c(0, 0, 0),
Omega=diag(c(.1,.1,.1)),
prec=diag(c(1.0E-6,1.0E-6,1.0E-6)))
# initial MCMC values (maybe these can be improved)
b.inits <- list(B=matrix(c(-7.69,.350,-.0058),
nrow=10,ncol=3,byrow=TRUE),
mu.beta=c(-7.69,.350,-.0058),
Tau.B=diag(c(.1,.1,.1)))
library(rjags)
Loading required package: coda
Loading required package: lattice
Linked to JAGS 3.4.0
Loaded modules: basemod,bugs
# fit model using JAGS
jags <- jags.model('hoftraj.bug',
data = b.data,
inits = b.inits,
n.chains = 1,
n.adapt = 100)
Compiling model graph
Resolving undeclared variables
Allocating nodes
Graph Size: 1626
Initializing model
update(jags, 1000)
# Collect the simulated draws of the proportion for 10,000 iterations
posterior <- coda.samples(jags, c("B"),
n.iter=10000)
# posterior means of regression vectors
post.means <- t(matrix(apply(posterior[[1]], 2, mean), 10, 3))
# compare individual and smoothed trajectories
# player is number of slugger (1 through 10)
sluggerdata$player <- as.integer(sluggerdata$Player)
# get individual estimates
ind.est <- function(d){
d$rate <- with(d, H / AB)
c(lm(rate ~ Age + I(Age^2), weights=AB, data=d)$coef,
min(d$Age), max(d$Age))
}
library(plyr)
OUT <- ddply(sluggerdata, .(player), ind.est)
# graph individual trajectories in blue
graph.est <- function(j, OUT, ...)
curve(OUT[j, 2] + OUT[j, 3] * x + OUT[j, 4] * x ^ 2,
from=OUT[j, 5], to=OUT[j, 6], add=TRUE, ...)
plot(30, .05, xlim=c(17, 43), ylim=c(.2, 0.37),
type="n", xlab="AGE", ylab="BATTING AVERAGE",
main="Career Trajectory Estimates")
sapply(1:10, graph.est, OUT, col="blue")
[,1] [,2] [,3] [,4] [,5] [,6]
x Numeric,101 Numeric,101 Numeric,101 Numeric,101 Numeric,101 Numeric,101
y Numeric,101 Numeric,101 Numeric,101 Numeric,101 Numeric,101 Numeric,101
[,7] [,8] [,9] [,10]
x Numeric,101 Numeric,101 Numeric,101 Numeric,101
y Numeric,101 Numeric,101 Numeric,101 Numeric,101
# overlay posterior means from exchangeable model in red
POST.MEANS <- cbind(1:10, t(post.means), OUT[, 5:6])
sapply(1:10, graph.est, POST.MEANS, col="red")
[,1] [,2] [,3] [,4] [,5] [,6]
x Numeric,101 Numeric,101 Numeric,101 Numeric,101 Numeric,101 Numeric,101
y Numeric,101 Numeric,101 Numeric,101 Numeric,101 Numeric,101 Numeric,101
[,7] [,8] [,9] [,10]
x Numeric,101 Numeric,101 Numeric,101 Numeric,101
y Numeric,101 Numeric,101 Numeric,101 Numeric,101
# add legend
legend(37, .36, legend=c("Individual", "Exchangeable"),
lty=1, col=c("blue", "red"))
# compare estimates of peak age
Individual.Peak.Age <- -OUT[, 3] / OUT[, 4] / 2
Exchangeable.Peak.Age <- -POST.MEANS[, 3] / POST.MEANS[, 4] / 2
par(yaxt="n")
plot(25, 0.5, xlim=c(18, 34), ylim=c(0, 1.5),
type="n", xlab="AGE", ylab="",
main="Estimates of Peak Age")
for(j in 1:10)
lines(c(Individual.Peak.Age[j], Exchangeable.Peak.Age[j]),
c(1.0, 0.5), type="b")
text(25, 1.2, "INDIVIDUAL ESTIMATE", cex=2)
text(25, 0.3, "EXCHANGEABLE ESTIMATE", cex=2)
