## Run height-earnings model
## source("http://web.as.uky.edu/statistics/users/pbreheny/701/S13/notes/4-2c.R")

## Older hispanic comparison
j <- 2
k <- 3
ind <- eid==j & aid==k
fit.mle <- lm(y~Height, subset=ind)
exp(5*confint(fit.mle, 2))

## CI plot
A <- fit$BUGSoutput$summary[1:12,c(5,3,7)]
B <- fit$BUGSoutput$summary[13:24,c(5,3,7)]
rownames(A) <- rownames(B) <- paste(rep(levels(Data$Ethnicity), 3), rep(levels(Age), rep(4,3)))
ord <- c(1,3,2,4) + rep(c(0,4,8), c(4,4,4))
CIplot(exp(A)[ord,], xlab="Average earnings for an average-height (5'7) person")
CIplot(exp(5*B[ord,]), xlab="Earnings ratio for a 5-inch difference")
abline(v=1, col="gray75")

## Weighted average
w <- prop.table(table(Data$Ethnicity))
b.weighted <- sweep(b, 2, w, "*")
b.avg <- apply(b.weighted, c(1,3), sum)
apply(b.avg, 2, mean)
apply(b.avg, 2, sd)
exp(5*psm(b.avg))

## Weighted average, incorporating uncertainty in w
E <- table(Data$Ethnicity)
pi <- rdirichlet(nrow(b), E+1)
b.avg <- apply(b, 3, function(b) apply(b*pi, 1, sum))
exp(5*psm(b.avg))