A distribuição Chi-quadrado

Ola, vamos dar uma olhada aqui na distribuição chi-quadrado e tentar entender melhor porque da para usar ela para selecionar modelos.

A gente já falou da distribuição f aqui e da distribuição t aqui, e agora adicionando a distribuição chi-quadrado ao repertório.

A distribuição chi-quadrado sempre tem um papel importante no teste de hipóteses sobre frequências.

Basicamente, se a gente tem um conjunto de n elementos \{Z_1,Z_2,\cdots,Z_n\} independentes e normalizados, a soma dos quadrados deles vai ter uma distribuição aleatória chi-quadrado com n graus de liberdades.

Então:

\chi^2=\sum\limits_{i=1}^n\{Z_i^2\}

Mas o que quer dizer normalizado? No post anterior, aqui, a gente realizou essa operação, é só diminuir o valor pela média e dividir pelo desvio.

Na distribuição normal a gente faz:

{x-\alpha} \over \sigma

Agora no teste de qui-quadrado a gente faz o famoso:

{Observado\ -\ Esperado} \over {Esperado}

Certo, agora vamos contextualizar isso.
Na graduação, muita gente deve ter participado em algum momento de um trabalho, pelo menos eu participei, de reproduzir moscas Drosophila melanogaster. Famosas nos trabalhinhos de genética.

Chega uma hora que a gente deixa elas cruzarem e depois ia ver se havia realmente uma segregação independente dos cromossomos. No caso do cromossomos ligados ao sexo, temos fêmeas XX e machos XY.

Ai lembrando as leis do nosso colega Gregor Mendel , podemos fazer a quadrado de punnet.

Fêmea X X M X XX XX a c Y XY XY h O

Ou seja esperamos 50% de machos e 50% de fêmeas, uma razão de 1:1.

Então se nasceu 100 moscas, dessas 50 machos e 50 fêmeas, isso é perfeitamente o que esperamos.
Mas se nascer 49 macho e 51 fêmeas, não é exatamente o esperado, mas esta bem próximo, agora se a gente ver 5 machos e 95 fêmeas, isso é bem longe do que esperamos. Mas como transformar esse próximo ou longe em um número?
Bang, chegamos ao chi-quadrado. Vamos fazer uma função de chi-quadrado no R para o nosso caso e fazer alguns teste.

chi2_moscas<-function(observado) {
    if(is.vector(observado) && length(observado)==2) {
        esperado<-rep(sum(observado)/2,2)
        estatistica<-((observado[1]-esperado[1])^2/esperado[1])+((observado[2]-esperado[2])^2/esperado[2])
        return(estatistica)
    } else {
        print("Entrada errada. Os dados devem ser entrados como vetor numericos de tamanho 2!")
    }
}

então essa função recebe como entrada um vetor de dois números, a nossa contagem de machos e fêmeas, e vai calcular a estatística chi-quadrado, dado que nossa hipótese para as mosquinhas realmente seja uma razão sexual de 1:1, como o senhor Mendel propôs, não estamos tirando da cartola esse valor de 50% para cada sexo.

Veja que na função, a gente adicionou um teste para ver que a gente está entrando com os dados corretamente, caso alguém mande os dados errados, é comum a gente ver verificações em funções, fizemos uma aqui.
Usando nossa função a gente vê.

> chi2_moscas(c(200)) [1] "Entrada errada. Os dados devem ser entrados como vetor numericos de tamanho 2!"

Opa, entramos com um dado que só tem uma contagem, isso não faz o menor sentido, então não calculamos nada, a gente só faz conta quando faz sentido.

Agora vamos supor que pegamos nossa garrafinha la, e contamos 221 moscas macho e 475 moscas fêmeas, nossa isso ta bem fora da casinha, ta muito diferente de 1:1, se a gente calcula o chi-quadrado a gente vê…

> chi2_moscas(c(221,475)) [1] 92.6954

Um número grandão. Mas é grande mesmo? Vamos pensar que em outra garrafa, a gente pegou 332 machos e 375 fêmeas, isso é bem mais razoável segundo a hipótese não? Ta mais com cara de 1:1, e quando a gente calcula o valor de chi-quadrado a gente vê…

> chi2_moscas(c(332,375)) [1] 2.615276

Um valor baixinho. Legal, então um valor grandão quer dizer uma discrepância muito grande da nossa hipótese, e um valor baixo quer dizer que estamos diante de algo bem próximo da hipótese.
Isso não é difícil de enxergar com a formula que estamos fazendo.

Estamos calculando na função chi-quadrado o seguinte

{\chi^2}={\sum\limits_{i=1}^2} { {(Observado\ -\ Esperado)^2} \over {Esperado} }

Veja que vai até 2 a somatória, porque temos machos e fêmeas aqui nesse problema. O quadrado é para não ter problemas com números negativos, mas o que acontece é que se o observado é próximo do número esperado, vai dar um número pequenininho ali no numerador, se for exatamente o mesmo número, da zero, e zero dividido por qualquer coisa da zero, se for um número pequeno, vai ser algo pequeno dividido por um número que vai dar um número pequeno, agora conforme essa diferença aumenta, o numerador aumenta, e cada vez o número fica maior.

Continuando com esse número total de 100 moscas, para ficar fácil fazer contas, podemos calcular o chi-quadrado para os casos de 49-51, 48-52 … até 1-99 e ver o que acontece.

01

Olha ai, conforme a discrepância em relação a nossa hipótese aumenta, o valor de chi-quadrado aumenta também, exponencialmente.

Mas num é qualquer discrepância que a gente vai para e assumir que o senhor Mendel tava errado. A gente tem que lembrar como as coisas estão acontecendo.

Cada mosca tem 50% de chance de ser macho o fêmea, independente da contagem atual, pode já ter nascido 99 moscas fêmeas, a próxima, se Mendel está certo, tem 50% de ser macho. Nos, seres humanos, somo péssimos geradores de números aleatórias por isso, se você contar 99 fêmeas, você vai primeiro pensar que algo acontece, as vezes sem levar em conta que esse resultado é perfeitamente possível dentro da teoria de Mendel, ele é raro, mas perfeitamente possível.

Agora a gente pode fazer um simulação no R e ver o quão raro é um caso desse, alias, a gente pode ao invés de guardar quantas moscas de qual sexo nasceram, simular os nascimentos e guardar somente o valor de chi-quadrado para o resultado.

Então é so a gente ir simulando os nascimento usando a função sample

> sample(c("Macho","Femea"),100,replace=T) [1] "Femea" "Macho" "Femea" "Femea" "Macho" "Femea" "Macho" "Macho" "Macho" "Femea" "Femea" [12] "Femea" "Femea" "Femea" "Femea" "Femea" "Femea" "Femea" "Macho" "Macho" "Femea" "Femea" [23] "Femea" "Femea" "Macho" "Macho" "Femea" "Macho" "Macho" "Macho" "Macho" "Macho" "Macho" [34] "Macho" "Macho" "Macho" "Macho" "Macho" "Macho" "Femea" "Macho" "Macho" "Femea" "Macho" [45] "Femea" "Macho" "Macho" "Femea" "Femea" "Macho" "Femea" "Macho" "Macho" "Macho" "Femea" [56] "Macho" "Femea" "Femea" "Femea" "Macho" "Femea" "Femea" "Femea" "Macho" "Macho" "Macho" [67] "Macho" "Femea" "Macho" "Macho" "Macho" "Femea" "Macho" "Femea" "Macho" "Macho" "Femea" [78] "Femea" "Femea" "Macho" "Macho" "Femea" "Macho" "Femea" "Macho" "Macho" "Femea" "Macho" [89] "Macho" "Femea" "Femea" "Macho" "Macho" "Femea" "Macho" "Femea" "Femea" "Femea" "Macho" [100] "Macho"

E ir guardando o resultado num vetor, já que é somente um número.

for(i in 1:10000) {
    experimento<-sample(c("Macho","Femea"),100,replace=T)
    macho<-sum(experimento=="Macho")
    femea<-sum(experimento=="Femea")
    estatistica[i]<-chi2_moscas(c(macho,femea))
}

Agora para facilitar a, a gente olha o resultado num histograma.

02

Humm, então existem poucas chances de valores extremos, mas valores próximos de zero são bem comuns.
Bem quando alguém descreve uma distribuição, o cara faz uma probability density function ou PDF, que nada mais é que uma função que descreve esse resultado. O R ja tem pronta um monte desses PDF, normalmente é d grudado no nome da distribuição, aqui é o dchisq.

Vamos dar uma olhada, se é parecido com o que encontramos.

03

Olha ai que legal, fica bem parecido com o que encontramos, ou seja antes da gente fazer a simulação, ja tinha como saber como seria o resultado com essa distribuição. Eu marquei o quantile de 95% também. O que esse risco azul ta mostrando, que 95% dos casos estão a esquerda desse número, e 5% a direita.

É fantástico pensar que isso funciona, podemos comparar esse numero que calculamos com o quantile com os valores de chi-quadrado que calculamos e ver se ele funciona de verdadinha.

> sum(estatistica<qchisq(0.95,1))/length(estatistica)
[1] 0.9456

E num é que funciona, 0.9456 são menor que ele, isso é praticamente 95%, muito boa predição.

Agora se o nosso valor que encontramos na garrifinha for muito extremo, podemos começar a suspeitar da teoria do Mendel, não que isso implique que ela é errada, mas podemos suspeitar, isso é o tal do p ser significativo, é a gente começar a suspeitar que algo estranho ta acontecendo, esse algo é permitido acontecer, mas incomum, muito incomum.

Vamos pensar em dois casos aqui, vamos pensar que contamos 43 machos e 57 femeas, se a gente calcular o valor de qui-quadrado

> experimento1<-c(43,57) > chi2_moscas(experimento1) [1] 1.96

Da 1.96, isso é bem baixinho, se a gente compara com os resultados da nossa simulação

> sum(chi2_moscas(experimento1)<estatistica)/length(estatistica)
[1] 0.1277

tem pelo menos 13% dos valores tão baixos quanto 1.96, não tem muito porque levantar suspeitas sobre a teoria de Mendel, e podemos comparar o procedimento que estamos fazendo com o teste de chi-quadrado já pronto do R, usando a função chisq.test

> chisq.test(experimento1,p=c(0.5,0.5)) Chi-squared test for given probabilities data: experimento1 X-squared = 1.96, df = 1, p-value = 0.1615

E a gente tem o mesmo valor da estatística calculada, 1.96 (estamos fazendo a mesma continha), e um valor de p que traz a mesma conclusão, 0.16.

Agora vamos pensar que em outra experimento vemos um resultado bem diferente, 21 machos e 79 fêmeas.
Aqui começa a ter cara de um desvio na razão sexual nervoso. A gente pode calcular o valor de chi-quadrado e vemos o seguinte.

> experimento2<-c(21,79) > chi2_moscas(experimento2) [1] 33.64

Nossa, 33.64, se a gente comprar com os valores que obtemos ao acaso, fazendo a simulação onde machos e fêmeas tem 50% de chance cada de nascer.

> sum(chi2_moscas(experimento2)<estatistica)/length(estatistica)
[1] 0

Nenhum, absolutamente nenhuma das 10 mil simulações foi tão extrema assim. Ou seja, algo deve estar acontecendo, ou a teoria do Mendel esta errada aqui, ou alguém esta pregando uma peça no laboratório, não sabemos o que, varias coisas podem influenciar esse resultado, mas sabemos que se for para confiar no azar para ver um resultado assim, esse resultado é muitoooo raro.

Podemos fazer o teste de chi-quadrado com a função do R aqui também.

> chisq.test(experimento2,p=c(0.5,0.5)) Chi-squared test for given probabilities data: experimento2 X-squared = 33.64, df = 1, p-value = 6.631e-09

E olha ai, um valor de p de 6.631e-09, lembrando que isso está em notação cientifica, esse e-09 significa que esse número é…

> format(6.631e-09,scientific=F) [1] "0.000000006631"

Muito legal né, então podemos usar essa distribuição para compara se algo está próximo da hipótese ou longe. Agora outro coisa, além de razão sexual, que tem uma distribuição que segue uma distribuição chi-quadrado são valores de verosimilhança, ou likelihood, e AIC também (AIC é loglikelihood corrigido meu, claro que tem que ser igual o likelihood quanto a distribuição).

E onde entra o chi-quadrado aqui?

Imagine que temos duas variáveis, um preditor e uma resposta, mas que nao tem relação nenhuma, os dois são amostras ao acaso.

> preditor<-rnorm(30) > resposta<-rnorm(30)

A gente pode tentar modelar eles tanto como variáveis que tem uma relação entre si, como um modelo onde a resposta simplesmente vem de uma distribuição normal, ou seja, assumir que ela simplesmente vem de uma distribuição normal é mais simples que assumir que ela tem relação com algo.

> modelo1<-lm(resposta~preditor) > modelo2<-lm(resposta~1)

A gente pode avaliar esses dois modelos e seus parâmetros.

> summary(modelo1) Call: lm(formula = resposta ~ preditor) Residuals: Min 1Q Median 3Q Max -1.7239 -0.6622 -0.1236 0.5255 2.1173 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 0.06622 0.17136 0.386 0.702 preditor 0.00845 0.19107 0.044 0.965 Residual standard error: 0.9379 on 28 degrees of freedom Multiple R-squared: 6.984e-05, Adjusted R-squared: -0.03564 F-statistic: 0.001956 on 1 and 28 DF, p-value: 0.965 > summary(modelo2) Call: lm(formula = resposta ~ 1) Residuals: Min 1Q Median 3Q Max -1.7339 -0.6600 -0.1162 0.5272 2.1214 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 0.06592 0.16826 0.392 0.698 Residual standard error: 0.9216 on 29 degrees of freedom

Mas e ai, qual é melhor? Ue o que tiver o melhor likelihood, o menor custo, aquele que acumula menos resíduo. Mas essa diferença pode ser ao acaso, a diferença de likelihood, e podemos testar isso usando a distribuição chi-quadrado.

> anova(modelo1,modelo2,test="Chisq") Analysis of Variance Table Model 1: resposta ~ preditor Model 2: resposta ~ 1 Res.Df RSS Df Sum of Sq Pr(>Chi) 1 28 24.628 2 29 24.630 -1 -0.0017202 0.9647

Olha ai, comparamos esses modelos e vemos que eles são iguais, tem a mesma capacidade de predição. Então por parcimônia, é melhor eu escolher o mais simples uai, porque ficar com algo mais complexo, se a complexidade nao ajuda em nada, essa é a idea da navalha de occan, e vamos seguir ela, mas tendo um teste para sustentar essa escolha, pode diminuir a discussão entre nos e nossos pares, afinal se todo mundo concorda com esses métodos, vamos todos chegar as mesmas conclusões.

Ficamos por aqui, agora com a distribuição chi-quadrado melhor compreendida, espero eu, os scripts sempre no repositório recologia e é isso, até o próximo post.

set.seed(123)
chi2_moscas<-function(observado) {
 
    if(is.vector(observado) && length(observado)==2) {
 
        esperado<-rep(sum(observado)/2,2)
        estatistica<-((observado[1]-esperado[1])^2/esperado[1])+((observado[2]-esperado[2])^2/esperado[2])
        return(estatistica)
 
    } else {
        print("Entrada errada. Os dados devem ser entrados como vetor numericos de tamanho 2!")
    }
 
}
 
 
chi2_moscas(c(200))
chi2_moscas(c(221,475))
chi2_moscas(c(332,375))
 
 
dados<-vector()
 
for(i in 1:50) {
    dados[i]<-chi2_moscas(c(50-i,50+i))
}
 
#figura 1
plot(dados,frame=F,type="p",pch=19,cex=0.5,xaxt="n",ylab="Estatisthica Chi-quadrado",xlab="Razão de Machos para Femeas")
axis(1,at=c(1,13,25,38,49),label=c("1:1","1:13","1:25","1:38","1:49"))
 
 
sample(c("Macho","Femea"),100,replace=T)
 
estatistica<-vector()
 
for(i in 1:10000) {
    experimento<-sample(c("Macho","Femea"),100,replace=T)
    macho<-sum(experimento=="Macho")
    femea<-sum(experimento=="Femea")
    estatistica[i]<-chi2_moscas(c(macho,femea))
}
 
#figura 2
hist(estatistica,prob=T,main="Estatística Chi-Quadrado",xlab="Estatística")
 
 
 
#figura 3
hist(estatistica,prob=T,main="Estatística Chi-Quadrado",xlab="Estatística")
curve(dchisq(x,1),0,15,add=T,lwd=3,lty=2,col="red")
abline(v=qchisq(0.95,1),col="blue",lty=3)
legend("topright",lwd=c(3,1),col=c("red","blue"),lty=c(2,3),legend=c("Distribuição Chi-Quadrado","Quantile de 95%"),
       bty="n")
 
 
sum(estatistica<qchisq(0.95,1))/length(estatistica)
 
 
experimento1<-c(43,57)
chi2_moscas(experimento1)
sum(chi2_moscas(experimento1)<estatistica)/length(estatistica)
chisq.test(experimento1,p=c(0.5,0.5))
 
 
experimento2<-c(21,79)
chi2_moscas(experimento2)
sum(chi2_moscas(experimento2)<estatistica)/length(estatistica)
chisq.test(experimento2,p=c(0.5,0.5))
 
format(6.631e-09,scientific=F)
 
preditor<-rnorm(30)
resposta<-rnorm(30)
 
modelo1<-lm(resposta~preditor)
modelo2<-lm(resposta~1)
 
summary(modelo1)
summary(modelo2)
 
anova(modelo1,modelo2,test="Chisq")

Construindo árvores filogenéticas do R com o pacote ape e PhyML

Para quem teve pouco contato e sempre quis saber como se fazem filogenias ou árvores filogenéticas, vamos fazer uma aqui somente pela curiosidade. Particularmente, vamos olhar algumas espécies de gatinhos, todo mundo gosta de gatinhos, para maiores detalhes consulte Johnson & O’Brien (1997), mas vamos olhar os genes 16s das mitocôndrias de gatinhos.

Mitocôndrias são particularmente mais simples de pensar, porque você só recebe mitocôndrias da sua mãe, além disso ele não tem duas cópias como os nossos cromossomos, é so comparar diretamente, não temos recombinação para complicar tudo. Mas vamos fazer apenas um exemplo trivial pela curiosidade mesmo.

Então ta, da onde diabos vem as árvores filogenéticas?

Ela vem da comparação de características, como características morfológicas, conjuntos de proteínas ou sequências de dna, e como conseguimos essas ultimas?

A gente pega algo que tenha dna, algum tecido, célula sei la, faz uma PCR para amplificar e ter muitas e muitas cópias do dna que nos interessa e mandamos para um empresa sequenciar, eu já ouvi falar de duas grandes tecnologias para sequenciar dna, a do Sanger, que inclusive ganhou um nobel por isso e a Pyrosequencing, outros laureados diga-se de passagem, realmente ambas são sacadas muito boas e que inundou o mundo com muito mais informação que conseguimos processar hoje em dia.

Mas então para pensar de forma mais simples, você tem um tubo de ensaio, faz pcr, ai posta no correio pra uma empresa e recebe um e-mail com um arquivo cheio de letrinhas, ai você guarda em lugares como o genbank essa informação, claro que acho que todo mundo segura essa informação enquanto pode, mas para publicar todo mundo tem que tornar público a informação, ainda bem, senão esse post não existiria porque não teriamos sequencias para testar como fazer árvores.

Ok, mas a partir daqui a gente começa. Como a gente olhou no artigo e descobriu os números pelos quais as sequências foram guardadas, a gente pode baixar da internet essa sequências e olhar. Particularmente o pacote ape deixa uma função pronta para fazer essa atividade.

Então fazemos uma lista das sequências que queremos baixar:

> lista [1] "AF006387" "AF006389" "AF006391" "AF006393" "AF006395" "AF006397" "AF006399" "AF006401" [9] "AF006403" "AF006405" "AF006407" "AF006409" "AF006411" "AF006413" "AF006415" "AF006417" [17] "AF006419" "AF006421" "AF006423" "AF006425" "AF006427" "AF006429" "AF006431" "AF006433" [25] "AF006435" "AF006437" "AF006439" "AF006441" "AF006443" "AF006445" "AF006447" "AF006449" [33] "AF006451" "AF006453" "AF006455" "AF006457" "AF006459"

E baixamos tudo usando a função read.GenBank(lista) do ape.

> felidseq16S 37 DNA sequences in binary format stored in a list. Mean sequence length: 373.892 Shortest sequence: 372 Longest sequence: 376 Labels: AF006387 AF006389 AF006391 AF006393 AF006395 AF006397 ... Base composition: a c g t 0.338 0.226 0.199 0.237

Olha ai, 37 sequências, mas elas tem tamanhos diferentes, a menor tem 372, a maior 376, os nomezinhos delas e olha ai a frequência de ocorrência das bases, fazer uma função para contar isso é o primeiro exercício do Rosalind.

Mas olhando então o tamanho das sequências vemos que:

> table(unlist(lapply(felidseq16S, length))) 372 373 374 375 376 3 9 15 9 1

Uma tabelinha com o tamanho delas.

Para desenhar uma árvore, precisamos que todas tenham o mesmo tamanho, ou seja, que estejam alinhadas, então temos que alinhar elas.

Se alguém tentar ler as sequências, vai ver que as bases estão representadas por números, mas ainda são os mesmo ACTG.
A função que le as sequências do GenBank tem uma opção para pegarmos tudo como letras, mas por algum motivo deve ser melhor ou mais eficiente pro pacote ape usar daquela forma, eu ainda num sei porque.

> read.GenBank(lista[1],as.character = T) $AF006387 [1] "t" "t" "t" "g" "t" "t" "c" "c" "c" "t" "a" "a" "a" "t" "a" "g" "g" "g" "a" "c" "t" "t" "g" [24] "t" "a" "t" "g" "a" "a" "c" "g" "g" "c" "c" "a" "c" "a" "c" "g" "a" "g" "g" "g" "c" "t" "t" [47] "t" "a" "c" "t" "g" "t" "c" "t" "c" "t" "t" "a" "c" "t" "t" "c" "c" "a" "a" "t" "c" "c" "g" [70] "t" "g" "a" "a" "a" "t" "t" "g" "a" "c" "c" "t" "t" "c" "c" "c" "g" "t" "g" "a" "a" "g" "a" [93] "g" "g" "c" "g" "g" "g" "a" "a" "t" "a" "t" "a" "a" "t" "a" "a" "t" "a" "a" "g" "a" "c" "g" [116] "a" "g" "a" "a" "g" "a" "c" "c" "c" "t" "a" "t" "g" "g" "a" "g" "c" "t" "t" "t" "a" "a" "t" [139] "t" "a" "a" "t" "t" "g" "a" "c" "c" "c" "a" "a" "a" "g" "a" "g" "a" "c" "c" "c" "a" "t" "t" [162] "a" "a" "t" "t" "t" "c" "a" "a" "c" "c" "g" "a" "c" "a" "g" "g" "a" "a" "t" "a" "a" "c" "a" [185] "a" "a" "c" "c" "t" "c" "t" "a" "t" "a" "t" "g" "g" "g" "c" "c" "a" "a" "c" "a" "a" "t" "t" [208] "t" "a" "g" "g" "t" "t" "g" "g" "g" "g" "t" "g" "a" "c" "c" "t" "c" "g" "g" "a" "g" "a" "a" [231] "c" "a" "g" "a" "a" "c" "a" "a" "c" "c" "t" "c" "c" "g" "a" "g" "t" "g" "a" "t" "t" "t" "a" [254] "a" "a" "t" "c" "a" "a" "g" "a" "c" "t" "a" "a" "c" "c" "a" "g" "t" "c" "a" "a" "a" "a" "g" [277] "t" "a" "t" "t" "a" "c" "a" "t" "c" "a" "c" "t" "a" "a" "t" "t" "g" "a" "t" "c" "c" "a" "a" [300] "a" "a" "a" "c" "t" "t" "g" "a" "t" "c" "a" "a" "c" "g" "g" "a" "a" "c" "a" "a" "g" "t" "t" [323] "a" "c" "c" "c" "t" "a" "g" "g" "g" "a" "t" "a" "a" "c" "a" "g" "c" "g" "c" "a" "a" "t" "c" [346] "c" "t" "a" "t" "t" "t" "t" "a" "g" "a" "g" "t" "c" "c" "a" "t" "a" "t" "c" "g" "a" "c" "a" [369] "a" "t" "a" "g" "g" "g" attr(,"species") [1] "Acinonyx_jubatus"

Mas como cada sequência tem um tamanho diferente, vamos alinhar elas. Vamos usar um serviço web porque minha maquina é uma uma droga e isso demora eras, quem não tem nem ideia do que é um alinhamento pode dar uma olhada no wikipedia aqui, mas vamos colocar pauzinhos e aumentar as sequências menores até que todas tenham o mesmo tamanho, é claro que tem uma lógica para se fazer isso, mas é assim que todas vão ficar do mesmo tamanho, então exportamos as sequencias em um arquivo fasta.

> write.dna(felidseq16S, "felidseq16S.fas", format = "fasta")

E alinhamos no site do Clustal.
So fazer ulpload do arquivo felidseq16S.fas que acabamos de salvar e importante, vamos fazer a arvore usando o programa PhylML, então a entrada tem que estar no formato Phylip, que é como deve estar organizado os dados, senão não conseguiremos progredir, mas esse formato é comum e basta mudarmos a opção no site para a saída, as sequencias alinhas serem retornadas no formato phylip.

Apos o alinhamento completo, que deve demorar alguns minutos, menos até, ja que são poucas sequências e curtas. Mas salvamos o resultado e vamos olhar como ficou.

> felidseq16Sali <- read.dna("felidae16s.phylip") > table(unlist(lapply(felidseq16Sali, length))) 379 37

Agora todo mundo ta alinhado, tem o mesmo tamanho, e agora além dos “ACTG”, temos uns “-” no meio das sequências. O DNA além de mutar, sofre inserções, deleções, o sequenciamento pode errar, comer uma base, etc, então por isso fica diferente, mas alinhamos levando isso em conta.

Feito isso, vamos desenhar a árvore usando Maximum Likelihood, da pra desenhar com MCMC também, mas aqui vamos usar um programa desenvolvido pelo Joe Felsenstein chamado PhyML, é so baixar gratuitamente na pagina do projeto. Usamos ele diretamente do R, a partir de uma função do ape chamada phymltest.

Aqui vão ser feitas muitas contas, que depois a gente vai estudando devagarinho, mas basicamente são combinações de uma tabela de chance de mutações, uma matriz na verdade dizendo qual a chance de um A ser trocado por um T na replicação do DNA, um C po um G, e não o complemento reverso, mas ser trocado mesmo. Outra parte são os modelos de evolução, as chances de acontecerem adições, deleções, etc.

Para cada combinação é desenhada uma árvore e calculada um ajuste para essa árvore com likelihood, o quão bom a árvore ficou, você pode escolher os modelos que quer usar, as combinações a excluir ou adicionar, suas próprias se quiser, mas vamos mandar tudo no default porque pelo menos eu não entendo muito bem dessas coisas mesmo, então deixa rodar pra ver o que acontece.

>phymltest.felid <- phymltest(seqfile="felidae16s.phylip", + execname="/home/augusto/git/recologia/PhyML-3.1/PhyML-3.1_linux32",append = F) > phymltest.felid nb.free.para loglik AIC JC69 1 -2622.287 5246.573 JC69+I 2 -2461.377 4926.753 JC69+G 2 -2427.008 4858.016 JC69+I+G 3 -2418.223 4842.445 K80 2 -2517.004 5038.008 K80+I 3 -2355.692 4717.384 K80+G 3 -2313.295 4632.590 K80+I+G 4 -2296.573 4601.146 F81 4 -2622.821 5253.641 F81+I 5 -2462.858 4935.716 F81+G 5 -2429.195 4868.390 F81+I+G 6 -2415.210 4842.420 F84 5 -2511.634 5033.267 F84+I 6 -2341.970 4695.940 F84+G 6 -2295.159 4602.319 F84+I+G 7 -2279.239 4572.477 HKY85 5 -2504.279 5018.559 HKY85+I 6 -2336.727 4685.454 HKY85+G 6 -2292.286 4596.573 HKY85+I+G 7 -2275.920 4565.840 TN93 6 -2477.334 4966.669 TN93+I 7 -2320.064 4654.129 TN93+G 7 -2264.062 4542.123 TN93+I+G 8 -2256.190 4528.379 GTR 9 -2461.696 4941.392 GTR+I 10 -2307.736 4635.472 GTR+G 10 -2255.337 4530.674 GTR+I+G 11 -2250.072 4522.144

Vai demorar um pouco calculando tudo, mas olha ai o resultado temos muitas combinações, cada linha dessa é uma árvore, e cada uma tem valores de ajustes, a gente ja viu AIC aqui, a ideia é a mesma, qualidade de ajuste, se a arvore é uma boa representação da filogenia, das diferenças e igualdades nas sequências que alinhamos.

Podemos ainda ver um súmario, assim como nas regressões usávamos um drop1, pra se os ajustes estão representando diferenças.

> summary(phymltest.felid) model1 model2 chi2 df P.val 1 JC69 JC69+I 321.82008 1 0.0000 2 JC69 JC69+G 390.55744 1 0.0000 3 JC69 JC69+I+G 408.12800 2 0.0000 4 JC69 K80 210.56508 1 0.0000 . . . 207 GTR GTR+I 307.92014 1 0.0000 208 GTR GTR+G 412.71776 1 0.0000 209 GTR GTR+I+G 423.24866 2 0.0000 210 GTR+I GTR+I+G 115.32852 1 0.0000 211 GTR+G GTR+I+G 10.53090 1 0.0012 >

São muitos casos, então eu suprimi um pouco, mas a ideia é que, se um modelo mais complexo representa melhor que um modelo mais simples as sequências, nessa caso árvores, ficamos com o mais complexo para não correr o risco de jogar informação valiosa fora, mas se um modelo simples e um complexo tem o mesmo ajuste, ficamos com o mais simples uai, porque complicar o que pode ser simples.

Podemos ver esses resultado, nessa forma gráfica aqui, esse gráfico ja vem pronto, é so dar plot no phymltest.felid que calculamos.

01

Bem agora podemos matar a curiosidade e ver como ficou, abrimos a árvore que ta salva em outro arquivo…

tr <- read.tree("felidae16s.phylip_phyml_tree.txt")

Escolhemos a GTR+I+G, que teve o melhor AIC (menor valor) e foi diferentes dos outros.

mltree.felid <- tr[[28]]

Trocamos os números de acesso do Genbank pelos nomes das especies

mltree.felid$tip.label <- taxa.felid[mltree.felid$tip.label]

Digite Galidia elegans no google e veja quem é esse carinha, ele não é um gatinho, ele esta aqui apenas para ser a raiz da árvore, um parente bem longe, então usamos ele como raiz mas depois tiramos ele fora porque não queremos ver ele na filogenia.

mltree.felid <- root(mltree.felid, "Galidia_elegans") mltree.felid <- drop.tip(mltree.felid, c("Crocuta_crocuta","Galidia_elegans"))

A agora é olhar o resultado do trabalho.

Comumente vemos representações assim.

02

Mas existem muitas outras formas de representar as arvores, aqui vão as possibilidades disponíveis no ape.

03

04

05

06

Bem, atropelamos bastante o processo todo, mas acho que isso aqui é igual a primeira vez que a gente faz qualquer coisa, a gente quer ver os finalmentes, depois vamos aprendendo os pormenores devagarinho, quando eu vi o arduino a primeira vez eu queria ver luzinhas piscando, depois eu vejo como é possível. Mas isso é bem complicado e tem muitos pormenores, o Joe Felsenstein disponibiliza aqui um livro, ebook, gratuitamente, é bem legal, bem complexo também, mas uma leitura interessante, para ser feita bem devagarinho. Existem também a lista de Phylogenia do R, a r-sig-phylo, onde o próprio autor principal e mantedor do pacote ape esta quase sempre respondendo a galera, o Emmanuel Paradis, alias o Felsenstein ja respondeu pessoas nessa lista também. Bem é isso, existe uma lista gigantesca de softwares para fazer exatamente esse mesmo processo que fizemos aqui, no wikipedia aqui tem uma lista dessas. Muito mais que ver somente o parentesco, a filogenia podem ser muito interessantes para perguntas como essa aqui.

Bem até a próxima, segue o script, além do repositório Recologia no github.

Referencias:
Paradis, E. 2012 Analysis of Phylogenetics and Evolution with R. 2nd ed. 386p. springer

Johnson W. E. & O’Brien S. J. 1997. Phylogenetic reconstruction of the Felidae using 16S rRNA and NADH-5 mitochondrial genes. Journal of Molecular Evolution 44: S98–S116.

#Instalando e abrindo o pacote para Phylogenia no R
#install.packages("ape")
library(ape)
 
#Criando uma lista com números de acesso ao Genbank
lista <- paste("AF006", seq(387, 459, 2), sep = "")
lista
 
#Fazendo Download pela internet
felidseq16S <-read.GenBank(lista)
 
#Olhando os dados
felidseq16S
 
table(unlist(lapply(felidseq16S, length)))
cat(felidseq16S[[1]], "\n", sep = "")
 
#Com letrinhas
read.GenBank(lista[1],as.character = T)
 
#Salvando o arquivo para alinhas
write.dna(felidseq16S, "felidseq16S.fas", format = "fasta")
 
#Apos o alinhamento, abra novamente e de uma olhada
felidseq16Sali <- read.dna("felidae16s.phylip")
table(unlist(lapply(felidseq16Sali, length)))
 
#Salvando os nomes das especies também
taxa.felid <- attr(felidseq16S, "species")
names(taxa.felid) <- names(felidseq16S)
 
#Desenhando as arvores, lembre-se de indicar corretamente o diretorio e o executavel do PhyML
getwd()
phymltest.felid <- phymltest(seqfile="felidae16s.phylip",
                             execname="/home/augusto/git/recologia/PhyML-3.1/PhyML-3.1_linux32",append = F)
 
#O resultado
phymltest.felid
 
#Testando diferenças nos modelos usados
summary(phymltest.felid)
 
#Visualizando os ajustes
plot(phymltest.felid)
 
#Abrindo a arvore que queremos olhar, definindo a raiz e tirando ela do plot
tr <- read.tree("felidae16s.phylip_phyml_tree.txt")
mltree.felid <- tr[[28]]
mltree.felid$tip.label <- taxa.felid[mltree.felid$tip.label]
mltree.felid <- root(mltree.felid, "Galidia_elegans")
mltree.felid <- drop.tip(mltree.felid, c("Crocuta_crocuta",
"Galidia_elegans"))
 
#Plots das arvores
plot(mltree.felid)
axisPhylo()
 
plot(mltree.felid,type="cladogram")
 
plot(mltree.felid,type="fan")
 
plot(mltree.felid,type="radial")
 
plot(mltree.felid,type="unrooted",cex=0.6)

Encontrando a máxima verosimilhança de um conjunto de dados usando um grid para a busca.

A boa e velha verosimilhança, ou likelihood em inglês, para quem preferir. Eu prefiro likelihood porque tem menos letras e não tem acento, eu sempre erro o acento da palavra verosimilhança, aliás eu erro a acentuação de um monte de palavras, talvez erre até um pouco na vida, infelizmente.

Mas a gente já falou a muito tempo atras de verossimilhança aqui. Aliás o que vamos fazer nesse post a gente já fez aqui também, quando falamos de aproximação por grid em estatística bayesiana, só que para a distribuição binomial.

Eu sempre imagino se sou só eu que demoro entender as coisas, ou tem mais gente lenta como eu.

Por exemplo, uma coisa que me perturbou muito tempo foi a diferenças entre  \mu (letra grega mu minúscula), que representa a média e porque as vezes tinha um  \bar{\mu} assim, com esse chapéuzinho. Se era a mesma coisa, se eram coisas diferentes.

Ai tinha a explicação que com chapéuzinho, com a barrinha em cima era a média das amostras, e sem chapeuzinho era a média da distribuição, e isso era estranho porque amostra vem dos dados, não deveria ser a mesma coisa?

Depois que eu vi nos livros do Mark Kéry a idéia de simular dados, para depois realizar uma análise de dados para ver se recuperávamos o padrão simulado que isso começou a fazer mais sentido, quem deu uma olhada nesse post aqui e nesse outro aqui viu que eu fiz um gráfico tentando resumir tudo, vou colocar ele denovo aqui em baixo.

02
03

A diferença é que quando eu simulo dados, eu tenho como colocar no gráfico as medidas “reais”, agora no segundo caso, quando usamos a análise em dados reais, não tem linha “real” nem nunca vai ter, podemos pensar que a natureza fez os dados, mas ela não mostra os parâmetros que usou, a gente estima e acredita que funcionou e para por ai. A gente nunca vai saber como diabos são os dados reais, e essa é a diferença entre  \mu e  \bar{\mu} , ou qualquer coisa com e sem chapeuzinho, qualquer parâmetro, a gente estima eles a partir de amostras, então a gente tem a medida com chapeuzinho, a medida sem chapeuzinho nunca saberemos ao lidar com dados reais, so que como as vezes a gente simula dados, simula amostras a partir de um  \mu e consegue ter um  \bar{\mu} parecido, a gente acredita que esse procedimento funciona e usa ele nos dados reais, mas não temos como ter certeza, nada garante. No caso da estatística bayesiana a gente ainda combina esse  \bar{\mu} com um prior, com algo que a gente acha para gerar uma distribuição posterior, que é um intervalo de valores que achamos que o  \mu ta la dentro, mas nada é exato.

Certo, um bocado de blablabla, mas vamos fazer algo legal, vamos estimar o likelihood usando um grid de valores de média ( \bar{\mu} ) e desvio ( \bar{\sigma} ).

Então a gente gera uma amostra a partir de uma distribuição usando a função rnorm do R que gera números ao acaso seguindo os parâmetros que estipulamos.

media <- 15
desvio <- 2.5
n <- 100
x<-rnorm(n,media,desvio)

Certo, agora imagine que esse passo foi feito por outro pessoa, que a gente não tem contato, a única coisa que a gente ve são os valores em x, ou seja os 100 valores que estão dentro de x.

Como sempre e todo o sempre, a primeira coisa que a gente faz é um gráfico para ver como são os dados que temos em mãos

01

Certo, so de bater o olho nos parece uma distribuição normal, mas porque?
Oras, porque tem esse jeitão em forma de sino, com bordinhas, ai um levantando no meio, uma barriguinha, acredito todo mundo deve ter uma imagem de sino na cabeça.

Mas então o primeiro passo é imaginar uma distribuição adequada. Uma idéia aqui também seria fazer uma tabelinha com as forma das distribuições mais conhecidas, e ir no olhômetro.

Mas supondo que estamos satisfeitos com a idéia de usar uma distribuição normal, precisamos ajustar os parâmetros da distribuição a usar, a distribuição normal é baseada em dois parâmetros, a média e o desvio (também conhecido como variância, precisão, depende do caso e do uso do freguês). Esses parâmetros, como falamos ali em cima, podem ser descritos comumente com as letras gregas  \mu e  \sigma (isso é um sigma minúsculo). Mas quando eu falo ajustar parâmetros eu falo que temos que dar valores para eles.

Mas vejamos, temos uma forma, a distribuição normal, agora precisamos encaixar ela naquele histograma la em cima.
Vamos ver o que acontece para diferentes valores desses parâmetros:

02

Olha ai, vejam que alguns casos tudo fica muito bonitinho, um encaixe bom, todas as barrinhas bem perto da linha, mas alguns casos, por exemplo média 18 e desvio 5, o encaixe é horrível, esses valores são muito ruins, mas com média 14 e desvio 3 até que a curvinha, a distribuição normal encaixa bonitinho nos dados. Mas aqui usamos 4 valores de média, 12, 14, 16, 18 e usamos 4 valores de desvio, 2, 3, 4 e 5. Ai fomos no olhômetro e ficamos tentando perceber qual a melhor combinação que ficava melhor.

Basicamente é isso que a gente quer, no entanto a gente poderia aumentar a precisão aumentando a quantidade de valores possíveis, ao invés de 12 a 18 pulando de 2 em 2, a gente podia pular de 1 em 1, de 0.5 em 0.5, ou menos, mas ai o problema é que conforme aumentamos esse grid de busca, aumentamos a quantidade de gráficos que temos que olhar, ou seja usando a técnica do olhômetro nos estamos limitados a quantos gráficos podemos olhar se o encaixe ficou bom ou não. E se aumentarmos a precisão então ferrou, porque as possibilidades crescem exponencialmente conforme aumentamos o número de desvios e médias que queremos olhar, ja que o número de casos é o número de médias vezes o número de desvios usado.

Agora a gente chegou ao ponto de o que diabos o likelihood ou verosimilhança faz, ela substitui o nosso “olhômetro”, ele ve se o encaixe ficou legal, a gente calcula um número que quanto maior, melhor o encaixe.

Algebrismos a parte, agora já deu para imaginar o procedimento, pegar o feeling da coisa, a gente vai escolher valores arbitrários, bastantes valores, claro que sendo razoável nessa escolha, não vou testar o caso de média um zilhão, porque pelo gráfico da para ver que não faria sentido um valor assim, então nem precisa desperdiçar esforço num número assim, mas escolhendo números razoáveis, numa precisão que de tempo para o computador processar, a gente vai fazer a curva, ver se ela encaixa legal, e registrar o valor do encaixe, ou o likelihood, depois é so ver qual é melhor.

Como a gente viu na figura anterior, como eu testo todas as combinações, eu posso registrar numa matriz todos os valores de todas as combinações possíveis de likelihood, e depois vejo o maior.

Notem que os valores sempre ficam negativos, mas o menos negativo (mais próximo de zero) é maior que um número muito muito grande mas negativo.

Depois de tudo isso eu posso olhar o resultado do maior valor na matrix, ou procurar direto o maior valor.

> max_lh<-max(log_lh,na.rm =T) > max_lh [1] -233.3122

E temos ai nosso melhor encaixe, é so olhar onde ele esta na matriz que temos o valor de média e desvio que deveríamos estar mais inclinados a acreditar que deva ser o valor real dos dados.

Mas eu sou sempre horrível para ver valores e matrizes, então eu faço gráficos.

Nesse caso a gente pode até usar uma analogia, a média e desvios reais são uma flecha que acertou o alvo em cheio, so que a gente não ve a flecha, a gente so sabe que ela acerta em cheio, o que a gente ta tentando construir é o alvo, a gente tenta desenha o alvo e espera que a flecha tenha acertado em cheio. Mas o alvo que a gente faz nunca é perfeito. As vezes pior ou melhor, dependendo das nossas escolhas para desenhar ele, mas olhando esse caso aqui específico ficou assim.

03

Aqui, olha ai como a gente sabe onde a flecha está, a gente pinta um pontinho cinza para ver se ta dando certo esse esquema, ai pintamos de vermelho a area mais provavel, e conforme a gente muda as cores do vermelho para o azul, a chance é menor (os likelihood são valores menores).

E a flechada acertou mesmo no alvo que a gente pintou, então a gente fez um bom trabalho, esse procedimento parece funcionar bem.

Se a gente compara o valor de média que ficou no melhor encaixe, o maior valor de likelihood, ou maximum likelihood, ou máxima verosimilhança, vemos o seguinte, em relação as possibilidades que testamos:

> mu_mle <- muhat[maxrow] > mu_mle [1] 14.37779 > mean(x) [1] 14.67122

Esse é nosso melhor chute de valor para média, ou o  \bar{\mu} , se fazemos um gráfico com todos os valores que usamos no grid para média, e colocamos esse como o separador de águas, menor que ele pintando de verde e maior pintando de azul vemos que:

04

Olha ai a média real ta quase junto do “divisor de águas”. Muito legal né.

Agora qual a relação disso com os modelos que fazemos todos os dias, a gente faz esse processo toda hora que chama o comando lm, ou glm para outras distribuições, ou muitos outros comandos do R.

Veja os valores que estimamos para média e desvio com nosso método, e compare com a boa e velha formulinha da média e do desvio padrão.

> mu_mle [1] 14.37779 > sd_mle [1] 2.58246

Agora note que são quase a mesma coisa que o valor de média e desvio usando as funções mean e sd.

> mean(x) [1] 14.67122 > sd(x) [1] 2.507243

Média e desvio são o corte de caminho para todo esse procedimento? Ai você pensa, então tudo isso é inútil?

Não, lembra la atras quando usamos um procedimento de grid para estatística bayesina? Isso é o procedimento geral, da para usar com qualquer distribuição, no mundo não existe somente a distribuição normal, existe muitas delas, na página do John Cook tem uma diagrama só para termos uma idéia das possibilidades de distribuições, ache a normal nesse bolo ai.

Agora a relação com quando a gente usa o comando lm do R é o seguinte. Como só temos um nível, o x, a gente pode usar ele para estimar a média falando na formula que so existe um nível (usando x~1).

Então usamos o comando:

summary(lm(x~1))

E recebemos a resposta:

Call: lm(formula = x ~ 1) Residuals: Min 1Q Median 3Q Max -6.2405 -1.7020 -0.5832 1.7584 5.8910 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 14.6712 0.2507 58.52 <2e-16 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 2.507 on 99 degrees of freedom

Agora olha que fascinante, notem que aqui:

Estimate Std. Error t value Pr(>|t|) (Intercept) 14.6712 0.2507 58.52 <2e-16 ***

Estimate é a média, não ta idêntico ao nosso valor porque aqui é usado cálculo infinitesimal, que é o mesmo que um grid infinitamente preciso, com infinitos números, mas como cálculo é complicado, começamos de um jeito mais fácil, com um grid, e temos um jeito agora de conferir se quando tentarmos isso com calculo estaremos indo no caminho certo, ou não, viu é muito útil seja para um fim como para outro.
Mas então, o estimate é o meio da bola vermelha do gráfico la em cima, e o Std. Error, ou erro, é a precisão com que podemos falar que esse valor ta certo. Lembre-se que na figura la em cima, o vermelho tem uma considerável área, a gente ta falando que qualquer valor ali ta valendo, ou seja, esse valor é relativo ao tamanho da bolota vermelha la naquele gráfico, se ele fosse menor, a bola seria menor, se ele for grande a bola é grande do local onde assumimos que a flecha acertou.

Ta certo, agora olha esse número aqui

Residual standard error: 2.507 on 99 degrees of freedom

Olha ai é o nosso desvio do modelo. Esses são os parâmetros do modelo que dão a forma aquela curva para encaixar nos dados.

Por último a gente ve se ficou centralizadinho, porque se não ficou centralizado, não encaixou legal, certinho, o que pode ocorrer se la em cima escolhemos o modelo errado, a distribuição errada.
Então diminuímos todos os valores de média dos dados e fazemos um histograma.

05

E se quiser fazer isso com o ajuste que o comando lm fez, é so usar o comando resid no modelo gerado pelo lm e fazer um histograma, dessa forma:

hist(resid(lm(x~1)))

Aqui eu to fazendo um histograma do resíduo do modelo produzido pelo lm, se achou complicado é so fazer parte a parte.

modelo<-lm(x~1)
residuos<-resid(modelo)
hist(residuos)

Ambos os jeitos vão fazer exatamente a mesma coisa, mas as vezes usar um comando dentro do outro pode ser confuso logo de primeira, até aprender as coisas.

Que legal né, agora a gente não precisa olhar os valores juntos, média e desvio juntos, a gente pode olhar de um em um, ou olhar o que nos interessa mais. Basta a gente pegar o Maximum likelihood para os outros valores (nesse caso o desvio) e olhar como fica a distribuição dos valores para o que nos interessa, nesse caso a média.

06

Isso que as pessoas chamam de profile, é olhar so o que interessa, você mantem todos os outros parâmetros no melhor encaixe e olha o que te interessa e tira suas conclusões.

Bem, eu acho legal esse tópico, primeiro que ajuda a fazer uma relação com qual a linha de raciocínio que as análises usam para produzir conclusões, além de ajudar a interpretar mais valores que somente o p quando fazemos alguma análise.

E outra coisa, é que, quando começamos a perceber esse esquema do likelihood, da para perceber melhor a relação entre diferentes análises, e ver como regressão linear, anova, regressão de poisson, regressão binomial ou logística, é tudo o mesmo esquema, variando alguns passos apenas.

Bem até o próximo post para quem teve saco de ler até aqui ^^

#################################################################
#Máxima verossimilhança da média e desvio utilizando
#grid para a busca
#veja o script original em:
#http://www.sortie-nd.org/lme/Course_Schedule_2012/Day_1/Lab%201%20-%20Section%201.R
##################################################################
set.seed(51)
 
#Simulando uma amostra com distribuição normal
media <- 15
desvio <- 2.5
n <- 100
x<-rnorm(n,media,desvio)
#A partir daqui a gente tem que fingir que não sabe os valores
#reais de média e desvio
 
#Histograma
hist(x,main="Histograma",ylab="Frequência",ylim=c(0,50))
 
 
#Grid testando alguns valores de parametros
par(mfrow=c(4,4))
for(i in seq(12,18,by=2)) {
  for(j in seq(2,5,by=1)) {
    hist(x,prob=T,main=paste("Média=",i,"Desvio=",j),ylim=c(0,0.3),xlim=c(6,24))
    curve(dnorm(x,i,j),add=T,col="red",lty=3,lwd=4)
  }
}
 
#Agora vamos fazer um grid para calcular a verosimilhança de diferentes
#combinações de estimativas
 
#Defina as dimensões do grid
range <- 1
increment <- 0.01
 
#Crie um vetor com a mudança dos valores do grid
muhat <- seq(mean(x)-(range*mean(x)),mean(x)+(range*mean(x)),by=mean(x)*increment)
sdhat <- seq(sd(x)-(range*sd(x)),sd(x)+(range*sd(x)),by=sd(x)*increment)
 
#Iniciando uma matrix para salvar os resultados
log_lh <- array(0,dim=c(length(sdhat),length(muhat)))
 
#Cuidado, na função original acho que o autor invertou o loop do j e do i
#Como a matrix original era quadrada não tinha problema, mas se vc alterar o
#grid vc tera um erro ao tentar olhar fora do limite da matrix.
 
for (i in 1:length(sdhat)) {
  for (j in 1:length(muhat)) {
    log_lh[i,j] <-  sum(dnorm(x,muhat[i],sdhat[j],log=T))
  }
}
 
#Essa linha:
#log_lh[i,j] <-  sum(dnorm(x,muhat[i],sdhat[j],log=T))
#É o mesmo que pensar nesse loop
#for (k in 1:length(x)) {
#  log_lh[i,j] <- log_lh[i,j] + dnorm(x[k],muhat[i],sdhat[j],log=T)
#}
#So para ir mais rapido :)
#Mas vale a pena olhar como fica o loop, pelo menos eu
#entendo melhor esses loops com for que operação com vetor
 
#Qual o Maximum likelihood
max_lh<-max(log_lh,na.rm =T)
max_lh
 
# Agora o gráfico do alvo de arco e flecha
xlimit <-range(muhat)
ylimit <-range(sdhat)
clevels <- seq(max_lh - 10,max_lh,by=0.5)
cores<-colorRampPalette(colors=c("blue","red"))
 
#Grafico Alvo
par(mfrow=c(1,1))
filled.contour(x=muhat,y=sdhat,z=log_lh,xlim=c(13,17),ylim=c(1,4),
               levels=clevels,color.palette=cores,frame=F,
               plot.axes = { axis(1); axis(2);
                             contour(x = muhat, y = sdhat, z = log_lh,
                                     xlim = c(13,17), ylim =c(1,4),
                                     levels=seq(min(clevels),max(clevels),
                                       length=11)[-11],drawlabels=F,add=T,lwd=2,
                                     lty=3)
                             points(media,desvio,pch=19,cex=1.5,col="darkgray")
                           }
               )
 
 
#Mas ainda não esta pronto, agora extraimos os indices dos maximos valores da
#maxima verosemelhança
mx<-max.col(log_lh)
 
#Agora vamos linha por linha e achamos qual linha tem o maior valor
#iniciando um valor para ser substituido, so colocar um valor gigante
max <- -1000000000
for (i in 1:nrow(log_lh)) {
  if (log_lh[i,mx[i]] > max)
    {  maxrow <- i
       maxcol <- mx[i]
       max <- log_lh[i,mx[i]]
     }
}
 
#Estimativa da maxima verosimilhança da média
mu_mle <- muhat[maxrow]
mu_mle
mean(x)
 
#Figura 4
plot(muhat,rep(1,length(muhat)),col=ifelse(muhat<=mu_mle,"green","blue"),
     yaxt="n",frame=F,xlim=c(12,16),ylab="",pch=19)
abline(v=mean(x),lwd=2,lty=2,col="red")
legend("topleft",bty="n",col=c("green","blue"),pch=19,
       legend=c("Valores de muhat<=mu_mle",
         "Valores de muhat>mu_mle"))
legend("topright",bty="n",lwd=2,lty=2,col="red",legend="Média Real")
 
 
#Estimativa da maxima verosimilhança do desvio padrão
sd_mle <- sdhat[maxcol]
sd_mle
sd(x)
 
#Usando o comando lm para estimar os valores
summary(lm(x~1))
 
#Agora calculamos os residuos (erros) e vemos se ele são
#normalmente distribuidos com média zero
res <- x - muhat[maxcol]
 
#Histograma com curva estimada
hist(res,freq=F,xlab="Residual",main="Residual standard error")
 
#Nos podemos sobrepor uma PDF normal nesse grafico, dados os MLEs
pdfx <- seq(from = -2*desvio, to = 2*desvio, by = 0.5)
pdfy <- dnorm(pdfx,0,sd_mle)
points(pdfx,pdfy,type="p",pch=19,col="Blue",cex=2)
lines(pdfx,pdfy,type="l",col="Blue",lwd=2)
 
 
#Profile likelihoods são graficos da mudança do likelihood de uma parametro
#em especial enquanto o outro parametro é mantido no maximo.
#Por exemplo o Profile da média
plot(muhat,log_lh[,maxcol],xlim=c(0,30),frame=F)
abline(max_lh-2,0)
abline(v=media,lwd=2,col="red",lty=2)
legend(0,-300,legend="Média Real",lty=2,lwd=2,col="red",bty="n")

Biscoitinhos e o Teorema de Bayes

biscoito

La vamos nos para outro exemplo usando o Teorema de Bayes.

Pequenos exemplo são sempre legais para ajudar a internalizar alguns conceitos. E aqui vamos em mais um relacionado ao Teorema de Bayes, referente ao uso dele.

Eu vi esse exemplo um dia espiando videos no youtube. Alias o video que estou falando é esse aqui. É meio longo mas esse parece ser um tiozinho legal que sabe o que esta falando. Mas vamos ao exemplo:

Suponha que existe dois potes cheios de biscoitos. O pote #1 tem 10 biscoitos de chocolate e 30 de morango, enquanto o pote #2 tem 20 biscoitos de cada um desses sabores. Nossa amiga Zildamara pega um desses potes ao acaso e depois pega um biscoito ao acaso dentro desse pote. Nos vamos assumir que a Zildamara não tem motivos para tratar um pote de forma diferente do outro, assim como os biscoitos. Dado que o biscoito que a Zildamara pegou é um biscoito de morango. Qual a probabilidade que o biscoito que tenha vindo do pote #1?

Certo, temos nossa pergunta, agora vamos usar a idéia do Teorema de Bayes para atacar esse problema.

Primeiro, a pergunta é qual a probabilidade, qual a chance de o biscoito ter vindo do pote #1?
Podemos começar pensando o seguinte, existem dois potes, o pote #1 e o pote #2, e como a Zildamara não liga para qual pote pegar, então a chance de cada pote ter sido escolhido é meio a meio, 50% para cada pote, então estabelecemos nossas chances:

p(H1)={1/2}
p(H2)={1/2}

Estamos representando aqui a chance de cada hipótese, eu estou chamando de p(H1) a hipótese do pote escolhido ter sido o #1, e p(H2) a chance do pote #2 ter sido escolhido. Então não tem porque se intimidar com essa nomenclatura, p é de probabilidade e o que esta dentro do parênteses é o que estamos tratando.

Até aqui bem fácil. Agora vamos para o segundo passo, vamos calcular a verosimilhança de cada uma dessas hipóteses. Ja falamos de verosimilhança antes, ou likelihood, nesse post, mas é só pensar no seguinte. Qual a chance desse biscoito de morango ter saído do pote #1?

Bom, estamos falando de um biscoito de morango, no pote #1 tem 10 biscoitos de chocolate e 30 biscoito de morango, então é meio obvio que com mais biscoitos de morango, a chance de sair um biscoito de morango é maior, mas quanto maior? Uai são 30 biscoitos de morango num pote que tem 40 biscoitos, que da 30/40 a chance de sair um biscoito de morango, que é o número de biscoitos de morango sobre o total de biscoitos do pote, que é o mesmo que 75% ou 3/4.
Agora para a verossimilhança para o pote #2 seguimos a mesma idéia, são 20 biscoitos de morango, num pote de 40 biscoitos, o que nos da 20/40, que nada mais é que 50% ou 2/4 cortando os zeros, que ainda podemos simplificar para 1/2.

Vamos anotar esses valores aqui para não esquecer

p(E|H1)=3/4
p(E|H2)=1/2

Agora temos um prior e a verosimilhança, já podemos usar o teorema de bayes para obter obter a probabilidade que queremos saber, que é a chance que o biscoito tenha vindo do pote #1, lembre-se que acima temos qual a chance do biscoito ter saído do pote #1 dado que só existe o pote #1, o p(E|H1) é o que esperamos dado que o pote #1 tenha sido escolhido (probabilidade condicional nesse post), mas não temos certeza que ele foi escolhido, ele tem somente 50% de chance de ter sido escolhido (p(H1)={1/2}).

Mas vamos relembrar o teorema de bayes que é:

p(A|B)={P(A) \cdot {P(B|A)}\over{P(B)}}

Mas que usando para a idéia de probabilidade temos:

p(H|E)={P(H) \cdot {P(E|H)}\over{P(E)}}

Aqui H é a hipótese em questão e E é a evidência que temos.
E olha só que loucura, a gente calculou o p(H) que chamamos de prior e calculamos o p(E|H) que é a verossimilhança, se não fosse esse p(E) ai faltando, era só substituir os valores e ver o milagre acontecer.

Mas p(E) vem da chance dada todas as hipóteses possíveis, o E é de evidência, então p(E) é a soma das verossimilhanças sobre todas as hipóteses possíveis, que nesse caso, so temos duas hipóteses possíveis, quais são elas?
Que o biscoito veio do pote #1 que é p(H1) \cdot {p(E|H1)} e que o biscoito veio do pote #2, que é, p(H2) \cdot p(E|H2), agora sim ficou fácil não?

Então vamos substituir os valores para responder o problema:

Lembrando os valores

p(H1)=1/2
p(E|H1)=3/4
p(E|H2)=1/2

Queremos:

p(H1|E)={(P(H1) \cdot {P(E|H1))} \over{P(E)}}

Que é a mesma coisa então que:

p(H1|E)={(P(H1) \cdot {P(E|H1))} \over {(p(H1) \cdot p(E|H1)+p(H2) \cdot p(E|H2))}}

Essa última ficou grande, mas é so olhar para ela que não tem nada que não vimos acima 🙂

Aqui eu só substitui o p(E) por p(H1) \cdot p(E|H1) + p(H2) \cdot p(E|H2), a evidência dado todas as possibilidades

Mas então colocando valores:

p(H1|E)=(1/2 \cdot 3/4)/(1/2 \cdot 3/4 + 1/2 \cdot 1/2)

Daqui para frente é só álgebra chata hehe

Primeiro multiplica as fração:

p(H1|E)=(3/8)/(3/8+1/4)

Calculo o mínimo divisor comum para poder somar as frações do denominador

p(H1|E)=(3/8)/(3/8+2/8)

Depois eu somo

p(H1|E)=(3/8)/(5/8)

Ai uma fração dividida pela outra, é só inverter elas e multiplicar, isso é regra de álgebra do segundo grau, eu num faço a mínima idéia como demonstra que isso funciona, mas acredito que um dia saberei, e também acredito que funciona.
p(H1|E)=(3/8).(8/5)

Ai eu multiplico
p(H1|E)=(24/40)

E simplifico para um número que da para entender
p(H1|E)=(3/5)

E 3/5 é o mesmo que 60%.

Ou seja, dado que o primeiro biscoito que saiu é de morango, a chance é de 60% que o pote escolhido tenha sido o pote #1. E isso não é difícil de aceitar intuitivamente, poxa tudo bem que cada pote tem 50% de chance de ser escolhido, mas no pote #1 tem muito mais biscoito de morango, faz sentido ele ter sido escolhido não ja que foi um biscoito de morando que a Zildamara pegou e não que o pote #2. Não possa produzir esse resultado também, mas a chance de produzir esse resultado é menor, a verosimilhança é menor. Assim é mais verossímil chutar que a Zildamara pegou o pote #1. Verossímil é uma palavra chique não?

Agora olha que legal, quem olhou esses post usando um conjugate prior e esse post usando a aproximação por grid, ambos sobre análise sem usar o MCMC, da para lembrar que eu falei que a distribuição posterior sempre vai estar entre o prior e o likelihood.

Aqui isso faz todo o sentido, o prior nesse caso é 50%, como assumimos que a Zildamara não ligava para qual pote pegava, assumimos que cada pote tem 50% de chance, mas a verosimilhança era de 3/4 para o pote #1, 3/4 é 75%, e a distribuição posterior é sempre algo entre o prior e a verosimilhança, que no nosso caso foi 3/5, ou 60%.

Nos dois exemplos anteriores usando conjugate prior e aproximação por grid, a gente sempre via isso, e fazendo exemplos la a gente sempre vai sempre ver isso, a distribuição posterior sempre vai ser algo entre o prior e a verosimilhança (likelihood).

Outra coisa é ver o que é o denominador, no caso da aproximação por grid, era uma somatória. E olhemos denovo o que fizemos.

p(E)=(p(H1) \cdot p(E|H1) + p(H2) \cdot p(E|H2))

Não deve ser difícil visualizar que existem duas hipóteses, e somamos as duas, então isso é uma generalização de:

p(E)=\Sigma (p(H) \cdot p(E|H))

Essa é a letra muito grega muito louca para somatório, ola senhor  \Sigma (chama-se sigma maiúsculo), a diferença para o conjugate prior é que somar coisas contínuas é cálculo, é uma integral, p(E)= \int {(p(H) \cdot p(E|H))}.

Bem esse é um exemplo bem legal, bobinho mas legal, mas que pode ajudar bastante a entender melhor como algumas coisas funcionam.

Aqui tem mais algum material do tiozinho do vídeo. E tiozinho é maneira de falar, não encarem isso de forma pejorativa, o nome do cara é Allen Downey e ele é professor de Ciências da Computação e aparenta ser bem inteligente :).

Parece que é um seminário sobre estatística bayesiana usando python. Mas isso é um tópico a ser explorado mais para frente, mas vale a pena uma olhadinha, na parte 1 existem mais exemplos e probleminhas legais para ajudar a entender melhor o teorema de bayes, que é tão simples a formula mas ao mesmo tempo tão difícil de internalizar.

Aproximação por grid em estatística bayesiana.

A gente já viu nesse post como usar estatística bayesiana sem o mcmc. Onde definimos um prior usando uma distribuição (a beta), calculávamos o likelihood e aplicávamos o teorema de bayes.

Acontece que existe outra forma ainda de fazer as coisas, naquele exemplo trabalhamos com uma distribuição continua, ou seja que que acomodava qualquer valor, seja 0.5 ou 0.55. Mas podemos definir as coisas sem usar uma distribuição, definindo um “grid” de possibilidades.

Ou seja eu defino um grid que represente todas as possibilidades, mas eu posso fechar ele de forma a não existir o 0.55, ou a chance de colonizar um arbusto é 50% ou é 60%, sem meio termo como na distribuição beta.

Isso é legal que agora o divisor do teorema de bayes vai ser uma soma, a soma de todas as possibilidades, e não uma integral como era no outro caso.

Então definimos nossa probabilidades do prior:

01

Uma coisa, a gente usa essas barras e não um histograma para ser claro na precisão. Ou seja eu to falando de 10 em 10%, eu poderia ser mais preciso, mas ai seriam mais barras e se quisesse ser infinitamente preciso e cobrir todos os valores possíveis eu usaria a distribuição beta horas.

Outra coisa legal é que agora vemos como é possível estabelecer um prior que não respeita as forma possíveis da distribuição beta como vimos anteriormente.
Nesse caso podemos pensar no prior como duas vertentes, uma que acha que a ocupação é baixa, está em torno dos 20% e outra que é alta, esta em torno dos 60%.

02

E vamos la, coletamos nossos dados, olhamos vários arbustos novamente, calculamos o likelehood, que de 10 arbustos, 7 tinham aranhas e começamos a achar que realmente é algo em torno de 60% a ocupação desses arbustos. Ou seja a parte do prior que estava ocupando uma baixa taxa de ocupação é não la muito verossímil com o que observamos, afinal quase todos os arbustos estão ocupado, não há como a ocupação não ser alta.

Bem acho que assim como o método exato falado no post anterior, a descrição das probabilidades por grid não vai ser muito útil, mas é legal para fazer alguns experimentos e entender melhor a aplicação do teorema de bayes, ainda mais que nem integral precisamos nesse caso, além criar alguns insights de qual o efeito da precisão usada no grid (usar poucos pontos ou muitos para descrever o mundo). Criar prior esquisitos e ver o que acontece ou ainda tentar com quantidades diferentes de dados para entender melhor o que vai acontecer. Daqui é so mudar o script, alterando a precisão do grid e mudar a coleta de dados (a forma como se gera ela com a função rbinom).

####################################################################
#Aproximação por grid
#Exemplo Baseado no livro Doing Bayesian Data Analysis
#http://www.indiana.edu/~kruschke/DoingBayesianDataAnalysis/
#################################################################
#Prior
Theta<-seq(0.01,0.99,0.05)
pTheta<-c(0,0.025,0.05,0.075,0.2,0.075,0.05,0.025,0,
0.025,0.05,0.075,0.2,0.075,0.05,0.025,0,0,0,0)
sum(pTheta)
 
#grafico do prior
plot(pTheta,type="n",frame.plot=F,xaxt="n",xlim=c(0,21))
segments(1:length(pTheta),0,1:length(pTheta),pTheta,lwd=3)
axis(1,at=1:20,labels =Theta,cex.axis=0.5)
 
#Dados
set.seed(123)
Data<-rbinom(10,1,0.8)
 
#Como fazer os gráficos
credib<-0.95
nToPlot<-length(Theta)
 
z = sum( Data==1 ) # Número de arbustos Ocupados
N = length( Data ) # Número de arbustos Observados
 
# likelihood dos dados para cada valor de Theta que usamos
pDataGivenTheta = Theta^z * (1-Theta)^(N-z)
 
# Computando a evidencia para a distribuição posterios
pData = sum( pDataGivenTheta * pTheta )
 
#Multiplicando o Prior vezes o Likelihood dividido pela Evidencia
#(Teorema Bayes)
pThetaGivenData = pDataGivenTheta * pTheta / pData
 
# Gráfico dos Resultados
layout( matrix( c( 1,2,3 ) ,nrow=3 ,ncol=1 ,byrow=FALSE ) ) # Painel 3×1
par( mar=c(3,3,1,0) , mgp=c(2,1,0) , mai=c(0.5,0.5,0.3,0.1) ) # Margens
dotsize = 1 #Tamanho das bolinhas
 
#Se escolhermos 1 milhão de valores para Theta, num vai dar certo o plot
#pois vão ser muitos pontos, mas da pra fazer um plot menor com media entre
#thetas, como definimos la em cima quantos ponto plotar
 
nteeth = length(Theta)
if ( nteeth > nToPlot ) {
thinIdx = seq( 1, nteeth , round( nteeth / nToPlot ) )
if ( length(thinIdx) < length(Theta) ) {
thinIdx = c( thinIdx , nteeth ) # makes sure last tooth is included
}
} else { thinIdx = 1:nteeth }
 
# Grafico do prior.
 
meanTheta = sum( Theta * pTheta ) # mean of prior, for plotting
plot( Theta[thinIdx] , pTheta[thinIdx] , type="p" , pch=19 , cex=dotsize ,
xlim=c(0,1) , ylim=c(0,1.1*max(pThetaGivenData)) , cex.axis=1.2 ,
xlab=bquote(theta) , ylab=bquote(p(theta)) , cex.lab=1.5 ,
main="Prior" , cex.main=1.5 )
if ( meanTheta > .5 ) {
textx = 0 ; textadj = c(0,1)
} else {
textx = 1 ; textadj = c(1,1)
}
text( textx , 1.0*max(pThetaGivenData) ,
bquote( "média(" * theta * ")=" * .(signif(meanTheta,3)) ) ,
cex=2.0 , adj=textadj )
# Grafico do likelihood: p(Data|Theta)
 
plot(Theta[thinIdx],pDataGivenTheta[thinIdx],type="p",pch=19,cex=dotsize
,xlim=c(0,1) ,cex.axis=1.2 ,xlab=bquote(theta)
,ylim=c(0,1.1*max(pDataGivenTheta))
,ylab=bquote( "p(D|" * theta * ")" )
,cex.lab=1.5 ,main="Likelihood" ,cex.main=1.5 )
if ( z > .5*N ) {
textx = 0 ; textadj = c(0,1)
} else {
textx = 1 ; textadj = c(1,1)
}
 
text(textx,1.0*max(pDataGivenTheta),cex=2.0,
bquote( "Data: z=" * .(z) * ",N=" * .(N) ) ,adj=textadj )
 
# Grafico do posterior.
meanThetaGivenData = sum( Theta * pThetaGivenData )
plot(Theta[thinIdx],pThetaGivenData[thinIdx],type="p" ,pch=19,cex=dotsize
,xlim=c(0,1) ,ylim=c(0,1.1*max(pThetaGivenData)) ,cex.axis=1.2
,xlab=bquote(theta) ,ylab=bquote( "p(" * theta * "|D)" )
,cex.lab=1.5 ,main="Posterior" ,cex.main=1.5 )
 
if ( meanThetaGivenData > .5 ) {
textx = 0 ; textadj = c(0,1)
} else {
textx = 1 ; textadj = c(1,1)
}
 
text(textx ,1.00*max(pThetaGivenData),cex=2.0,
bquote("média(" * theta * "|D)=" * .(signif(meanThetaGivenData,3)) ),adj=textadj )
text(textx,0.75*max(pThetaGivenData),cex=2.0,
bquote( "p(D)=" * .(signif(pData,3)) ) ,adj=textadj )
 
# Marcando o intervalo de maxima densidade, decidimos o intervalo la em cima.
#função para calcular o HDI
HDIofGrid = function( probMassVec , credMass=0.95 ) {
sortedProbMass = sort( probMassVec , decreasing=T )
HDIheightIdx = min( which( cumsum( sortedProbMass ) >= credMass ) )
HDIheight = sortedProbMass[ HDIheightIdx ]
HDImass = sum( probMassVec[ probMassVec >= HDIheight ] )
return( list( indices = which( probMassVec >= HDIheight ) ,
mass = HDImass , height = HDIheight ) )
}
 
HDIinfo = HDIofGrid(pThetaGivenData,credMass=credib)
points(Theta[ HDIinfo$indices ],rep(HDIinfo$height,length( HDIinfo$indices ))
,pch=19,cex=0.5,col="red" )
text(mean(Theta[HDIinfo$indices]), HDIinfo$height,
bquote( .(100*signif(HDIinfo$mass,3)) * "% HDI" ),
adj=c(0.5,-1.5),cex=1.5 )
 
#Marcando uma linha no Intervado de 95%
lowLim = Theta[min(HDIinfo$indices)]
highLim = Theta[max(HDIinfo$indices)]
lines(c(lowLim,lowLim),c(-0.5,HDIinfo$height),type="l",lty=2,lwd=1.5)
lines(c(highLim,highLim),c(-0.5,HDIinfo$height),type="l",lty=2,lwd=1.5)
segments(lowLim,HDIinfo$height,highLim,HDIinfo$height,lty=2,lwd=1.5)
text(lowLim,HDIinfo$height,bquote(.(round(lowLim,3))),
adj=c(1.1,-0.1),cex=1.2 )
text(highLim,HDIinfo$height,bquote(.(round(highLim,3))),
adj=c(-0.1,-0.1),cex=1.2 )