Nothing
R/misclassification.cov.R
In episensr: Basic Sensitivity Analysis of Epidemiological Results
Defines functions
misclassification.cov
Documented in
misclassification.cov
#' Sensitivity analysis for covariate misclassification.
#'
#' Simple sensitivity analysis to correct for a misclassified covariate (a potential
#' confounder or effect measure modifier).
#'
#' @param case Outcome variable. If a variable, this variable is tabulated against.
#' @param exposed Exposure variable.
#' @param covariate Covariate to stratify on.
#' @param bias_parms Vector defining the bias parameters. This vector has 4 elements
#' between 0 and 1, in the following order:
#' \enumerate{
#' \item Sensitivity of confounder classification among those with the outcome,
#' \item Sensitivity of confounder classification among those without the outcome,
#' \item Specificity of confounder classification among those with the outcome,and
#' \item Specificity of confounder classification among those without the outcome.
#' }
#' @param alpha Significance level.
#'
#' @return A list with elements:
#' \item{obs.data}{The analyzed stratified 2 x 2 tables from the observed data.}
#' \item{corr.data}{The expected stratified observed data given the true data assuming
#' misclassification.}
#' \item{obs.measures}{A table of observed relative risk and odds ratio with
#' confidence intervals.}
#' \item{adj.measures}{A table of adjusted relative risk and odds ratio.}
#' \item{bias.parms}{Input bias parameters.}
#'
#' @references Lash, T.L., Fox, M.P, Fink, A.K., 2009 \emph{Applying Quantitative
#' Bias Analysis to Epidemiologic Data}, pp.79--108, Springer.
#'
#' @examples
#' # The data for this example come from:
#' # Berry, R.J., Kihlberg, R., and Devine, O. Impact of misclassification of in vitro
#' # fertilisation in studies of folic acid and twinning: modelling using population
#' # based Swedish vital records.
#' # BMJ, doi:10.1136/bmj.38369.437789.82 (published 17 March 2004)
#' misclassification.cov(array(c(1319, 38054, 5641, 405546,
#' 565, 3583, 781, 21958,
#' 754, 34471, 4860, 383588),
#' dimnames = list(c("Twins+", "Twins-"),
#' c("Folic acid+", "Folic acid-"), c("Total", "IVF+", "IVF-")),
#' dim = c(2, 2, 3)),
#' bias_parms = c(.6, .6, .95, .95))
#' @export
#' @importFrom stats qnorm
misclassification.cov <- function(case,
exposed,
covariate,
bias_parms = NULL,
alpha = 0.05){
if(is.null(bias_parms))
bias_parms <- c(1, 1, 1, 1)
else bias_parms <- bias_parms
if(length(bias_parms) != 4)
stop('The argument bias_parms should be made of the following components: (1) Sensitivity of classification among those with the outcome, (2) Sensitivity of classification among those without the outcome, (3) Specificity of classification among those with the outcome, and (4) Specificity of classification among those without the outcome.')
if(!all(bias_parms >= 0 & bias_parms <=1))
stop('Bias parameters should be between 0 and 1.')
if(inherits(case, c("table", "array")))
tab <- case
else {
tab.df <- table(case, exposed, covariate)
tab <- tab.df[2:1, 2:1, ]
}
a <- as.numeric(tab[1, 1, 1])
b <- as.numeric(tab[1, 2, 1])
c <- as.numeric(tab[2, 1, 1])
d <- as.numeric(tab[2, 2, 1])
Ac1 <- as.numeric(tab[1, 1, 2])
Bc1 <- as.numeric(tab[1, 2, 2])
Cc1 <- as.numeric(tab[2, 1, 2])
Dc1 <- as.numeric(tab[2, 2, 2])
Ac0 <- as.numeric(tab[1, 1, 3])
Bc0 <- as.numeric(tab[1, 2, 3])
Cc0 <- as.numeric(tab[2, 1, 3])
Dc0 <- as.numeric(tab[2, 2, 3])
obs.rr <- (a/(a + c)) / (b/(b + d))
se.log.obs.rr <- sqrt((c/a) / (a+c) + (d/b) / (b+d))
lci.obs.rr <- exp(log(obs.rr) - qnorm(1 - alpha/2) * se.log.obs.rr)
uci.obs.rr <- exp(log(obs.rr) + qnorm(1 - alpha/2) * se.log.obs.rr)
obs.or <- (a/c) / (b/d)
se.log.obs.or <- sqrt(1/a + 1/b + 1/c + 1/d)
lci.obs.or <- exp(log(obs.or) - qnorm(1 - alpha/2) * se.log.obs.or)
uci.obs.or <- exp(log(obs.or) + qnorm(1 - alpha/2) * se.log.obs.or)
SMR_RR <- a / (((Ac1+Cc1) * Bc1/(Bc1+Dc1)) + ((Ac0+Cc0) * Bc0/(Bc0+Dc0)))
SMR_OR <- a / ((Cc1 * Bc1/Dc1) + (Cc0 * Bc0/Dc0))
SMR_RR_C <- obs.rr / SMR_RR
SMR_OR_C <- obs.or / SMR_OR
MH_RR <- (Ac1 * (Bc1+Dc1) / sum(tab[,,2]) + Ac0 * (Bc0+Dc0) / sum(tab[,,3])) /
(Bc1 * (Ac1+Cc1) / sum(tab[,,2]) + Bc0 * (Ac0+Cc0) / sum(tab[,,3]))
MH_OR <- (Ac1 * Dc1 / sum(tab[,,2]) + Ac0 * Dc0 / sum(tab[,,3])) /
(Bc1 * Cc1 / sum(tab[,,2]) + Bc0 * Cc0 / sum(tab[,,3]))
MH_RR_C <- obs.rr / MH_RR
MH_OR_C <- obs.or / MH_OR
Ac1_Cr <- (Ac1 - (1 - bias_parms[3]) * a) / (bias_parms[1] - (1 - bias_parms[3]))
Bc1_Cr <- (Bc1 - (1 - bias_parms[3]) * b) / (bias_parms[1] - (1 - bias_parms[3]))
Cc1_Cr <- (Cc1 - (1 - bias_parms[4]) * c) / (bias_parms[2] - (1 - bias_parms[4]))
Dc1_Cr <- (Dc1 - (1 - bias_parms[4]) * d) / (bias_parms[2] - (1 - bias_parms[4]))
Ac0_Cr <- a - Ac1_Cr
Bc0_Cr <- b - Bc1_Cr
Cc0_Cr <- c - Cc1_Cr
Dc0_Cr <- d - Dc1_Cr
if(Ac1_Cr < 1 | Ac0_Cr < 1 | Bc1_Cr < 1 | Bc0_Cr < 1 | Cc1_Cr < 1 | Cc0_Cr < 1 |
Dc1_Cr < 1 | Dc0_Cr < 1)
stop('Parameters chosen lead to negative cell(s) in adjusted stratified 2x2 table.')
corr.tab1 <- matrix(c(Ac1_Cr, Bc1_Cr, Cc1_Cr, Dc1_Cr), nrow = 2, byrow = TRUE)
corr.tab0 <- matrix(c(Ac0_Cr, Bc0_Cr, Cc0_Cr, Dc0_Cr), nrow = 2, byrow = TRUE)
corr.SMR_RR <- a / (((Ac1_Cr+Cc1_Cr) * Bc1_Cr / (Bc1_Cr+Dc1_Cr)) +
((Ac0_Cr+Cc0_Cr) * Bc0_Cr / (Bc0_Cr+Dc0_Cr)))
corr.SMR_OR <- a / ((Cc1_Cr * Bc1_Cr / Dc1_Cr) + (Cc0_Cr * Bc0_Cr / Dc0_Cr))
corr.RR_C <- obs.rr / corr.SMR_RR
corr.OR_C <- obs.or / corr.SMR_OR
corr.MH_RR <- (Ac1_Cr * (Bc1_Cr+Dc1_Cr) / sum(corr.tab1) + Ac0_Cr *
(Bc0_Cr+Dc0_Cr) / sum(corr.tab0)) /
(Bc1_Cr * (Ac1_Cr+Cc1_Cr) / sum(corr.tab1) + Bc0_Cr * (Ac0_Cr+Cc0_Cr) /
sum(corr.tab0))
corr.MH_OR <- (Ac1_Cr * Dc1_Cr / sum(corr.tab1) + Ac0_Cr * Dc0_Cr / sum(corr.tab0)) /
(Bc1_Cr * Cc1_Cr / sum(corr.tab1) + Bc0_Cr * Cc0_Cr / sum(corr.tab0))
corr.MH_RR_C <- obs.rr / corr.MH_RR
corr.MH_OR_C <- obs.or / corr.MH_OR
if (is.null(rownames(tab)))
rownames(tab) <- paste("Row", 1:2)
if (is.null(colnames(tab)))
colnames(tab) <- paste("Col", 1:2)
if (is.null(rownames(tab))){
rownames(corr.tab1) <- paste("Row", 1:2)
rownames(corr.tab0) <- paste("Row", 1:2)
} else {
rownames(corr.tab1) <- row.names(tab)
rownames(corr.tab0) <- row.names(tab)
}
if (is.null(colnames(tab))){
colnames(corr.tab1) <- paste("Col", 1:2)
colnames(corr.tab0) <- paste("Col", 1:2)
} else {
colnames(corr.tab1) <- colnames(tab)
colnames(corr.tab0) <- colnames(tab)
}
rmat <- rbind(c(obs.rr, lci.obs.rr, uci.obs.rr),
c(obs.or, lci.obs.or, uci.obs.or))
rownames(rmat) <- c("Observed Relative Risk:",
" Observed Odds Ratio:")
colnames(rmat) <- c(" ",
paste(100 * (alpha/2), "%", sep = ""),
paste(100 * (1 - alpha/2), "%", sep = ""))
rmata <- rbind(SMR_RR, SMR_RR_C, MH_RR, MH_RR_C,
SMR_OR, SMR_OR_C, MH_OR, MH_OR_C)
rmatar <- rbind(corr.SMR_RR, corr.RR_C, corr.MH_RR, corr.MH_RR_C,
corr.SMR_OR, corr.OR_C, corr.MH_OR, corr.MH_OR_C)
rmatc <- cbind(rmata, rmatar)
rownames(rmatc) <- c(" SMR RR adjusted for confounder:",
" RR due to confounding by misclassified confounder:",
" Mantel-Haenszel RR adjusted for confounder:",
"MH RR due to confounding by misclassified confounder:",
" SMR OR adjusted for confounder:",
" OR due to confounding by misclassified confounder:",
" Mantel-Haenszel OR adjusted for confounder:",
"MH OR due to confounding by misclassified confounder:")
colnames(rmatc) <- c("Observed", "Corrected")
res <- list(model = "misclassification.cov",
obs.data = tab,
corr.data = c(corr.tab1, corr.tab0),
obs.measures = rmat,
adj.measures = rmatc,
bias.parms = bias_parms)
class(res) <- c("episensr", "list")
res
}
Try the episensr package in your browser
Any scripts or data that you put into this service are public.
episensr documentation built on Aug. 30, 2023, 5:09 p.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 misclassification.cov
Documented in misclassification.cov
#' Sensitivity analysis for covariate misclassification.
#'
#' Simple sensitivity analysis to correct for a misclassified covariate (a potential
#' confounder or effect measure modifier).
#'
#' @param case Outcome variable. If a variable, this variable is tabulated against.
#' @param exposed Exposure variable.
#' @param covariate Covariate to stratify on.
#' @param bias_parms Vector defining the bias parameters. This vector has 4 elements
#' between 0 and 1, in the following order:
#' \enumerate{
#' \item Sensitivity of confounder classification among those with the outcome,
#' \item Sensitivity of confounder classification among those without the outcome,
#' \item Specificity of confounder classification among those with the outcome,and
#' \item Specificity of confounder classification among those without the outcome.
#' }
#' @param alpha Significance level.
#'
#' @return A list with elements:
#' \item{obs.data}{The analyzed stratified 2 x 2 tables from the observed data.}
#' \item{corr.data}{The expected stratified observed data given the true data assuming
#' misclassification.}
#' \item{obs.measures}{A table of observed relative risk and odds ratio with
#' confidence intervals.}
#' \item{adj.measures}{A table of adjusted relative risk and odds ratio.}
#' \item{bias.parms}{Input bias parameters.}
#'
#' @references Lash, T.L., Fox, M.P, Fink, A.K., 2009 \emph{Applying Quantitative
#' Bias Analysis to Epidemiologic Data}, pp.79--108, Springer.
#'
#' @examples
#' # The data for this example come from:
#' # Berry, R.J., Kihlberg, R., and Devine, O. Impact of misclassification of in vitro
#' # fertilisation in studies of folic acid and twinning: modelling using population
#' # based Swedish vital records.
#' # BMJ, doi:10.1136/bmj.38369.437789.82 (published 17 March 2004)
#' misclassification.cov(array(c(1319, 38054, 5641, 405546,
#' 565, 3583, 781, 21958,
#' 754, 34471, 4860, 383588),
#' dimnames = list(c("Twins+", "Twins-"),
#' c("Folic acid+", "Folic acid-"), c("Total", "IVF+", "IVF-")),
#' dim = c(2, 2, 3)),
#' bias_parms = c(.6, .6, .95, .95))
#' @export
#' @importFrom stats qnorm
misclassification.cov <- function(case,
exposed,
covariate,
bias_parms = NULL,
alpha = 0.05){
if(is.null(bias_parms))
bias_parms <- c(1, 1, 1, 1)
else bias_parms <- bias_parms
if(length(bias_parms) != 4)
stop('The argument bias_parms should be made of the following components: (1) Sensitivity of classification among those with the outcome, (2) Sensitivity of classification among those without the outcome, (3) Specificity of classification among those with the outcome, and (4) Specificity of classification among those without the outcome.')
if(!all(bias_parms >= 0 & bias_parms <=1))
stop('Bias parameters should be between 0 and 1.')
if(inherits(case, c("table", "array")))
tab <- case
else {
tab.df <- table(case, exposed, covariate)
tab <- tab.df[2:1, 2:1, ]
}
a <- as.numeric(tab[1, 1, 1])
b <- as.numeric(tab[1, 2, 1])
c <- as.numeric(tab[2, 1, 1])
d <- as.numeric(tab[2, 2, 1])
Ac1 <- as.numeric(tab[1, 1, 2])
Bc1 <- as.numeric(tab[1, 2, 2])
Cc1 <- as.numeric(tab[2, 1, 2])
Dc1 <- as.numeric(tab[2, 2, 2])
Ac0 <- as.numeric(tab[1, 1, 3])
Bc0 <- as.numeric(tab[1, 2, 3])
Cc0 <- as.numeric(tab[2, 1, 3])
Dc0 <- as.numeric(tab[2, 2, 3])
obs.rr <- (a/(a + c)) / (b/(b + d))
se.log.obs.rr <- sqrt((c/a) / (a+c) + (d/b) / (b+d))
lci.obs.rr <- exp(log(obs.rr) - qnorm(1 - alpha/2) * se.log.obs.rr)
uci.obs.rr <- exp(log(obs.rr) + qnorm(1 - alpha/2) * se.log.obs.rr)
obs.or <- (a/c) / (b/d)
se.log.obs.or <- sqrt(1/a + 1/b + 1/c + 1/d)
lci.obs.or <- exp(log(obs.or) - qnorm(1 - alpha/2) * se.log.obs.or)
uci.obs.or <- exp(log(obs.or) + qnorm(1 - alpha/2) * se.log.obs.or)
SMR_RR <- a / (((Ac1+Cc1) * Bc1/(Bc1+Dc1)) + ((Ac0+Cc0) * Bc0/(Bc0+Dc0)))
SMR_OR <- a / ((Cc1 * Bc1/Dc1) + (Cc0 * Bc0/Dc0))
SMR_RR_C <- obs.rr / SMR_RR
SMR_OR_C <- obs.or / SMR_OR
MH_RR <- (Ac1 * (Bc1+Dc1) / sum(tab[,,2]) + Ac0 * (Bc0+Dc0) / sum(tab[,,3])) /
(Bc1 * (Ac1+Cc1) / sum(tab[,,2]) + Bc0 * (Ac0+Cc0) / sum(tab[,,3]))
MH_OR <- (Ac1 * Dc1 / sum(tab[,,2]) + Ac0 * Dc0 / sum(tab[,,3])) /
(Bc1 * Cc1 / sum(tab[,,2]) + Bc0 * Cc0 / sum(tab[,,3]))
MH_RR_C <- obs.rr / MH_RR
MH_OR_C <- obs.or / MH_OR
Ac1_Cr <- (Ac1 - (1 - bias_parms[3]) * a) / (bias_parms[1] - (1 - bias_parms[3]))
Bc1_Cr <- (Bc1 - (1 - bias_parms[3]) * b) / (bias_parms[1] - (1 - bias_parms[3]))
Cc1_Cr <- (Cc1 - (1 - bias_parms[4]) * c) / (bias_parms[2] - (1 - bias_parms[4]))
Dc1_Cr <- (Dc1 - (1 - bias_parms[4]) * d) / (bias_parms[2] - (1 - bias_parms[4]))
Ac0_Cr <- a - Ac1_Cr
Bc0_Cr <- b - Bc1_Cr
Cc0_Cr <- c - Cc1_Cr
Dc0_Cr <- d - Dc1_Cr
if(Ac1_Cr < 1 | Ac0_Cr < 1 | Bc1_Cr < 1 | Bc0_Cr < 1 | Cc1_Cr < 1 | Cc0_Cr < 1 |
Dc1_Cr < 1 | Dc0_Cr < 1)
stop('Parameters chosen lead to negative cell(s) in adjusted stratified 2x2 table.')
corr.tab1 <- matrix(c(Ac1_Cr, Bc1_Cr, Cc1_Cr, Dc1_Cr), nrow = 2, byrow = TRUE)
corr.tab0 <- matrix(c(Ac0_Cr, Bc0_Cr, Cc0_Cr, Dc0_Cr), nrow = 2, byrow = TRUE)
corr.SMR_RR <- a / (((Ac1_Cr+Cc1_Cr) * Bc1_Cr / (Bc1_Cr+Dc1_Cr)) +
((Ac0_Cr+Cc0_Cr) * Bc0_Cr / (Bc0_Cr+Dc0_Cr)))
corr.SMR_OR <- a / ((Cc1_Cr * Bc1_Cr / Dc1_Cr) + (Cc0_Cr * Bc0_Cr / Dc0_Cr))
corr.RR_C <- obs.rr / corr.SMR_RR
corr.OR_C <- obs.or / corr.SMR_OR
corr.MH_RR <- (Ac1_Cr * (Bc1_Cr+Dc1_Cr) / sum(corr.tab1) + Ac0_Cr *
(Bc0_Cr+Dc0_Cr) / sum(corr.tab0)) /
(Bc1_Cr * (Ac1_Cr+Cc1_Cr) / sum(corr.tab1) + Bc0_Cr * (Ac0_Cr+Cc0_Cr) /
sum(corr.tab0))
corr.MH_OR <- (Ac1_Cr * Dc1_Cr / sum(corr.tab1) + Ac0_Cr * Dc0_Cr / sum(corr.tab0)) /
(Bc1_Cr * Cc1_Cr / sum(corr.tab1) + Bc0_Cr * Cc0_Cr / sum(corr.tab0))
corr.MH_RR_C <- obs.rr / corr.MH_RR
corr.MH_OR_C <- obs.or / corr.MH_OR
if (is.null(rownames(tab)))
rownames(tab) <- paste("Row", 1:2)
if (is.null(colnames(tab)))
colnames(tab) <- paste("Col", 1:2)
if (is.null(rownames(tab))){
rownames(corr.tab1) <- paste("Row", 1:2)
rownames(corr.tab0) <- paste("Row", 1:2)
} else {
rownames(corr.tab1) <- row.names(tab)
rownames(corr.tab0) <- row.names(tab)
}
if (is.null(colnames(tab))){
colnames(corr.tab1) <- paste("Col", 1:2)
colnames(corr.tab0) <- paste("Col", 1:2)
} else {
colnames(corr.tab1) <- colnames(tab)
colnames(corr.tab0) <- colnames(tab)
}
rmat <- rbind(c(obs.rr, lci.obs.rr, uci.obs.rr),
c(obs.or, lci.obs.or, uci.obs.or))
rownames(rmat) <- c("Observed Relative Risk:",
" Observed Odds Ratio:")
colnames(rmat) <- c(" ",
paste(100 * (alpha/2), "%", sep = ""),
paste(100 * (1 - alpha/2), "%", sep = ""))
rmata <- rbind(SMR_RR, SMR_RR_C, MH_RR, MH_RR_C,
SMR_OR, SMR_OR_C, MH_OR, MH_OR_C)
rmatar <- rbind(corr.SMR_RR, corr.RR_C, corr.MH_RR, corr.MH_RR_C,
corr.SMR_OR, corr.OR_C, corr.MH_OR, corr.MH_OR_C)
rmatc <- cbind(rmata, rmatar)
rownames(rmatc) <- c(" SMR RR adjusted for confounder:",
" RR due to confounding by misclassified confounder:",
" Mantel-Haenszel RR adjusted for confounder:",
"MH RR due to confounding by misclassified confounder:",
" SMR OR adjusted for confounder:",
" OR due to confounding by misclassified confounder:",
" Mantel-Haenszel OR adjusted for confounder:",
"MH OR due to confounding by misclassified confounder:")
colnames(rmatc) <- c("Observed", "Corrected")
res <- list(model = "misclassification.cov",
obs.data = tab,
corr.data = c(corr.tab1, corr.tab0),
obs.measures = rmat,
adj.measures = rmatc,
bias.parms = bias_parms)
class(res) <- c("episensr", "list")
res
}
Try the episensr 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.