New in episensr 0.9.3 - Multiple bias analysis

New version of episensr available on CRAN! Now you can more easily run multiple bias analysis. Very often, epidemiologic studies suffer from multiple biases. Up to now, it took a bit of “manipulation” to apply a sequence of bias function from episensr. So here’s the function multiple.bias that allows to pipe corrected 2-by-2 table from one bias function to an other.

Let’s take for example the study from Chien et al, 2006 looking at the association between antidepressant use (AD) and breast cancer (BC):

chien <- matrix(c(118, 832, 103, 884),
                dimnames = list(c("BC+", "BC-"), c("AD+", "AD-")),
                nrow = 2, byrow = TRUE)
library(knitr)
kable(chien)
AD+ AD-
BC+ 118 832
BC- 103 884

Antidepressant use was established based on pharmacy records. However this information was not available for every participant, self-report of medication was also used leading to a possible misclassification. The cross-classification frequencies of interview responses and pharmacy records allowed to determine sensitivity (Se) and specificity (Sp) for cases and controls, used in the misclassification bias sensitivity analysis:

library(episensr)
## Registered S3 methods overwritten by 'ggplot2':
##   method         from 
##   [.quosures     rlang
##   c.quosures     rlang
##   print.quosures rlang
library(magrittr)
chien %>% misclassification(., type = "exposure", bias_parms = c(.56, .58, .99, .97))
## --Observed data-- 
##          Outcome: BC+ 
##        Comparing: AD+ vs. AD- 
## 
##     AD+ AD-
## BC+ 118 832
## BC- 103 884
## 
##                                        2.5%     97.5%
## Observed Relative Risk: 1.1012443 0.9646019 1.2572431
##    Observed Odds Ratio: 1.2172330 0.9192874 1.6117443
## ---
##                                                              2.5%    97.5%
## Misclassification Bias Corrected Relative Risk: 1.272939                  
##    Misclassification Bias Corrected Odds Ratio: 1.676452 1.150577 2.442679

Let’s follow this bias by a selection bias:

chien %>% 
    misclassification(., type = "exposure", bias_parms = c(.56, .58, .99, .97)) %>%
    multiple.bias(., bias_function = "selection", bias_parms = c(.73, .61, .82, .76))
## 
## Multiple bias analysis
## ---
##                                                 
## Selection Bias Corrected Relative Risk: 1.192461
##    Selection Bias Corrected Odds Ratio: 1.512206

The misclassification output is piped (%>%) into the selection bias function thanks to multiple.bias. We can also add an unmeasured confounder the same way:

chien %>%
    misclassification(., type = "exposure", bias_parms = c(.56, .58, .99, .97)) %>%
    multiple.bias(., bias_function = "selection",
                  bias_parms = c(.73, .61, .82, .76)) %>%
    multiple.bias(., bias_function = "confounders",
                  type = "OR", bias_parms = c(.92, .3, .44))
## 
## Multiple bias analysis
## ---
##                                        Adjusted OR
## Standardized Morbidity Ratio: 1.494938    1.011609
##              Mantel-Haenszel: 1.494938    1.011609

Probabilistic methods

The same can be realized within a probabilistic framework with the probsens series of functions. Let’s look first at each single bias analysis, then by combining them.

Misclassification

set.seed(123)
mod1 <- chien %>%
    probsens(., type = "exposure", reps = 100000,
             seca.parms = list("trapezoidal", c(.45, .5, .6, .65)),
             seexp.parms = list("trapezoidal", c(.4, .48, .58, .63)),
             spca.parms = list("trapezoidal", c(.95, .97, .99, 1)),
             spexp.parms = list("trapezoidal", c(.96, .98, .99, 1)),
             corr.se = .8, corr.sp = .8)
library(ggplot2)
ggplot(mod1$sim.df, aes(x = corr.OR)) + 
    geom_histogram(aes(y = ..density..), colour = "black", fill = "white") +
    geom_density(alpha = .2, fill = "#FF6666") +
    labs(title = "1. Misclassification bias",
         x = "OR without random error incorporated") +
    coord_cartesian(xlim = c(0, 3))
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Misclassification bias.

Figure 1: Misclassification bias.

Selection bias

set.seed(123)
mod2 <- chien %>%
    probsens.sel(., reps = 100000,
                 case.exp = list("logit-normal", c(-1.1, 0, 0, 1)),
                 case.nexp = list("trapezoidal", c(.75, .85, .95, 1)),
                 ncase.exp = list("logit-normal", c(-1.2, 0, 0, 1)),
                 ncase.nexp = list("trapezoidal", c(0.7, 0.8, 0.9, 1)))
ggplot(mod2$sim.df, aes(x = corr.or)) + 
    geom_histogram(aes(y = ..density..), colour = "black", fill = "white") +
    geom_density(alpha = .2, fill = "#FF6666") +
    labs(title = "2. Selection bias",
         x = "OR without random error incorporated") +
    coord_cartesian(xlim = c(0, 3))
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Selection bias

Figure 2: Selection bias

Confounding

set.seed(123)
mod3 <- chien %>%
    probsens.conf(., reps = 100000,
                  prev.exp = list("logit-normal", c(-0.75, 0.8, 0, 1)),
                  prev.nexp = list("logit-normal", c(-0.4, 0.8, 0, 1)),
                  risk = list("trapezoidal", c(.2, .58, 1.01, 1.24)))
ggplot(mod3$sim.df, aes(x = OR.SMR.or)) + 
    geom_histogram(aes(y = ..density..), colour = "black", fill = "white") +
    geom_density(alpha = .2, fill = "#FF6666") +
    labs(title = "3. Confounding bias",
         x = "OR without random error incorporated") +
    coord_cartesian(xlim = c(0, 3))
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Confounding.

Figure 3: Confounding.

Model 2 then model 3

set.seed(123)
chien %>%
    probsens(., type = "exposure", reps = 100000,
             seca.parms = list("trapezoidal", c(.45, .5, .6, .65)),
             seexp.parms = list("trapezoidal", c(.4, .48, .58, .63)),
             spca.parms = list("trapezoidal", c(.95, .97, .99, 1)),
             spexp.parms = list("trapezoidal", c(.96, .98, .99, 1)),
             corr.se = .8, corr.sp = .8) %>%
    multiple.bias(., bias_function = "probsens.sel",
                  case.exp = list("logit-normal", c(-1.1, 0, 0, 1)),
                  case.nexp = list("trapezoidal", c(.75, .85, .95, 1)),
                  ncase.exp = list("logit-normal", c(-1.2, 0, 0, 1)),
                  ncase.nexp = list("trapezoidal", c(0.7, 0.8, 0.9, 1)))
## 
## Multiple bias analysis
## ---
##                                               Median 2.5th percentile
##            Odds Ratio -- systematic error: 1.1074217        0.9202678
## Odds Ratio -- systematic and random error: 1.1086260        0.8189285
##                                            97.5th percentile
##            Odds Ratio -- systematic error:         1.3517598
## Odds Ratio -- systematic and random error:         1.5020760

Model 2 then model 3 then model 4

set.seed(123)
mod6 <- chien %>%
    probsens(., type = "exposure", reps = 100000,
             seca.parms = list("trapezoidal", c(.45, .5, .6, .65)),
             seexp.parms = list("trapezoidal", c(.4, .48, .58, .63)),
             spca.parms = list("trapezoidal", c(.95, .97, .99, 1)),
             spexp.parms = list("trapezoidal", c(.96, .98, .99, 1)),
             corr.se = .8, corr.sp = .8) %>%
    multiple.bias(., bias_function = "probsens.sel",
                  case.exp = list("logit-normal", c(-1.1, 0, 0, 1)),
                  case.nexp = list("trapezoidal", c(.75, .85, .95, 1)),
                  ncase.exp = list("logit-normal", c(-1.2, 0, 0, 1)),
                  ncase.nexp = list("trapezoidal", c(0.7, 0.8, 0.9, 1))) %>% 
    multiple.bias(., bias_function = "probsens.conf",
                  prev.exp = list("logit-normal", c(-0.75, 0.8, 0, 1)),
                  prev.nexp = list("logit-normal", c(-0.4, 0.8, 0, 1)),
                  risk = list("trapezoidal", c(.2, .58, 1.01, 1.24)))
ggplot(mod6$sim.df, aes(x = OR.SMR.or)) + 
    geom_histogram(aes(y = ..density..), colour = "black", fill = "white") +
    geom_density(alpha = .2, fill = "#FF6666") +
    labs(title = "6. Misclassification, then\nselection, then confounding bias",
         x = "OR without random error incorporated") +
    coord_cartesian(xlim = c(0, 3))
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Misclassification followed by selection bias and then confounding bias.

Figure 4: Misclassification followed by selection bias and then confounding bias.

ggplot(mod6$sim.df, aes(x = tot.ORadj.smr)) + 
    geom_histogram(aes(y = ..density..), colour = "black", fill = "white") +
    geom_density(alpha = .2, fill = "#FF6666") +
    labs(title = "6. Misclassification, then\nselection, then confounding bias",
         x = "OR with random error incorporated") +
    coord_cartesian(xlim = c(0, 3))
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Same as previous figure but with random error incorporated.

Figure 5: Same as previous figure but with random error incorporated.