data(race)
race.south <- race.nonsouth <- race
race.south[, "south"] <- 1
race.nonsouth[, "south"] <- 0
if (FALSE) {
# Fit EM algorithm ML model with constraint
ml.constrained.results <- ictreg(y ~ south + age + male + college,
data = race, treat = "treat", J=3, method = "ml",
overdispersed = FALSE, constrained = TRUE)
# Calculate average predictions for respondents in the South
# and the the North of the US for the MLE model, replicating the
# estimates presented in Figure 1, Imai (2011)
avg.pred.south.mle <- predict(ml.constrained.results,
newdata = race.south, avg = TRUE, interval = "confidence")
avg.pred.nonsouth.mle <- predict(ml.constrained.results,
newdata = race.nonsouth, avg = TRUE, interval = "confidence")
# A plot of a single estimate and its confidence interval
plot(avg.pred.south.mle, labels = "South")
# A plot of the two estimates and their confidence intervals
# use c() to combine more than one predict object for plotting
plot(c(avg.pred.south.mle, avg.pred.nonsouth.mle), labels = c("South", "Non-South"))
# The difference option can also be used to simultaneously
# calculate separate estimates of the two sub-groups
# and the estimated difference. This can also be plotted.
avg.pred.diff.mle <- predict(ml.constrained.results,
newdata = race.south, newdata.diff = race.nonsouth,
se.fit = TRUE, avg = TRUE, interval="confidence")
plot(avg.pred.diff.mle, labels = c("South", "Non-South", "Difference"))
# Social desirability bias plots
# Estimate logit for direct sensitive question
data(mis)
mis.list <- subset(mis, list.data == 1)
mis.sens <- subset(mis, sens.data == 1)
# Fit EM algorithm ML model
fit.list <- ictreg(y ~ age + college + male + south,
J = 4, data = mis.list, method = "ml")
# Fit logistic regression with directly-asked sensitive question
fit.sens <- glm(sensitive ~ age + college + male + south,
data = mis.sens, family = binomial("logit"))
# Predict difference between response to sensitive item
# under the direct and indirect questions (the list experiment).
# This is an estimate of the revealed social desirability bias
# of respondents. See Blair and Imai (2010).
avg.pred.social.desirability <- predict(fit.list,
direct.glm = fit.sens, se.fit = TRUE)
plot(avg.pred.social.desirability)
}
Run the code above in your browser using DataLab