ANOVA para comparar modelos

Estou procurando neste site um workshop sobre o GAM em R: http://qcbs.ca/wiki/r_workshop8

No final da seção, 2. Multiple smooth termseles mostram um exemplo, onde eles usam anovapara comparar três modelos diferentes para determinar o melhor modelo de ajuste. A saída é

  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

Com base nisso, eles concluem que o modelo 3 é o melhor. Minha pergunta é como eles vêem isso?

Meu entendimento atual é: O Pr(>Chi)valor-é pequeno para os modelos 2 e 3; portanto, eles são melhores que o modelo 1. No entanto, que outra variável eles estão usando para determinar que 3 é melhor que 2?

A saída de anova()é uma série de testes de razão de verossimilhança. As linhas na saída são:

1. A primeira linha na saída corresponde ao modelo mais simples, com apenas uma suavização de x1(estou ignorando o fator, x0pois não está em consideração no seu exemplo) - isso não é testado contra nada mais simples, portanto, as últimas entradas da coluna são esvaziar.

2. A segunda linha é um teste de razão de verossimilhança entre o modelo na linha 1 e o modelo na linha 2. Ao custo de 0.97695graus extras de liberdade, o desvio residual é diminuído em 1180.2. Essa redução no desvio (ou, inversamente, o aumento do desvio explicado), ao custo de <1 grau de liberdade, é altamente improvável se o verdadeiro efeito de x2for 0.

   Por que os 0.97695graus de liberdade aumentam? Bem, a função linear de x2adicionaria 1 df ao modelo, mas o mais suave x1será penalizado um pouco mais do que antes e, portanto, usará um pouco menos graus efetivos de liberdade, daí a mudança <1 nos graus gerais de liberdade.

3. A terceira linha é exatamente a mesma que eu descrevi acima, mas para uma comparação entre o modelo na segunda linha e o modelo na terceira linha: ou seja, a terceira linha está avaliando a melhoria na mudança de modelagem x2como um termo linear para modelagem x2como uma função suave. Novamente, essa melhoria no ajuste do modelo (a mudança no desvio está agora 2211.8ao custo de 7.37288mais graus de liberdade) é improvável se os parâmetros extras associados a s(x2)todos forem iguais a 0.

Em resumo, a linha 2 diz que o Modelo 2 se encaixa melhor que o Modelo 1, portanto, uma função linear de x2é melhor do que nenhum efeito x1. Mas a linha 3 diz que o Modelo 3 ajusta os dados melhor que o Modelo 2, portanto, uma função suave de x2é preferida a uma função linear de x2. Esta é uma análise seqüencial de modelos, não uma série de comparações com o modelo mais simples.

Contudo…

O que eles estão mostrando não é a melhor maneira de fazer isso - a teoria recente sugere que a saída summary(m3)teria as propriedades de cobertura mais "corretas". Além disso, para selecionar entre os modelos, provavelmente, deve-se usar select = TRUEao ajustar o modelo completo (aquele com dois suaves), o que permitiria o encolhimento dos termos que incluiriam o modelo com efeito linear x2ou mesmo sem efeito dessa variável. Eles também não são adequados usando a seleção de suavidade REML ou ML, que muitos de nós usuários do mgcv considerariam a opção padrão (mesmo que não seja o padrão real gam())

O que eu faria é:

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

A linha final produz o seguinte:

> 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 os dois termos suaves são significativamente diferentes das funções nulas.

O que select = TRUEestá fazendo é aplicar uma penalidade extra no espaço nulo da penalidade (essa é a parte do spline que é perfeitamente suave). Se você não tiver isso, a seleção de suavidade poderá penalizar apenas uma volta suave para uma função linear (porque a penalidade que está fazendo a seleção de suavidade funciona apenas nas partes não suaves (as onduladas) da base). Para realizar a seleção, precisamos ser capazes de penalizar o espaço nulo (as partes suaves da base).

select = TRUEconsegue isso através do uso de uma segunda penalidade adicionada a todos os termos suaves do modelo (Marra e Wood, 2011). Isso age como uma espécie de retração, puxando todos os termos suaves de alguma maneira para 0, mas puxa termos supérfluos para 0 muito mais rapidamente, portanto, selecioná-los fora do modelo, se não tiverem nenhum poder explicativo. Pagamos um preço por isso ao avaliar a importância dos smooths; observe a Ref.dfcoluna acima (o 9 vem do valor padrão de k = 10, que para splines de chapas finas com restrições de centralização significa 9 funções básicas), em vez de pagar algo como 2,5 e 7,7 graus de liberdade pelos splines, pagamos 9 graus de liberdade cada. Isso reflete o fato de termos feito a seleção, de não termos certeza de quais termos deveriam estar no modelo.

Nota: é importante que você não use anova(m1, m2, m3)chamadas de tipo em modelos usando select = TRUE. Como observado em ?mgcv:::anova.gam, a aproximação usada pode ser muito ruim para suavizações com penalidades em seus espaços nulos.

Nos comentários, @BillyJean mencionou o uso da AIC para a seleção. Trabalhos recentes de Simon Wood e colegas (Wood et al, 2016) derivaram uma AIC que explica a incerteza extra devido a termos estimado os parâmetros de suavidade no modelo. Essa AIC funciona razoavelmente bem, mas há alguma discussão sobre o comportamento de sua derivação da AIC quando as suavizações do IIRC estão próximas das funções lineares. De qualquer forma, a AIC nos daria:

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

Observe que reajustei tudo isso com a seleção de suavidade de ML, pois não tenho certeza do que o AIC faz quando select = TRUEe você deve ter cuidado ao comparar modelos com diferentes efeitos fixos, que não são totalmente penalizados, usando REML.

Novamente a inferência é clara; o modelo com suavidades x1e x2ajuste substancialmente melhor do que qualquer um dos outros dois modelos.

Marra, G. & Wood, SN Seleção prática de variáveis ​​para modelos de aditivos generalizados. Comput. Estado. Data Anal. 55, 2372- 2387 (2011).

Wood, SN, Pya, N. & Säfken, B. Parâmetro de suavização e seleção de modelos para modelos suaves gerais. Geléia. Estado. Assoc. 111, 1548-1563 (2016).

Você pode querer testar os dois modelos com 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 ***

Embora a adição de uma função suave a ambos os termos complique o modelo, a melhoria na probabilidade de log é significativa. Isso não deveria surpreender, porque os dados foram gerados por um simulador GAM.

Você também pode imprimir as estatí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

e

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

Observe a diferença de desvio explicada (é enorme). O modelo mais complicado também possui melhor R-sq. (Adj). O segundo termo de suavização é altamente significativo e se ajusta bem aos dados.