ANOVA zum Vergleichen von Modellen

Ich suche auf dieser Seite nach einem Workshop zu GAM in R: http://qcbs.ca/wiki/r_workshop8

Am Ende des Abschnitts 2. Multiple smooth termszeigen sie ein Beispiel, in anovadem drei verschiedene Modelle verglichen werden, um das am besten passende Modell zu ermitteln. Die Ausgabe ist

  Analysis of Deviance Table
  Model 1: y ~ x0 + s(x1)
  Model 2: y ~ x0 + s(x1) + x2
  Model 3: y ~ x0 + s(x1) + s(x2)
    Resid. Df Resid. Dev      Df Deviance  Pr(>Chi)    
  1    394.08     5231.6                               
  2    393.10     4051.3 0.97695   1180.2 < 2.2e-16 ***
  3    385.73     1839.5 7.37288   2211.8 < 2.2e-16 ***
  ---
  Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Daraus schließen sie, dass Modell 3 am besten ist. Meine Frage ist, wie sie das sehen?

Mein derzeitiges Verständnis ist: Der Pr(>Chi)-Wert ist sowohl für Modell 2 als auch für Modell 3 klein, daher sind diese besser als Modell 1. Welche andere Variable verwenden sie jedoch, um zu bestimmen, dass 3 besser als 2 ist?

Die Ausgabe von anova()ist eine Reihe von Likelihood-Ratio-Tests. Die Zeilen in der Ausgabe sind:

1. Die erste Zeile in der Ausgabe entspricht dem einfachsten Modell mit nur einer Glättung von x1(ich ignoriere den Faktor, x0da er in Ihrem Beispiel nicht berücksichtigt werden muss) - dies wird nicht gegen etwas Einfacheres getestet, daher sind die letzten Spalteneinträge leeren.

2. Die zweite Zeile ist ein Likelihood-Ratio-Test zwischen dem Modell in Zeile 1 und dem Modell in Zeile 2. Auf Kosten 0.97695zusätzlicher Freiheitsgrade wird die Restabweichung um verringert 1180.2. Diese Verringerung der Abweichung (oder umgekehrt die Erhöhung der Abweichung erklärt) auf Kosten von <1 Freiheitsgrad ist höchst unwahrscheinlich, wenn der wahre Effekt von x20 wäre.

   Warum nehmen die 0.97695Freiheitsgrade zu? Nun, die lineare Funktion von x2würde dem Modell 1 df hinzufügen, aber die glattere Funktion x1wird etwas mehr als zuvor bestraft und verwendet daher etwas weniger effektive Freiheitsgrade, daher die <1-Änderung der Gesamtfreiheitsgrade.

3. Die dritte Zeile ist genau die gleiche wie oben beschrieben, jedoch für einen Vergleich zwischen dem Modell in der zweiten Zeile und dem Modell in der dritten Zeile: dh die dritte Zeile bewertet die Verbesserung beim Übergang von der Modellierung x2als linearem Begriff zur Modellierung x2als glatte Funktion. Auch diese Verbesserung der Modellanpassung (Änderung der Abweichung geht jetzt 2211.8zu Lasten von 7.37288mehr Freiheitsgraden) ist unwahrscheinlich, wenn die damit verbundenen zusätzlichen Parameter s(x2)alle gleich 0 wären .

Zusammenfassend sagt Zeile 2, dass Modell 2 besser passt als Modell 1, sodass eine lineare Funktion von x2besser ist als kein Effekt von x1. In Zeile 3 heißt es jedoch, dass Modell 3 besser zu den Daten passt als Modell 2, sodass eine glatte Funktion von x2einer linearen Funktion von vorgezogen wird x2. Dies ist eine sequentielle Analyse von Modellen, keine Reihe von Vergleichen mit dem einfachsten Modell.

Jedoch…

Was sie zeigen, ist nicht der beste Weg, dies zu tun - die jüngste Theorie würde vorschlagen, dass die Ausgabe von summary(m3)die "korrektesten" Abdeckungseigenschaften aufweist. Um zwischen Modellen zu wählen, sollte man wahrscheinlich select = TRUEbeim Anpassen des vollständigen Modells (das mit zwei Glättungen) verwenden, was eine Schrumpfung von Begriffen ermöglichen würde, die das Modell mit linearer x2oder sogar keiner Auswirkung dieser Variablen einschließen würden . Sie passen auch nicht zur REML- oder ML-Glättungsauswahl, die viele von uns mgcv- Benutzern als Standardoption betrachten würden (obwohl dies nicht die tatsächliche Standardeinstellung ist gam()).

Was ich tun würde ist:

library("mgcv")
gam_data <- gamSim(eg=5)
m3 <- gam(y ~ x0 + s(x1) + s(x2), data = gam_data, select = TRUE,
          method = "REML")
summary(m3)

Die letzte Zeile erzeugt Folgendes:

> summary(m3)

Family: gaussian 
Link function: identity 

Formula:
y ~ x0 + s(x1) + s(x2)

Parametric coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   8.4097     0.2153  39.053  < 2e-16 ***
x02           1.9311     0.3073   6.284 8.93e-10 ***
x03           4.4241     0.3052  14.493  < 2e-16 ***
x04           5.7639     0.3042  18.948  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Approximate significance of smooth terms:
        edf Ref.df     F p-value    
s(x1) 2.487      9 25.85  <2e-16 ***
s(x2) 7.627      9 76.03  <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

R-sq.(adj) =  0.769   Deviance explained = 77.7%
-REML = 892.61  Scale est. = 4.5057    n = 400

Wir können sehen, dass sich beide glatten Terme signifikant von Nullfunktionen unterscheiden.

Was select = TRUEgetan wird, ist eine zusätzliche Strafe auf den Nullraum der Strafe zu setzen (dies ist der Teil des Splines, der vollkommen glatt ist). Wenn Sie dies nicht haben, kann die Glättungsauswahl nur einen glatten Rücken zu einer linearen Funktion bestrafen (da die Strafe für die Glättungsauswahl nur für die nicht glatten (wackeligen) Teile der Basis gilt). Um eine Auswahl durchführen zu können, müssen wir auch den Nullraum (die glatten Teile der Basis) bestrafen können.

select = TRUEDies wird durch die Verwendung einer zweiten Strafe erreicht, die zu allen glatten Begriffen im Modell hinzugefügt wird (Marra und Wood, 2011). Dies wirkt als eine Art Schrumpfung, bei der alle glatten Terme etwas gegen 0 gezogen werden, aber überflüssige Terme werden viel schneller gegen 0 gezogen, sodass sie aus dem Modell ausgewählt werden, wenn sie keine Erklärungskraft haben. Wir zahlen dafür einen Preis, wenn wir die Bedeutung der Glättungen bewerten. Beachten Sie die Ref.dfobige Spalte (die 9 stammt aus dem Standardwert von k = 10, was für dünne Platten-Splines mit Zentrierungsbeschränkungen 9 Basisfunktionen bedeutet), anstatt etwa 2,5 und 7,7 Freiheitsgrade für die Splines zu zahlen, zahlen wir 9 Grad von Freiheit jeder. Dies spiegelt die Tatsache wider, dass wir die Auswahl getroffen haben und nicht sicher waren, welche Begriffe im Modell enthalten sein sollten.

Hinweis: Es ist wichtig, dass Sie keine anova(m1, m2, m3)Typaufrufe für Modelle verwenden, die verwenden select = TRUE. Wie in erwähnt ?mgcv:::anova.gam, kann die verwendete Annäherung für Glättungen mit Strafen für ihre Nullräume sehr schlecht sein.

In den Kommentaren erwähnte @BillyJean die Verwendung von AIC zur Auswahl. Jüngste Arbeiten von Simon Wood und Kollegen (Wood et al., 2016) haben einen AIC abgeleitet, der die zusätzliche Unsicherheit berücksichtigt, da wir die Glättungsparameter im Modell geschätzt haben. Dieser AIC funktioniert ziemlich gut, aber es gibt einige Diskussionen über das Verhalten ihrer Ableitung von AIC, wenn IIRC-Glättungen nahe an linearen Funktionen liegen. Wie auch immer, AIC würde uns geben:

m1 <- gam(y ~ x0 + s(x1), data = gam_data, method = "ML")
m2 <- gam(y ~ x0 + s(x1) + x2, data = gam_data, method = "ML")
m3 <- gam(y ~ x0 + s(x1) + s(x2), data = gam_data, method = "ML")
AIC(m1, m2, m3)

> AIC(m1, m2, m3)
          df      AIC
m1  7.307712 2149.046
m2  8.608444 2055.651
m3 16.589330 1756.890

Hinweis: Ich habe all dies mit der ML-Glättungsauswahl ausgestattet, da ich nicht sicher bin, was der AIC wann tut, select = TRUEund Sie müssen vorsichtig sein, wenn Sie Modelle mit verschiedenen festen Effekten, die nicht vollständig bestraft werden, mit REML vergleichen.

Wieder ist die Folgerung klar; Das Modell mit Glättungen x1und x2hat eine wesentlich bessere Passform als jedes der beiden anderen Modelle.

Marra, G. & Wood, SN Praktische Variablenauswahl für verallgemeinerte additive Modelle. Comput. Stat. Daten Anal. 55, 2372–2387 (2011).

Wood, SN, Pya, N. & Säfken, B. Glättungsparameter und Modellauswahl für allgemeine glatte Modelle. Marmelade. Stat. Assoc. 111, 1548–1563 (2016).

Möglicherweise möchten Sie die beiden Modelle mit testen lrest.

lrtest(two_term_model, two_smooth_model)

Model 1: y ~ x0 + s(x1) + x2
Model 2: y ~ x0 + s(x1) + s(x2)
      #Df  LogLik    Df  Chisq Pr(>Chisq)    
1  8.1107 -995.22                            
2 15.0658 -848.95 6.955 292.55  < 2.2e-16 ***

Während das Hinzufügen einer glatten Funktion zu beiden Begriffen das Modell tatsächlich kompliziert, ist die Verbesserung der Log-Wahrscheinlichkeit signifikant. Dies sollte nicht überraschen, da die Daten von einem GAM-Simulator generiert wurden.

Möglicherweise möchten Sie auch die zusammenfassende Statistik ausdrucken:

Link function: identity 

Formula:
y ~ x0 + s(x1) + x2

Parametric coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  11.6234     0.3950  29.429  < 2e-16 ***
x02           2.1147     0.4180   5.059 6.48e-07 ***
x03           4.3813     0.4172  10.501  < 2e-16 ***
x04           6.2644     0.4173  15.010  < 2e-16 ***
x2           -6.4110     0.5212 -12.300  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Approximate significance of smooth terms:
        edf Ref.df     F p-value    
s(x1) 2.111  2.626 64.92  <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

R-sq.(adj) =  0.583   Deviance explained = 58.9%
GCV = 8.7944  Scale est. = 8.6381    n = 400

und

Family: gaussian 
Link function: identity 

Formula:
y ~ x0 + s(x1) + s(x2)

Parametric coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   8.3328     0.2074  40.185  < 2e-16 ***
x02           2.1057     0.2955   7.125 5.15e-12 ***
x03           4.3715     0.2934  14.901  < 2e-16 ***
x04           6.1197     0.2935  20.853  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Approximate significance of smooth terms:
        edf Ref.df     F p-value    
s(x1) 2.691  3.343 95.00  <2e-16 ***
s(x2) 7.375  8.356 85.07  <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

R-sq.(adj) =  0.796   Deviance explained = 80.2%
GCV = 4.3862  Scale est. = 4.232     n = 400

Beachten Sie den Unterschied in der Abweichung erklärt (es ist riesig). Das kompliziertere Modell hat auch ein besseres R-Quadrat (adj). Der zweite Glättungsterm ist von großer Bedeutung und passt gut zu den Daten.