Robustez de um teste estatístico, testando o teste t

Vamos supor que eu tenho duas populações amostrais diferentes x_{1},x_{2},x_{3}, \dots ,x_{m} e y_{1},y_{2},y_{3}, \dots ,y_{n}.

A partir do teste t de student ou teste de diferença de médias eu posso perguntar se as médias são iguais.

H_{0}: \mu_x = \mu_y

Bem \mu_x é a média da distribuição normal que saíram os x_{1},x_{2},x_{3}, \dots ,x_{m}, se tivéssemos todos os infinitos valores de x, x_{1},x_{2},x_{3}, \dots ,x_{\infty}, a média seria \mu_x, mas como temos um número discreto de amostras, a gente denota que tem a média da amostra \bar{X}

Ta então a gente tem \bar{X} e \bar{Y} e seus respectivos desvios padrões  s_x e s_y

Com isso a gente pode fazer o teste t de student e verificar se as médias \bar{X} e \bar{Y} são iguais.

Realizamos a seguinte continha:

T = { { \bar{X}-\bar{Y} } \over{ S_{p}  \cdot  \sqrt{ {1 \over m }+ {1 \over n} } } }

E aquele S_{p} é o seguinte:

S_{p} = \sqrt{ (m-1) \cdot s_x^2 + (n-1) \cdot s_y^2 \over{m+n-2} }

E assim calculamos a estatística T, e sabemos que sobre a hipótese H_{0}: \mu_x = \mu_y a estatística T tem uma distribuição t com m+n-2 graus de liberdade. Claro que assumindo o seguinte:

Ambos x_{m} e y_{n} são amostras independentes e aleatorias de distribuições normais.
Os desvios padrões de x_{m} e y_{n}, os  \sigma_x e \sigma_y são iguais.

Certo então quando a gente faz um teste t, a gente ta assumindo essas duas coisas são verdades.

Ai o teste t vai consistir de uma comparação, a gente vai escolher um nível de significância \alpha que a gente quiser e vai comparar o valor que calculamos da estatística t

\vert T \vert \geq t_{n+m-2,{\alpha \over 2}}

onde t_{df,{\alpha \over 2}} é o quantile de  1 - {\alpha \over 2} de uma variável aleatória com df graus de liberdade. (df é de Degree of Freedom, grau de liberdade em inglês).

Bem se você leu até aqui parabéns, isso é meio chato né, parece muto abstrato, mas é exatamente assim que o teste t é feito. No R ou qualquer outro software. O teste t para diferenças de médias com variâncias iguais. Essa formula que mostramos ai em cima esta no wikipedia também aqui. Embaixo dela tem uma outra forma para não ter que assumir variâncias iguais, mas vamos testar esse caso aqui.

Vamos primeiro escrever na linguagem do R aquelas formulas que colocamos ali em cima.

Essa aqui:

S_{p} = \sqrt{ (m-1) \cdot s_x^2 + (n-1) \cdot s_y^2 \over{m+n-2} }

Que é o desvio padrão agrupado (pooled standard deviation) ou sei la como chama em português fica assim:

sp<-sqrt(((m-1)*sd(x)^2+(n-1)*sd(y)^2)/(m+n-2))

É so transcrevera formula para a linguagem do R, sqrt é a função que calcula raiz quadrado, ai como na formula ta tudo dentro dela, ai é uma coisa em cima e uma coisa em baixo, em cima temos m-1 onde m é o tamanho do vetor x, daria pra colocar length(x), mas salvamos o tamanho na variável m pra ficar mais curtinho e fácil de ver, ^2 eleva ao quadrado, menos pro y e o / faz a divisão. e pronto temos calculado toda aquela parte que multiplica em baixo a estatística T, agora é so fazer o mesmo esqueminha pro T inteiro.

Pegamos a formula:

T = { { \bar{X}-\bar{Y} } \over{ S_{p}  \cdot  \sqrt{ {1 \over m }+ {1 \over n} } } }

E passamos para a linguagem do R.

t<-(mean(x)-mean(y))/(sp*sqrt(1/m+1/n))

Essa foi mais fácil não? diminui as média e divide pelo sp multiplicando tudo aquilo la e pimba, temos a estatística t.

Então vamos fazer um exemplo com essa funções e ver se é isso mesmo que o R faz.

> set.seed(123456) > x<-rnorm(10,mean=50,sd=10) > y<-rnorm(10,mean=50,sd=10) > m<-length(x) > n<-length(y) > sp<-sqrt(((m-1)*sd(x)^2+(n-1)*sd(y)^2)/(m+n-2)) > t<-(mean(x)-mean(y))/(sp*sqrt(1/m+1/n)) > t [1] 1.189262

Então geramos 10 valores ao acaso, com média 50 e desvio 10, contamos o tamanho deles e salvamos em m e n e então usamos nossas formulinhas, e pimba, temos o valor t para esse caso, que é de 1.189262.

Mas para confirmar que estamos no caminho certo, o r tem um comando para fazer o teste t, o comando é t.test(), se usarmos ele e declararmos que as variações são iguais, que é o teste que estamos tentado realizar podemos conferir o valor t que ele vai calcular.

> t.test(x,y,var.equal =T,alternative = c("two.sided")) Two Sample t-test data: x and y t = 1.1893, df = 18, p-value = 0.2498 alternative hypothesis: true difference in means is not equal to 0 95 percent confidence interval: -4.132755 14.915122 sample estimates: mean of x mean of y 57.93401 52.54283

E olha la em cima, t=1.1893, ele apenas arredondou o valor para 4 casas depois da virgula, mas é exatamente o mesmo valor, então implementamos corretamente a formula no R.

Por curiosidade, quando a gente faz um modelo linear, usando lm por exemplo:

> summary(lm(c(x,y)~rep(c("a","b"),c(m,n)))) Call: lm(formula = c(x, y) ~ rep(c("a", "b"), c(m, n))) Residuals: Min 1Q Median 3Q Max -13.682 -10.094 0.407 7.069 17.092 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 57.934 3.205 18.073 5.49e-13 *** rep(c("a", "b"), c(m, n))b -5.391 4.533 -1.189 0.25 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 10.14 on 18 degrees of freedom Multiple R-squared: 0.07285, Adjusted R-squared: 0.02134 F-statistic: 1.414 on 1 and 18 DF, p-value: 0.2498

Notem que calculamos o valor t também, é a mesma conta, so que ele saiu negativo, porque nao esta apresentado em modulo, mas o valor é o mesmo, somente arredondado para três casas decimais.

Então a gente pode colocar todos aqueles passos dentro de uma função, assim:

estatistica.t<-function(x,y) {
    m<-length(x)
    n<-length(y)
    sp<-sqrt(((m-1)*sd(x)^2+(n-1)*sd(y)^2)/(m+n-2))
    t<-(mean(x)-mean(y))/(sp*sqrt(1/m+1/n))
    return(t)
    }

E temos nossa própria função de calcular a estatística t, podemos até dar uma testada, é so entrar com dois vetores cheio de números e temos a estatística t.

> dados.x<-c(1,4,3,6,5) > dados.y<-c(5,4,7,6,10) > estatistica.t(dados.x,dados.y) [1] -1.937926

Hehe, muito divertido. Certo agora a gente pode fazer uma pequena simulação para entender como o teste funciona. Basicamente, a gente vai estabelecer m como 10 e n como 10 também, ai vamos gerar 10 valores ao acaso para um vetor x e um y, preenchendo então 10 amostras em x_{1},x_{2},x_{3}, \dots ,x_{m} e y_{1},y_{2},y_{3}, \dots ,y_{n}
vamos pegar, gerar 10 valores de uma distribuição normal com média 0 e desvio padrão de um, depois vamos calcular o valor da estatistica t com nossa função para esse exemplo e fazer essa comparação:

\vert T \vert \geq t_{n+m-2,{\alpha \over 2}}

Ai vamos registrar quantas vezes essa comparação vai ser verdadeira. Na comparação com a distribuição t, para estabelecer o valor para comparar, vamos usar um valor de alpha de 0.05, ja que é uma convenção usar esse nível de significância, ou essa chance para assumir que H_{0}: \mu_x = \mu_y esta errada,é falsa. Aqui ta meio confuso, mas o que a gente faz é basicamente isso.

#Definir um alpha > alpha<-0.05 #Tamano das amostras > m<-10 > n<-10 #Número de simulações > N<-10000 #Uma variável para contar se há diferença ou não entre as médias > n.rejeitar<-0 #Simulando N vezes > for (i in 1:N) { #Criando as duas populações com média 0, desvio 1 e m ou n amostras + x<-rnorm(m,mean=0,sd=1) + y<-rnorm(n,mean=0,sd=1) #Calculando estatistica t + t<-estatistica.t(x,y) #Se T calculado é maior que o quantile da distribuição t, conto 1 rejeição + if (abs(t)>qt(1-alpha/2,n+m-2)) { + n.rejeitar<-n.rejeitar+1 + } #E começamos denovo, até 10000 vezes + } #Depois é so ver quantas vezes isso aconteceu dado o número de simulações > nivel.sig.verdadeiro<-n.rejeitar/N > nivel.sig.verdadeiro [1] 0.0484

Certo, então rejeitamos a hipótese que as populações são iguais 0.0484 vezes, ou 04.84% das vezes, praticamente o alpha de 0.05, ou de 5% ou 1/20. Então so existe 5% de chance de a gente falar que duas populações são diferentes no caso de elas serem iguais. Quando a gente ve isso nos nossos dados, o valor T calculado é muito grande, é mais fácil procurar um motivo para a diferença nas médias que achar que elas foram diferentes ao acaso, a gente so tem 5% de chance de ta fazendo cagada.

Agora vamos fazer uma pequena modificação, ao invés de descartar o valor T que calculamos a cada teste, vamos guardar eles e fazer um gráfico para ver como é a distribuição deles para compararmos com a distribuição t, que estamos usando para tomar essa decisão.

alpha<-0.05
m<-10
n<-10
N<-10000
n.rejeitar<-0
t.caso01<-vector()
 
for (i in 1:N) {
    x<-rnorm(m,mean=0,sd=1)
    y<-rnorm(n,mean=0,sd=1)
    t.caso01[i]<-estatistica.t(x,y)
    if (abs(t.caso01[i])>qt(1-alpha/2,n+m-2)) {
        n.rejeitar<-n.rejeitar+1
        }
    }

Veja que agora eu to criando um vetor vazio t.caso01, para guardar os resultados, e toda hora que eu calculo a estatística t, eu vou guardar numa posição desse vetor.
Ai fazemos a simulação, agora guardando as valores da estatística t e colocamos isso num gráfico.

01

E olha que legal, ficam praticamente idênticos, a estatística t que calculamos e a distribuição t. Isso porque essa é exatamente a distribuição t.

Agora até aqui, a gente esta respeitando todas as premissas, mas e se a gente quebrar alguma premissa, não respeitar elas.

Vamos colocar desvios muito diferentes e ver o que acontece:

> #Duas populações normais com média zero e desvios muito diferentes (σx=1, σy=10). > nivel.sig.verdadeiro<-n.rejeitar/N > nivel.sig.verdadeiro [1] 0.065

Humm, aumentamos um pouco o número de rejeições, mas bem pouco, quer dizer que não estaremos assumindo que existe diferença onde não existe de verdade com muito maior frequência, mas se fosse o juros da minha conta do banco, eu gostaria de 1% a mais, depende então do prejuízo que esse 1% pode causar, mas a gente começar a procurar algo a mais onde não existe. Mas variâncias diferentes pode ser contornado, so ver a próxima formula no wikipedia, na pagina de teste t.
Vamos fazer o gráfico também para ver o que esse 1% a mais dessa simulação significou em termos de diferença na distribuição da estatística t.

02

Bem, ainda está bem parecido, um pouquinho desviado se você desejar ser muito detalhista, mas pra mim parece ok.
Vamos testar outro caso, duas populações vinda de uma distribuição t com 4 graus de liberdade, ela vai ser um pouco mais aberta que a distribui normal.

> #Populações T com 4 graus de liberdade e desvios iguais > nivel.sig.verdadeiro<-n.rejeitar/N > nivel.sig.verdadeiro [1] 0.046

Agora rejeitamos menos vezes? É claro que cada vez que simulamos, esse número vai mudar um pouquinho, mas são dez mil simulações, a tendencia é sempre sair algo perto disso dae. Mas estamos na verdade considerando populações iguais mais que deveríamos. Olhando o gráfico vemos…

03

O pico ficou encolhido, sei la, mas ainda ficou bem bom, e acho que menos de 1% de diferença no nível de significância é um erro aceitável, podia ser pior.

Vamos testar agora se os dados viessem de outra distribuição, uma distribuição exponencial, o que aconteceria.

> #Populações exponenciais com μx = μy = 1. > nivel.sig.verdadeiro<-n.rejeitar/N > nivel.sig.verdadeiro [1] 0.0442

Um nível de significância um pouco menor denovo, menor até que o caso anterior, olhando amostras de uma distribuição t, mas 1% até que da para levar, se isso fosse o teste para ver a eficiência de um remédio, acho que daria para usar e levar a frente as pesquisas. Se vimos uma diferença quer dizer até que o remédio pode ser um pouquinho melhor, ja que a rejeição é menor que o alpha básico de 0.05 que achamos que estamos comparando, mas vamos olhar o gráfico nesse caso e ver como fica.

04

Novamente, a gente come um pedaço do cume da montainha que é a distribuição aqui, mas ainda ta parecido, mas não é a mesma coisa, mas esta parecido.

Agora vamos testar caso de cada população vem de uma distribuição diferente.

> #Uma população normal (μx = 10, σx = 2) e uma população exponencial (μy = 10). > nivel.sig.verdadeiro<-n.rejeitar/N > nivel.sig.verdadeiro [1] 0.1042

Rapaz, agora aqui ferrou tudo, aqui esculhambou geral a estatística, o nível de significância que deveria ta pertinho do 0.05, ta em 0.1042 na verdade. Isso quer dizer que coisas que são iguais, a gente vai ta falando que são diferente, a gente vai fica procurando diferença em um monte de coisa que não tem diferença, em que as médias são iguais, se a gente olhar o gráfico desse caso…

05

Aqui não tem nem o que discutir muito, as curvas estão muito diferentes. Mas muito mesmo, nem estão no mesmo lugar. Se a gente usar o teste t quebrando as premissas desse jeito vai dar problema, um grande problema.

Mas fazendo uma revisão então:

06

Da para ver que quando as pessoas falam que um teste é robusto, eu acho que é algo no sentido que dependendo de como você quebra as premissas, ele ainda tem uma boa confiabilidade, não vamos ter conclusões muito ruins. Mas dependendo de como quebramos, como visto no ultimo exemplo, nossas conclusões vão ser péssimas, os valores que estamos lidando não tem quase nada a ver com a distribuição que estamos usando para comparar, e ficaremos tentando ver efeito onde não existe.
É claro que existe muitas outras situações ficar pensando em todas as possibilidades para simular e ver que tipo de erro vai estar acontecendo não vale muito a pena, basta ter sempre em mente que, se você esta quebrando premissas, ou não esta levando elas em conta, pelo menos olhe com cuidado as estatísticas que esta calculado e pense bem nas conclusões que esta tomando baseado nessas estatísticas.

Bem, como sempre o script fica guardado no repositório Recologia além de aqui em baixo.

set.seed(123456)
x<-rnorm(10,mean=50,sd=10)
y<-rnorm(10,mean=50,sd=10)
 
m<-length(x)
n<-length(y)
 
sp<-sqrt(((m-1)*sd(x)^2+(n-1)*sd(y)^2)/(m+n-2))
sp
 
t<-(mean(x)-mean(y))/(sp*sqrt(1/m+1/n))
t
 
t.test(x,y,var.equal =T,alternative = c("two.sided"))
summary(lm(c(x,y)~rep(c("a","b"),c(m,n))))
 
estatistica.t<-function(x,y) {
    m<-length(x)
    n<-length(y)
    sp<-sqrt(((m-1)*sd(x)^2+(n-1)*sd(y)^2)/(m+n-2))
    t<-(mean(x)-mean(y))/(sp*sqrt(1/m+1/n))
    return(t)
    }
 
dados.x<-c(1,4,3,6,5)
dados.y<-c(5,4,7,6,10)
estatistica.t(dados.x,dados.y)
 
alpha<-0.05
m<-10
n<-10
N<-10000
n.rejeitar<-0
 
for (i in 1:N) {
    x<-rnorm(m,mean=0,sd=1)
    y<-rnorm(n,mean=0,sd=1)
    t<-estatistica.t(x,y)
    if (abs(t)>qt(1-alpha/2,n+m-2)) {
        n.rejeitar<-n.rejeitar+1
        }
    }
 
nivel.sig.verdadeiro<-n.rejeitar/N
nivel.sig.verdadeiro
 
 
 
#Duas populações normais com médias zero e desvios iguais (σx = σy = 1).
alpha<-0.05
m<-10
n<-10
N<-10000
n.rejeitar<-0
t.caso01<-vector()
 
for (i in 1:N) {
    x<-rnorm(m,mean=0,sd=1)
    y<-rnorm(n,mean=0,sd=1)
    t.caso01[i]<-estatistica.t(x,y)
    if (abs(t.caso01[i])>qt(1-alpha/2,n+m-2)) {
        n.rejeitar<-n.rejeitar+1
        }
    }
 
nivel.sig.verdadeiro<-n.rejeitar/N
nivel.sig.verdadeiro
 
plot(density(t.caso01),xlim=c(min(t.caso01),max(t.caso01)),ylim=c(0,0.5),frame=F,main="",lwd=2)
lines(seq(min(t.caso01),max(t.caso01),by=0.1),dt(seq(min(t.caso01),max(t.caso01),by=0.1),df=18),lty=3,col="red",lwd=3)
legend("topright",c("Simulação","t(18)"),lwd=c(2,3),lty=c(2,3),col=c("black","red"),bty="n")
 
#Duas populações normais com média zero e desvios muito diferentes (σx=1, σy=10).
alpha<-0.05
m<-10
n<-10
N<-10000
n.rejeitar<-0
t.caso02<-vector()
 
for (i in 1:N) {
    x=rnorm(m,mean=0,sd=1)
    y=rnorm(n,mean=0,sd=10)
    t.caso02[i]<-estatistica.t(x,y)
    if (abs(t.caso02[i])>qt(1-alpha/2,n+m-2)) {
        n.rejeitar<-n.rejeitar+1
        }
    }
 
nivel.sig.verdadeiro<-n.rejeitar/N
nivel.sig.verdadeiro
 
plot(density(t.caso02),xlim=c(min(t.caso02),max(t.caso02)),ylim=c(0,0.5),frame=F,main="",lwd=2)
lines(seq(min(t.caso02),max(t.caso02),by=0.1),dt(seq(min(t.caso02),max(t.caso02),by=0.1),df=18),lty=3,col="red",lwd=3)
legend("topright",c("Simulação","t(18)"),lwd=c(2,3),lty=c(2,3),col=c("black","red"),bty="n")
 
#Populações T com 4 graus de liberdade e desvios iguais
alpha<-0.05
m<-10
n<-10
N<-10000
n.rejeitar<-0
t.caso03<-vector()
 
for (i in 1:N) {
    x=rt(m,df=4)
    y=rt(n,df=4)
    t.caso03[i]<-estatistica.t(x,y)
    if (abs(t.caso03[i])>qt(1-alpha/2,n+m-2)) {
        n.rejeitar<-n.rejeitar+1
        }
    }
 
nivel.sig.verdadeiro<-n.rejeitar/N
nivel.sig.verdadeiro
 
plot(density(t.caso03),xlim=c(min(t.caso03),max(t.caso03)),ylim=c(0,0.5),frame=F,main="",lwd=2)
lines(seq(min(t.caso03),max(t.caso03),by=0.1),dt(seq(min(t.caso03),max(t.caso03),by=0.1),df=18),lty=3,col="red",lwd=3)
legend("topright",c("Simulação","t(18)"),lwd=c(2,3),lty=c(2,3),col=c("black","red"),bty="n")
 
#Populações exponenciais com μx = μy = 1.
alpha<-0.05
m<-10
n<-10
N<-10000
n.rejeitar<-0
t.caso04<-vector()
 
for (i in 1:N) {
    x=rexp(m,rate=1)
    y=rexp(n,rate=1)
    t.caso04[i]<-estatistica.t(x,y)
    if (abs(t.caso04[i])>qt(1-alpha/2,n+m-2)) {
        n.rejeitar<-n.rejeitar+1
        }
    }
 
nivel.sig.verdadeiro<-n.rejeitar/N
nivel.sig.verdadeiro
 
plot(density(t.caso04),xlim=c(min(t.caso04),max(t.caso04)),ylim=c(0,0.5),frame=F,main="",lwd=2)
lines(seq(min(t.caso04),max(t.caso04),by=0.1),dt(seq(min(t.caso04),max(t.caso04),by=0.1),df=18),lty=3,col="red",lwd=3)
legend("topright",c("Simulação","t(18)"),lwd=c(2,3),lty=c(2,3),col=c("black","red"),bty="n")
 
#Uma população normal (μx = 10, σx = 2) e uma população exponencial (μy = 10).
alpha<-0.05
m<-10
n<-10
N<-10000
n.rejeitar<-0
t.caso05<-vector()
 
for (i in 1:N) {
    x=rnorm(m,mean=10,sd=2)
    y=rexp(n,rate=1/10)
    t.caso05[i]<-estatistica.t(x,y)
    if (abs(t.caso05[i])>qt(1-alpha/2,n+m-2)) {
        n.rejeitar<-n.rejeitar+1
        }
    }
 
nivel.sig.verdadeiro<-n.rejeitar/N
nivel.sig.verdadeiro
 
plot(density(t.caso05),xlim=c(min(t.caso05),max(t.caso05)),ylim=c(0,0.5),frame=F,main="",lwd=2)
lines(seq(min(t.caso05),max(t.caso05),by=0.1),dt(seq(min(t.caso05),max(t.caso05),by=0.1),df=18),lty=3,col="red",lwd=3)
legend("topright",c("Simulação","t(18)"),lwd=c(2,3),lty=c(2,3),col=c("black","red"),bty="n")
 
#Graficão
par(mfrow=c(3,2))
 
plot(density(t.caso01),xlim=c(min(t.caso01),max(t.caso01)),ylim=c(0,0.5),frame=F,
     main="Duas populações normais com médias zero e desvios iguais",lwd=2)
lines(seq(min(t.caso01),max(t.caso01),by=0.1),dt(seq(min(t.caso01),max(t.caso01),by=0.1),df=18),lty=3,col="red",lwd=3)
legend("topright",c("Simulação","t(18)"),lwd=c(2,3),lty=c(2,3),col=c("black","red"),bty="n")
 
plot(density(t.caso02),xlim=c(min(t.caso02),max(t.caso02)),ylim=c(0,0.5),frame=F,
     main="Duas populações normais com média zero e desvios muito diferentes",lwd=2)
lines(seq(min(t.caso02),max(t.caso02),by=0.1),dt(seq(min(t.caso02),max(t.caso02),by=0.1),df=18),lty=3,col="red",lwd=3)
legend("topright",c("Simulação","t(18)"),lwd=c(2,3),lty=c(2,3),col=c("black","red"),bty="n")
 
plot(density(t.caso03),xlim=c(min(t.caso03),max(t.caso03)),ylim=c(0,0.5),frame=F,
     main="Populações T com 4 graus de liberdade e desvios iguais",lwd=2)
lines(seq(min(t.caso03),max(t.caso03),by=0.1),dt(seq(min(t.caso03),max(t.caso03),by=0.1),df=18),lty=3,col="red",lwd=3)
legend("topright",c("Simulação","t(18)"),lwd=c(2,3),lty=c(2,3),col=c("black","red"),bty="n")
 
plot(density(t.caso04),xlim=c(min(t.caso04),max(t.caso04)),ylim=c(0,0.5),frame=F,
     main="Populações exponenciais",lwd=2)
lines(seq(min(t.caso04),max(t.caso04),by=0.1),dt(seq(min(t.caso04),max(t.caso04),by=0.1),df=18),lty=3,col="red",lwd=3)
legend("topright",c("Simulação","t(18)"),lwd=c(2,3),lty=c(2,3),col=c("black","red"),bty="n")
 
plot(density(t.caso05),xlim=c(min(t.caso05),max(t.caso05)),ylim=c(0,0.5),frame=F,
     main="Uma população normal e uma população exponencial",lwd=2)
lines(seq(min(t.caso05),max(t.caso05),by=0.1),dt(seq(min(t.caso05),max(t.caso05),by=0.1),df=18),lty=3,col="red",lwd=3)
legend("topright",c("Simulação","t(18)"),lwd=c(2,3),lty=c(2,3),col=c("black","red"),bty="n")

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)

Ancova Binomial Bayesiana avaliando Prevalências

Esses nomes de posts são de certa forma uma sacanagem. Mas normalmente é sempre assim que a gente vê as coisas.
Análises normalmente estão “empacotadas” em nomes, teste t, diferença de médias, com mesma variação, sem mesma variação, anova, ancova, anova de medidas repetidas , regressão linear e etc.
Eu entendo que capítulos de livros ou posts precisam de títulos, índices, mas a gente tem que tomar sempre cuidado para não focar demais num nome e não olhar o todo, se você focar apenas em um ponto você pode perder todo o esplendor do por do Sol, como diria Bruce Lee em Operação Dragão (Isso que é citação heim).

Após o blablabla inicial, vamos voltar ao GLM binomial que estávamos fazendo aqui.

Certo, no caso anterior estávamos contando e vendo a proporção, ou a palavra mais usada para parasitismos que é prevalência, como ela muda entre duas espécies.

Apesar de estarmos simulando dados, vamos pensar em algo legal. Um exemplo que eu acho legal é o Pardal ou Passer domesticus, esse que infesta as cidades e a gente vê um bocado. A distribuição geográfica era originalmente, a maior parte da Europa e da Ásia, mas ele foi introduzido na América, África, Austrália e Nova Zelândia, e olha que esses dois últimos são ilhas e eles não chegaram voando eu acho. Mas por algum motivo sua população tem decaído na europa. Inclusive o jornal The Independent oferece 5 mil euros ou libras, não sei reconhecer o simbolo do dinheiro que eles usam, mas então aqui tem a reportagem que eles abrem essa oportunidade. Ja houveram três entradas sérias para ganhar o prêmio, onde alguns caras publicaram artigos em revistas científicas e colocaram suas teorias para explicar o que esta acontecendo. Mas até agora nada. Bem apesar de a gente não ligar muito pro Pardal, são 5 mil euros no bolso, pode ser interessante tentar entender o que ta acontecendo. Mas como a gente estava pensando em parasitismos e prevalências, vamos continuar pensando em exemplos assim, uma ideia interessante é ver se o mesmo parasita observado anteriormente, o Plasmodium que causa a Malaria Aviaria não pode ter um dedo nesse caso. Como assim?

Uma teoria interessante é essa:

immunityinvasion

Espécies tem sucesso ao invadir novos ambientes, porque elas trazem seus parasitinhas, e ja vivem bem com eles, mas ai esses parasitas podem atacar novas espécies ao chegarem, ai os parasitas rebentam todo mundo do local enquanto a espécie nova procria e aumenta sua população tranquila, lembra portugueses chegando ao Brasil e os índios morrendo não? Mas essa é a idea.

Acontece que os pardais estão se ferrando no seu local de origem. Tudo é bem complicado de se encaixar. Mas suponha que tenhamos os seguintes dados:

> dados status local peso 1 1 Mata Atlantica 79 2 1 Mata Atlantica 25 3 1 Mata Atlantica 68 . . . 277 1 Mata Atlantica 23 278 1 Mata Atlantica 19 279 1 Amazonia 62 280 1 Amazonia 28 281 1 Amazonia 38 . . . 496 1 Amazonia 20 497 1 Amazonia 50 498 1 Amazonia 15 499 1 Cerrado 50 500 1 Cerrado 34 . . . 726 1 Cerrado 53 727 1 Cerrado 26 728 1 Cerrado 78

A gente tem pego alguns pardais aqui no Brasil, cada linha é um passarinho diferente, e vendo toda essa discussão, a gente resolveu dar uma olhadinha aqui, será que a prevalência dos pardais varia aqui dentro do Brasil, o que a gente ve no cerrado é igual na amazônia ou na mata atlântica?

Então temos três locais, e nesses locais a gente foi colocando redes pra capturar passarinhos, ai pegávamos eles, tirávamos sangue de um por um, pesavamos e soltavamos, depois a gente pega o sangue, joga num frasco, poe um primer que amplifica a sequência do plasmodium, faz uma pcr e olha se amplificou, se amplificou ta doente, status é 1, ta doente, tem o parasita, caso contrario é zero.

Esses são os dados resultantes desse processo todo. A primeira coisa a notar, é que não temos muito como controlar o número de amostras, o número de passarinhos em cada caso, claro que precisamos de um minimo, se pegar apenas um ou nenhum pardal num da para querer ver muita coisa, mas as populações variam, a chance de captura depende de muitas coisas, então o número de pardais que teremos vai variar entre locais, não tem muito jeito, daria para homogeniezar isso jogando fora amostras, mas ai nos teríamos que se basear no pior caso, no local que temos menos amostras para ver quantas amostras usar, e convenhamos, faz pouco sentido ou nenhum jogar informação fora. Mas para nossa satisfação é possível usar a distribuição binomial para descrever a prevalência, ou chance e não teremos problemas com a variação do número de passarinhos por amostra.

Agora aqui temos uma diferença forte em relação ao que estávamos fazendo no post anterior. Cada individuo tem uma característica somente sua, o peso, então não da para simplesmente ver a contagem total de passarinhos com o parasita e sem de acordo com o local.

Bem antes de começar a pensar em como fazer o modelo, vamos olhar melhor os dados.

01

Primeiro vamos olhar quantos passarinhos estão parasitados e não parasitados de forma geral, usando todos os passarinhos.
Aqui parece que temos meio a meio, mais ou menos.

Mas vamos olhar se isso continua do mesmo jeito por local.

02

Nossa, agora virou uma zona, na Amazônia parece que é meio a meio, agora no cerrado parece que tem muito mais indivíduos parasitados que não parasitados, e pra piorar tudo parece que a mata atlântica apresenta um padrão oposto.
Mas ainda tem o peso dos passarinhos na jogada.

Colocando num gráfico vemos

03

Bem como so tem 0 ou 1, os pontinhos ficam no teto do gráfico ou no chão. Mas a gente pode fazer uma retinha chamada Loess, que segue a seguinte ideia, pensa que a gente pega os 5 primeiro pontinhos, sendo que eles estão organizados do menor pro maior peso, que podem ser 0 ou 1, somamos e dividimos por 5, ai temos um pontinho, uma prevalência desses 5 pontinhos, ai a gente faz isso com os próximos 5 pontinhos e fazemos mais um pontinho, e depois denovo e denovo, e juntamos uma linha nesses pontinhos, essa é a linha vermelha, que nesse caso ajuda a visualizar como a prevalência esta mudando de acordo com o peso. E vemos uma mudança nos três locais, mais mais forte no cerrado aparentemente, menos nos outros dois casos, mas mudança de acordo com o peso. Não é difícil imaginar motivos que ocasionariam tal tendência, maior peso pode significar indivíduos mais velhos, então que ficaram mais tempo expostos ao ambiente com o parasita, existem ainda a possibilidade do sistema imuno piorar com a idade, muitas possíveis causas, mas antes de cogitar isso, temos que verificar que isso esta acontecendo, e parece que é o caso.

Então parece que podemos tentar pensar da seguinte forma:

04

O status, o estado relativo a presença do parasita depende do peso e do local, mas essas variáveis interagem, uma influência a outra, como sugerido no segundo e terceiro gráfico, que mostra padrões diferentes de acordo com o local e vice versa.

Então fazemos nosso modelinho para mandar para o OpenBugs.

model {
# Priors
 for (i in 1:n.local) {
    alpha[i] ~ dnorm(0, 0.01)  # Interceptos (Um pra cada local)
    beta[i] ~ dnorm(0, 0.01)   # Inclinações (Um pra cada local)
 }
 
# Likelihood
 for (i in 1:n) {
    status[i] ~ dbin(p[i], 1)
    logit(p[i]) <- alpha[local[i]] + beta[local[i]]* peso[i]
 }

Olha, primeiro estabelecemos nossos prior, vamos assumir nenhum efeito e zero parasitismos com esses prior. Veja que essas distribuições do jeito que estão são um pico sobre o zero.
Outra coisa aqui é que precisamos de três equações da reta, uma para cada local, então temos três interceptos e três inclinações, três a+bx.

Agora o likelihood fica assim, o status vem de uma distribuição binomial, dai o uso do status[i] ~ dbin(p[i], 1), como é 0 ou 1, temos o dbin com a probabilidade do indivíduo, e para 1 indivíduo, os dois parâmetros dentro do dbin.
Dai a segunda linha é a parte determinística do modelo, que é a+bx, aqui alpha[local[i]] + beta[local[i]]* peso[i], agora o alpha e o beta, o a e o b dependem do local que o passarinho esta e do seu peso. E acabou o modelo. Como estamos trabalhando com a distribuição binomial, ou seja um GLM, temos que transformar essa probabilidade com a função logística, veja que tem uma função pronta no bugs pra isso para facilitar a vida, isso é a parte escrita logit(p[i]).

Ai antes de fechar a chave do modelo calculamos mais alguns parametros de interesse.

# Recuperando os efeitos baseados no caso base "Amazonia"
 a.effe2 <- alpha[2] - alpha[1]		# Intercepto Amazonia vs Cerrado
 a.effe3 <- alpha[3] - alpha[1]		# Intercepto Amazonia vs Mata Atlantica
 b.effe2 <- beta[2] - beta[1]		# Inclinação Amazonia vs Cerrado
 b.effe3 <- beta[3] - beta[1]		# Inclinação Amazonia vs Mata Atlantica
 
# Teste Personalizado (Respondendo a minha pergunta)
 test1 <- alpha[2] - alpha[1]		# Diferença entre Cerrado a Amazonia
 test2 <- alpha[2] - alpha[3]		# Diferença entre Cerrado e Mata Atlantica
}

Antes de fechar a chave do modelo, recuperamos os efeitos baseados num intercepto geral, para comparar com os valores usando estatística frequentista, usando o comando glm do R e a última parte eu faço meus testes que tenho interesse, que é ver se a prevalência muda do cerrado para a Amazônia e Mata Atlântica, afinal antes de querer saber porque, eu tenho que saber que ha mudança, se não houver mudança a pergunta porque não faz sentido, que é o meu interesse a principio. Isso seria algo como comparações entre grupos, mas eu não vou sair comparando tudo com tudo que nem louco, eu quero saber isso a principio e não tentar todas as comparações possíveis.

Pronto, agora juntamos os dados numa lista.

bugs.data <- list(status = dados$status, local = as.numeric(dados$local), n.local = 3, n = nrow(dados),peso=dados$peso)

Lembrando que mandamos os locais como números, então temos que lembrar depois quem é o 1, quem é 2 e quem é três, como isso vem de um fator do R, a não ser que declaremos uma ordem, ela vai ser a alfabética ou lexografica, então é Amazônia, que começa com A é 1, Cerrado é 2 e Mata atlântica é 3. Simples.

Outro detalhe é que o peso foi Uniformizado, meramente isso é diminuir cada valor da média e dividir pelo desvio padrão, isso nos deixa com os valores variando em unidades de desvio padrão e com média zero. As tendencias vão continuar as mesmas com ou sem essa transformação, mas podem haver problemas se o bugs tiver que lidar com números muitos grandes, por isso ta assim, talvez não haja diferença aqui, ja que os pesos so tem duas casas decimais, mas transformamos de qualquer forma.

Fazemos uma função para gerar valores para o inicio da cadeia ao acaso:

> # Função para iniciar a cadeia > inits <- function(){ list(alpha = rlnorm(3, 3, 1), beta = rlnorm(3, 2, 1))} # Note log-normal inits > #Ela so gera valores ao acaso para iniciar as cadeias, veja so > inits() $alpha [1] 6.996297 22.373664 6.808293 $beta [1] 12.366757 6.475181 5.738183

Olha o que isso faz, gera 3 interceptos e 3 inclinações para iniciar a cadeia, esses valores tem que ser razoáveis, mas como vamos queimar o início da cadeia não são muito relevantes, mas como no caso em que calculamos o pi, o começo sempre é ruim por isso é bom descartar, mas temos que começar de algum lugar, senão não da.

Depois disso falamos também quais parâmetros vamos guardar e configuramos o MCMC.

Mandamos tudo do R para o OpenBugs e recebemos de volta resultado.

> print(out,digits=3) Inference for Bugs model at "ancovabin.txt", Current: 3 chains, each with 1500 iterations (first 500 discarded), n.thin = 5 Cumulative: n.sims = 3000 iterations saved mean sd 2.5% 25% 50% 75% 97.5% Rhat n.eff alpha[1] -0.054 0.136 -0.325 -0.141 -0.051 0.036 0.216 1.001 3000 alpha[2] 2.205 0.274 1.731 2.011 2.194 2.372 2.789 1.001 3000 alpha[3] -0.521 0.128 -0.765 -0.609 -0.522 -0.435 -0.272 1.001 3000 beta[1] 0.286 0.137 0.017 0.191 0.283 0.378 0.553 1.003 870 beta[2] 1.526 0.288 0.982 1.330 1.519 1.712 2.130 1.001 3000 beta[3] 0.392 0.126 0.146 0.304 0.392 0.478 0.638 1.001 3000 a.effe2 2.258 0.303 1.695 2.053 2.242 2.456 2.878 1.001 3000 a.effe3 -0.468 0.187 -0.845 -0.595 -0.466 -0.340 -0.101 1.001 3000 b.effe2 1.240 0.321 0.637 1.024 1.234 1.452 1.888 1.003 3000 b.effe3 0.107 0.187 -0.258 -0.018 0.108 0.234 0.469 1.003 860 test1 2.258 0.303 1.695 2.053 2.242 2.456 2.878 1.001 3000 test2 2.726 0.302 2.164 2.523 2.713 2.913 3.363 1.001 2400 deviance 817.706 3.492 812.900 815.100 817.100 819.600 826.202 1.001 3000 For each parameter, n.eff is a crude measure of effective sample size, and Rhat is the potential scale reduction factor (at convergence, Rhat=1). DIC info (using the rule, pD = Dbar-Dhat) pD = 6.0 and DIC = 823.7 DIC is an estimate of expected predictive error (lower deviance is better).

E temos ai nosso resultado, o intercepto e inclinação de cada caso, cada local, parâmetros para desenhar uma função logística para cada quadro da figura 3, so pensar em uma função logística desenhada em cada quadradinho ali.
Outra coisa é que recuperamos os coeficientes baseados no intercepto geral, o caso base.

Se verificarmos, para ter certeza que a análise funciona, vemos que recuperamos os valores originais usados para fazer a simulação dos dados.

> beta.vec [1] 0.0 2.0 -0.5 0.5 1.0 0.0

alpha[1]=-0.054 é o intercepto básico, o primeiro valor ali em cima, e a.effe2=2.258 e a.effe3=-0.468 que são as diferença de intercepto, estimativas muito boas, não estamos concluindo nada errado até aqui.
Dai temos as inclinações beta[1]=0.286 e as interações b.effe2=1.240 e b.effe3=0.107, que são razoáveis, apesar de termos subestimado um pouco talvez a inclinação base, ja que a média deu 0.286, mas o valor real usado na simulação de 0.5 está dentro do intervalo de confiança de 95% que vai de 0.017 até 0.553, so olhar ali o beta[1].

Usando o comando glm do R podemos ainda fazer a mesma analise num framework frequentista e temos:

> summary(glm(status ~ local * peso, family = binomial,data=dados)) Call: glm(formula = status ~ local * peso, family = binomial, data = dados) Deviance Residuals: Min 1Q Median 3Q Max -3.0142 -1.0129 0.2111 1.0985 1.7024 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) -0.04842 0.13629 -0.355 0.722375 localCerrado 2.21404 0.30130 7.348 2.01e-13 *** localMata Atlantica -0.46942 0.18620 -2.521 0.011699 * peso 0.27947 0.13866 2.015 0.043858 * localCerrado:peso 1.20495 0.31480 3.828 0.000129 *** localMata Atlantica:peso 0.11294 0.18557 0.609 0.542798 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 (Dispersion parameter for binomial family taken to be 1) Null deviance: 999.97 on 727 degrees of freedom Residual deviance: 811.67 on 722 degrees of freedom AIC: 823.67 Number of Fisher Scoring iterations: 6

Estimativas muito semelhantes dos parâmetros, aqui teríamos que fazer as contas para extrair os parâmetros de cada local, mas como a gente viu no modelo bayesiano, são contas de mais e menos, então não é muito difícil. E temos que chegar a mesma conclusão de qualquer forma, senão algo estaria errado.

Mas após verificar os resíduos e validar o modelo, damos uma olhada no resultado dos nossos testes específicos.

05

E vemos uma diferença significativa, ou seja é pouco provável que não haja diferença, veja como a distribuição das diferenças esta muito longe do zero, ou seja podemos afirmar que existe uma diferença entre a prevalência do cerrado e os outros locais. Agora porque? Porque será que isso deveria acontecer? Eu agora começaria a sequenciar os Plasmodiuns e veria se eles são muito diferentes entre os locais, isso poderia ajudar a pensar se o que esta acontecendo tem haver com o ambiente, ou com a adaptação do parasita as espécies, algo assim, mas são dados simulados, a gente só queria fazer a análise pra ver como era 🙂

Bem é isso, analise de covariância binomial, covariância porque tem fatores contínuos (peso) e fatores discretos (locais). Binomial que os dados que tinha era contagem de 0 e 1, e é um caso de modelo geral linar (glm) ja que eu lianerizo uma relação que não é linear com a função logística, pra estimar as coisas como se fossem uma reta, usando a distribuição binomial, então glm binomial.

O script inteiro vai estar no repositório Recologia no Github e aqui em baixo. Até a próxima.

###############################################################
#    ANCOVA Binomial                                          #
###############################################################
#Gerando dados para um exemplo
set.seed(51)
quantidade<-round(runif(3, 200, 300))
peso<-round(runif(sum(quantidade),10,80))
hist(peso)
peso<-(peso-mean(peso))/sd(peso)
hist(peso)
 
#Preditores
dados<-data.frame(local=factor(rep(c("Mata Atlantica", "Amazonia", "Cerrado"),quantidade)),peso)
dados
 
Xmat <- model.matrix(~ local*peso,data=dados)   #Matriz de Efeitos
print(Xmat, dig = 2)
beta.vec <- c(0, 2, -0.5, 0.5, 1, 0)            #Valor dos parametros
 
lin.pred <- Xmat[,] %*% beta.vec	        # Valor do preditor Linear
exp.p <- exp(lin.pred) / (1 + exp(lin.pred))    # Chance esperada
exp.p
 
#Resposta
status <- rbinom(n = nrow(dados), size = 1, prob = exp.p) # Adicionando erro binomial
status				                          # Checando o que foi criado
 
#Juntando tudo num data.frame
dados<-data.frame(status,dados)
 
##################################################################
 
#Olhando os dados
dados
 
#Grafico dos Diagnósticos
#jpeg("01.jpg")
barplot(table(dados$status),ylim=c(0,500),xlab="Diagnóstico",ylab="Número de passarinhos")
#dev.off()
#Grafico dos diagnósticos por locais
#jpeg("02.jpg")
barplot(table(dados$status,dados$local),ylim=c(0,300),xlab="Diagnóstico",ylab="Número de passarinhos",beside=T,
        legend.text = c("Não parasitado","Parasitado"))
#dev.off()
#Grafico dos diagnósticos por locais e peso.
library(lattice)
 
panel.smoother <- function(x, y) {
  panel.xyplot(x, y,pch=19,col="black")             # Plota os pontos
  panel.loess(x, y,span=0.6,lwd=2,lty=2,col="red")  # Mostra a linha que muda de acordo com a média
}
 
#jpeg("03.jpg")
xyplot(status~peso|local,data=dados,panel=panel.smoother,layout = c(1, 3))
#dev.off()
 
#Grafo do modelo
library(igraph)
#jpeg("04.jpg")
plot(graph.formula( Local -+ Status , Peso -+ Status , Local ++ Peso ),vertex.size=25,edge.color="black",
     vertex.color="white")
#dev.off()
 
### Analise usando OpenBugs
# Definindo um modelo
sink("ancovabin.txt")
cat("
model {
 
# Priors
 for (i in 1:n.local) {
    alpha[i] ~ dnorm(0, 0.01)  # Interceptos (Um pra cada local)
    beta[i] ~ dnorm(0, 0.01)   # Inclinações (Um pra cada local)
 }
 
# Likelihood
 for (i in 1:n) {
    status[i] ~ dbin(p[i], 1)
    logit(p[i]) <- alpha[local[i]] + beta[local[i]]* peso[i]
 }
 
# Recuperando os efeitos baseados no caso base "Amazonia"
 a.effe2 <- alpha[2] - alpha[1]		# Intercepto Amazonia vs Cerrado
 a.effe3 <- alpha[3] - alpha[1]		# Intercept Amazonia vs Mata Atlantica
 b.effe2 <- beta[2] - beta[1]		# Inclinação Amazonia vs Cerrado
 b.effe3 <- beta[3] - beta[1]		# Inclinação Amazonia vs Mata Atlantica
 
# Teste Personalizado (Respondendo a minha pergunta)
 test1 <- alpha[2] - alpha[1]		# Diferença entre Cerrado a Amazonia
 test2 <- alpha[2] - alpha[3]		# Diferença entre Cerrado e Mata Atlantica
}
",fill=TRUE)
sink()
 
# Juntando os dados numa lista (Veja que o local a gente manda como um número apenas,
#vc tem que lembrar que número é o que)
bugs.data <- list(status = dados$status, local = as.numeric(dados$local), n.local = 3, n = nrow(dados),peso=dados$peso)
 
# Função para iniciar a cadeia
inits <- function(){ list(alpha = rlnorm(3, 3, 1), beta = rlnorm(3, 2, 1))} # Note log-normal inits
#Ela so gera valores ao acaso para iniciar as cadeias, veja so
inits()
 
# Parametros que queremos salvar os resultados
params <- c("alpha", "beta","a.effe2","a.effe3","b.effe2","b.effe3","test1","test2")
 
# Configurações do MCMC
ni <- 1500
nb <- 500
nt <- 5
nc <- 3
 
# Iniciando o Gibbs Sampler
library(R2OpenBUGS)
out <-bugs(data = bugs.data, inits = inits, parameters.to.save = params, model.file = "ancovabin.txt",
           n.thin = nt, n.chains = nc, n.burnin = nb, n.iter = ni)
 
print(out,digits=3)
beta.vec
summary(glm(status ~ local * peso, family = binomial,data=dados))
 
str(out)
 
#jpeg("05.jpg")
par(mfrow=c(2,1))
hist(out$sims.list$test1,xlim=c(0,max(out$sims.list$test2)),main="Diferença entre Cerrado a Amazonia",
     xlab="Tamanho da Diferença")
lines(c(0,0),c(0,400),lty=2,lwd=2,col="red")
 
hist(out$sims.list$test2,xlim=c(0,max(out$sims.list$test2)),main="Diferença entre Cerrado e Mata Atlantica",
     xlab="Tamanho da Diferença")
lines(c(0,0),c(0,400),lty=2,lwd=2,col="red")
#dev.off()

Estimando o valor de pi usando o método de Monte Carlo

Ae, denovo com algum tempo livre para mais alguns posts 🙂

Eu já fiz alguns posts sobre Markov Chain, para pelo menos termos uma ideia do que se tratava, mas tudo que a gente faz em estatística Bayesiana é sempre usando o Monte Carlo Markov Chain (MCMC), a gente ja deve ter alguma noção do que se trata o Markov Chain, mas e o tal do Monte Carlo?

Bem eu ja comentei em algum post que não me lembro, que Monte Carlo é um bairro famoso de Mônaco por seus cassinos, e cassino é o lugar da probabilidade. Basicamente, o método de Monte Carlo é um algorítimo para se obter soluções numéricas de optimização, integral e amostrar distribuições probabilísticas.

O método de Monte Carlo mais ou menos é o seguinte:

1. Defina um dominio de todas as possibilidades 2. Gere aleatoriamente entradas de uma distribuição de probabilidade que cobre o dominio. 3. Faça uma computação deterministica das entradas 4. Reuna os resultados

Certo, num entendi porra nenhuma. Eu também não, mas é tentando com algo fácil que sabemos o resultado que a gente aprende.

Então vamos calcular a area de um circulo:

Então desenhamos um circulo no R

01

Certo, podemos ver que o raio desse circulo é 1. E sabemos que a área do círculo é

 Area = \pi \cdot r^{2}

Então a área desse exemplo tem que ser:

> #Area do circulo > pi*1^2 [1] 3.141593

Certo, mas como aplicar aquela ideia do Monte Carlo?
Bem o primeiro passo é definir um domínio de todas as possibilidades, vamos supor que a gente não soube-se o valor de \pi(pi), a gente tava ferrado porque não tinha como calcular a área, mas uma figura que é mais simples de deduzir qual a área é o quadrado. Então a gente desenha um quadrado em volta desse círculo.

02

Sabemos aqui que o lado desse quadrado é 2 e a área do quadrado é:

 Area = lado^{2}

Então a área do quadrado é:

> #Area do Quadrado > 2^2 [1] 4

Pode parecer idiota, mas se o nosso objetivo aqui é saber a área do circulo dentro do quadrado, podemos dizer que ele contem toda a área do circulo e mais um pouquinho.

O primeiro passo ta pronto. Definimos algo que contém todas as possibilidades.

Agora vamos para o passo dois, gere entradas de uma distribuição que cobre todo o domínio.

Vamos la, pense que se quisermos colocar um pontinho nesse quadrado num lugar qualquer, se sortearmos um ponto no eixo x ao acaso de uma distribuição uniforme dos valores -1 até o valor 1, e o mesmo para o eixo y, um ponto ao acaso de uma distribuição uniforme dos valores -1 até o valor 1, teremos um ponto x e y , que com certeza vai cair dentro desse quadrado.

Podemos simular isso no r. Então sorteamos um ponto que são duas coordenadas, x e y, de uma distribuição uniforme que vai de -1 a 1.

Então basicamente estamos fazendo é isso.

03

Bang, jogamos um pontinho ali dentro, e ele caiu dentro do circulo legal. Se funcionou, vamos gerar uns 500 pontos desses dai e ver o que vai acontecer

04

Olha que legal, os pontinhos vermelhos cairam dentro do círculo e os azuis cairam fora do circulo. Tem muito mais pontos vermelhos que azuis certo, isso porque a área do circulo é grande, alias, ela é uma proporção da área do quadrado. Proporção é a palavra, se o quadrado é o total, então a parte que queremos dividido pelo total é uma proporção que nos permite saber da area total quanto é área do circulo, ou seja, essa proporção vezes a área do quadrado é a área do circulo. Mas como a gente vê essa proporção a partir desse monte de ponto sorteados ao acaso, vamos fazer uma simulação maior, cinco mil pontos, então o número de pontos que ficaram dentro do circulo dividido pelo total de pontos que simulamos no quadrado é a proporção da área entre os dois.

Então fazemos a simulação da seguinte forma:

#Fazendo essa simulação 5 mil vezes e olhando a proporção de pontos dentro do circulo
n<-5000
x.pos <- runif(n, min=-1, max=1)
y.pos <- runif(n, min=-1, max=1)
local.pos <- ifelse(x.pos^2 + y.pos^2 <= 1, TRUE, FALSE)
dentro <- length(which(local.pos == TRUE))

Escolhemos um tamanho n para a simulação, cinco mil, dai sorteamos 5 mil pontos para x e para y, veja que é exatamente o que tamo fazendo no gráfico ali visualmente, só que 5 mil pontos ia deixar poluído o gráfico.

Dai o passo três é fazermos uma computação determinística, olhamos se caiu dentro ou fora do circulo. Se esta difícil pegar o que esta quarta linha esta fazendo, so lembrar pra que serve o famoso  c^2 = a^2 + b^2 , conhecido como teorema de pitágoras, mas eu também apanhei para entender essa parte, mas aqui tem um tutorialzinho legal para entender como fazer um gráfico de um circulo.
Por ultimo, agregamos os resultados, vemos a razão de pontos dentro do circulo do total de pontos que simulamos.

Temos nesse caso que:

> #Propoção > dentro/n [1] 0.786

Mas a gente sabe qual deve ser essa proporção, com infinitos pontos, a area do circulo dividido pela area do quadrado que é o total.
E para o estamos, essa proporção da.

> #Exato > (pi*1^2) / (2^2) [1] 0.7853982

Puxa vida, como diabos esse treco funciona heheheh, é impressionante. Se antes tava complicado de entender o algorítimo de Metropolis-hasting, olha ai como é bem parecido, ficamos gerando alguma coisa, temos um função pra aceitar ou não essas coisas que geramos e depois vemos uma distribuição.

Agora se multiplicarmos essa proporção pela área do quadrado, que no exemplo é 4, redescobriremos a área do circulo, que é o mesmo valor de \pi nesse caso, estimamos o valor de \pi, que legal eeeee.

> 4*(dentro/n) [1] 3.144

Bem mas podemos ainda ver o seguinte, com 1 pontinho apenas, a estimativa seria péssima, porque seria ou 1/1 ou 0/1, que daria 1 ou 0, que multiplicando por 4, daria ou 4 ou 0 a nossa estimativa de pi, mas com 2 pontos, ja temos 5 possibilidades 4/4, 3/4, 2/4, 1/4 e 0/4, note que quanto mais pontos sorteamos, mais possibilidades de valores, mas o que vai acontecendo de acordo com que sorteamos mais pontos, de acordo com que iteramos mais nessa ideia, nesse algorítimo?

pi

Olha ai, no começo o valor oscila muito, e vai de ruim a pior, mas conforme usamos mais pontos, mais iterações, a estimativa fica melhor e melhor, e começa a variar cada vez menos, ficamos mais e mais próximos do valor de pi, mas vamos chegar no valor de pi?
Ai eu ja não entendo, pi não é um número real, então num tem como saber exatamente, poxa mas a estimativa é muito boa. O valor de pi mesmo que temos no R tem um número finito de casas decimais, enquanto pi realmente tem um número infinito de casas decimais. Mas isso é uma discussão complexa e recorrente em fóruns, mas tem mais haver com computação que matemática nesse caso, mas tem um comentário interessante no faq do R aqui.

Mas beleza, agora a gente ja tem uma ideia melhor do que Monte Carlo, do que é Markov Chain, então da para tentar progredir pro MCMC talvez entendo melhor o que ta acontecendo.

Bem segue o script, ele esta no repositório do github também e até o próximo post 🙂

set.seed(171)
t <- seq(0.2*pi,length=1000)
coords <- t(rbind(0+sin(t)*1,0+cos(t)*1))
 
#desenhando um circulo
#jpeg("01.jpg")
plot(0,0,type="n",xlim=c(-1.5,1.5),ylim=c(-1.5,1.5),frame=F,xlab="Eixo x",ylab="Eixo y")
points(coords,cex=0.1)
abline(h=c(0),v=c(0),lty=3,lwd=0.5)
lines(c(0.01,0.99),c(0.01,0.01),lwd=2,lty=2,col="red")
text(0.5,0.1,"Raio=1")
#dev.off()
 
#Desenhando um quadrado sobre o circulo
#jpeg("02.jpg")
plot(0,0,type="n",xlim=c(-1.5,1.5),ylim=c(-1.5,1.5),frame=F,xlab="Eixo x",ylab="Eixo y")
points(coords,cex=0.1)
abline(h=c(0),v=c(0),lty=3,lwd=0.5)
lines(c(-1,-1,-1,1,1,1,1,-1),c(-1,1,1,1,1,-1,-1,-1),lwd=2,col="darkgray")
lines(c(-1.01,1.01),c(-1.01,-1.01),lwd=2,lty=2,col="red")
text(0,-1.08,"Lado=2")
#dev.off()
 
#Area do circulo
pi*1^2
 
#Area do Quadrado
2^2
 
#proporção
(pi*1^2) / (2^2)
 
#Jogando um ponto ao acaso dentro do quadrado
#jpeg("03.jpg")
plot(0,0,type="n",xlim=c(-1.5,1.5),ylim=c(-1.5,1.5),frame=F,xlab="Eixo x",ylab="Eixo y")
points(coords,cex=0.1)
abline(h=c(0),v=c(0),lty=3,lwd=0.5)
lines(c(-1,-1,-1,1,1,1,1,-1),c(-1,1,1,1,1,-1,-1,-1),lwd=2,col="darkgray")
 
x.pos <- runif(1, min=-1, max=1)
y.pos <- runif(1, min=-1, max=1)
local.pos <- ifelse(x.pos^2 + y.pos^2 <= 1, TRUE, FALSE)
points(x.pos,y.pos,pch=19,cex=0.5,col=ifelse(local.pos,"red","blue"))
#dev.off()
 
#Jogando 500 pontos ao acaso dentro do quadrado
#jpeg("04.jpg")
plot(0,0,type="n",xlim=c(-1.5,1.5),ylim=c(-1.5,1.5),frame=F,xlab="Eixo x",ylab="Eixo y")
points(coords,cex=0.1)
abline(h=c(0),v=c(0),lty=3,lwd=0.5)
lines(c(-1,-1,-1,1,1,1,1,-1),c(-1,1,1,1,1,-1,-1,-1),lwd=2,col="darkgray")
 
for(i in 1:500) {
x.pos <- runif(1, min=-1, max=1)
y.pos <- runif(1, min=-1, max=1)
local.pos <- ifelse(x.pos^2 + y.pos^2 <= 1, TRUE, FALSE)
points(x.pos,y.pos,pch=19,cex=0.5,col=ifelse(local.pos,"red","blue"))
}
#dev.off()
 
#Fazendo essa simulação 5 mil vezes e olhando a proporção de pontos dentro do circulo
n<-5000
x.pos <- runif(n, min=-1, max=1)
y.pos <- runif(n, min=-1, max=1)
local.pos <- ifelse(x.pos^2 + y.pos^2 <= 1, TRUE, FALSE)
dentro <- length(which(local.pos == TRUE))
 
#Propoção
dentro/n
#Exato
(pi*1^2) / (2^2)
 
#Multiplicando pera area do quadrado para termos o valor de pi
4*(dentro/n)
#Valor de pi
pi
 
#Guardando as estimativas num vetor
estimativa<-vector()
 
for(i in 1:5000) {
    estimativa[i] <- length(which(local.pos[1:i] == TRUE)) / (i)
}
 
#Grafico animado
 
passos<-round(seq(1,5000,length=200))
estimativa<-estimativa*4
 
for(i in passos) {
jpeg(sprintf("pi%05d.jpg",i), width = 300, height = 300)
plot(c(1:5000),estimativa[1:5000]*4,type="n",ylim=c(min(estimativa),max(estimativa)),xlab="Iterações",
     ylab="Valor de Pi",main=paste("Iteração",i),frame=F)
points(c(1:i),estimativa[1:i],type="l",col="blue",lwd=2)
abline(h=pi,lty=3,col="red")
dev.off()
}
system("convert pi*.jpg -delay 10 pi.gif")
system("rm pi*.jpg")

Outro exemplo do paradoxo de Simpson.

Exemplos são sempre muito legais para ajudar a entender melhorar as coisas e esse aqui eu acho bem legal.

Suponha que você é o reitor ou reitora de uma universidade e existe uma suspeita que existe um desvio na razão sexual entre alunos que ingressam na universidade.

Para verificar se essa suspeita é verdadeira, você pega os dados referentes ao número de candidatos e aprovados de alguns cursos.

Começando com engenharia, a gente vê o seguinte:

, , curso = Engenharia sexo status Feminino Masculino Aprovados 80 450 Candidatos 100 900

Aqui temos para engenharia uma tabelinha resumindo os dados, então temos o número de candidatos e sexo de cada um.
Bem, o indice de aprovação das meninas é bem alto, 80%, se pegarmos o numero de aprovados e dividirmos pelo número de candidatos obtemos essa razão. Ja não podemos dizer o mesmo para os meninos, que nesse caso, 450 alunos de 900 foram aprovados. Uma razão 0.5, ou seja apenas 50% dos candidatos estão ingressando, bem inferior a frequência com que meninas são aprovadas. Tudo bem mas engenharia tem um número de candidatos maior do sexo masculino do que do sexo feminino, vamos olhar um curso de humanas para ver se temos algo diferente.

, , curso = Humanas sexo status Feminino Masculino Aprovados 180 10 Candidatos 900 100

Certo, aqui temos muito mais candidatas que candidatos, e vamos ver como estão as frequências de aprovação, aqui temos 900 candidatas e 180 aprovadas, isso da 20% de aprovação, dividindo 180 por 900 da 0.2 e multiplicando por 100 temos 20%. Ja para os meninos temos apenas 10% de aprovação, dos 100 candidatos apenas 10 entraram no curso.

Observamos dois cursos aqui e aparentemente em ambos as meninas tem um indice de aprovação muito maior que dos meninos, quase o dobro dos meninos sempre. Parece que existe sim um desvio na razão sexual, mas mulheres tem uma taxa de aprovação muito superior a dos homens.

No entanto essa conclusão pode não expressar toda a verdade.

Caso realizemos a mesma conta ignorando o curso, temos a seguinte tabela:

sexo status Feminino Masculino Aprovados 260 460 Candidatos 1000 1000

Agora se fizermos a mesma conta aqui, temos que a frequência de aprovação das meninas é de apenas 26% no total, contra 46% dos meninos, o sexo masculino tem praticamente o dobro de Aprovados que mulheres, mesmo que estas tenha uma aprovação superior em todos os cursos que observamos. Muito estranho não? Mas isso é um efeito conhecido como paradoxo de Simpson que ja vimos aqui anteriormente.

O que acontece é que existe uma interação entre o sexo e o curso, então apesar de em todos os cursos as mulheres obterem maior indice de aprovação, no geral existe mais homens entrando nessa universidade que mulheres, quase o dobro, e o desvio da razão sexual esta beneficiando pessoas do sexo masculino na verdade.

E é possível verificar isso apenas apresentando de forma diferente a mesma tabela, os mesmo dados e a mesma contagem.

Isso traz a tona um fato interessante, que é, nem sempre tudo que parece é, nem sempre a conclusão obvia é a que realmente vai trazer a decisão certa. Se houve-se uma politica para melhorar os indices de aprovação masculino nos cursos que estão deficientes, com certeza diminuiríamos a números muitos baixos a população feminina na sua universidade. Mas sem uma observação mais minuciosa dos dados, somos passiveis de sermos enganados facilmente para tomar essa decisão. Talvez incentivar um maior numero de mulheres a se inscreverem e tentar ingressar na universidade traria a razão sexual da universidade a 1:1.

Esse exemplo eu vi no curso de introdução a estatística do udacity, nesse link aqui. Um cursinho bem legal onde o professor assume um conhecimento básico das quatro operações básicas da matemática, então qualquer pessoa pode acompanhar sem problemas, uma excelente introdução para quem não sabe absolutamente nada e odeio cursos onde as pessoas só te ensinam distribuições e valores p.

Até a próxima. E como não tenho nenhum gráfico aqui para esse post, ficamos com a tirinha sobre o teorema de bayes do xkcd publicado essa semana 🙂

seashell

dados<-data.frame(
curso=c("Engenharia","Engenharia","Engenharia","Engenharia","Humanas","Humanas","Humanas",
"Humanas"),
status=c("Candidatos","Candidatos","Aprovados","Aprovados","Candidatos","Candidatos","Aprovados","Aprovados"),
sexo=c("Masculino","Feminino","Masculino","Feminino","Masculino","Feminino","Masculino",
"Feminino"),
n=c(900,100,450,80,100,900,10,180))
 
#Tabela por curso
xtabs(n~status+sexo+curso,data=dados)
 
#Tabela independente do curso
xtabs(n~status+sexo,data=dados)