• Entropia e informação
• Discrepância por aproximação e por estimativa
• Distância de Kullback-Leibler
• Critério de informação de Akaike
• Delta-AIC, peso de evidência, tabela de seleção de modelos

O Critério de Informação de Akaike é um estimador da distância relativa esperada entre dois modelos probabilísticos.

Se o modelo verdadeiro $$f(x)$$ é conhecido, a discrepância de um outro modelo $$g(x)$$ em relação a ele é dada pela distância de Kullback-Leibler:

$$I(f,g)=\int f(x)\ ln \ frac{f(x)}{g(x | \theta)}$$

onde $$g(x | \theta)$$ representa o modelo $$g$$ com os valores de seus parâmetros $$\theta$$ que dão a melhor aproximação para o modelo verdadeiro $$f$$.

Para distribuições discretas esta distância é:

$$I(f,g)=\sum f(x_i)\ ln \ frac{f(x_i)}{g(x_i | \theta)}$$

Nos dois casos esta medida expressa a discrepância em relação a um modelo fixo de referência. Por isso podemos deduzir uma medida de distância relativa, que para distribuições discretas é

$$I(f,g)=\C - sum f(x_i)\ ln \ g(x_i | \theta)\ = \ C - E_x\ [ln \ g(x_i| \theta)]$$

Onde $$C \ = \ sum f(x_i) \ ln \ f(x_i) \ = \ E_x\ [ln\ f(x_i)]$$

E para distribuições contínuas:

$$I(f,g)=\C - int f(x_i)\ ln \ g(x_i | \theta) \ = \ C - E_x\ [ln \ g(x| \theta)]$$

Onde $$C \ = \ int f(x_i) \ ln \ f(x_i) \ = \ E_x\ [ln \ f(x_i)]$$

Como o modelo de referência é o mesmo, podemos usar $$E_x\ [ln \ g(x| \theta)]$$ como uma medida de distância relativa entre modelos.

Vamos imaginar que nosso modelo real é uma variável binomial com parâmetros $$N=10$$ e $$p=0,05$$, o que resulta em uma média de $$Np=10 \times 0,05 = 0,5$$. A Poisson que melhor aproxima este modelo é a que tem a mesma média, portanto com parâmetro $$\lambda=0,5$$. A sequência de comandos abaixo faz os cálculos da distância de Kullback-Leibler entre este modelo Poisson e o verdadeiro:

x <- 0:10 ## Todos os valores de x
f.x <- dbinom(x,size=10,p=0.05) ## valores de f(x) para cada x
g.x <- dpois(x,lambda=0.5) ## valores de g(x) para cada x
log.g.x <- log(g.x) ## ln (g(x))

Uma tabela com estes valores pode ser obtida com:

> round(data.frame(x,f.x,g.x,log.g.x,E.log.x=-f.x*log.g.x),3)
    x   f.x   g.x log.g.x E.log.x
1   0 0.599 0.607  -0.500   0.299
2   1 0.315 0.303  -1.193   0.376
3   2 0.075 0.076  -2.579   0.193
4   3 0.010 0.013  -4.371   0.046
5   4 0.001 0.002  -6.451   0.006
6   5 0.000 0.000  -8.753   0.001
7   6 0.000 0.000 -11.238   0.000
8   7 0.000 0.000 -13.877   0.000
9   8 0.000 0.000 -16.650   0.000
10  9 0.000 0.000 -19.540   0.000
11 10 0.000 0.000 -22.536   0.000

E o obtemos valor de distância relativa com o negativo da soma da última coluna:

> -sum(f.x*log.g.x)
[1] 0.9204515

No mundo real nem o modelo verdadeiro nem os parâmetros de qualquer modelo que melhor aproximam o verdadeiro são conhecidos. Tudo o que temos são modelos simples definidos a partir de nosso conhecimento do sistema, e estimativas dos parâmetros destes modelos, obtidos a partir de uma amostra.

Por exemplo, tomando uma amostra de tamanho 5 da binomial

> set.seed(42)
> obs <- rbinom(5,size=10,prob=0.05)
> obs
[1] 2 2 0 1 1
> mean(obs)
[1] 1.2

A melhor estimativa do parâmetro de um modelo Poisson para aproximá-la será a média amostral, $$\bar x = 1.2$$.

Os comandos abaixo produzem uma figura com modelo verdadeiro binomial, o modelo Poisson que melhor o aproxima, e um modelo Poisson com parâmetro $$\lambda=1.2$$, estimado da amostra tomada acima.

plot(0:10,dbinom(0:10,size=10,p=0.05),,type="h",
     lwd=3,xlab="N de sucessos", ylab="p", ylim=c(0,0.7),
     main="X ~ Bin(N=10,p=0,05), Aproximação Poisson")
points(0:10,dpois(0:10,lambda=0.5),col="red", pch=19,
       cex=0.75, type="b")
points(0:10,dpois(0:10,lambda=mean(obs)),col="blue",
       pch=19, cex=0.75, type="b")
legend("topright",
       legend=c(bquote(lambda==.(mean(obs))),
         expression(paste(lambda," = 0.5"))),pch=19,
       col=c("blue","red"),cex=1.5)

Você deve obter uma figura como esta:

Nestes casos, as estimativas dos parâmetros obtidas da amostra, $$\hat \theta (y)$$ são variáveis aleatórias. Portanto, a distância relativa do modelo, que já é uma esperança, agora depende de outra esperança, que é o valor médio de $$g(x|\theta(y))$$, definido pela distribuição de valores possíveis das estimativas dos parâmetros. Com isto, temos uma esperança dupla, que expressa a distância relativa média, condicionada a um certo conjunto de dados:

$$T \ = \E_y\ E_x \ [ln \ g(x|\hat \theta(y))]$$

o Critério de Informação de Akaike é usado para estimar esta dupla esperança a partir da log-verossimilhança negativa:

$$AIC \ = \ -2\ ln L( \theta | y) + 2 K$$

Onde $$L( \theta | y)$$ é a log-verossimilhança, e $$K$$ é o número de parâmetros do modelo.

Agora vamos simular várias amostras da mesma binomial, calcular o AIC e a distância relativa esperada de Kullback-Leibler para os modelos Poisson com parâmetros estimados de cada uma destas amostras.

Primeiro definimos uma função para o cálculo da distância estimada,¹⁾ dada uma amostra y e os parâmetros do modelo verdadeiro (parâmetros da binomial, argumentos N e p):

KLD.bp <- function(y,N,p){
  x <- 0:N
  log.gx <- -dpois(x,lambda=mean(y),log=T)
  length(y)*sum(dbinom(x,size=N,prob=p)*log.gx)
}

E uma função para o cálculo do AIC de uma amostra y para a Poisson com parâmetro $$\lambda$$ estimado pela média amostral:

AIC.pois <- function(x){
  LL <- sum(dpois(x,lambda=mean(x),log=T))
  -2*LL + 2
}

Agora criamos 1000 amostras de tamanho 100, tomadas da binomial com parâmetros $$N=10$$ e $$p=0.05$$. Criamos um objeto matriz para guardar estes valores, sendo cada linha uma amostra:

sim <- matrix(rbinom(100000,size=10,p=0.05),nrow=1000)

Aplicando as funções criadas acima a cada linha da matriz teremos os valores de AIC e de Distância relativa estimada para cada amostra:

AIC.p <- apply(sim,1,AIC.pois)
KLD.p <- apply(sim,1,KLD.bp,N=10,p=0.05)

Faça um gráfico de densidade empírica dos valores de AIC com a função density, e acrescente uma linha vertical marcando o dobro da média dos valores da distância de Kullback-Leibler²⁾:

plot(density(AIC.p))
abline(v=2*mean(KLD.p),lty=2)

Qual sua avaliação?

Aqui retomamos o exercício da unidade de verossimilhança sobre cipós e sobrevivência de árvores.

A tabela abaixo apresenta o número de árvores em diferentes classes de incidência de cipó (incidência crescente). Os dados no ano de 1997 se referem a árvores com as quais se iniciou o estudo e os dados no ano de 2000 se referem às árvores que sobreviveram.

Do exercício anterior já temos os perfis de verossimilhança das taxas de sobrevivência de cada grupo, que indicam que há diferenças entre elas:

É realista aplicar uma mesma taxa de sobrevivência para todas as árvores independentemente da classe? Esta pergunta é respondida comparando-se dois modelos:

1. Todas as classes têm a mesma sobrevivência.
2. Cada classe tem uma taxa de sobrevivência diferente.
3. As classes de baixa e média infestação têm a mesma taxa de sobrevivência, que difere daquela da classe de alta infestação.
4. As classes de alta e média infestação têm a mesma taxa de sobrevivência, que difere daquela da classe de baixa infestação.

O número de sobreviventes pode ser descrito como uma variável binomial, cujo parâmetro $$N$$ é fixado pelo número de árvores em 1997. A estimativa de máxima verossimilhança da taxa de sobrevivência é a razão entre sobreviventes e este total.

Vamos definir uma função para calcular a log-verossimilhança negativa:

LL.bin <- function(x,p,N){
  -sum(dbinom(x,prob=p,size=N,log=T))
}

E a utilizamos para calcular a log-verossimilhança negativa de cada modelo. Para o modelo 1, ela será a soma dos logaritmos das probabilidades atribuídas pela binomial com parâmetro $$p=2689/3024 = 0.889$$ às contagens de sobreviventes em cada grupo:

 
N1 <- 1146
N2 <- 797
N3 <- 1081
n1 <- 1044
n2 <- 711
n3 <- 934
modelo1 <- LL.bin(p=(n1+n2+n3)/(N1+N2+N3),
                  N=N1+N2+N3,x=c(n1,n2,n3))

Para o modelo 2, calculamos as log-verossimilhanças para cada classe, usando como parâmetros o total de árvores e a proporção de sobreviventes em cada classe. Como as log-verossimilhanças são aditivas, podemos calcular cada uma separadamente e depois somá-las, para obter a log-verossimilhança do modelo completo:

modelo2 <- LL.bin(p=n1/N1,N=N1,x=n1)+
  LL.bin(p=n2/N2,N=N2,x=n2)+
  LL.bin(p=n3/N3,N=N3,x=n3)

Para o modelo 3, calculamos a log-verossimilhança das contagens nas duas primeiras classes para parâmetros $$N = 1146+797 = 1943$$ e $$p=(1044+711)/1943 = 0.903$$, e com parâmetros $$N=1031$$ e $$p=934/1081= 0.864$$ para a classe de infestação alta. A soma destas log-verossimilhanças será a log-verossimilhança total do modelo:

modelo3 <- LL.bin(p=(n1+n2)/(N1+N2),N=N1+N2,x=c(n1,n2))+
  LL.bin(p=n3/N3,N=N3,x=n3)

Usando a mesma lógica, calculamos a log-verossimlhança para o modelo 4:

modelo4 <- LL.bin(p=(n3+n2)/(N3+N2),N=N3+N2,x=c(n3,n2))+
  LL.bin(p=n1/N1,N=N1,x=n1)

O primeiro modelo tem um parâmetro estimado, que é uma taxa de sobrevivência para todas as classes. Seu AIC é:

> 2*modelo1 + 2
[1] 17811.13

E como o modelo 2 tem três parâmetros estimados, seu AIC será:

> 2*modelo2 + 6
[1] 25.23459

Por fim, calculamos os AICs dos modelos 3 e 4, que têm dois parâmetros estimados:

> 2*modelo3 + 4
[1] 5103.06
> 2*modelo4 + 4
[1] 4178.069

O modelo 2 é o mais plausível.

O primeiro exemplo de Teste t de Zar(1989)³⁾ são medidas de tempo de coagulação em pacientes que receberam dois tipos de medicamentos.

drugB <- c(8.8,8.4,7.9,8.7,9.1,9.6)
drugG <- c(9.9,9,11.1,9.6,8.7,10.4,9.5)

O experimento foi conduzido para avaliar se estas drogas afetam de maneira diferente a coagulação. As médias e desvios-padrão amostrais são:

media.B <- mean(drugB)
media.G <- mean(drugG)
sd.B <- sd(drugB)
sd.G <- sd(drugG)

Sendo a média e desvio-padrão das duas amostras combinadas (pooled mean e pooled standard deviation):

media.pool <- mean(c(drugB,drugG))
sd.pool <- (var(drugB)*(length(drugB)-1)+
            var(drugG)*(length(drugG)-1))/
            (length(drugB)+length(drugG)-2)

Estas quantidades são usadas para calcular a estatística $$t$$, cuja probabilidade de assumir o valor observado sob a hipótese nula é usado para se decidir se há diferenças nas médias das populações de onde vieram as amostras.

As duas hipóteses deste teste são os seguintes modelos

• Modelo 1: as medidas são realizações de uma só variável normal.
• Modelo 2: As medidas tomadas de cada grupo são realizações de duas variáveis normais, com parâmetro $$\mu$$ diferente, mas com o mesmo parâmetro $$\sigma$$.

Outro modelo, não contemplado neste teste é

• Modelo 3: As medidas tomadas de cada grupo são realizações de duas variáveis normais, com parâmetros $$\mu$$ e $$\sigma$$ diferentes.

Para comparar estes modelos, calculamos suas log-verossimilhanças negativas. Para cada modelo usaremos como valor dos parâmetros as estimativas de médias e desvios-padrão amostrais. Note que os estimadores de desvio-padrão usados no teste t não são MLEs.

LL1 <- -sum(dnorm(c(drugB,drugG),mean=media.pool,
            sd=sd.pool,log=T))
LL2 <- -sum(dnorm(drugB,mean=media.B,sd=sd.pool,log=T))+ 
       (-sum(dnorm(drugG,mean=media.G,sd=sd.pool,log=T)))
LL3 <- -sum(dnorm(drugB,mean=media.B,sd=sd.B,log=T))+ 
       (-sum(dnorm(drugG,mean=media.G,sd=sd.G,log=T)))

E com estes valores calculamos os valores de AIC. Como a amostra é pequena, é preciso usar o fator de correção

$$c = frac{n}{n-K-1}$$

Onde $$n$$ é o tamanho da amostra e $$k$$ o número de parâmetros do modelo. Criamos uma função para este cálculo:

fc <- function(k,n)(n/(n-k-1))

E calculamos os valores de AICc⁴⁾:

> 2*LL1+4*fc(2,13)
[1] 45.04798
> 2*LL2+6*fc(3,13)
[1] 36.70419
> 2*LL3+8*fc(4,13)
[1] 38.59355

O que podemos sintetizar em uma tabela de seleção de modelos

O modelo 3, não considerado pelo teste, é tão plausível quanto o modelo 2, que corresponde à hipótese alternativa, na abordagem de teste de significância.

No entanto, as estimativas de desvio-padrão nos modelos não são MLEs. Podemos fazer isto com os comandos:

sd.mle = function(x)  sqrt( (1/length(x)) * sum( (x - mean(x))^2 ) )
sd.m <- sd.mle(c(drugB,drugG))
sd.B.m <- sd.mle(drugB)
sd.G.m <- sd.mle(drugG)

LL1.mle <- -sum(dnorm(c(drugB,drugG),mean=media.pool,
            sd=sd.m,log=T))
LL2.mle <- -sum(dnorm(drugB,mean=media.B,sd=sd.m,log=T))+ 
       (-sum(dnorm(drugG,mean=media.G,sd=sd.m,log=T)))
LL3.mle <- -sum(dnorm(drugB,mean=media.B,sd=sd.B.m,log=T))+ 
       (-sum(dnorm(drugG,mean=media.G,sd=sd.G.m,log=T)))

O que resulta em valores diferentes, mas que levam à mesma conclusão:

Aula 5: Seleção de Modelos

• Página de David Anderson, um dos principais divulgadores da seleção de modelos com o AIC, autor com Ken Burnham da cartilha mais conhecida sobre o assunto (Burham & Anderson 2002). Vários links relevantes, como revisões críticas de testes de hipótese nula, respostas às críticas feitas ao uso do AIC (AIC Myths and misunderstandings), e pdfs de artigos selecionados sobre o assunto.

• Comentário de H. Akaike sobre a descoberta do AIC, publicado na coleção Citation Classics.

• Akaike without tears, um artigo de Ken Burnham & Elliot Sober. Veja vários outros artigos relacionados na página de E. Sober.

• Key Concepts in Model Selection, uma revisão crítica da seleção de modelos ( que não se esgota em AIC). Por Malcom Foster, filósofo que considera a seleção de modelos uma das operações básicas das ciências quantiativas.

Indique os tutoriais e partes da leitura básica que merecem mais atenção na discussão aqui.