Nothing
R/SMR.R
In ems: Epimed Solutions Collection for Data Editing, Analysis, and Benchmark of Health Units
Defines functions
forest.SMR
SMR.table
SMR
Documented in
forest.SMR
SMR
SMR.table
#' Standardized Mortality Ratio (SMR)
#'
#' @name SMR
#'
#' @description Calculates the standardized mortality ratio and its confidence interval. SMR, for a group, is defined as the ratio of the observed deaths in this group and the sum of the predicted individual probabilities of death by any model (expected deaths).
#'
#' \code{SMR.table} estimates at once the overall SMR and the SMR across several groups, e.g. ICU units or clinical characteristics. The \code{SMR.table} can be ordered by the SMR estimate or its confidence intervals, facilitating the comparison of the units ranks.
#'
#' \code{forest.SMR} shows the \code{SMR.table} output as a forest plot. The plot opens two windows and plot at the left side the values from the \code{SMR.table} and at the right side the points and lines graphically representing each SMR and its confidence interval.
#'
#' @param obs.var Observed death. Accepted values are 0 (absence) or 1 (presence) in a vector. For \code{SMR.table} it must be a character indicating the name of the variable in the data.
#'
#' @param pred.var Death individual predictions (ranging from 0 to 1) in a vector. For \code{SMR.table} it must be a character indicating the name of the variable in the data.
#'
#' @param digits Number of digits for rounding the output.
#'
#' @param ci.method Method to estimate the confidence interval. "Hosmer" (default) or "Byar" are acceptable values.
#'
#' @param ci.level Level of the confidence interval. Default is 0.95.
#'
#' @param data For \code{SMR.table}, a dataset where pred.var, obs.var and group.var are in.
#'
#' @param group.var For \code{SMR.table}, this is a character vector indicating the name(s) of the variable(s) in the data that will form the groups where SMR will be calculated. The variables must be factors.
#'
#' @param use.label Logical. Default is FALSE. For \code{SMR.table} this option will replace the variables names by its labels in var.labels argument.
#'
#' @param var.labels A character vector with variables labels. The default is to replace the variable name by the label stored at attr(data, "var.labels"). But one may specify labels directly.
#'
#' @param reorder Default is "no". Possible values are: "no", "SMR","lower.Cl", and "upper.Cl". It will make the \code{SMR.table} to be ordered within each varibale by its original order, or by SMR order, or by lower.Cl order, or by upper.Cl.
#'
#' @param decreasing Logical. When 'reorderd' is TRUE, should the order be decreasing or incresing? See \code{\link[base]{order}}
#'
#' @param x For the \code{forest.SMR} this is the output of \code{SMR.table}.
#'
#' @param mar1,mar.SMR Values to set the margins (mar parameter) of left and right windows. See \code{\link[graphics]{par}}
#'
#' @param overall.arg A list of arguments passed to \code{\link[graphics]{text}} for plotting the overall label. Internally, 'y' coordinate is replaced.
#'
#' @param NOE.overall.args A list of arguments passed to \code{\link[graphics]{text}} for plotting the overall N (number of observations), O (observed deaths) and E (expected deaths). Internally, 'labels' and 'y' arguments are replaced.
#'
#' @param var.labels.arg A list of arguments passed to \code{\link[graphics]{text}} for plotting the variables labels. Internally, 'y' coordinate is replaced.
#'
#' @param cat.labels.arg A list of arguments passed to \code{\link[graphics]{text}} for plotting the categories labels. Internally, 'y' coordinate is replaced.
#'
#' @param N.values.arg A list of arguments passed to \code{\link[graphics]{text}} for plotting the values of N (number of observations) of each subgroup. Internally, the arguments 'label' and 'y' coordinate are replaced.
#'
#' @param O.values.arg A list of arguments passed to \code{\link[graphics]{text}} for plotting the values of Observed deaths of each subgroup. Internally, the arguments 'label' and 'y' coordinate are replaced.
#'
#' @param E.values.arg A list of arguments passed to \code{\link[graphics]{text}} for plotting the values of Expected deaths of each subgroup. Internally, the arguments 'label' and 'y' coordinate are replaced.
#'
#' @param NOE.head.arg A list of arguments passed to \code{\link[graphics]{text}} for plotting the labels of the columns N, E and O on the top of the graph. Internally, the 'x' and 'y' coordinates are replaced. The x coordinates are taken from the x in \code{N.values.arg}, \code{O.values.arg} and \code{E.values.arg}.
#'
#' @param Overall.seg.arg A list of arguments passed to \code{\link[graphics]{segments}} for plotting the lines corresponding to overall SMR confidence intervals. Internally, 'x' and 'y' coordinates are replaced.
#'
#' @param Overall.p.arg A list of arguments passed to \code{\link[graphics]{points}} for plotting the points corresponding to overall SMR. Internally, 'x' and 'y' coordinates are replaced.
#'
#' @param Overall.est.arg A list of arguments passed to \code{\link[graphics]{text}} for plotting the overall SMR beside the graph. Internally, 'y' coordinate and 'label' argument are replaced.
#'
#' @param cat.seg.arg A list of arguments passed to \code{\link[graphics]{segments}} for plotting the lines corresponding to SMR confidence intervals for all groups. Internally, 'x' and 'y' coordinates are replaced.
#'
#' @param cat.p.arg A list of arguments passed to \code{\link[graphics]{points}} for plotting the points corresponding to all categoreis SMR. Internally, 'x' and 'y' coordinates are replaced.
#'
#' @param cat.est.arg A list of arguments passed to \code{\link[graphics]{text}} for plotting the categories SMR beside the graph. Internally, 'y' coordinate and 'label' arguments are replaced.
#'
#' @param SMR.head.arg A list of arguments passed to \code{\link[graphics]{text}} for plotting the label of the SMR column on the top of the graph. Internally, the 'y' coordinate is replaced.
#'
#' @param smr.xlab Label of the x axis. Default is "Standardized Mortality Ratio".
#'
#' @param smr.xlim Limits of x axis of the \code{forest.SMR} plot. Default is "auto", which internally will pick the highest values of all upper.Cl and the lowest lower.Cl. Besides "auto", only a vector of 2 numbers is valid, and will be passed to \code{\link[graphics]{plot.default}}.
#'
#' @param grid Logical. If TRUE (default), it will draw a grid with the \code{\link[graphics]{grid}} default arguments.
#'
#' @return
#' If SMR, then:
#' \itemize{
#' \item \code{N} Number of subjects analyzed.
#' \item\code{O} Observed number of deaths.
#' \item \code{E} Expected number of deaths.
#' \item \code{SMR} Standardized mortality ratio.
#' \item \code{lower.Cl} lower confidence limit.
#' \item \code{upper.Cl} upper confidence limit.
#' }
#'
#' If SMR.table, then a data.frame with the same information as above, and the
#' additional information is returned: "Variables" (variables names), "Levels" (variables levels).
#'
#' If forest.SMR, then a plot is returned.
#'
#' @author Lunna Borges and Pedro Brasil
#'
#' @seealso \code{\link{SRU}}, \code{\link{reclass}}, \code{\link{funnel}}
#'
#' @examples
#' # Loading a example data
#' data(icu)
#'
#' # Setting variable labels to data
#' attr(icu, "var.labels")[match(c("Unit", "IsMechanicalVentilation1h",
#' "AdmissionTypeName_pri","Vasopressors_D1"), names(icu))] <-
#' c("ICU unit","Mechanichal ventilation","Admission type","Vasopressors at admission")
#'
#' # Some editing
#' icu$Saps3DeathProbabilityStandardEquation <- icu$Saps3DeathProbabilityStandardEquation /100
#' icu$IsMechanicalVentilation1h <- as.factor(ifelse(icu$IsMechanicalVentilation1h == 1, "Yes", "No"))
#' icu$AdmissionTypeName_pri <- as.factor(icu$AdmissionTypeName_pri)
#' levels(icu$AdmissionTypeName_pri) <- c("Clinical","Elective surgery", "Urgent surgery")
#' icu$Vasopressors_D1 <- as.factor(ifelse(icu$Vasopressors_D1 == 1, "Yes", "No"))
#'
#' # The overall SMR for the whole sample
#' SMR(icu$UnitDischargeName, icu$Saps3DeathProbabilityStandardEquation)
#'
#' # The overall SMR and for some subgroups
#' x <- SMR.table(data = icu, obs.var = "UnitDischargeName",
#' pred.var = "Saps3DeathProbabilityStandardEquation",
#' group.var = c( "IsMechanicalVentilation1h",
#' "AdmissionTypeName_pri","Vasopressors_D1"),
#' reorder = "no",
#' decreasing = TRUE,
#' use.label = TRUE)
#' x
#'
#' # A forest plot for all groups SMR (resize the window may be required)
#' forest.SMR(x, digits = 2)
#'
#' # The same thing but reordering the categories
#' x <- SMR.table(data = icu, obs.var = "UnitDischargeName",
#' pred.var = "Saps3DeathProbabilityStandardEquation",
#' group.var = c( "IsMechanicalVentilation1h",
#' "AdmissionTypeName_pri", "Vasopressors_D1"),
#' reorder = "SMR",
#' decreasing = TRUE,
#' use.label = TRUE)
#' forest.SMR(x, digits = 2)
#'
#' # The overall SMR and for all Units
#' x <- SMR.table(data = icu, obs.var = "UnitDischargeName",
#' pred.var = "Saps3DeathProbabilityStandardEquation",
#' group.var = "Unit",
#' reorder = "no",
#' decreasing = TRUE,
#' use.label = TRUE)
#' x
#'
#' # A forest plot for all Units
#' forest.SMR(x, digits = 2)
#'
#' # The same thing but reordering the categories
#' x <- SMR.table(data = icu, obs.var = "UnitDischargeName",
#' pred.var = "Saps3DeathProbabilityStandardEquation",
#' group.var = "Unit",
#' reorder = "SMR",
#' decreasing = TRUE,
#' use.label = TRUE)
#' forest.SMR(x, digits = 2)
#'
#' rm(x, icu)
#' @references
#'
#' David W. Hosmer and Stanley Lemeshow. Confidence intervals estimates of an index of quality performance basend on logistic regression models. Statistics in Medicine , vol. 14, 2161-2172 (1995)
#'
#' @import stats
#' @export
SMR <- function(obs.var, pred.var, digits = 5, ci.method = c("Hosmer", "Byar"), ci.level = 0.95) {
if (length(obs.var) != length(pred.var)){
stop("Length of pred.var and obs.var differ.")
}
if (any(is.na(obs.var) | is.na(pred.var))) {
stop("There are NAs either at pred.var or obs.var. Either remove the NAs or impute!")
}
if (!is.numeric(pred.var)) {
stop("Predicted death variable must be numeric.")
}
if (any(min(pred.var) < 0 | max(pred.var) > 1)) {
stop("The individual predicted death must range from 0 to 1.")
}
if (!is.numeric(obs.var)) {
stop("Observed death variable must be numeric.")
}
if (any(obs.var != 0 & obs.var != 1)) {
stop("Observed death variable must be coded as 0 and 1.")
}
if (ci.method[1] != "Byar" && ci.method[1] != "Hosmer") {
stop("ci.method must be either 'Byar' or 'Hosmer'.")
}
if (!is.numeric(ci.level) || ci.level > 1 || ci.level < 0) {
stop("ci.level must be numeric between 0 and 1.")
}
N <- length(pred.var)
O <- length(obs.var[which(obs.var == 1)])
E <- sum(pred.var)
SMR <- O / E
if (ci.method[1] == "Byar") {
lowerCL <- O / E * (1 - 1 / (9 * O) - qnorm(1 - (1 - ci.level) / 2) / (3 * sqrt(O))) ^ 3
upperCL <- (O + 1) / E * (1 - (1 / (9 * (O + 1))) + qnorm(1 - (1 - ci.level) / 2) /
(3 * sqrt(O + 1))) ^ 3
}
if (ci.method[1] == "Hosmer") {
V2 <- sum(pred.var * (1 - pred.var))
lowerCL <- exp(log(SMR) - qnorm(1 - (1 - ci.level) / 2) * sqrt(V2 / O ^ 2))
upperCL <- exp(log(SMR) + qnorm(1 - (1 - ci.level) / 2) * sqrt(V2 / O ^ 2))
}
output <- c(N = N, Observed = O, Expected = E, SMR = SMR, lower.Cl = lowerCL,
upper.Cl = upperCL)
output <- round(output, digits = digits)
output
}
#' @rdname SMR
#' @export
SMR.table <- function(data, group.var, obs.var, pred.var, digits = 5, use.label = FALSE, var.labels = attr(data, "var.labels")[match(group.var, names(data))], ci.method = c("Hosmer", "Byar"), ci.level = 0.95, reorder = c("no","SMR","lower.Cl","upper.Cl"), decreasing = FALSE) {
if (any(is.na(match(group.var, names(data))))) {
stop("One or more variables in group var is not a variable of data.")
}
if (!all(sapply(data[ , group.var], is.factor))) {
stop("All variables in group.var must be factors.")
}
if (any(is.na(match(c(pred.var, obs.var), names(data))))) {
stop("Either 'pred.var' or 'obs.var' is not a variable in data.")
}
if (use.label) {
if (is.null(var.labels)) {
stop("Either there is no label in attr(data,'var.labels') or 'var.labels argument was not set.")
}
if (length(group.var) != length(var.labels)){
stop("The number of variables and labels are different.")
}
if (!is.character(var.labels)) {
stop("'var.labels ' must be a character vector.")
}
}
if (reorder[1] != "no" & reorder[1] != "SMR" & reorder[1] != "lower.Cl" & reorder[1] != "upper.Cl") {
stop("'reorder' must be one of 'no','SMR','lower.Cl' or 'upper.Cl'")
}
cont.table <- sapply(1:length(group.var), function(i) table(data[, group.var[i]]))
if (any(unlist(cont.table) == 0)) {
# warning("In SMR analysis, levels were deleted before analysis due to zero observations.")
# which levels were deleted
if(is.matrix(unlist(cont.table))){
deleted_levels <- rownames(unlist(cont.table))[which(unlist(cont.table) == 0)]
} else {
deleted_levels <- names(unlist(cont.table))[which(unlist(cont.table) == 0)]
}
warning(paste("In SMR analysis, level",deleted_levels,"were deleted before analysis due to zero observations.\n"))
data[ , group.var] <- droplevels(data[ , group.var])
}
Variables <- c("Overall", unlist(sapply(1:length(group.var), function(i) rep(group.var[i], nlevels(data[ , group.var[i]])))))
if (length(group.var) == 1) {
Levels <- c(NA, levels(data[ , group.var]))
} else {
Levels <- c(NA, unname(unlist(sapply(data[ , group.var], levels))))
}
x <- data.frame(Variables, Levels, stringsAsFactors = TRUE)
x$N <- NA
x$Observed <- NA
x$Expected <- NA
x$SMR <- NA
x$lower.Cl <- NA
x$upper.Cl <- NA
x[1,3:8] <- SMR(data[ , obs.var], data[, pred.var], digits = digits, ci.method = ci.method, ci.level = ci.level)
for(i in 2:nrow(x)){
cond <- which(data[ , as.character(x[i , "Variables"])] == as.character(x[i , 2]))
x[i,3:8] <- SMR(data[cond , obs.var], data[cond , pred.var], digits = digits, ci.method = ci.method, ci.level = ci.level)
}
x$Variables <- as.character(x$Variables)
if (reorder[1] == "SMR") {
for (i in group.var) {
# i = 1
cond <- grep(i, x$Variables)
y <- x[ cond, ]
y <- y[ order(y$SMR, decreasing = decreasing), ]
x[ cond, ] <- y
}
}
if (reorder[1] == "lower.Cl") {
# i = 2
for (i in group.var) {
cond <- grep(i, x$Variables)
y <- x[ cond, ]
y <- y[ order(y$lower.Cl, decreasing = decreasing), ]
x[ cond, ] <- y
}
}
if (reorder[1] == "upper.Cl") {
for (i in group.var) {
cond <- grep(i, x$Variables)
y <- x[ cond, ]
y <- y[ order(y$upper.Cl, decreasing = decreasing), ]
x[ cond, ] <- y
}
}
if (use.label) {
names(group.var) <- var.labels
for (i in 2:length(x$Variables)) {
x$Variables[i] <- names(group.var)[which(x$Variables[i] == group.var)]
}
}
for (i in length(x$Variables):2) {
if (x$Variables[i - 1] == x$Variables[i]) {x$Variables[i] <- NA}
}
x
}
#' @rdname SMR
#' @import graphics
#' @import grDevices
#' @export
forest.SMR <- function(x,
mar1 = c(5.1, 1, 4.1, 1),
mar.SMR = c(5.1, 7, 4.1, 1),
overall.arg = list(x = .01, font = 2, las = 1, labels = var.labels[1], xpd = NA, adj = 0),
NOE.overall.args = list(x = c(N.values.arg$x, O.values.arg$x, E.values.arg$x), font = 2, las = 1, xpd = NA),
var.labels.arg = list(x = .01, font = 2, las = 1, cex = 1, xpd = NA, adj = 0),
cat.labels.arg = list(x = .1, font = 3, las = 1, cex = .95, col = gray(.4), xpd = NA , adj = 0),
N.values.arg = list(x = .5, col = gray(.4), xpd = NA),
O.values.arg = list(x = .675, col = gray(.4), xpd = NA),
E.values.arg = list(x = .85, col = gray(.4), xpd = NA),
NOE.head.arg= list(font = 2, labels = c("N","O","E"), xpd = NA),
Overall.seg.arg = list(col = "navyblue", xpd = NA, lwd = 2),
Overall.p.arg = list(pch = 23, cex = 2, col = "black", bg = gray(.4), xpd = NA),
Overall.est.arg = list(x = smr.xlim[1] - .06, las = 1, font = 2, xpd = NA, adj = 1),
cat.seg.arg = list(col = "navyblue", xpd = NA, lwd = 2),
cat.p.arg = list(pch = 22 ,cex= 1, col = "black", bg = gray(.4), xpd = NA),
cat.est.arg = list(x = smr.xlim[1] - .06, las = 1, col = gray(.4), xpd = NA, adj = 1),
SMR.head.arg = list(smr.xlim[1] - .06, font = 2, labels = "SMR [95% CIs]", xpd = NA, adj = 1),
smr.xlab = "Standardized Mortality Ratio",
smr.xlim = "auto",
grid = TRUE,
digits = 3){
opar <- par(mfrow = c(1,1), mar = c(5.1, 4.1, 4.1, 2.1))
on.exit(par(opar))
# Separando os valores da medida preincipal
main.pos <- 1
main.smr <- unlist(x[which(x$Variables == "Overall"), c("SMR", "lower.Cl", "upper.Cl")])
main.n <- unlist(x[which(x$Variables == "Overall"), "N"])
main.observed <- round(unlist(x[which(x$Variables == "Overall"), "Observed"]), digits = digits)
main.expected <- round(unlist(x[which(x$Variables == "Overall"), "Expected"]), digits = digits)
# Separando as medidas do SMR para montar o grafico
smr.estimates <- unlist(x$SMR[-1])
smr.ll <- unlist(x$lower.Cl[-1])
smr.ul <- unlist(x$upper.Cl[-1])
# Separando os valores para serem impressos na direita do grafico
N <- unlist(x$N[-1])
observed <- round(unlist(x$Observed[-1]), digits = digits)
expected <- round(unlist(x$Expected[-1]), digits = digits)
# Achandos os limites do grafico no eixo horizontal
if (smr.xlim[1] == "auto"){
smr.xlim <- c(min(unlist(c(main.smr[2], smr.ll)), na.rm = TRUE), max(unlist(c(main.smr[3], smr.ul)), na.rm = TRUE))
}
# Separando os valores das vari?veis e categorias
var.labels <- x$Variables[-which(is.na(x$Variables))]
cat.labels <- x$Levels[-which(is.na(x$Levels))]
# Achando as posi??es no eixo vertical
nLevels <- c(which(!is.na(x$Variables))[-1],(nrow(x) + 1))
for(i in 1:length(nLevels)){
nLevels[i] <- nLevels[i+1] - nLevels[i]
}; rm(i)
nLevels <- nLevels[which(!is.na(nLevels))]
# x$Levels
# i = 3
cat.pos <- seq(1:nLevels[1])
if (length(nLevels) == 1) {
cat.pos <- cat.pos + 2
} else {
for(i in 2:length(nLevels)){
ini <- cat.pos[length(cat.pos)]
cat.pos <- c(cat.pos, (ini + 3):(ini + 2 + nLevels[i]))
}
cat.pos <- cat.pos + 2
}
if (length(nLevels) == 1) {
var.pos <- cat.pos[length(cat.pos)] + 1
} else {
var.pos = nLevels[1] + 1
# i = 2
for(i in 2:length(nLevels)){
var.pos <- c(var.pos, var.pos[i-1] + nLevels[i] + 2)
}
var.pos <- var.pos + 2
}
ylim = c(1, var.pos[length(var.pos)])
# 5.1 4.1 4.1 2.1
par(mfrow = c(1,2))
# Janela com textos com variaveis e categorias
par(mar = mar1)
plot(NA, NA, xlim = c(0,1), ylim = ylim , xlab = "", ylab = "",yaxt = "n", xaxt = "n", xaxs = "i", bty = "n", xpd = NA)
# Efeito principal
overall.arg$y <- main.pos
do.call(text, overall.arg)
# N O E do overall
NOE.overall.args$labels <- c(main.n, main.observed, main.expected)
NOE.overall.args$y <- main.pos
do.call(text, NOE.overall.args)
# Demais vari?veis
var.labels.arg$labels <- var.labels[2:length(var.labels)]
var.labels.arg$y <- var.pos
do.call(text, var.labels.arg)
# Categorias
cat.labels.arg$labels <- cat.labels
cat.labels.arg$y <- cat.pos
do.call(text, cat.labels.arg)
# Colunas com os valores de N O e E
N.values.arg$labels <- N
N.values.arg$y <- cat.pos
do.call(text, N.values.arg)
O.values.arg$labels <- observed
O.values.arg$y <- cat.pos
do.call(text, O.values.arg)
E.values.arg$labels <- expected
E.values.arg$y <- cat.pos
do.call(text, E.values.arg)
# Cabecalho das colunas N O E
NOE.head.arg$y <- ylim[2] + 1
NOE.head.arg$x <- c(N.values.arg$x, O.values.arg$x, E.values.arg$x)
do.call(text, NOE.head.arg)
# SMR
par(mar = mar.SMR)
plot(NA, NA, xlim = smr.xlim, ylim = ylim , xlab = smr.xlab, ylab = "", yaxt = "n", xaxt = "n", xaxs = "i", bty = "n", xpd = NA); axis(1)
if(grid) grid()
# Segmento do IC Efeito principal
Overall.seg.arg$x0 <- main.smr[2]
Overall.seg.arg$y0 <- main.pos
Overall.seg.arg$x1 <- main.smr[3]
Overall.seg.arg$y1 <- main.pos
do.call(segments, Overall.seg.arg)
# Ponto do Efeito principal
Overall.p.arg$x <- main.smr[1]
Overall.p.arg$y <- main.pos
Overall.p.arg$type <- "p"
do.call(points, Overall.p.arg)
# Estimativa do efeito principal
Overall.est.arg$y <- main.pos
Overall.est.arg$labels <- sprintf(paste0("%.",digits,"f [%.",digits,"f ; %.",digits,"f]"),main.smr[1],main.smr[2],main.smr[3])
do.call(text, Overall.est.arg)
# Segmentos das categorias
cat.seg.arg$x0 <- smr.ll
cat.seg.arg$y0 <- cat.pos
cat.seg.arg$x1 <- smr.ul
cat.seg.arg$y1 <- cat.pos
do.call(segments, cat.seg.arg)
# Ponto das estimativas das categorias
cat.p.arg$x <- smr.estimates
cat.p.arg$y <- cat.pos
cat.p.arg$type <- "p"
do.call(points, cat.p.arg)
# Estimativas das categorias
cat.est.arg$y <- cat.pos
cat.est.arg$labels <- sprintf(paste0("%.",digits,"f [%.",digits,"f ; %.",digits,"f]"),smr.estimates,smr.ll,smr.ul)
do.call(text, cat.est.arg)
# Cabecalho dos SMRs
SMR.head.arg$y <- ylim[2] + 1
do.call(text, SMR.head.arg)
}
Try the ems package in your browser
Any scripts or data that you put into this service are public.
ems documentation built on Dec. 11, 2021, 9:15 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions forest.SMR SMR.table SMR
Documented in forest.SMR SMR SMR.table
#' Standardized Mortality Ratio (SMR)
#'
#' @name SMR
#'
#' @description Calculates the standardized mortality ratio and its confidence interval. SMR, for a group, is defined as the ratio of the observed deaths in this group and the sum of the predicted individual probabilities of death by any model (expected deaths).
#'
#' \code{SMR.table} estimates at once the overall SMR and the SMR across several groups, e.g. ICU units or clinical characteristics. The \code{SMR.table} can be ordered by the SMR estimate or its confidence intervals, facilitating the comparison of the units ranks.
#'
#' \code{forest.SMR} shows the \code{SMR.table} output as a forest plot. The plot opens two windows and plot at the left side the values from the \code{SMR.table} and at the right side the points and lines graphically representing each SMR and its confidence interval.
#'
#' @param obs.var Observed death. Accepted values are 0 (absence) or 1 (presence) in a vector. For \code{SMR.table} it must be a character indicating the name of the variable in the data.
#'
#' @param pred.var Death individual predictions (ranging from 0 to 1) in a vector. For \code{SMR.table} it must be a character indicating the name of the variable in the data.
#'
#' @param digits Number of digits for rounding the output.
#'
#' @param ci.method Method to estimate the confidence interval. "Hosmer" (default) or "Byar" are acceptable values.
#'
#' @param ci.level Level of the confidence interval. Default is 0.95.
#'
#' @param data For \code{SMR.table}, a dataset where pred.var, obs.var and group.var are in.
#'
#' @param group.var For \code{SMR.table}, this is a character vector indicating the name(s) of the variable(s) in the data that will form the groups where SMR will be calculated. The variables must be factors.
#'
#' @param use.label Logical. Default is FALSE. For \code{SMR.table} this option will replace the variables names by its labels in var.labels argument.
#'
#' @param var.labels A character vector with variables labels. The default is to replace the variable name by the label stored at attr(data, "var.labels"). But one may specify labels directly.
#'
#' @param reorder Default is "no". Possible values are: "no", "SMR","lower.Cl", and "upper.Cl". It will make the \code{SMR.table} to be ordered within each varibale by its original order, or by SMR order, or by lower.Cl order, or by upper.Cl.
#'
#' @param decreasing Logical. When 'reorderd' is TRUE, should the order be decreasing or incresing? See \code{\link[base]{order}}
#'
#' @param x For the \code{forest.SMR} this is the output of \code{SMR.table}.
#'
#' @param mar1,mar.SMR Values to set the margins (mar parameter) of left and right windows. See \code{\link[graphics]{par}}
#'
#' @param overall.arg A list of arguments passed to \code{\link[graphics]{text}} for plotting the overall label. Internally, 'y' coordinate is replaced.
#'
#' @param NOE.overall.args A list of arguments passed to \code{\link[graphics]{text}} for plotting the overall N (number of observations), O (observed deaths) and E (expected deaths). Internally, 'labels' and 'y' arguments are replaced.
#'
#' @param var.labels.arg A list of arguments passed to \code{\link[graphics]{text}} for plotting the variables labels. Internally, 'y' coordinate is replaced.
#'
#' @param cat.labels.arg A list of arguments passed to \code{\link[graphics]{text}} for plotting the categories labels. Internally, 'y' coordinate is replaced.
#'
#' @param N.values.arg A list of arguments passed to \code{\link[graphics]{text}} for plotting the values of N (number of observations) of each subgroup. Internally, the arguments 'label' and 'y' coordinate are replaced.
#'
#' @param O.values.arg A list of arguments passed to \code{\link[graphics]{text}} for plotting the values of Observed deaths of each subgroup. Internally, the arguments 'label' and 'y' coordinate are replaced.
#'
#' @param E.values.arg A list of arguments passed to \code{\link[graphics]{text}} for plotting the values of Expected deaths of each subgroup. Internally, the arguments 'label' and 'y' coordinate are replaced.
#'
#' @param NOE.head.arg A list of arguments passed to \code{\link[graphics]{text}} for plotting the labels of the columns N, E and O on the top of the graph. Internally, the 'x' and 'y' coordinates are replaced. The x coordinates are taken from the x in \code{N.values.arg}, \code{O.values.arg} and \code{E.values.arg}.
#'
#' @param Overall.seg.arg A list of arguments passed to \code{\link[graphics]{segments}} for plotting the lines corresponding to overall SMR confidence intervals. Internally, 'x' and 'y' coordinates are replaced.
#'
#' @param Overall.p.arg A list of arguments passed to \code{\link[graphics]{points}} for plotting the points corresponding to overall SMR. Internally, 'x' and 'y' coordinates are replaced.
#'
#' @param Overall.est.arg A list of arguments passed to \code{\link[graphics]{text}} for plotting the overall SMR beside the graph. Internally, 'y' coordinate and 'label' argument are replaced.
#'
#' @param cat.seg.arg A list of arguments passed to \code{\link[graphics]{segments}} for plotting the lines corresponding to SMR confidence intervals for all groups. Internally, 'x' and 'y' coordinates are replaced.
#'
#' @param cat.p.arg A list of arguments passed to \code{\link[graphics]{points}} for plotting the points corresponding to all categoreis SMR. Internally, 'x' and 'y' coordinates are replaced.
#'
#' @param cat.est.arg A list of arguments passed to \code{\link[graphics]{text}} for plotting the categories SMR beside the graph. Internally, 'y' coordinate and 'label' arguments are replaced.
#'
#' @param SMR.head.arg A list of arguments passed to \code{\link[graphics]{text}} for plotting the label of the SMR column on the top of the graph. Internally, the 'y' coordinate is replaced.
#'
#' @param smr.xlab Label of the x axis. Default is "Standardized Mortality Ratio".
#'
#' @param smr.xlim Limits of x axis of the \code{forest.SMR} plot. Default is "auto", which internally will pick the highest values of all upper.Cl and the lowest lower.Cl. Besides "auto", only a vector of 2 numbers is valid, and will be passed to \code{\link[graphics]{plot.default}}.
#'
#' @param grid Logical. If TRUE (default), it will draw a grid with the \code{\link[graphics]{grid}} default arguments.
#'
#' @return
#' If SMR, then:
#' \itemize{
#' \item \code{N} Number of subjects analyzed.
#' \item\code{O} Observed number of deaths.
#' \item \code{E} Expected number of deaths.
#' \item \code{SMR} Standardized mortality ratio.
#' \item \code{lower.Cl} lower confidence limit.
#' \item \code{upper.Cl} upper confidence limit.
#' }
#'
#' If SMR.table, then a data.frame with the same information as above, and the
#' additional information is returned: "Variables" (variables names), "Levels" (variables levels).
#'
#' If forest.SMR, then a plot is returned.
#'
#' @author Lunna Borges and Pedro Brasil
#'
#' @seealso \code{\link{SRU}}, \code{\link{reclass}}, \code{\link{funnel}}
#'
#' @examples
#' # Loading a example data
#' data(icu)
#'
#' # Setting variable labels to data
#' attr(icu, "var.labels")[match(c("Unit", "IsMechanicalVentilation1h",
#' "AdmissionTypeName_pri","Vasopressors_D1"), names(icu))] <-
#' c("ICU unit","Mechanichal ventilation","Admission type","Vasopressors at admission")
#'
#' # Some editing
#' icu$Saps3DeathProbabilityStandardEquation <- icu$Saps3DeathProbabilityStandardEquation /100
#' icu$IsMechanicalVentilation1h <- as.factor(ifelse(icu$IsMechanicalVentilation1h == 1, "Yes", "No"))
#' icu$AdmissionTypeName_pri <- as.factor(icu$AdmissionTypeName_pri)
#' levels(icu$AdmissionTypeName_pri) <- c("Clinical","Elective surgery", "Urgent surgery")
#' icu$Vasopressors_D1 <- as.factor(ifelse(icu$Vasopressors_D1 == 1, "Yes", "No"))
#'
#' # The overall SMR for the whole sample
#' SMR(icu$UnitDischargeName, icu$Saps3DeathProbabilityStandardEquation)
#'
#' # The overall SMR and for some subgroups
#' x <- SMR.table(data = icu, obs.var = "UnitDischargeName",
#' pred.var = "Saps3DeathProbabilityStandardEquation",
#' group.var = c( "IsMechanicalVentilation1h",
#' "AdmissionTypeName_pri","Vasopressors_D1"),
#' reorder = "no",
#' decreasing = TRUE,
#' use.label = TRUE)
#' x
#'
#' # A forest plot for all groups SMR (resize the window may be required)
#' forest.SMR(x, digits = 2)
#'
#' # The same thing but reordering the categories
#' x <- SMR.table(data = icu, obs.var = "UnitDischargeName",
#' pred.var = "Saps3DeathProbabilityStandardEquation",
#' group.var = c( "IsMechanicalVentilation1h",
#' "AdmissionTypeName_pri", "Vasopressors_D1"),
#' reorder = "SMR",
#' decreasing = TRUE,
#' use.label = TRUE)
#' forest.SMR(x, digits = 2)
#'
#' # The overall SMR and for all Units
#' x <- SMR.table(data = icu, obs.var = "UnitDischargeName",
#' pred.var = "Saps3DeathProbabilityStandardEquation",
#' group.var = "Unit",
#' reorder = "no",
#' decreasing = TRUE,
#' use.label = TRUE)
#' x
#'
#' # A forest plot for all Units
#' forest.SMR(x, digits = 2)
#'
#' # The same thing but reordering the categories
#' x <- SMR.table(data = icu, obs.var = "UnitDischargeName",
#' pred.var = "Saps3DeathProbabilityStandardEquation",
#' group.var = "Unit",
#' reorder = "SMR",
#' decreasing = TRUE,
#' use.label = TRUE)
#' forest.SMR(x, digits = 2)
#'
#' rm(x, icu)
#' @references
#'
#' David W. Hosmer and Stanley Lemeshow. Confidence intervals estimates of an index of quality performance basend on logistic regression models. Statistics in Medicine , vol. 14, 2161-2172 (1995)
#'
#' @import stats
#' @export
SMR <- function(obs.var, pred.var, digits = 5, ci.method = c("Hosmer", "Byar"), ci.level = 0.95) {
if (length(obs.var) != length(pred.var)){
stop("Length of pred.var and obs.var differ.")
}
if (any(is.na(obs.var) | is.na(pred.var))) {
stop("There are NAs either at pred.var or obs.var. Either remove the NAs or impute!")
}
if (!is.numeric(pred.var)) {
stop("Predicted death variable must be numeric.")
}
if (any(min(pred.var) < 0 | max(pred.var) > 1)) {
stop("The individual predicted death must range from 0 to 1.")
}
if (!is.numeric(obs.var)) {
stop("Observed death variable must be numeric.")
}
if (any(obs.var != 0 & obs.var != 1)) {
stop("Observed death variable must be coded as 0 and 1.")
}
if (ci.method[1] != "Byar" && ci.method[1] != "Hosmer") {
stop("ci.method must be either 'Byar' or 'Hosmer'.")
}
if (!is.numeric(ci.level) || ci.level > 1 || ci.level < 0) {
stop("ci.level must be numeric between 0 and 1.")
}
N <- length(pred.var)
O <- length(obs.var[which(obs.var == 1)])
E <- sum(pred.var)
SMR <- O / E
if (ci.method[1] == "Byar") {
lowerCL <- O / E * (1 - 1 / (9 * O) - qnorm(1 - (1 - ci.level) / 2) / (3 * sqrt(O))) ^ 3
upperCL <- (O + 1) / E * (1 - (1 / (9 * (O + 1))) + qnorm(1 - (1 - ci.level) / 2) /
(3 * sqrt(O + 1))) ^ 3
}
if (ci.method[1] == "Hosmer") {
V2 <- sum(pred.var * (1 - pred.var))
lowerCL <- exp(log(SMR) - qnorm(1 - (1 - ci.level) / 2) * sqrt(V2 / O ^ 2))
upperCL <- exp(log(SMR) + qnorm(1 - (1 - ci.level) / 2) * sqrt(V2 / O ^ 2))
}
output <- c(N = N, Observed = O, Expected = E, SMR = SMR, lower.Cl = lowerCL,
upper.Cl = upperCL)
output <- round(output, digits = digits)
output
}
#' @rdname SMR
#' @export
SMR.table <- function(data, group.var, obs.var, pred.var, digits = 5, use.label = FALSE, var.labels = attr(data, "var.labels")[match(group.var, names(data))], ci.method = c("Hosmer", "Byar"), ci.level = 0.95, reorder = c("no","SMR","lower.Cl","upper.Cl"), decreasing = FALSE) {
if (any(is.na(match(group.var, names(data))))) {
stop("One or more variables in group var is not a variable of data.")
}
if (!all(sapply(data[ , group.var], is.factor))) {
stop("All variables in group.var must be factors.")
}
if (any(is.na(match(c(pred.var, obs.var), names(data))))) {
stop("Either 'pred.var' or 'obs.var' is not a variable in data.")
}
if (use.label) {
if (is.null(var.labels)) {
stop("Either there is no label in attr(data,'var.labels') or 'var.labels argument was not set.")
}
if (length(group.var) != length(var.labels)){
stop("The number of variables and labels are different.")
}
if (!is.character(var.labels)) {
stop("'var.labels ' must be a character vector.")
}
}
if (reorder[1] != "no" & reorder[1] != "SMR" & reorder[1] != "lower.Cl" & reorder[1] != "upper.Cl") {
stop("'reorder' must be one of 'no','SMR','lower.Cl' or 'upper.Cl'")
}
cont.table <- sapply(1:length(group.var), function(i) table(data[, group.var[i]]))
if (any(unlist(cont.table) == 0)) {
# warning("In SMR analysis, levels were deleted before analysis due to zero observations.")
# which levels were deleted
if(is.matrix(unlist(cont.table))){
deleted_levels <- rownames(unlist(cont.table))[which(unlist(cont.table) == 0)]
} else {
deleted_levels <- names(unlist(cont.table))[which(unlist(cont.table) == 0)]
}
warning(paste("In SMR analysis, level",deleted_levels,"were deleted before analysis due to zero observations.\n"))
data[ , group.var] <- droplevels(data[ , group.var])
}
Variables <- c("Overall", unlist(sapply(1:length(group.var), function(i) rep(group.var[i], nlevels(data[ , group.var[i]])))))
if (length(group.var) == 1) {
Levels <- c(NA, levels(data[ , group.var]))
} else {
Levels <- c(NA, unname(unlist(sapply(data[ , group.var], levels))))
}
x <- data.frame(Variables, Levels, stringsAsFactors = TRUE)
x$N <- NA
x$Observed <- NA
x$Expected <- NA
x$SMR <- NA
x$lower.Cl <- NA
x$upper.Cl <- NA
x[1,3:8] <- SMR(data[ , obs.var], data[, pred.var], digits = digits, ci.method = ci.method, ci.level = ci.level)
for(i in 2:nrow(x)){
cond <- which(data[ , as.character(x[i , "Variables"])] == as.character(x[i , 2]))
x[i,3:8] <- SMR(data[cond , obs.var], data[cond , pred.var], digits = digits, ci.method = ci.method, ci.level = ci.level)
}
x$Variables <- as.character(x$Variables)
if (reorder[1] == "SMR") {
for (i in group.var) {
# i = 1
cond <- grep(i, x$Variables)
y <- x[ cond, ]
y <- y[ order(y$SMR, decreasing = decreasing), ]
x[ cond, ] <- y
}
}
if (reorder[1] == "lower.Cl") {
# i = 2
for (i in group.var) {
cond <- grep(i, x$Variables)
y <- x[ cond, ]
y <- y[ order(y$lower.Cl, decreasing = decreasing), ]
x[ cond, ] <- y
}
}
if (reorder[1] == "upper.Cl") {
for (i in group.var) {
cond <- grep(i, x$Variables)
y <- x[ cond, ]
y <- y[ order(y$upper.Cl, decreasing = decreasing), ]
x[ cond, ] <- y
}
}
if (use.label) {
names(group.var) <- var.labels
for (i in 2:length(x$Variables)) {
x$Variables[i] <- names(group.var)[which(x$Variables[i] == group.var)]
}
}
for (i in length(x$Variables):2) {
if (x$Variables[i - 1] == x$Variables[i]) {x$Variables[i] <- NA}
}
x
}
#' @rdname SMR
#' @import graphics
#' @import grDevices
#' @export
forest.SMR <- function(x,
mar1 = c(5.1, 1, 4.1, 1),
mar.SMR = c(5.1, 7, 4.1, 1),
overall.arg = list(x = .01, font = 2, las = 1, labels = var.labels[1], xpd = NA, adj = 0),
NOE.overall.args = list(x = c(N.values.arg$x, O.values.arg$x, E.values.arg$x), font = 2, las = 1, xpd = NA),
var.labels.arg = list(x = .01, font = 2, las = 1, cex = 1, xpd = NA, adj = 0),
cat.labels.arg = list(x = .1, font = 3, las = 1, cex = .95, col = gray(.4), xpd = NA , adj = 0),
N.values.arg = list(x = .5, col = gray(.4), xpd = NA),
O.values.arg = list(x = .675, col = gray(.4), xpd = NA),
E.values.arg = list(x = .85, col = gray(.4), xpd = NA),
NOE.head.arg= list(font = 2, labels = c("N","O","E"), xpd = NA),
Overall.seg.arg = list(col = "navyblue", xpd = NA, lwd = 2),
Overall.p.arg = list(pch = 23, cex = 2, col = "black", bg = gray(.4), xpd = NA),
Overall.est.arg = list(x = smr.xlim[1] - .06, las = 1, font = 2, xpd = NA, adj = 1),
cat.seg.arg = list(col = "navyblue", xpd = NA, lwd = 2),
cat.p.arg = list(pch = 22 ,cex= 1, col = "black", bg = gray(.4), xpd = NA),
cat.est.arg = list(x = smr.xlim[1] - .06, las = 1, col = gray(.4), xpd = NA, adj = 1),
SMR.head.arg = list(smr.xlim[1] - .06, font = 2, labels = "SMR [95% CIs]", xpd = NA, adj = 1),
smr.xlab = "Standardized Mortality Ratio",
smr.xlim = "auto",
grid = TRUE,
digits = 3){
opar <- par(mfrow = c(1,1), mar = c(5.1, 4.1, 4.1, 2.1))
on.exit(par(opar))
# Separando os valores da medida preincipal
main.pos <- 1
main.smr <- unlist(x[which(x$Variables == "Overall"), c("SMR", "lower.Cl", "upper.Cl")])
main.n <- unlist(x[which(x$Variables == "Overall"), "N"])
main.observed <- round(unlist(x[which(x$Variables == "Overall"), "Observed"]), digits = digits)
main.expected <- round(unlist(x[which(x$Variables == "Overall"), "Expected"]), digits = digits)
# Separando as medidas do SMR para montar o grafico
smr.estimates <- unlist(x$SMR[-1])
smr.ll <- unlist(x$lower.Cl[-1])
smr.ul <- unlist(x$upper.Cl[-1])
# Separando os valores para serem impressos na direita do grafico
N <- unlist(x$N[-1])
observed <- round(unlist(x$Observed[-1]), digits = digits)
expected <- round(unlist(x$Expected[-1]), digits = digits)
# Achandos os limites do grafico no eixo horizontal
if (smr.xlim[1] == "auto"){
smr.xlim <- c(min(unlist(c(main.smr[2], smr.ll)), na.rm = TRUE), max(unlist(c(main.smr[3], smr.ul)), na.rm = TRUE))
}
# Separando os valores das vari?veis e categorias
var.labels <- x$Variables[-which(is.na(x$Variables))]
cat.labels <- x$Levels[-which(is.na(x$Levels))]
# Achando as posi??es no eixo vertical
nLevels <- c(which(!is.na(x$Variables))[-1],(nrow(x) + 1))
for(i in 1:length(nLevels)){
nLevels[i] <- nLevels[i+1] - nLevels[i]
}; rm(i)
nLevels <- nLevels[which(!is.na(nLevels))]
# x$Levels
# i = 3
cat.pos <- seq(1:nLevels[1])
if (length(nLevels) == 1) {
cat.pos <- cat.pos + 2
} else {
for(i in 2:length(nLevels)){
ini <- cat.pos[length(cat.pos)]
cat.pos <- c(cat.pos, (ini + 3):(ini + 2 + nLevels[i]))
}
cat.pos <- cat.pos + 2
}
if (length(nLevels) == 1) {
var.pos <- cat.pos[length(cat.pos)] + 1
} else {
var.pos = nLevels[1] + 1
# i = 2
for(i in 2:length(nLevels)){
var.pos <- c(var.pos, var.pos[i-1] + nLevels[i] + 2)
}
var.pos <- var.pos + 2
}
ylim = c(1, var.pos[length(var.pos)])
# 5.1 4.1 4.1 2.1
par(mfrow = c(1,2))
# Janela com textos com variaveis e categorias
par(mar = mar1)
plot(NA, NA, xlim = c(0,1), ylim = ylim , xlab = "", ylab = "",yaxt = "n", xaxt = "n", xaxs = "i", bty = "n", xpd = NA)
# Efeito principal
overall.arg$y <- main.pos
do.call(text, overall.arg)
# N O E do overall
NOE.overall.args$labels <- c(main.n, main.observed, main.expected)
NOE.overall.args$y <- main.pos
do.call(text, NOE.overall.args)
# Demais vari?veis
var.labels.arg$labels <- var.labels[2:length(var.labels)]
var.labels.arg$y <- var.pos
do.call(text, var.labels.arg)
# Categorias
cat.labels.arg$labels <- cat.labels
cat.labels.arg$y <- cat.pos
do.call(text, cat.labels.arg)
# Colunas com os valores de N O e E
N.values.arg$labels <- N
N.values.arg$y <- cat.pos
do.call(text, N.values.arg)
O.values.arg$labels <- observed
O.values.arg$y <- cat.pos
do.call(text, O.values.arg)
E.values.arg$labels <- expected
E.values.arg$y <- cat.pos
do.call(text, E.values.arg)
# Cabecalho das colunas N O E
NOE.head.arg$y <- ylim[2] + 1
NOE.head.arg$x <- c(N.values.arg$x, O.values.arg$x, E.values.arg$x)
do.call(text, NOE.head.arg)
# SMR
par(mar = mar.SMR)
plot(NA, NA, xlim = smr.xlim, ylim = ylim , xlab = smr.xlab, ylab = "", yaxt = "n", xaxt = "n", xaxs = "i", bty = "n", xpd = NA); axis(1)
if(grid) grid()
# Segmento do IC Efeito principal
Overall.seg.arg$x0 <- main.smr[2]
Overall.seg.arg$y0 <- main.pos
Overall.seg.arg$x1 <- main.smr[3]
Overall.seg.arg$y1 <- main.pos
do.call(segments, Overall.seg.arg)
# Ponto do Efeito principal
Overall.p.arg$x <- main.smr[1]
Overall.p.arg$y <- main.pos
Overall.p.arg$type <- "p"
do.call(points, Overall.p.arg)
# Estimativa do efeito principal
Overall.est.arg$y <- main.pos
Overall.est.arg$labels <- sprintf(paste0("%.",digits,"f [%.",digits,"f ; %.",digits,"f]"),main.smr[1],main.smr[2],main.smr[3])
do.call(text, Overall.est.arg)
# Segmentos das categorias
cat.seg.arg$x0 <- smr.ll
cat.seg.arg$y0 <- cat.pos
cat.seg.arg$x1 <- smr.ul
cat.seg.arg$y1 <- cat.pos
do.call(segments, cat.seg.arg)
# Ponto das estimativas das categorias
cat.p.arg$x <- smr.estimates
cat.p.arg$y <- cat.pos
cat.p.arg$type <- "p"
do.call(points, cat.p.arg)
# Estimativas das categorias
cat.est.arg$y <- cat.pos
cat.est.arg$labels <- sprintf(paste0("%.",digits,"f [%.",digits,"f ; %.",digits,"f]"),smr.estimates,smr.ll,smr.ul)
do.call(text, cat.est.arg)
# Cabecalho dos SMRs
SMR.head.arg$y <- ylim[2] + 1
do.call(text, SMR.head.arg)
}
Try the ems package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.