Warum liefert heteroskedastische AICc/n niedrigere Werte als Einkompartiment-Modelle? Probleme bei der Parameterzählung ⇐ Python

[Anonymous]

Ich passe Kompartimentmodelle an biologische Daten an und stoße auf ein kontraintuitives Ergebnis: Ein Modell, das die Anpassung dramatisch verbessert (18 % NRMSE-Reduktion), hat einen schlechteren AICc als ein Modell mit minimaler Verbesserung (0,63 % NRMSE-Reduktion), obwohl beide die gleiche Anzahl von Parametern haben. Ich glaube, dass dies mit unterschiedlichen Varianzskalen in den verschiedenen Kompartimenten zusammenhängt, aber ich suche nach einer Klärung des richtigen statistischen Ansatzes. Ich habe Burnham & Anderson (2002) gelesen, konnte aber keine spezifische Hilfe zu meinem Fall finden.
Hintergrund
Ich habe drei Datensätze aus demselben System:
Datensatz 1 (Sagen wir A): Einzelnes Kompartiment (C1), n=19 (19 experimentelle Datenpunkte)
Datensatz 2 (B): Einzelnes Kompartiment (C1), n=19
Datensatz 3 (C): Zwei Kompartimente (C1 + P1), n=38 (19 + 19)
Ich verwende die Standardkorrektur für kleine Stichproben:
AICc = -2*log(L) + 2K + (2K(K+1)/n-k-1)
Das Skalierungsverhältnis ist ungefähr 130:1 zwischen C1 und P1. Ich möchte dasselbe Modell in allen drei Datensätzen vergleichen, um festzustellen, welcher Datensatz am besten passt. Um Vergleiche aussagekräftig zu machen, berechne ich AICc pro Beobachtung (AICc/n) als normalisierte Metrik. Ich vergleiche zwei Modelle mit jeweils 1 freien Parameter (k=1):



Modell
NRMSE-Reduktion vs. Basislinie
Erwartet Ergebnis
Tatsächlicher AICc




Modell A
-0,63 %
Höher AICc
411,90


Modell B
-18,06 %
Unterer AICc
412,87



Modell B passt viel besser, hat aber einen schlechteren AICc. Warum?
Ich habe zunächst einen homoskedastischen Ansatz (gepoolte Varianz) verwendet:

Code: Select all

all_residuals = concatenate([compartment1_residuals, compartment2_residuals])
n_total = 38

# Single pooled variance
rss_total = sum(all_residuals**2)
sigma2_pooled = rss_total / n_total

# Log-likelihood
logL = -n_total/2 * (log(2π * sigma2_pooled) + 1)

# AICc
k = 1 + 1  # parameter + variance
AICc = 2*k - 2*logL + 2*k*(k+1)/(n_total - k - 1)

Ein Blick auf die aktuellen RSS-Beiträge:
Modell A:

Code: Select all

C1 RSS:     76,750  (99.98% of total)
P1 RSS:       14  (0.02% of total)
Total RSS:   76,764

Modell B:

Code: Select all

C1 RSS:     79,153  (99.98% of total)  [+2,403 worse!]
P1 RSS:       10  (0.01% of total)  [-4 better]
Total RSS:   79,163  [Net: +2,399 worse]

Ich vermutete also, dass die P1-Verbesserung statistisch gesehen wertlos ist, da P1 0,02 % zum gesamten RSS beiträgt. Deshalb dachte ich, dass ich für jedes Kompartiment separate Varianzen benötige, da sie sehr unterschiedliche Skalen haben und ich einen heteroskedastischen Ansatz verwendet habe.

Code: Select all

# Separate variances for each compartment
sigma2_csf = rss_csf / n_csf
sigma2_plasma = rss_plasma / n_plasma

# Separate log-likelihoods
logL_csf = -n_csf/2 * (log(2π * sigma2_csf) + 1)
logL_plasma = -n_plasma/2 * (log(2π * sigma2_plasma) + 1)

# Sum log-likelihoods
logL_total = logL_csf + logL_plasma

# AICc (but what is k?)
k = 1
AICc = 2*k - 2*logL_total + 2*k*(k+1)/(n_total - k - 1)

Beim heteroskedastischen Ansatz zeigt Modell B jetzt einen niedrigeren AICc.



Modell
NRMSE-Reduktion vs. Basislinie
Erwartetes Ergebnis
Tatsächlicher AICc




Modell A
-0,63 %
Höherer AICc
255,12


Modell B
-18,06 %
Niedriger AICc
239,64



Der heteroskedastische Ansatz stuft Modell B korrekt als besser ein (entspricht NRMSE). Aber jetzt bricht der datensatzübergreifende Vergleich ab!
Ich muss das dasselbe Modell über mehrere Datensätze hinweg vergleichen, um herauszufinden, wo es am besten passt. Hier sind AICc/n-Werte für Modell A:



Datensatz
Struktur
Ansatz
Roh-AICc
n
AICc/n
k




Datensatz 1
1 Fach
Standard
212.11
19
11.16
1


Datensatz 2
1 Fach
Standard
209.21
19
11.01
1


Datensatz 3
2 Kompartimente
Homoskedastisch
411,90
38
10,84
1


Datensatz 3
2 Kompartimente
Heteroskedastisch
255,02
38
6,71
1



Problem: Beim heteroskedastischen Ansatz sinkt der AICc/n von Datensatz 3 auf den Bereich 6–7, während die Datensätze 1 und 2 im Bereich 10–11 liegen. Beim homoskedastischen Ansatz liefern alle drei Datensätze vergleichbare AICc/n-Werte (10-11), aber dies führt zu einem falschen Modellranking für Datensatz 3.
Die Standardreferenz (Burnham & Anderson 2002) befasst sich nicht explizit mit diesem Fall, und ich bitte die Statistik-Community um Rat zu Best Practices. Ich wäre wirklich dankbar, wenn ich Hilfe und Rat bekommen könnte, um herauszufinden, ob ich es richtig mache oder nicht? Und wenn der heteroskedastische Ansatz richtig ist, wie geht man dann mit dem Problem des unteren Bereichs um?