^[[start|Modelos Estatísticos: Abordagem da Verossimilhança]]^ \\ 5. Seleção de Modelos \\ ====== Conceitos ====== * 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 ====== Tutoriais ====== ===== AIC como uma Estimativa ===== 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: {{:biometria:verossim:bin_pois.png|}} 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,((note que a distância é multiplicada pelo tamanho da amostra, veja Burnham & Anderson 2002 p. 59)) 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 [[http://finzi.psych.upenn.edu/R/library/stats/html/density.html|density]], e acrescente uma linha vertical marcando o dobro da média dos valores da distância de Kullback-Leibler((O AIC estima o dobro da distância KL)): plot(density(AIC.p)) abline(v=2*mean(KLD.p),lty=2) Qual sua avaliação? ===== Cipós e Sobrevivência de Árvores ===== Aqui retomamos o exercício da unidade de verossimilhança sobre [[biometria:verossim:02-veros#cipos_e_sobrevivencia_de_arvores|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: {{:biometria:verossim:ex_ver_2.png|Perfis de verossimilhança das taxas de sobrevivência para cada grupo}} É realista aplicar uma mesma taxa de sobrevivência para todas as árvores independentemente da classe? Esta pergunta é respondida comparando-se dois modelos: - Todas as classes têm a mesma sobrevivência. - Cada classe tem uma taxa de sobrevivência diferente. - 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. - 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. ===== Modelos no Teste t ===== O primeiro exemplo de Teste t de Zar(1989)((Zar,J. 1999. Biostatistical Analysis. 4th Ed. Prentice)) 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 [[http://en.wikipedia.org/wiki/Maximum_likelihood_estimation|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((usamos esta sigla para o AIC com correção para pequenas amostras)): > 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 | ====== Recursos para Estudo ====== ===== Resumo da Aula ===== {{:biometria:verossim:aula5_selecao_modelos.pdf|Aula 5: Seleção de Modelos}} ===== Na Internet ===== * Página de [[http://welcome.warnercnr.colostate.edu/~anderson/|David Anderson]], um dos principais divulgadores da seleção de modelos com o AIC, autor com [[http://welcome.warnercnr.colostate.edu/~kenb/|Ken Burnham]] da cartilha mais conhecida sobre o assunto [[http://www.springer.com/statistics/statistical+theory+and+methods/book/978-0-387-95364-9|(Burham & Anderson 2002)]]. Vários //links// relevantes, como [[http://welcome.warnercnr.colostate.edu/~anderson/null.html|revisões críticas de testes de hipótese nula]], respostas às críticas feitas ao uso do AIC [[http://welcome.warnercnr.colostate.edu/~anderson/PDF_files/AIC%20Myths%20and%20Misunderstandings.pdf|(AIC Myths and misunderstandings)]], e //pdfs// de [[http://welcome.warnercnr.colostate.edu/~anderson/sel_reprints.html|artigos selecionados]] sobre o assunto. * [[http://www.garfield.library.upenn.edu/classics1981/A1981MS54100001.pdf|Comentário de H. Akaike]] sobre a descoberta do AIC, publicado na coleção [[http://garfield.library.upenn.edu/classics.html|Citation Classics]]. * [[http://philosophy.wisc.edu/sober/akaike%20-%20b&s%203.PDF|Akaike without tears]], um artigo de Ken Burnham & Elliot Sober. Veja vários outros artigos relacionados na página de [[http://philosophy.wisc.edu/sober/index.html|E. Sober]]. * [[http://philosophy.wisc.edu/forster/papers/model.htm|Key Concepts in Model Selection]], uma revisão crítica da seleção de modelos ( que não se esgota em AIC). Por [[http://philosophy.wisc.edu/forster/|Malcom Foster]], filósofo que considera a seleção de modelos uma das operações básicas das ciências quantiativas. ====== Pesquisa ====== Indique os tutoriais e partes da leitura básica que merecem mais atenção na discussão [[http://ecologia.ib.usp.br/bie5782/doku.php?id=bie5782:verossim:pesquisa1#modelos_com_par%C3%A2metros_constantes_e_sele%C3%A7%C3%A3o_de_modelos|aqui]].