Usando o AIC e Maximum likelihood, um pequeno exemplo.

Ha muitas armadilhas para interpretar resultados estatísticos, ainda mais se tratando de valores p como ja visto em um post que escrevi aqui.
Mas valores p não são a única forma de se obter conclusões, ou pensar em relação a análise de dados.
Como a gente já viu em outro post anterior, podemos pensar no ajuste do modelo aos dados.

Isso tem raízes na teoria da informação eu acredito, assim a máxima verosimilhança e o critério de informação de Akaike podem ser usadas como medidas de de o quão bom é o ajuste do modelo aos dados.

Então podemos pensar ao invés do famoso “Existe diferença significativa”, qual o modelo que melhor explica nossos dados, dentro de uma gama de vários modelos que podemos produzir e então utilizar o mais adequado que veio de todos os modelos que produzimos.

Vou tentar ilustrar com um exemplo de como usar esse tipo de abordagem funciona, pelo menos para esses simples exemplos.

Então vamos supor uma uma regressão, onde temos o modelo “real” que é resposta=a+b*preditor. Simulando isso para um a=5 e b=2 temos:

Certo, vamos então ajustar 2 modelos.
Um modelo vai ser a regressão normal, resposta = a + b * preditor e o segundo modelo vai ser resposta = a , ou seja o preditor não prevê nada (o que sabemos que não é verdade).

Gerado os dados, podemos usar o comando logLik() para calcular um valor pro “quão bom esse modelo é”. Minha terminologia deve ser meio ruim, o primeiro lugar a verificar definições mais corretas é na documentação da própria função usando o comando ?logLik() no R.
Mas continuando. Faremos isso para os 2 modelos que ajustamos:

> logLik(modelo1) 'log Lik.' -46.87914 (df=3) > logLik(modelo2) 'log Lik.' -78.03321 (df=2)

Então temos que o valor do modelo que sabemos que é o correto (já que simulamos ele, mas a natureza nunca nos conta isso…) é de -46, enquanto o modelo mais simples (e errado) foi de -78, vemos então que quanto melhor o ajuste (e teoricamente melhor o modelo) temos um valor maior de logLik.

Agora o AIC é mais ou menos esse valor, só que ele da uma penalidade por parâmetro estimado.
Como assim parâmetro? Assim, o modelo 1 tem a+bx, então precisamos de um parâmetro, o “a”, dois parâmetros, o “b” e a quantidade de desvio do modelo (deviance), variância. Olhem la em cima que o df=3, ou seja o modelo tem esses 3 ingredientes.

Agora o modelo 2 tem menos ingredientes, ele só tem um “a” e desvio, a variação dos dados, olhem como o df=2. Ou seja, se ele explicar bem, mesmo que um pouco pior que o modelo 1, ele deveria ganhar uma “colher de chá”, já que ele tenta explicar os mesmo dados com menos parâmetros, de forma mais simples.

Esse é o aic, que podemos calcular usando o comando extractAIC(), para ver certinho a forma como é calculado, entre outras informações, use o comando ?extractAIC().

Mas utilizando ele no nosso exemplo temos:

> extractAIC(modelo1) [1] 2.00000 12.62196 > extractAIC(modelo2) [1] 1.0000 72.9301

Agora o modelo 1, com 2 parâmetros(a e b, variação todo modelo tem que ter senão não seria probabilístico então não contamos aqui) tem um valor baixo de 12 enquanto o modelo 2 (que sabemos que é o errado) tem 1 parâmetro e o valor de 72 de AIC, então valores menores devem ser bons. Somos inclinados a usar o modelo 1 que tem os 2 parâmetros, intercepto e inclinação.

Mas podemos ainda fazer um teste F, vendo a razão entre os logLik para saber se os modelos são diferentes.

> anova(modelo1,modelo2) Analysis of Variance Table Model 1: resposta ~ preditor Model 2: resposta ~ 1 Res.Df RSS Df Sum of Sq F Pr(>F) 1 28 39.99 2 29 319.11 -1 -279.12 195.44 3.745e-14 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

E vemos o que? Que existe diferença entre os modelos, então ficamos com aquele com mais parâmetros, o modelo 1, o mais complexo, não podemos abandonar ele pelo mais simples, ja que ele explica muita coisa que o modelo 2 mais simples não deu conta de explicar.

Podemos ainda fazer isso usando o comando step(), que vai pegar o modelo inicial que você fornecer, e se vc mandar ele simplificar com o argumento direction =”backward”, ele vai tirando um por um dos parâmetros, e ve se o modelo fica diferente, se continuar explicando bem os dados, ele descarta esse parâmetro, so deixar o trace=1 para ver o processo, como o modelo mais complexo aqui, o modelo 1 tem apenas 2 parâmetros, ele só faz isso uma vez, mas se você tiver muitos parâmetros ele vai fazendo isso pra você.
Então temos:

ps(trace=1 e nesse caso o direction=”backward” são os default da função)

> mod.final1<-step(modelo1) Start: AIC=12.62 resposta ~ preditor Df Sum of Sq RSS AIC 39.99 12.622 - preditor 1 279.12 319.11 72.930

ele tira um parâmetro, o modelo cresce muito o valor de AIC ai ele fala, não da pra tirar esse parâmetro senão o modelo vai ficar ruim, então ele deixa o modelo como está.

> mod.final1 Call: lm(formula = resposta ~ preditor) Coefficients: (Intercept) preditor 4.696 2.069

Até agora nos vimos um caso onde havia uma relação, ou seja havia uma inclinação, o b tinha um valor de 2 nos dados que simulamos, mas e o contrario será que funciona?
Então simulamos um caso onde não existe relação:

Podemos começar como no exemplo anterior, vamos chamar agora o modelo de regressão de modelo 3, e o modelo que é somente uma média de modelo 4, sabemos (porque simulamos) que o modelo 4 é o correto, e o 3 é errado, então olhando os valores de logLik.

> logLik(modelo3) 'log Lik.' -73.34481 (df=3) > logLik(modelo4) 'log Lik.' -73.34932 (df=2)

Agora eles são praticamente os mesmo, os dois modelos parecem explicar bem os dados, a diferença é que o modelo 3 intercepto, inclinação e variação e o modelo 4 usa somente intercepto e variação, ele explica o mesmo sendo mais “económico”.

Utilizando o AIC vemos:

> extractAIC(modelo3) [1] 2.00000 65.55331 > extractAIC(modelo4) [1] 1.00000 63.56234

Aha, agora o modelo 3 com 2 parâmetro tem o valor de AIC maior. Ou seja apesar de eles explicarem bem os dados, ele é penalizado por ter um parâmetro a mais e estar explicando a mesma coisa do modelo 4.

Se fizermos um teste F entre os modelos agora nos vemos que:

> anova(modelo3,modelo4) Analysis of Variance Table Model 1: resposta2 ~ preditor2 Model 2: resposta2 ~ 1 Res.Df RSS Df Sum of Sq F Pr(>F) 1 28 233.45 2 29 233.52 -1 -0.070269 0.0084 0.9275

Eles são iguais, logo a gente escolhe o mais simples, se ambos explicam igualmente os dados. Qual a lógica disso? Uma introdução rápida a ideia da navalha de Occlam, imagine que cada parâmetro vai explodir, mas ele também é a chave para a verdade, se você explodir não vai dar para ver a verdade, mas parâmetros podem valer muito, então a gente pesa, será que eu fico com 2 parâmetros e tenho 2 chances de explodir ou como eu sei que um não serve para muita coisa jogo ele fora, e fico com uma chance de ver a verdade e uma chance somente de explodir. Algo mais ou menos assim.

Podemos fazer esse processo com o comando step() também e ficamos com:

> mod.final2<-step(modelo3) Start: AIC=65.55 resposta2 ~ preditor2 Df Sum of Sq RSS AIC - preditor2 1 0.070269 233.52 63.562 233.45 65.553 Step: AIC=63.56 resposta2 ~ 1 > mod.final2 Call: lm(formula = resposta2 ~ 1) Coefficients: (Intercept) 19.69

Agora ele viu que o modelo fica igual, quase não aumenta o AIC então ele tira a inclinação do modelo e vai para o próximo passo, como já só tem um parâmetro no próximo passo ele para, e nosso modelo final selecionado pelo step() é somente intercepto, igual ao modelo 4.

Bem essa é somente uma ideia de como funciona o AIC, mas é preciso cuidado já que os modelos usados nos exemplos são extremamente simplistas. Com modelos mais complexos as coisas precisam de mais pensamento e reflexão, e nem step() ainda é melhor que um pouco de raciocínio, sempre temos que olhar se o resultado ou respostas que estamos tendo fazem sentido com a reliadade. Mas a partir daqui já da para ter uma ideia do que se trata AIC, logLik e compania, já que esses valores saem junto com modelos que ajustamos na maioria das funções no R e em outros programas estatísticos.
Outra coisa a notar, é que a conclusão sobre o que é importante aqui é a mesma usando valores p ou AIC/logLik, so dar um summary() nos modelos e veremos como no final a nossa conclusão sobre o que é importante é a mesma. Temos muitas possibilidades hoje de chegar a conclusões, mas é importante ter pelo menos uma ideia de todas para facilitar a leitura de artigos e saber as possibilidades existentes para resolver nossos próprios problemas:

set.seed(5)
#exemplo1
#simulação
preditor<-runif(30,5,10)
resposta<-rnorm(30,5+2*preditor)
 
#grafico
plot(resposta~preditor,pch=19,cex=1.6,xlab="Preditor",ylab="Resposta")
 
#ajuste dos modelos
modelo1<-lm(resposta~preditor)
modelo2<-lm(resposta~1)
 
#Qualidade do ajuste
logLik(modelo1)
logLik(modelo2)
 
extractAIC(modelo1)
extractAIC(modelo2)
 
anova(modelo1,modelo2)
 
#seleção automatica de modelo
mod.final1<-step(modelo1)
mod.final1
 
#compare com os valores p
summary(modelo1)
 
#exemplo2
#simulação, note que eu uso a media e desvio da resposta anterior
preditor2<-runif(30,5,10)
resposta2<-rnorm(30,mean(resposta),sd(resposta))
 
#grafico
plot(resposta2~preditor2,pch=19,cex=1.6,xlab="Preditor",ylab="Resposta")
 
#ajuste dos modelos
modelo3<-lm(resposta2~preditor2)
modelo4<-lm(resposta2~1)
 
#avaliando a qualidade do ajuste
logLik(modelo3)
logLik(modelo4)
 
extractAIC(modelo3)
extractAIC(modelo4)
 
anova(modelo3,modelo4)
 
#seleção automatica de modelo
mod.final2<-step(modelo3)
mod.final2
 
#compare com os valores p
summary(modelo3)

2 thoughts on “Usando o AIC e Maximum likelihood, um pequeno exemplo.

  1. Mas e aí, o que é melhor? Penalizar o número de parâmetros ou a likelihood? Porque pensando no que é o ajuste (likelihood) se você penalizar muito o número de parâmetros você se distancia do likelihood, mas ao mesmo tempo você não inclui o efeito de muitos parâmetros nos modelos. Com isso, qual seria o meio termo, ou como as pessoas tem feito com esse tipo de análise?

    • Ai acho que entra na parte de ponderar em relação a realidade, tem que cruzar o seu conhecimento de teoria de biologia com os resultados da modelagem que você esta fazendo, não existe solução ótima, tanto que no exemplo com AIC como com logLik a conclusão é a mesma.
      Então é realizar os modelos, de preferência até usando as duas abordagens e depois comparar os resultados.
      So que se tiver um milhão de parâmetros, usando o step tudo seja mais rápido para ter uma primeira ideia de como estão as coisas, ai depois ir confirmando dom logLik. Eu so tive uma experiência de usar AIC em publicação até hoje, e foi pra simplificar um modelinho ja bem simples hehe.

Leave a Reply

Your email address will not be published. Required fields are marked *