2
$\begingroup$

I have created a machine learning model (let's call it Model B) which I hope is better at classifying disease-states than the currently used model (Model A).

I have a dataset of 600 cases (roughly split equally between disease and control) and I've split this into a Training dataset (n=450) and a Validation dataset (n=150).

[The Validation dataset was completely withheld until the end of Model B's creation, so that Model B could be tested against an unseen dataset, thus testing the Model's ability to generalise]

Case 1 : CV on Training dataset

I created Model B using the Training dataset for feature selection. I then performed a 10-fold CV on Training with 1000 repeats. This involved Training being randomly split into a Train (80%) and Test (20%) datasets. Both Model A and B trained on the Train and tested on the Test, and I obtained results on both Models' sensitivity, specificity, AUC and kappa values. As there were 1000 repeats, and both Models were tested on the same Train/Test splits each time, I believe a paired Student T test is sufficient to determine if there is a statistically significant difference between the metrics I tested in Model A vs Model B.

Case 2 : Testing the model on an unseen dataset

I then ran Model B, which I trained on the whole Training dataset (n=450) on the previously unseen Validation dataset (n=150). I did the same with Model A.

The results of the confusion matrix I obtained are below :

  • Model A

truth prediction case control case 71 8 control 17 54

  • Model B

    truth prediction case control case 82 9 control 6 53

So if Model A has Cohen's kappa of 0.66 and Model B has 0.79, what test can I run to check if there's a statistically significant difference? (and the same question for sensitivity, specificity etc.) Because I only run this test once (running it 1000 times would just result in 1000 of the same confusion matrices) I'm not sure if the T test would work as there's no variance. I only have a single value for each metric (i.e sensitivity = 81% vs 93%; specificity = 87% vs 85%).

1. Am I right in using a paired T test to compare results of the 10-fold CV when I ran both models 1000 times within the Training dataset?

2. How can I test if there's a statistically significant difference between the two models from a single confusion matrix?

$\endgroup$
3
  • 1
    $\begingroup$ I'd start reading here stats.stackexchange.com/questions/312780/… and considering whether accuracy is really the important figure of merit for your problem if you are actually more concerned with an alternative measure of model quality. $\endgroup$ – Sycorax Feb 13 '20 at 20:50
  • $\begingroup$ Thanks very much for this, you're quite right that accuracy isn't a useful measure but I'm not looking at accuracy in my comparisons : rather, I'm aiming to compare the two models based on their AUC-ROC (area under the receiver operator curve) and Cohen's kappa values. $\endgroup$ – schmitzhugen Feb 13 '20 at 22:36
  • $\begingroup$ Does this answer your question? stats.stackexchange.com/questions/358101/… $\endgroup$ – Sycorax Feb 15 '20 at 23:08
1
$\begingroup$

I suggest to adapt a forecast encompassing test, e.g. in the way this paper does it. They also have a formal test of significance of difference between two forecasts.

The idea of encompassing test is to say whether one model forecast encompasses all the relevant info in the other. One way to do it is to run a regression: $Y=\alpha+\beta_1\hat Y_1+\beta_2\hat Y_2+\varepsilon,$ where $\hat Y_i$ are model forecasts. In this case if $\beta_1$ is significant and $\beta_2$ not significant then Model 2 encompasses Model 1.

$\endgroup$
1
$\begingroup$

You could use the bootstrap to see if the sampling distribution of the AUC (or Kappa) of each model is different.

First fit some "AI" models using glmnet.

library("glmnet")
library("pROC")
library("boot")
library("reshape2")
library("ggplot2")

load(system.file("data/BinomialExample.RData", package="glmnet"))

ridge <- cv.glmnet(x, y, family = "binomial", alpha = 0)
lasso <- cv.glmnet(x, y, family = "binomial", alpha = 1)

We can look at confusion matrices and AUC for each.

yhat_ridge <- predict(ridge, newx = x, type = "class")
conf_ridge <- table(yhat_ridge, y)
conf_ridge
#>           y
#> yhat_ridge  0  1
#>          0 39  1
#>          1  5 55
roc_ridge <- roc(y, as.numeric(yhat_ridge))
#> Setting levels: control = 0, case = 1
#> Setting direction: controls < cases
auc_ridge <- auc(roc_ridge)
auc_ridge
#> Area under the curve: 0.9343

yhat_lasso <- predict(lasso, newx = x, type = "class")
conf_ridge <- table(yhat_lasso, y)
conf_ridge
#>           y
#> yhat_lasso  0  1
#>          0 40  2
#>          1  4 54
roc_lasso <- roc(y, as.numeric(yhat_lasso))
#> Setting levels: control = 0, case = 1
#> Setting direction: controls < cases
auc_lasso <- auc(roc_lasso)
auc_lasso
#> Area under the curve: 0.9367

Then we can bootstrap the AUC by resampling with replacement the predictions (and true labels). We can then compare the sampling distributions.

boot_lasso <- boot(yhat_lasso,
function(data, ind) {
    suppressMessages({
      roc <- roc(y[ind], as.numeric(yhat_lasso[ind]))
      auc(roc)
    })
  },
  R = 1000
)

boot_ridge <- boot(yhat_ridge,
  function(data, ind) {
    suppressMessages({
      roc <- roc(y[ind], as.numeric(yhat_ridge[ind]))
      auc(roc)
    })
  },
  R = 1000
)

df <- data.frame(lasso = boot_lasso$t, ridge = boot_ridge$t)
mdf <- melt(df)
#> No id variables; using all as measure variables
ggplot(mdf, aes(x = value, color = variable)) +
  geom_density()

Alternatively we can look at the sampling distribution of the difference in AUC.

boot_diff <- boot(yhat_ridge,
  function(data, ind) {
    suppressMessages({
      roc_ridge <- roc(y[ind], as.numeric(yhat_ridge[ind]))
      roc_lasso <- roc(y[ind], as.numeric(yhat_lasso[ind]))
      auc(roc_ridge) - auc(roc_lasso)
    })
  },
  R = 1000
)

plot(density(boot_diff$t))
abline(v = 0, lty = "dashed")

Somewhat unsurprisingly, fitting the same model with ridge and lasso doesn't give massively different results!

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.