The hardware and bandwidth for this mirror is donated by dogado GmbH, the Webhosting and Full Service-Cloud Provider. Check out our Wordpress Tutorial.
If you wish to report a bug, or if you are interested in having us mirror your free-software or open-source project, please feel free to contact us at mirror[@]dogado.de.

R package cases: overview

The goal of is this vignette is to illustrate the R package cases by some elementary code examples.

Preparation

Load the package:

library(cases)

Important functions

categorize()

Often, binary predictions are not readily available but rather need to be derived from continuous (risk) scores. This can be done via the categorize function.

# real data example from publication here
set.seed(123)
M <- as.data.frame(mvtnorm::rmvnorm(10, mean=rep(0, 3), sigma=2*diag(3)))
M
#>             V1         V2         V3
#> 1  -0.79263226 -0.3255201  2.2043464
#> 2   0.09971392  0.1828405  2.4254682
#> 3   0.65183395 -1.7890668 -0.9713566
#> 4  -0.63026120  1.7311131  0.5088536
#> 5   0.56677642  0.1565290 -0.7860781
#> 6   2.52707679  0.7040669 -2.7812167
#> 7   0.99186703 -0.6686280 -1.5101308
#> 8  -0.30826308 -1.4509894 -1.0308079
#> 9  -0.88393901 -2.3853446  1.1848098
#> 10  0.21690234 -1.6095687  1.7731621

## categorize at 0 by default
yhat <- categorize(M)
yhat
#>    V1_0 V2_0 V3_0
#> 1     0    0    1
#> 2     1    1    1
#> 3     1    0    0
#> 4     0    1    1
#> 5     1    1    0
#> 6     1    1    0
#> 7     1    0    0
#> 8     0    0    0
#> 9     0    0    1
#> 10    1    0    1

## define multiple cutpoints to define multiple decision rules per marker
C <- c(0, 1, 0, 1, 0, 1)
a <- c(1, 1, 2, 2, 3, 3)
categorize(M, C, a)
#>    V1_0 V1_1 V2_0 V2_1 V3_0 V3_1
#> 1     0    0    0    0    1    1
#> 2     1    0    1    0    1    1
#> 3     1    0    0    0    0    0
#> 4     0    0    1    1    1    0
#> 5     1    0    1    0    0    0
#> 6     1    1    1    0    0    0
#> 7     1    0    0    0    0    0
#> 8     0    0    0    0    0    0
#> 9     0    0    0    0    1    1
#> 10    1    0    0    0    1    1


## this can even be used to do multi-class classification, like this:
C <- matrix(rep(c(-1, 0, 1, -2, 0, 2), 3), ncol=3, byrow = TRUE)
C
#>      [,1] [,2] [,3]
#> [1,]   -1    0    1
#> [2,]   -2    0    2
#> [3,]   -1    0    1
#> [4,]   -2    0    2
#> [5,]   -1    0    1
#> [6,]   -2    0    2
categorize(M, C, a)
#>    V1_a V1_b V2_a V2_b V3_a V3_b
#> 1     1    1    1    1    3    3
#> 2     2    2    2    2    3    3
#> 3     2    2    0    1    1    1
#> 4     1    1    3    2    2    2
#> 5     2    2    2    2    1    1
#> 6     3    3    2    2    0    0
#> 7     2    2    1    1    0    1
#> 8     1    1    0    1    0    1
#> 9     1    1    0    0    3    2
#> 10    2    2    0    1    3    2

compare()

In supervised classification, it is assumed that we have a true set of labels. In medical testing, this is usually called the reference standard provided by an established diagnostic/prognostic tool. We need to compare model predictions against these labels in order to compute model accuracy.

## consider binary prediction from 3 models from previous r chunk
names(yhat) <- paste0("rule", 1:ncol(yhat))
yhat
#>    rule1 rule2 rule3
#> 1      0     0     1
#> 2      1     1     1
#> 3      1     0     0
#> 4      0     1     1
#> 5      1     1     0
#> 6      1     1     0
#> 7      1     0     0
#> 8      0     0     0
#> 9      0     0     1
#> 10     1     0     1

## assume true labels
y <- c(rep(1, 5), rep(0, 5))

## compare then results in 
compare(yhat, y)
#> $specificity
#>    rule1 rule2 rule3
#> 6      0     0     1
#> 7      0     1     1
#> 8      1     1     1
#> 9      1     1     0
#> 10     0     1     0
#> 
#> $sensitivity
#>   rule1 rule2 rule3
#> 1     0     0     1
#> 2     1     1     1
#> 3     1     0     0
#> 4     0     1     1
#> 5     1     1     0

evaluate()

Main function of the package

evaluate(compare(yhat, y))
#> [cases] evaluation results:
#> $specificity
#>   parameter hypothesis estimate  lower  upper   pval pval_all reject reject_all
#> 1     rule1    == 0.75      0.4 0.0080 0.7920 0.0401   0.2266  FALSE      FALSE
#> 2     rule2    == 0.75      0.8 0.4799 1.1201 0.3797   0.3797  FALSE      FALSE
#> 3     rule3    == 0.75      0.6 0.2080 0.9920 0.2266   0.2266  FALSE      FALSE
#> 
#> $sensitivity
#>   parameter hypothesis estimate lower upper   pval pval_all reject reject_all
#> 1     rule1    == 0.75      0.6 0.208 0.992 0.2266   0.2266  FALSE      FALSE
#> 2     rule2    == 0.75      0.6 0.208 0.992 0.2266   0.3797  FALSE      FALSE
#> 3     rule3    == 0.75      0.6 0.208 0.992 0.2266   0.2266  FALSE      FALSE

More details on the dta function are provided in the last section

draw_data()

cases includes a few functions for synthetic data generation

draw_data_lfc(n=20)
#> $specificity
#>       model1 model2 model3 model4 model5 model6 model7 model8 model9 model10
#>  [1,]      1      1      1      1      1      1      1      1      1       1
#>  [2,]      1      1      1      1      1      0      1      1      1       1
#>  [3,]      0      1      1      1      1      1      1      1      1       1
#>  [4,]      1      1      1      1      1      1      1      1      1       1
#>  [5,]      1      1      0      1      1      1      1      1      1       1
#>  [6,]      1      1      1      1      1      1      0      1      0       1
#>  [7,]      1      1      1      1      1      1      1      1      1       1
#>  [8,]      1      1      0      1      1      1      1      1      1       1
#>  [9,]      1      1      1      1      1      0      1      1      1       1
#> [10,]      1      1      1      1      1      1      1      1      1       1
#> 
#> $sensitivity
#>       model1 model2 model3 model4 model5 model6 model7 model8 model9 model10
#>  [1,]      1      1      1      1      1      1      1      1      1       1
#>  [2,]      1      1      1      1      1      1      1      1      1       1
#>  [3,]      1      1      1      1      1      1      1      1      1       0
#>  [4,]      1      1      1      1      1      1      1      1      1       0
#>  [5,]      1      1      1      1      0      1      1      1      1       1
#>  [6,]      1      1      1      1      0      1      1      1      1       1
#>  [7,]      1      1      1      1      1      1      1      0      1       0
#>  [8,]      1      1      1      0      1      1      1      1      1       1
#>  [9,]      1      1      1      1      1      1      1      0      1       1
#> [10,]      1      1      1      1      1      1      1      1      1       0
#> 
#> attr(,"info")
#>      model     b  se  sp
#> 1   model1  TRUE 1.0 0.8
#> 2   model2 FALSE 0.8 1.0
#> 3   model3  TRUE 1.0 0.8
#> 4   model4 FALSE 0.8 1.0
#> 5   model5 FALSE 0.8 1.0
#> 6   model6  TRUE 1.0 0.8
#> 7   model7  TRUE 1.0 0.8
#> 8   model8 FALSE 0.8 1.0
#> 9   model9  TRUE 1.0 0.8
#> 10 model10 FALSE 0.8 1.0
draw_data_roc(n=20)
#> $specificity
#>       model1 model2 model3 model4 model5 model6 model7 model8 model9 model10
#>  [1,]      1      0      1      0      1      1      1      1      1       1
#>  [2,]      1      1      1      1      1      1      1      1      1       1
#>  [3,]      0      1      1      1      1      1      1      1      1       1
#>  [4,]      0      0      1      0      1      1      1      1      1       1
#>  [5,]      1      1      1      0      0      1      1      1      1       1
#>  [6,]      1      1      0      1      1      1      1      1      1       1
#>  [7,]      0      1      1      1      1      1      1      1      1       1
#>  [8,]      1      1      0      0      1      1      1      1      1       1
#>  [9,]      1      1      1      1      1      1      1      1      1       1
#> [10,]      0      1      1      1      1      1      1      1      1       1
#> 
#> $sensitivity
#>       model1 model2 model3 model4 model5 model6 model7 model8 model9 model10
#>  [1,]      1      0      1      1      1      1      1      1      1       0
#>  [2,]      1      0      0      0      0      0      0      0      0       0
#>  [3,]      1      0      1      1      1      1      1      1      0       0
#>  [4,]      1      1      1      1      1      1      1      1      1       1
#>  [5,]      1      1      1      1      1      1      1      1      1       1
#>  [6,]      1      0      1      1      1      1      1      1      1       1
#>  [7,]      1      1      1      1      1      1      1      1      1       1
#>  [8,]      1      1      1      1      1      1      1      1      1       1
#>  [9,]      0      1      1      1      1      1      1      1      1       1
#> [10,]      0      1      1      1      1      1      1      1      1       1
#> 
#> attr(,"info")
#>      model   auc    cutoff        se        sp
#> 1   model1 0.850 0.7026076 0.7773072 0.7588499
#> 2   model2 0.875 0.9744355 0.7429298 0.8350798
#> 3   model3 0.900 0.7414402 0.8579035 0.7707867
#> 4   model4 0.925 0.6637751 0.9149729 0.7465829
#> 5   model5 0.925 0.8579379 0.8805752 0.8045366
#> 6   model6 0.950 0.5861100 0.9590761 0.7210992
#> 7   model7 0.950 0.9744355 0.9117706 0.8350798
#> 8   model8 0.950 1.0909332 0.8916296 0.8623489
#> 9   model9 0.950 1.1685983 0.8764814 0.8787172
#> 10 model10 0.950 1.4792588 0.8014789 0.9304644

Remark: Synthetic data comes at the ‘compared’ level meaning the labels 1 and 0 indicate correct and false predictions, respectively. No need to compare() in addition.

Common workflows

The pipe operator ‘%>%’ allows us to chain together subsequent operations in R. This is useful, as the dta function expects preprocessed data indicating correct (1) and false (0) predictions.

M %>%
  categorize() %>%
  compare(y) %>%
  evaluate()
#> [cases] evaluation results:
#> $specificity
#>   parameter hypothesis estimate  lower  upper   pval pval_all reject reject_all
#> 1      V1_0    == 0.75      0.4 0.0080 0.7920 0.0401   0.2266  FALSE      FALSE
#> 2      V2_0    == 0.75      0.8 0.4799 1.1201 0.3797   0.3797  FALSE      FALSE
#> 3      V3_0    == 0.75      0.6 0.2080 0.9920 0.2266   0.2266  FALSE      FALSE
#> 
#> $sensitivity
#>   parameter hypothesis estimate lower upper   pval pval_all reject reject_all
#> 1      V1_0    == 0.75      0.6 0.208 0.992 0.2266   0.2266  FALSE      FALSE
#> 2      V2_0    == 0.75      0.6 0.208 0.992 0.2266   0.3797  FALSE      FALSE
#> 3      V3_0    == 0.75      0.6 0.208 0.992 0.2266   0.2266  FALSE      FALSE

Multiple testing for co-primary endpoints

Specification of hypotheses

The R command

?evaluate 

gives an overview over the function arguments of the evaluate function.

Together this implies the hypotheses system that is considered, namely

\(H_0: \forall g \forall j: \theta_j^g \leq \theta_0^g\)

In the application of primary interest, diagnostic accuracy studies, this simplifies to \(G=2\) with \(\theta_1 = Se\) and \(\theta_2 =Sp\) indicating sensitivity and specificity of a medical test or classication rule. In this case we aim to reject the global null hypothesis

\(H_0: \forall j: Se_j \leq Se_0 \wedge Sp_j \leq Sp_0\)

Comparison vs. confidence regions

In the following, we highlight the difference between the “co-primary” analysis (comparison regions) and a “full” analysis (confidence regions).

set.seed(1337)

data <- draw_data_roc(n=120, prev=c(0.25, 0.75), m=4,
                      delta=0.05, e=10, auc=seq(0.90, 0.95, 0.025), rho=c(0.25, 0.25))

head(data)
#> $specificity
#>       model1 model2 model3 model4
#>  [1,]      1      1      1      1
#>  [2,]      0      0      1      1
#>  [3,]      1      1      1      1
#>  [4,]      1      1      1      1
#>  [5,]      1      1      1      1
#>  [6,]      1      1      1      1
#>  [7,]      1      1      1      1
#>  [8,]      1      1      1      1
#>  [9,]      1      1      1      1
#> [10,]      1      1      1      1
#> [11,]      1      1      1      1
#> [12,]      1      1      1      1
#> [13,]      1      1      1      1
#> [14,]      1      1      1      1
#> [15,]      1      1      0      0
#> [16,]      0      0      0      0
#> [17,]      0      1      1      1
#> [18,]      1      1      1      1
#> [19,]      1      1      1      1
#> [20,]      1      1      0      0
#> [21,]      1      1      1      1
#> [22,]      0      0      1      1
#> [23,]      1      1      1      1
#> [24,]      1      1      1      1
#> [25,]      1      1      1      1
#> [26,]      1      1      1      1
#> [27,]      0      0      1      1
#> [28,]      1      1      1      1
#> [29,]      1      1      1      1
#> [30,]      1      1      0      0
#> [31,]      1      1      1      1
#> [32,]      1      1      1      1
#> [33,]      1      1      1      1
#> [34,]      1      1      1      1
#> [35,]      0      0      1      1
#> [36,]      1      1      1      1
#> [37,]      0      0      1      1
#> [38,]      1      1      1      1
#> [39,]      1      1      1      1
#> [40,]      1      1      1      1
#> [41,]      1      1      1      1
#> [42,]      1      1      1      1
#> [43,]      1      1      1      1
#> [44,]      0      0      1      1
#> [45,]      1      1      1      1
#> [46,]      1      1      1      1
#> [47,]      1      1      1      1
#> [48,]      1      1      1      1
#> [49,]      1      1      1      1
#> [50,]      0      0      1      1
#> [51,]      1      1      1      1
#> [52,]      1      1      1      1
#> [53,]      1      1      1      1
#> [54,]      1      1      1      1
#> [55,]      1      1      1      1
#> [56,]      1      1      1      1
#> [57,]      1      1      1      1
#> [58,]      1      1      1      1
#> [59,]      1      1      1      1
#> [60,]      1      1      1      1
#> [61,]      1      1      1      1
#> [62,]      0      1      1      1
#> [63,]      1      1      1      1
#> [64,]      1      1      1      1
#> [65,]      0      1      1      1
#> [66,]      1      1      1      1
#> [67,]      0      0      1      1
#> [68,]      1      1      1      1
#> [69,]      1      1      1      1
#> [70,]      1      1      1      1
#> [71,]      1      1      1      1
#> [72,]      1      1      1      1
#> [73,]      0      1      0      0
#> [74,]      1      1      1      1
#> [75,]      1      1      1      1
#> [76,]      0      0      1      1
#> [77,]      1      1      1      1
#> [78,]      0      1      1      1
#> [79,]      1      1      1      1
#> [80,]      1      1      1      1
#> [81,]      1      1      1      1
#> [82,]      0      0      1      1
#> [83,]      0      1      1      1
#> [84,]      1      1      1      1
#> [85,]      1      1      1      1
#> [86,]      1      1      1      1
#> [87,]      1      1      1      1
#> [88,]      1      1      1      1
#> [89,]      1      1      1      1
#> [90,]      1      1      1      1
#> 
#> $sensitivity
#>       model1 model2 model3 model4
#>  [1,]      0      0      1      1
#>  [2,]      1      1      1      1
#>  [3,]      1      1      1      1
#>  [4,]      1      1      0      0
#>  [5,]      1      1      1      1
#>  [6,]      1      1      1      1
#>  [7,]      1      1      1      1
#>  [8,]      1      1      1      1
#>  [9,]      1      1      1      1
#> [10,]      1      0      1      1
#> [11,]      1      1      1      0
#> [12,]      1      1      1      1
#> [13,]      1      1      0      0
#> [14,]      1      1      0      0
#> [15,]      1      1      1      1
#> [16,]      1      1      1      1
#> [17,]      1      1      1      1
#> [18,]      1      1      0      0
#> [19,]      1      1      0      0
#> [20,]      1      1      1      1
#> [21,]      1      1      1      1
#> [22,]      1      1      1      1
#> [23,]      1      1      1      1
#> [24,]      1      1      1      1
#> [25,]      1      1      1      1
#> [26,]      1      1      1      1
#> [27,]      1      1      1      1
#> [28,]      1      1      1      1
#> [29,]      1      0      0      0
#> [30,]      1      1      1      1
## comparison regions
results_comp <- data %>% evaluate(alternative ="greater",
                                  alpha=0.025,
                                  benchmark = c(0.7, 0.8),
                                  analysis = "co-primary",
                                  regu = TRUE,
                                  adj = "maxt")
visualize(results_comp)

## confidence regions
results_conf <- data %>% evaluate(alternative = "greater",
                                  alpha = 0.025,
                                  benchmark = c(0.7, 0.8),
                                  analysis = "full",
                                  regu = TRUE,
                                  adj = "maxt")
visualize(results_conf)

As we can see, the comparison regions are more liberal compared to the confidence regions.

Real data example

A second vignette shows an application of the cases package to the Breast Cancer Wisconsin Diagnostic (wdbc) data set.

vignette("example_wdbc", "cases")

References

These binaries (installable software) and packages are in development.
They may not be fully stable and should be used with caution. We make no claims about them.
Health stats visible at Monitor.