# \donttest{
# DATA
data(healthsurvey)
# the order of response levels decreases from the best health to
# the worst health; hence the hopit() parameter decreasing.levels
# is set to TRUE
levels(healthsurvey$health)
# fit a model
model1 <- hopit(latent.formula = health ~ hypertension + high_cholesterol +
heart_attack_or_stroke + poor_mobility + very_poor_grip +
depression + respiratory_problems +
IADL_problems + obese + diabetes + other_diseases,
thresh.formula = ~ sex + ageclass + country,
decreasing.levels = TRUE,
control = list(trace = FALSE),
data = healthsurvey)
# Example 1 ---------------------
# bootstrapping cut-points
# a function to be bootstrapped
cutpoints <- function(model) getCutPoints(model)$cutpoints
B <- boot_hopit(model = model1, func = cutpoints, nboot = 100)
# calculate lower and upper bounds using the percentile method
cutpoints.CI <- percentile_CI(B)
# print estimated cutpoints and their confidence intervals
cutpoints(model1)
cutpoints.CI
# Example 2 ---------------------
# bootstrapping differences in health levels
# a function to be bootstrapped
diff_BadHealth <- function(model) {
hl <- getLevels(model = model, formula=~ sex + ageclass, sep=' ')
hl$original[,1] + hl$original[,2] - hl$adjusted[,1]- hl$adjusted[,2]
}
# estimate the difference
est.org <- diff_BadHealth(model = model1)
# perform the bootstrap
B <- boot_hopit(model = model1, func = diff_BadHealth, nboot = 100)
# calculate lower and upper bounds using the percentile method
est.CI <- percentile_CI(B)
# plot the difference and its (asymmetrical) confidence intervals
pmar <- par('mar'); par(mar = c(9.5,pmar[2:4]))
m <- max(abs(est.CI))
pos <- barplot(est.org, names.arg = names(est.org), las = 3,
ylab = 'Original - Adjusted',
ylim=c(-m, m), density = 20, angle = c(45, -45),
col = c('blue', 'orange'))
for (k in seq_along(pos)) lines(c(pos[k,1],pos[k,1]),
est.CI[,k], lwd = 2, col = 2)
abline(h = 0); box(); par(mar = pmar)
# }
Run the code above in your browser using DataLab