Nothing
# INTERNAL : updating a model using new set ofvariable coefficients
#
# @details
# Gradient and Hessian of log-likelihood function as well as variance-covariance matrix are not re-calcualted in this function!
# @param model a fitted \code{hopit} model.
# @param newregcoef a new set of variable coefficients. A vector of length lower or equal to the length of \code{model$coef}.
# The coeficients are replaced starting from the first one.
# @param data a data used to fit original model.
# @author Maciej J. Danko
#' @noRd
#' @keywords internal
update.latent <- function(model, newregcoef, data){
coefnames <- names(model$coef)
thresh.names <- colnames(model$thresh.mm)
model$coef[seq_along(newregcoef)] <- newregcoef
class(model) <- 'hopit'
names(model$coef) <- coefnames
colnames(model$thresh.mm) <- thresh.names
p <- hopit_ExtractParameters(model)
model$alpha <- hopit_Threshold(p$thresh.lambda, p$thresh.gamma, model)
model$y_latent_i <- hopit_Latent(p$latent.params, model)
model$maxobservedlatentrange <- range(model$y_latent_i)
model$Ey_i <- factor(colSums(sapply(1L : model$N, function(k) model$alpha[k,]<model$y_latent_i[k])),levels=1L:model$J)
levels(model$Ey_i) <- levels(model$y_i)
model$coef.ls <- p
model$deviance <- -2 * model$LL
if (!length(model$design)) {
k <- 2
model$AIC <- model$deviance + k * (length(model$coef.ls$latent.params)+
length(model$coef.ls$thresh.lambda)+
length(model$coef.ls$thresh.gamma)+model$hasdisp)
} else model$AIC <- NA
return(model)
}
#' Bootstrapping hopit model
#'
#' \code{boot_hopit} performs the bootstrap of a function dependent on a fitted model.
#' In each of the bootstrap repetitions, a set of new model coefficients is drawn from the multivariate normal distribution,
#' assuming the originally estimated model coefficients (see \code{\link{coef.hopit}})
#' as a mean and using the model variance-covariance matrix (see \code{\link{vcov.hopit}}).
#' The drawn coefficients are then used to calculate the measure of interest using a function delivered by the \code{func} parameter.
#' @param model a fitted \code{hopit} model.
#' @param data data used to fit the model.
#' @param func a function to be bootstrapped of the form \code{func(model, ...)}.
#' @param nboot a number of bootstrap replicates.
#' @param unlist a logical indicating whether to unlist the boot object.
#' @param boot.only.latent a logical indicating whether to perform the bootstrap on latent variables only.
#' @param parallel.flag a logical if to use parallel computations.
#' @param parallel.nb_cores number of cores (<= number of CPU cores on the current host).
#' @param parallel.packages list of packages needed to run "func".
#' @param parallel.variables list of global variables and functions needed to run "func".
#' @param robust.vcov see \code{\link{vcov.hopit}}.
#' @param ... other parameters passed to the \code{func}.
#' @importFrom MASS mvrnorm
#' @import parallel
#' @author Maciej J. Danko
#' @return a list with bootstrapped elements.
#' @export
#' @seealso \code{\link{percentile_CI}}, \code{\link{getLevels}}, \code{\link{getCutPoints}}, \code{\link{latentIndex}}, \code{\link{standardiseCoef}}, \code{\link{hopit}}.
#' @examples
#' \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)
#' }
boot_hopit<-function(model, func, data=model$frame, nboot = 500, unlist = TRUE,
boot.only.latent = TRUE, parallel.flag=FALSE, parallel.nb_cores=NULL,
parallel.packages = NULL, parallel.variables = NULL,
robust.vcov, ...){
if (parallel.flag && !length(parallel.nb_cores)){
chk <- Sys.getenv("_R_CHECK_LIMIT_CORES_", "")
if (nzchar(chk) && chk == "TRUE") {
parallel.nb_cores <- 2L
} else {
parallel.nb_cores <- parallel::detectCores()
}
}
if (missing(robust.vcov)) {
if (length(model$design)) robust.vcov <- FALSE else robust.vcov <- TRUE
}
data <- model$na.action(data)
if (model$control$transform.latent != 'none')
data <- transform.data(model$latent.formula, data,
model$control$transform.latent)
if (model$control$transform.thresh != 'none')
data <- transform.data(model$thresh.formula, data,
model$control$transform.thresh)
VCOV <- vcov.hopit(model, robust.vcov)
if (boot.only.latent) N <- seq_len(model$parcount[1]) else N <- nrow(VCOV)
if (length(VCOV) < 2) stop(call. = NULL, hopit_msg(23))
bootsample <- MASS::mvrnorm(nboot, mu = model$coef[N], Sigma = VCOV[N,N])
if (parallel.flag){
cl<-parallel::makeCluster(parallel.nb_cores)
if (length(parallel.variables)) parallel::clusterExport(cl,parallel.variables, envir = .GlobalEnv)
parallel::clusterExport(cl, c("bootsample","nboot","N","model","data","func"), envir = environment())
parallel::clusterEvalQ(cl, {library("hopit")})
if (length(parallel.packages)) for (k in seq_along(length(parallel.packages)))
parallel::clusterEvalQ(cl, {eval(parse(text=paste0('library(',parallel.packages[k],')')))})
boots <- parallel::parLapply(cl, seq_len(nboot), function(k)
func(model = update.latent(model, bootsample[k,N],data = data), ...))
parallel::stopCluster(cl)
} else {
boots <- lapply(seq_len(nboot), function(k)
func(model = update.latent(model, bootsample[k,N],data = data), ...))
}
if (unlist[1]) {
boots <- sapply(boots,'[')
class(boots) <- 'hopit.boot'
} else class(boots) <- c('hopit.boot', 'list')
boots
}
# boot_hopit<-function(model, func, data=model$frame, nboot = 500, unlist = TRUE,
# boot.only.latent = TRUE, robust.vcov, ...){
# if (missing(robust.vcov)) {
# if (length(model$design)) robust.vcov <- FALSE else robust.vcov <- TRUE
# }
# data <- model$na.action(data)
# if (model$control$transform.latent != 'none')
# data <- transform.data(model$latent.formula, data,
# model$control$transform.latent)
# if (model$control$transform.thresh != 'none')
# data <- transform.data(model$thresh.formula, data,
# model$control$transform.thresh)
#
# VCOV <- vcov.hopit(model, robust.vcov)
# if (boot.only.latent) N <- seq_len(model$parcount[1]) else N <- nrow(VCOV)
# if (length(VCOV) < 2) stop(call. = NULL, hopit_msg(23))
# bootsample <- MASS::mvrnorm(nboot, mu = model$coef[N], Sigma = VCOV[N,N])
# boots <- lapply(seq_len(nboot), function(k)
# func(model = update.latent(model,
# bootsample[k,N],data = data), ...))
# if (unlist[1]) {
# boots <- sapply(boots,'[')
# class(boots) <- 'hopit.boot'
# } else class(boots) <- c('hopit.boot', 'list')
# boots
# }
#' Calculating the confidence intervals of the bootstrapped function using the percentile method
#'
#' Calculate the confidence intervals of the bootstrapped function using the percentile method.
#' @param boot a matrix or a list of vectors with bootstrapped elements. If it is list, then each element of the list is one replication.
#' @param alpha a significance level.
#' @param bounds which bounds to return; one of \code{"both"}, \code{"lo"}, \code{"up"}.
#' @author Maciej J. Danko
#' @seealso \code{\link{boot_hopit}}, \code{\link{getLevels}}, \code{\link{getCutPoints}}, \code{\link{latentIndex}}, \code{\link{standardiseCoef}}, \code{\link{hopit}}.
#' @export
#' @examples
#' # see examples in boot_hopit() function.
percentile_CI <- function(boot, alpha = 0.05, bounds = c('both', 'lo', 'up')){
bounds <- tolower(bounds[1])
if (inherits(boot,'list')) boot <- sapply(boot,'[')
probs <- switch(bounds,
up = 1-alpha/2,
lo = alpha/2,
both = c(alpha/2, 1-alpha/2))
apply(boot, 1, stats::quantile, probs = probs)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.