5. 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,1) 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-Leibler2):
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.
Incidência de Cipó | |||
---|---|---|---|
Ano | Baixa | Média | Alta |
1997 | 1146 | 797 | 1081 |
2000 | 1044 | 711 | 934 |
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:
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)3) 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
Outro modelo, não contemplado neste teste é
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 AICc4):
> 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
Modelo | Log-Verossimilhança | k | $$AIC_c$$ | $$\DeltaAIC_c$$ |
---|---|---|---|---|
1 | 19,9 | 2 | 45,1 | 8,4 |
2 | 14,1 | 3 | 36,7 | 0,0 |
3 | 12,8 | 4 | 38,6 | 1,9 |
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:
Modelo | Log-Verossimilhança | k | $$AIC_c$$ | $$\DeltaAIC_c$$ |
---|---|---|---|---|
1 | 16,0 | 2 | 37,2 | 1,2 |
2 | 13,7 | 3 | 36,0 | 0,0 |
3 | 12,7 | 4 | 38,4 | 3,6 |
Indique os tutoriais e partes da leitura básica que merecem mais atenção na discussão aqui.