Nothing
#' @keywords internal
.performance <- function(data, estvarname, se, true, empse_ref = NULL, rho = NULL, ci.limits, df, control) {
### Adjust 'true' to use internal column if not NULL
if (!is.null(true)) true <- ":true"
### Make object to return
obj <- list()
### Compute performance measures
# Number of non-missing point estimates, standard errors:
if (!is.null(se)) {
nsim <- sum(!is.na(data[[estvarname]]) & !is.na(data[[se]]), na.rm = control$na.rm)
} else {
nsim <- sum(!is.na(data[[estvarname]]), na.rm = control$na.rm)
}
# Mean, median and variance of betas
theta_mean <- mean(data[[estvarname]], na.rm = control$na.rm)
theta_median <- stats::median(data[[estvarname]], na.rm = control$na.rm)
# theta_var <- stats::var(data[[estvarname]], na.rm = control$na.rm)
# Mean, median and average of ses
if (!is.null(se)) {
se2_mean <- mean(data[[se]]^2, na.rm = control$na.rm)
se2_median <- stats::median(data[[se]]^2, na.rm = control$na.rm)
se2_var <- stats::var(data[[se]]^2, na.rm = control$na.rm)
}
# Bias
if (!is.null(true)) {
bias <- 1 / nsim * sum(data[[estvarname]] - data[[":true"]], na.rm = control$na.rm)
}
# Relative bias
if (!is.null(true)) {
rbias <- 1 / nsim * sum((data[[estvarname]] - data[[":true"]]) / data[[":true"]], na.rm = control$na.rm)
}
# Empirical standard error
empse <- sqrt(1 / (nsim - 1) * sum((data[[estvarname]] - mean(data[[estvarname]], na.rm = control$na.rm))^2, na.rm = control$na.rm))
# Mean squared error
if (!is.null(true)) {
mse <- 1 / nsim * sum((data[[estvarname]] - data[[":true"]])^2, na.rm = control$na.rm)
}
# Relative change in precision
if (!is.null(empse_ref) & !is.null(rho)) {
relprec <- 100 * ((empse_ref / empse)^2 - 1)
} else {
relprec <- NA
}
names(relprec) <- NULL
# Model-based standard error
if (!is.null(se)) modelse <- sqrt(1 / nsim * sum(data[[se]]^2, na.rm = control$na.rm))
# Relative error in model-based standard error
if (!is.null(se)) relerror <- 100 * (modelse / empse - 1)
# Compute critical values (for coverage) from either a normal or a t distribution
if (!is.null(df)) {
crit <- stats::qt(1 - (1 - control$level) / 2, df = data[[df]])
} else {
crit <- stats::qnorm(1 - (1 - control$level) / 2)
}
# Compute critical value (for power) from either a normal or a t distribution
crit_power <- ifelse(is.null(control$power_df), stats::qnorm(1 - (1 - control$level) / 2), stats::qt(1 - (1 - control$level) / 2, df = control$power_df))
# Coverage of a nominal (1 - level)% confidence interval
if (!is.null(true) & !is.null(se)) {
if (is.null(ci.limits)) {
cover <- 1 / nsim * sum(data[[":true"]] >= data[[estvarname]] - crit * data[[se]] & data[[":true"]] <= data[[estvarname]] + crit * data[[se]], na.rm = control$na.rm)
} else {
if (is.character(ci.limits)) {
cover <- 1 / nsim * sum(data[[":true"]] >= data[[ci.limits[1]]] & data[[":true"]] <= data[[ci.limits[2]]], na.rm = control$na.rm)
} else if (is.numeric(ci.limits)) {
data[["lower"]] <- ci.limits[1]
data[["upper"]] <- ci.limits[2]
cover <- 1 / nsim * sum(data[[":true"]] >= data[["lower"]] & data[[":true"]] <= data[["upper"]], na.rm = control$na.rm)
}
}
}
# Bias-corrected coverage of a nominal (1 - level)% confidence interval
if (!is.null(se)) {
if (is.null(ci.limits)) {
becover <- 1 / nsim * sum(mean(data[[estvarname]], na.rm = control$na.rm) >= data[[estvarname]] - crit * data[[se]] & mean(data[[estvarname]], na.rm = control$na.rm) <= data[[estvarname]] + crit * data[[se]], na.rm = control$na.rm)
} else {
if (is.character(ci.limits)) {
becover <- 1 / nsim * sum(mean(data[[estvarname]], na.rm = control$na.rm) >= data[[ci.limits[1]]] & mean(data[[estvarname]], na.rm = control$na.rm) <= data[[ci.limits[2]]], na.rm = control$na.rm)
} else if (is.numeric(ci.limits)) {
data[["lower"]] <- ci.limits[1]
data[["upper"]] <- ci.limits[2]
becover <- 1 / nsim * sum(mean(data[[estvarname]], na.rm = control$na.rm) >= data[[estvarname]] - data[["lower"]] & mean(data[[estvarname]], na.rm = control$na.rm) <= data[["upper"]], na.rm = control$na.rm)
}
}
}
# Power of a significance test at the `level` level
if (!is.null(se)) power <- 1 / nsim * sum(abs(data[[estvarname]]) >= crit_power * data[[se]], na.rm = control$na.rm)
### Compute Monte Carlo SEs if requested:
if (control$mcse) {
if (!is.null(true)) {
bias_mcse <- sqrt(1 / (nsim * (nsim - 1)) * sum((data[[estvarname]] - mean(data[[estvarname]], na.rm = control$na.rm))^2, na.rm = control$na.rm))
rbias_i <- (data[[estvarname]] - data[[":true"]]) / data[[":true"]]
rbias_mcse <- sd(rbias_i) / sqrt(nsim)
}
empse_mcse <- empse / sqrt(2 * (nsim - 1))
if (!is.null(true)) {
mse_mcse <- sqrt(sum(((data[[estvarname]] - data[[":true"]])^2 - mse)^2, na.rm = control$na.rm) / (nsim * (nsim - 1)))
}
if (!is.null(empse_ref) & !is.null(rho)) {
relprec_mcse <- 200 * (empse_ref / empse)^2 * sqrt((1 - rho^2) / (nsim - 1))
} else {
relprec_mcse <- NA
}
names(relprec_mcse) <- NULL
if (!is.null(se)) modelse_mcse <- sqrt(se2_var / (4 * nsim * modelse^2))
if (!is.null(se)) relerror_mcse <- 100 * (modelse / empse) * sqrt(se2_var / (4 * nsim * modelse^4) + 1 / (2 * (nsim - 1)))
if (!is.null(true) & !is.null(se)) cover_mcse <- sqrt(cover * (1 - cover) / nsim)
if (!is.null(se)) becover_mcse <- sqrt(becover * (1 - becover) / nsim)
if (!is.null(se)) power_mcse <- sqrt(power * (1 - power) / nsim)
}
# Easy thing to do is to add what can't be computed as null
if (is.null(true)) {
bias <- rbias <- cover <- mse <- NULL
bias_mcse <- rbias_mcse <- cover_mcse <- mse_mcse <- NULL
}
if (is.null(se)) {
se2_mean <- se2_median <- se2_var <- modelse <- relerror <- cover <- becover <- power <- NULL
modelse_mcse <- relerror_mcse <- cover_mcse <- becover_mcse <- power_mcse <- NULL
}
### Assemble object to return
obj$stat <- c("nsim", "thetamean", "thetamedian", "se2mean", "se2median", "bias", "rbias", "empse", "mse", "relprec", "modelse", "relerror", "cover", "becover", "power")
if (is.null(true)) obj$stat <- obj$stat[!(obj$stat %in% c("bias", "rbias", "cover", "mse"))]
if (is.null(se)) obj$stat <- obj$stat[!(obj$stat %in% c("se2mean", "se2median", "modelse", "relerror", "cover", "becover", "power"))]
obj$est <- c(nsim, theta_mean, theta_median, se2_mean, se2_median, bias, rbias, empse, mse, relprec, modelse, relerror, cover, becover, power)
if (control$mcse) {
obj$mcse <- c(rep(NA, ifelse(is.null(se), 3, 5)), bias_mcse, rbias_mcse, empse_mcse, mse_mcse, relprec_mcse, modelse_mcse, relerror_mcse, cover_mcse, becover_mcse, power_mcse)
}
obj <- as.data.frame(obj, stringsAsFactors = FALSE)
### Return object
return(obj)
}
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.