Abordagem estatística para comparar a calibração entre modelos

Parece um problema comum, mas não consigo encontrar uma solução.

Eu tenho um conjunto de observações binárias e dois modelos diferentes, cada um com previsões para cada observação. Eu quero comparar a calibração dos modelos.

Existem várias abordagens para comparar a discriminação desses modelos (ou seja, consulte o roc.test no pacote pROC em R), mas nenhuma abordagem para comparar a calibração. A maioria dos trabalhos empíricos apenas lista os valores p de dois testes de calibração diferentes que estão testando se a calibração de cada modelo está desativada (por exemplo, Hosmer-Lemeshow, pontuação Brier).

O que estou procurando é uma comparação estatística direta da calibração entre dois modelos.

Aqui está um conjunto de dados de teste extremo. Os valores do teste Brier, Spiegelhalter Z-test, etc, suportam que p2 seja melhor calibrado, e sabemos que é. Alguém pode transformar isso em um teste estatístico formal?

library("pROC")
y <- rbinom(100,1,1:100/100)
p1 <- 1:100/10001
p2 <- 1:100/101
val.prob(p1,y)
val.prob(p2,y)

Como você sabe, a pontuação Brier mede a calibração e é o erro quadrado médio, , entre as previsões, e as respostas, . Como a pontuação de Brier é uma média, a comparação de duas pontuações de Brier é basicamente uma comparação de médias e você pode usar a fantasia como quiser. Vou sugerir duas coisas e apontar para uma terceira:$\bar B = n^{-1} \sum (\hat y_i - y_i)^2$$\hat y,$$y$

Uma opção: faça um teste t

Minha resposta imediata quando ouço comparações de meios é fazer um teste t. Os erros ao quadrado provavelmente não são normalmente distribuídos em geral, portanto, é possível que este não seja o teste mais poderoso. Parece bom no seu exemplo extremo. Abaixo, testo a hipótese alternativa que p1tem maior MSE do que p2:

y <- rbinom(100,1,1:100/100)
p1 <- 1:100/10001
p2 <- 1:100/101

squares_1 <- (p1 - y)^2
squares_2 <- (p2 - y)^2

t.test(squares_1, squares_2, paired=T, alternative="greater")
#> 
#>  Paired t-test
#> 
#> data:  squares_1 and squares_2
#> t = 4.8826, df = 99, p-value = 2.01e-06
#> alternative hypothesis: true difference in means is greater than 0
#> 95 percent confidence interval:
#>  0.1769769       Inf
#> sample estimates:
#> mean of the differences 
#>               0.2681719

Temos um valor p super-baixo. Fiz um teste t emparelhado, pois, observação por observação, os dois conjuntos de previsões se comparam com o mesmo resultado.

Outra opção: teste de permutação

Se a distribuição dos erros ao quadrado o preocupa, talvez você não queira fazer suposições de um teste t. Você poderia, por exemplo, testar a mesma hipótese com um teste de permutação:

library(plyr)

observed <- mean(squares_1) - mean(squares_2)
permutations <- raply(500000, {
  swap <- sample(c(T, F), 100, replace=T)
  one <- squares_1
  one[swap] <- squares_2[swap]

  two <- squares_2
  two[swap] <- squares_1[swap]

  mean(one) - mean(two)
})

hist(permutations, prob=T, nclass=60, xlim=c(-.4, .4))
abline(v=observed, col="red")

# p-value. I add 1 so that the p-value doesn't come out 0
(sum(permutations > observed) + 1)/(length(permutations) + 1) 
#> [1] 1.999996e-06

Os dois testes parecem concordar estreitamente.

Algumas outras respostas

Uma pesquisa rápida deste site na comparação de MPEs aponta para o teste de Diebold-Mariano (veja a resposta aqui e um comentário aqui ). Parece que é simplesmente o teste de Wald e acho que ele terá um desempenho semelhante ao teste t acima.

Se eu entendi direito, você quer uma maneira de comparar dois modelos de regressão logística ou qualquer alternativa para modelar resultados binários.

Para mim, é importante ver que a maneira 'correta' de comparar modelos depende do objetivo de sua análise.

Se apenas a previsão binária (sim / não) é importante, um modelo que preveja p = 0,51 para cada caso que seja efetivamente verdadeiro e preveja p = 0,49 para cada caso que seja efetivamente falso, será perfeito, enquanto o escore Brier não será que bom. Nesse caso, eu compararia modelos com base em% de previsão binária correta.

Além disso, pode ser que um falso positivo seja pior que um falso negativo. Você pode definir uma função de pontuação que incorpore esse recurso (compare a previsão binária, mas com uma penalidade maior por um falso positivo).

Obviamente, se é importante prever a probabilidade tão boa quanto possível, medidas como o Brier-Score são melhores.

Finalmente, se a previsão é o objetivo (binário ou probabilidade), sempre seria considerado o uso de validação cruzada no cálculo das pontuações. É mais interessante avaliar como um modelo prevê 'novos' dados em vez dos dados de treinamento em si.

Para referência futura, a IMO, a primeira resposta, não trata do problema de calibração. Considere previsões$\hat{y}_1,\hat{y}_2 ..., \hat{y}_n$ feita por um modelo razoável e bem calibrado para valores de entrada $x_1, x_2,..., x_n$. Agora considere um segundo conjunto de previsões$\tilde{y}_1, \tilde{y}_2, ..., \tilde{y}_n$que são feitas por um modelo que simplesmente embaralha as previsões do primeiro modelo em cada uma das duas classes e as produz em ordem aleatória. É provável que o segundo modelo seja mal calibrado em comparação com o primeiro modelo bem calibrado, mas as pontuações dos brier dos dois modelos serão as mesmas.

Como afirmado na pergunta original, sugiro examinar o teste de Hosmer – Lemeshow e comparar as estatísticas de teste de HL calculadas para as previsões de cada um dos dois modelos (uma estatística de HL maior sugere uma calibração mais fraca).