ANOVA para comparar modelos

9

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?

BillyJean
fonte
3
Dica: quais variáveis ​​aparecem no Modelo 3 que não aparecem no Modelo 2? Para responder a isso, você precisa saber o que "s (x2)" significa - e isso depende de como a função sé definida. (Presumo que seja algum tipo de spline, mas reluto em presumir mais do que isso.) Podemos perceber pela saída que é bastante complicado - a mudança x2para s(x2)adiciona graus de liberdade - mas é tudo o que podemos determinar sobre isso a partir desta saída. 7.37288
whuber
11
Doing AIC(model1, model2, model3)revela que o modelo 3 tem um valor mais baixo AIC. Esta poderia ser mais uma prova de que é o modelo ideal entre os três
BillyJean
3
A enorme diferença de desvio entre os modelos (2) e (3) é convincente.
whuber

Respostas:

12

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).

Gavin Simpson
fonte
11
+1 para a resposta detalhada com código e referência. Vou ler e aprender esta resposta em detalhes mais tarde. Uma pergunta, você disse que é o teste da razão de verossimilhança. mas em ?anova.lmque não há essa possibilidade, pode ser F chisq ou CP
Haitao Du
11
@ hxd1011 este mgcv:::anova.gamnão está usando o método para lmmodelos. Estes são testes de análise de desvio, mas é a mesma coisa que os índices de probabilidade.
Gavin Simpson
11
Obrigado. Você poderia responder minha pergunta aqui com algum resumo de alto nível? Ou sua resposta está aqui já está coberta.
Haitao Du
11
@ hxd1011 Você precisaria ser mais específico. Que tipos de modelos? Há muitas suposições por trás, anova()mas quais dependem de qual é o modelo. Frequentemente, para modelos não gaussianos, eles estão realizando testes de razão de verossimilhança ou testes semelhantes, mas as suposições variam; eles variam mesmo para GLMs e GAMs. anova()é uma função de conveniência, mas não está executando o ANOVA sensu no modelo linear geral para nada além de um modelo linear geral (ajustado via lm()say).
Gavin Simpson
2
@ DeltaIV Desculpe, isso foi uma redação ruim da minha parte. O que select = TRUEfaz é penalizar totalmente todos os termos suaves, que o AFAIU faz comparações usando REML OK. Não observei os detalhes da nova AIC para o GAMS para ver o que faria com as multas extras adicionadas ao usar select = TRUE. Portanto, se deixamos de select = TRUElado o lado seguro, temos o problema de que a REML não é uma probabilidade verdadeira e não é usada nas comparações da AIC porque depende dos efeitos fixos no modelo. Contabilizar as duas preocupações significa que eu usei method = "ML"(não method = "REML") ao montar.
Gavin Simpson
3

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.

Olá Mundo
fonte
2
Isso não produz apenas outro exemplo como esse na questão? Você poderia indicar mais explicitamente como ele responde à pergunta?
whuber