ANOVA para comparar modelos

Question

Estoy mirando este sitio para un taller sobre GAM en R: http://qcbs.ca/wiki/r_workshop8

Al final de la sección 2. Multiple smooth terms muestran un ejemplo, en el que utilizan anova para comparar tres modelos diferentes y determinar el que mejor se ajusta. El resultado es

  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

Basándose en esto, concluyen que el modelo 3 es el mejor. Mi pregunta es ¿cómo lo ven?

Mi entendimiento actual es: El Pr(>Chi) -es pequeño tanto para el modelo 2 como para el 3, por lo que estos son mejores que el modelo 1. Sin embargo, ¿qué otra variable están utilizando para determinar que el 3 es mejor que el 2?

Answer 1

La salida de anova() es una serie de pruebas de razón de verosimilitud. Las líneas en la salida son:

1. La primera línea de la salida corresponde al modelo más sencillo, con sólo un liso de x1 (Estoy ignorando el factor x0 ya que no está en consideración en su ejemplo) - esto no está probado contra nada más sencillo de ahí que las últimas entradas de la columna estén vacías.

2. La segunda línea es una prueba de razón de verosimilitud entre el modelo de la línea 1 y el modelo de la línea 2. A costa de 0.97695 grados de libertad adicionales, la desviación residual disminuye en 1180.2 . Esta reducción de la desviación (o a la inversa, el aumento de la desviación explicada), a costa de <1 grado de libertad, es muy poco probable si el verdadero efecto de x2 eran 0.

   Por qué 0.97695 ¿aumentan los grados de libertad? Pues bien, la función lineal de x2 añadiría 1 df al modelo pero el suavizador para x1 se penalizará un poco más que antes y, por lo tanto, utilizará un poco menos de grados de libertad efectivos, de ahí el cambio de <1 en los grados de libertad totales.

3. La tercera línea es exactamente la misma que he descrito anteriormente pero para una comparación entre el modelo de la segunda línea y el modelo de la tercera línea: es decir, la tercera línea está evaluando la mejora al pasar de modelar x2 como un término lineal para modelar x2 como una función suave. De nuevo, esta mejora en el ajuste del modelo (el cambio en la desviación es ahora 2211.8 a costa de 7.37288 más grados de libertad) es poco probable si los parámetros adicionales asociados a s(x2) eran todos iguales a 0.

En resumen, la línea 2 dice que el Modelo 2 se ajusta mejor que el Modelo 1, por lo que una función lineal de x2 es mejor que ningún efecto de x1 . Pero la línea 3 dice que el Modelo 3 se ajusta a los datos mejor que el Modelo 2, por lo que una función suave de x2 es preferible a una función lineal de x2 . Se trata de un análisis secuencial de modelos, no de una serie de comparaciones con el modelo más sencillo.

Sin embargo

Lo que muestran no es la mejor manera de hacer esto - la teoría reciente sugeriría que la salida de summary(m3) tendría las propiedades de cobertura más "correctas". Además, para seleccionar entre los modelos, probablemente habría que utilizar select = TRUE al ajustar el modelo completo (el que tiene dos alisados), lo que permitiría reducir los términos que incluiría el modelo con lineal x2 o incluso ningún efecto de esta variable. Tampoco están ajustando utilizando REML o la selección de suavidad ML, que muchos de nosotros mgcv los usuarios considerarían la opción por defecto (aunque no lo sea en realidad en gam() )

Lo que yo haría es:

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

La última línea produce lo siguiente:

> 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

Podemos ver que ambos términos suaves son significativamente diferentes de las funciones nulas.

Qué select = TRUE está haciendo es poner una penalización extra en el espacio nulo de la penalización (esta es la parte de la spline que es perfectamente suave). Si no se tiene esto, la selección de suavidad sólo puede penalizar una función suave a una función lineal (porque la penalización que está haciendo la selección de suavidad sólo funciona en las partes no suaves (las onduladas) de la base). Para llevar a cabo la selección tenemos que ser capaces de penalizar el espacio nulo (las partes suaves de la base) también.

select = TRUE consigue esto mediante el uso de una segunda penalización añadida a todos los términos suaves del modelo (Marra y Wood, 2011). Esto actúa como un tipo de contracción, tirando de todos los términos suaves un poco hacia 0, pero tirará de los términos superfluos hacia 0 mucho más rápidamente, por lo que los seleccionará fuera del modelo si no tienen ningún poder explicativo. Pagamos un precio por esto cuando evaluamos la importancia de los términos suaves; observe el Ref.df de la columna anterior (el 9 proviene del valor por defecto de k = 10 (que para las splines de placa fina con restricciones de centrado significa 9 funciones base), en lugar de pagar algo así como 2,5 y 7,7 grados de libertad por las splines, estamos pagando 9 grados de libertad por cada una. Esto refleja el hecho de que hemos hecho la selección, que no estábamos seguros de qué términos deberían estar en el modelo.

Nota: es importante que no utilice anova(m1, m2, m3) en los modelos que utilizan select = TRUE . Como se señala en ?mgcv:::anova.gam la aproximación utilizada puede ser muy mala para los alisados con penalizaciones en sus espacios nulos.

En los comentarios, @BillyJean mencionó el uso de AIC para la selección. Un trabajo reciente de Simon Wood y sus colegas (Wood et al, 2016) derivó un AIC que tiene en cuenta la incertidumbre adicional debida a que hemos estimado los parámetros de suavidad en el modelo. Este AIC funciona razonablemente bien, pero hay alguna discusión en cuanto al comportamiento de su derivación del AIC cuando los suavizados del IIRC son cercanos a funciones lineales. De todos modos, el AIC nos daría

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

Tenga en cuenta que volví a ajustar todo esto con la selección de suavidad ML, ya que no estoy seguro de lo que hace el AIC cuando select = TRUE y hay que tener cuidado al comparar modelos con diferentes efectos fijos, que no están totalmente penalizados, utilizando REML.

De nuevo, la inferencia es clara; el modelo con suavidad de x1 y x2 tiene un ajuste sustancialmente mejor que cualquiera de los otros dos modelos.

---

Marra, G. & Wood, S. N. Practical variable selection for generalized additive models. Comput. Stat. Data Anal. 55, 2372-2387 (2011).

Wood, S. N., Pya, N. & Säfken, B. Smoothing Parameter and Model Selection for General Smooth Models. J. Am. Stat. Assoc. 111, 1548-1563 (2016).

Answer 2

Puede probar los dos modelos con 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 ***

Aunque la adición de una función suave a ambos términos complica el modelo, la mejora en la probabilidad logarítmica es significativa. Esto no debería sorprender porque los datos fueron generados por un simulador GAM.

También puede imprimir las estadísticas resumidas:

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

y

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

Obsérvese la diferencia en la desviación explicada (es enorme). El modelo más complicado también tiene mejor R-sq.(adj). El segundo término de suavización es muy significativo y se ajusta muy bien a los datos.