data(race)
if (FALSE) {
## Multiple chain MCMC list experiment regression
## starts with overdispersed MLE starting values
## Multiple item two level hierarchical model - varying intercepts
mle.estimates.multi <- ictreg(y ~ male + college, data = multi,
constrained = TRUE)
draws <- mvrnorm(n = 3, mu = coef(mle.estimates.multi),
Sigma = vcov(mle.estimates.multi) * 9)
bayesDraws.1 <- ictregBayesHier(y ~ male + college,
formula.level.2 = ~ 1,
delta.start.level.1 = list(draws[1, 8:9], draws[1, 2:3], draws[1, 5:6]),
data = multi, treat = "treat",
delta.tune = list(rep(0.005, 2), rep(0.05, 2), rep(0.05, 2)),
alpha.tune = rep(0.001, length(unique(multi$state))),
J = 3, group.level.2 = "state",
n.draws = 100000, burnin = 50000, thin = 100)
bayesDraws.2 <- ictregBayesHier(y ~ male + college,
formula.level.2 = ~ 1,
delta.start.level.1 = list(draws[2, 8:9], draws[2, 2:3], draws[2, 5:6]),
data = multi, treat = "treat",
delta.tune = list(rep(0.005, 2), rep(0.05, 2), rep(0.05, 2)),
alpha.tune = rep(0.001, length(unique(multi$state))),
J = 3, group.level.2 = "state",
n.draws = 100000, burnin = 50000, thin = 100)
bayesDraws.3 <- ictregBayesHier(y ~ male + college,
formula.level.2 = ~ 1,
delta.start.level.1 = list(draws[3, 8:9], draws[3, 2:3], draws[3, 5:6]),
data = multi, treat = "treat",
delta.tune = list(rep(0.005, 2), rep(0.05, 2), rep(0.05, 2)),
alpha.tune = rep(0.001, length(unique(multi$state))),
J = 3, group.level.2 = "state",
n.draws = 100000, burnin = 50000, thin = 100)
bayesHierTwoLevel <- as.list(bayesDraws.1, bayesDraws.2, bayesDraws.3)
summary(bayesHierTwoLevel)
## Multiple item two level hierarchical model - including covariates
mle.estimates.multi <- ictreg(y ~ male + college, data = multi,
constrained = TRUE)
draws <- mvrnorm(n = 3, mu = coef(mle.estimates.multi),
Sigma = vcov(mle.estimates.multi) * 9)
bayesDraws.1 <- ictregBayesHier(y ~ male + college,
formula.level.2 = ~ age,
delta.start.level.1 = list(draws[1, 8:9], draws[1, 2:3], draws[1, 5:6]),
data = multi, treat = "treat",
delta.tune = list(rep(0.005, 2), rep(0.05, 2), rep(0.05, 2)),
alpha.tune = rep(0.001, length(unique(multi$state))),
J = 3, group.level.2 = "state",
n.draws = 100000, burnin = 50000, thin = 100)
bayesDraws.2 <- ictregBayesHier(y ~ male + college,
formula.level.2 = ~ age,
delta.start.level.1 = list(draws[2, 8:9], draws[2, 2:3], draws[2, 5:6]),
data = multi, treat = "treat",
delta.tune = list(rep(0.005, 2), rep(0.05, 2), rep(0.05, 2)),
alpha.tune = rep(0.001, length(unique(multi$state))),
J = 3, group.level.2 = "state",
n.draws = 100000, burnin = 50000, thin = 100)
bayesDraws.3 <- ictregBayesHier(y ~ male + college,
formula.level.2 = ~ age,
delta.start.level.1 = list(draws[3, 8:9], draws[3, 2:3], draws[3, 5:6]),
data = multi, treat = "treat",
delta.tune = list(rep(0.005, 2), rep(0.05, 2), rep(0.05, 2)),
alpha.tune = rep(0.001, length(unique(multi$state))),
J = 3, group.level.2 = "state",
n.draws = 100000, burnin = 50000, thin = 100)
bayesHierTwoLevel <- as.list(bayesDraws.1, bayesDraws.2, bayesDraws.3)
summary(bayesHierTwoLevel)
}
Run the code above in your browser using DataLab