Statistischer Ansatz zum Vergleich der Kalibrierung zwischen Modellen

Dies scheint ein häufiges Problem zu sein, aber ich kann keine Lösung finden.

Ich habe eine Reihe von binären Beobachtungen und zwei verschiedene Modelle, jedes mit Vorhersagen für jede Beobachtung. Ich möchte die Kalibrierung der Modelle vergleichen.

Es gibt verschiedene Ansätze zum Vergleichen der Unterscheidung dieser Modelle (dh siehe den roc.test im pROC-Paket in R), aber keinen Ansatz zum Vergleichen der Kalibrierung. Die meisten empirischen Arbeiten listen nur die p-Werte von zwei verschiedenen Kalibrierungstests auf, die testen, ob die Kalibrierung jedes Modells deaktiviert ist (dh Hosmer-Lemeshow, Brier-Score).

Was ich suche, ist ein direkter statistischer Vergleich der Kalibrierung zwischen zwei Modellen.

Hier ist ein extremer Testdatensatz. Die Werte des Brier-Tests, des Spiegelhalter-Z-Tests usw. unterstützen alle, dass p2 besser kalibriert ist, und wir wissen, dass dies der Fall ist. Kann jemand daraus einen formalen statistischen Test machen?

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)

Wie Sie wissen, misst der Brier-Score die Kalibrierung und ist der mittlere quadratische Fehler zwischen den Vorhersagen und den Antworten . Da der Brier-Score ein Mittelwert ist, ist der Vergleich zweier Brier-Scores im Grunde ein Vergleich der Mittelwerte, und Sie können damit so ausgefallen sein, wie Sie möchten. Ich werde zwei Dinge vorschlagen und auf ein drittes hinweisen:$\bar B = n^{-1} \sum (\hat y_i - y_i)^2$$\hat y,$$y$

Eine Option: Machen Sie einen T-Test

Meine unmittelbare Antwort, wenn ich Vergleiche von Mitteln höre, ist ein T-Test. Quadratische Fehler werden wahrscheinlich im Allgemeinen nicht verteilt, daher ist dies möglicherweise nicht der leistungsstärkste Test. In Ihrem extremen Beispiel scheint es in Ordnung zu sein. Im Folgenden teste ich die alternative Hypothese, die p1eine größere MSE aufweist als 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

Wir bekommen einen super niedrigen p-Wert. Ich habe einen gepaarten T-Test durchgeführt, da die beiden Vorhersagesätze Beobachtung für Beobachtung mit demselben Ergebnis verglichen werden.

Eine weitere Option: Permutationstests

Wenn Sie sich über die Verteilung der quadratischen Fehler Sorgen machen, möchten Sie möglicherweise keine Annahmen über einen T-Test treffen. Sie könnten beispielsweise dieselbe Hypothese mit einem Permutationstest testen:

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

Die beiden Tests scheinen eng übereinzustimmen.

Einige andere Antworten

Eine schnelle Suche auf dieser Website zum Vergleich von MSEs weist auf den Diebold-Mariano-Test hin (siehe Antwort hier und Kommentar hier ). Das sieht so aus, als wäre es einfach Walds Test und ich denke, er wird ähnlich wie der obige T-Test funktionieren.

Wenn ich es richtig verstehe, möchten Sie zwei Modelle einer logistischen Regression oder eine Alternative zur Modellierung des binären Ergebnisses vergleichen.

Für mich ist es wichtig zu sehen, dass die „richtige“ Art, Modelle zu vergleichen, vom Ziel Ihrer Analyse abhängt.

Wenn nur die binäre Vorhersage (Ja / Nein) von Bedeutung ist, ist ein Modell, das p = 0,51 für jeden Fall vorhersagt, der effektiv wahr ist, und p = 0,49 für jeden Fall, der effektiv falsch ist, perfekt, während der Brier-Score nicht stimmt so gut. In diesem Fall würde ich Modelle vergleichen, die auf einer% korrekten binären Vorhersage basieren.

Außerdem kann es sein, dass ein falsches Positiv schlechter ist als ein falsches Negativ. Sie können eine Score-Funktion definieren, die diese Funktion enthält (vergleichen Sie die binäre Vorhersage, jedoch mit einer größeren Strafe für ein falsches Positiv).

Wenn es wichtig ist, die Wahrscheinlichkeit als solche so gut wie möglich vorherzusagen, sind Maßnahmen wie der Brier-Score natürlich besser.

Wenn die Vorhersage das Ziel ist (binär oder Wahrscheinlichkeit), wird bei der Berechnung der Punktzahlen immer die Verwendung einer Kreuzvalidierung in Betracht gezogen. Es ist interessanter zu beurteilen, wie ein Modell "neue" Daten anstelle der Trainingsdaten selbst vorhersagt.

Zum späteren Nachschlagen, IMO, wird in der ersten Antwort das Kalibrierungsproblem nicht behandelt. Betrachten Sie Vorhersagen die von einem vernünftigen, gut kalibrierten Modell für Eingabewerte . Betrachten Sie nun einen zweiten Satz von Vorhersagen , die von einem Modell erstellt wurden, das einfach die Vorhersagen des ersten Modells in jedem der Modelle verschlüsselt zwei Klassen und gibt sie in zufälliger Reihenfolge aus. Das zweite Modell ist im Vergleich zum ersten gut kalibrierten Modell wahrscheinlich schlecht kalibriert, aber die Brier-Scores der beiden Modelle sind gleich.$\hat{y}_1,\hat{y}_2 ..., \hat{y}_n$$x_1, x_2,..., x_n$$\tilde{y}_1, \tilde{y}_2, ..., \tilde{y}_n$

Wie in der ursprünglichen Frage angegeben, schlage ich vor, den Hosmer-Lemeshow-Test zu betrachten und die für die Vorhersagen jedes der beiden Modelle berechneten HL-Teststatistiken zu vergleichen (eine größere HL-Statistik deutet auf eine schlechtere Kalibrierung hin).