Intervalo de confiança e erro da média

Algo que sempre é bem confuso, pelo menos pra mim sempre foi, é esse tal de intervalo de confiança e erro da média. A gente pode dar uma olhadinha com mais cuidado nisso, até porque todo tipo de modelo que a gente faz, do lado de todo parâmetro tem o bendito erro padrão, que é disso que estamos falando aqui, mas são contas tão estranhas, com tantas raízes quadradas que fica meio difícil as vezes, e ja da desanimo so de olhar aquelas formulas.

Então vamos olhar com um exemplo simples, vamos supor que existam umas ervinhas no chão do campo, e que vamos olhar a altura delas, porque queremos saber o quanto elas crescem, temos uma hipótese evolutiva sobre isso, que se elas crescerem demais ou de menos seu fitness é ruim. Dai vamos começar a coletar as alturas delas, que é ir com uma régua e medir elas no chão.

Como nossas ervinhas são virtuais, vamos gerar elas no R.

1
universo<-rnorm(1000000,5,2)

E pimba, simples assim, existem um milhão de indivíduos dessa ervilha no chão da fazenda que estamos coletando, essas um milhão são todas as ervas que existem. Acontece que suas costas não aguenta se abaixar para medir um milhão de de indivíduos, isso levaria uma vida inteira, além de que da hora que você começar a coletar, essa é uma erva anual, entes de você terminar de coletar, elas vão ter morrido, e você não tem como pagar para muitas pessoas ficarem la medindo, e muitos e muitos outros motivos, você tira uma amostra.

Como criamos nossa criação virtual, sabemos que elas vem de uma distribuição com média 5 e desvio padrão 2, mas essa informação está reservada a mãe natureza, você mero cientista não sabe disso.

Mas como seu trabalho é muito legal, da pra fazer ele com uma amostra, então você vai la e coleta 30 indivíduos, mede eles, calcular média, desvio padrão e com essas duas informações você consegui calcular o intervalo de confiança. Vamos começar uma figura aqui para representar isso.

01

O intervalo de confiança de 95%, 95% porque eu decidi assim, mas você pode usar qualquer porcentagem, precisam que quiser, eu peguei 95% porque é um valor comumente usado, mas esse número é uma escolha sua.
Então nessa figura, eu peguei 30 amostras

1
amostra<-sample(universo,30)

Depois calculei a média e o desvio padrão

> mean(amostra) [1] 5.250596 > sd(amostra) [1] 2.033259

E intervalo de confiança é a média mais 1.96 vezes o desvio padrão para cima e para baixo, então na figura o ponto é a média e esses dois valores são essas habinhas, esses guardachuvinhas, chame como quiser.

Para o que diabos isso serve? Serve para mim predizer o que eu vou achar se continuar coletando.
Por exemplo, se depois dessas trinta ervinhas, eu for coletar mais 100, veja so.

02

Cada ponto é uma nova ervinha, e eu pintei de azul aquelas que estão dentro do meu intervalo de confiança dessa primeira amostra de 30 indivíduos, e vermelho que ficou fora. Incrivelmente, baseado no meu intervalo de confiança de 95% encobriu 95 plantinhas e só ficou fora do intervalo 5, caramba, tudo bem que nem sempre temos algo exato assim, mas a ideia do intervalo de confiança é isso. Poxa para que continuar coletando, se com uma amostra eu consigo já saber como são todas as plantas que estão ali no chão. Então a partir de uma amostra, eu calculei essas métricas, esse intervalo de confiança e sei como é a maioria das ervinhas ali. Tudo bem, deve ter alguns mutantes ali muito muito altos, ou bem baixinhos, mas em média, eu já consigo prever com certeza precisa, de 95% que eu escolhi, como são as coisas ali. Claro que temos que levar em conta o número de amostras, como tiramos essas amostras para não pseudo-replicar, etc, etc. Mas é isso que está fazendo o intervalo de confiança, está me dando uma informação de como são as ervinhas ali, baseado na minha amostra.

Agora você pode se perguntar, poxa, mas e se fosse outras 30 ervas, não as 30 que eu peguei, daria certo?

Vamos tentar, vamos ver como ficaria o intervalo de confiança, para umas 30 outras amostras diferentes

03

E tharam, usando 30 amostras, não fica tudo idêntico, mas de modo geral, essa ideia ia funcionar, a gente ia conseguir uma ideia geral de como estão essas plantas la.

Podemos ainda comparar com intervalo de confiança do nosso universo.

04

Aqui em verde está o intervalo de confiança real, veja que os intervalos de confiança não são perfeitos, as vezes maiores, as vezes menores, as médias, que usamos para definir eles, as vezes ficam um pouco mais pra cima, um pouco mais pra baixo, mas de modo geral eles vão funcionar.

Alias aqui entra uma ideia legal, meio recursiva até, algo inception, mas legal. Veja que as médias, as bolinhas, ali para todos os intervalos de confiança oscilam, as vezes mais pra cima, as vezes mais pra baixo, e isso vai levar o intervalo de confiança mais pra cima ou mais pra baixo. Mas podemos, calcular um intervalo de confiança para essas médias, assim como calculamos para nossas varias amostras.

Vamos fazer isso.

05

Veja que legal, agora temos qual a variação na média, qual a região que ela deve estar, e veja que dentro dessa região, está o valor de média real, essa é a parte legal, podemos prever com alguma certeza, aqui 95%, onde a média real deve estar. Mas nesse momento você deveria parar e pensar, mas poxa, você salvou a média de 30 amostras de 30 plantas. 900 ervinhas, 900 abaixadas para medir por esse número, eu me daria por satisfeito com o intervalo de confiança e o que consigo fazer com ele.

Concordo, mas o legal é que com aquelas primeiras 30 ervinhas, minha primeira amostra, eu consigo calcular esse intervalo, com alguma precisão, esse é o tal do erro padrão.

Chegamos no ponto onde eu queria. Primeiro vamos calcular o erro padrão para nossas médias de plantinhas somente com a primeira amostra.

1
erro_media<-sd(amostra_inicial)/sqrt(30)

É muito simples, é simplesmente o desvio padrão da amostras dividido pela raiz quadrada do número de indivíduos medidos, que foi 30.
Vamos comparar isso com o nosso intervalo de confiança de 30 médias.

07

Novamente, não é perfeito, mas é muito bom, e engloba a média real, que é a bolinha verde.
Então veja que em cinza escuro é o erro padrão a partir do intervalo de confiança de várias amostras de 30 ervinhas (900 ervinhas) e o cinza claro é o erro padrão a partir de apenas uma amostra de 30 ervinhas.

Agora a média é um parâmetro de um modelo usado para descrever ervinhas, estamos usando uma distribuição normal, que é um modelo para descrever a distribuição do tamanho das ervinhas no chão da fazenda.
Então o intervalo de confiança a partir da média e do desvio serviu para ver como devem ser todas as plantas, mesmo as que não coletamos, veja que sempre temos muitos pontos azuis e poucos vermelhos, as vezes mais de 5, as vezes menos de 5, mas o intervalo de confiança é de 95%, então quase sempre acertamos 95% de qual altura deve ser o próximo ponto, ou seja essa é a informação de como são as ervinhas.

O erro padrão, é o que a gente acha que é o parâmetro do nosso modelo, a média. A distribuição tem dois parâmetros, média de desvio padrão, e não temos certeza de qual vai ser a média real, mas podemos afirmar, com 95% de certeza, quanto que no mínimo ela é e quanto no máximo ela é. E vimos que isso funciona.

Agora esse número a gente vê direto, e muitas vezes não da moral pra ele nos modelos.

Veja que calculamos aqui

> erro_media [1] 0.3347009

E onde esses cara ta nos modelos que a gente vive calculando?

> summary(lm(amostra_inicial~1)) Call: lm(formula = amostra_inicial ~ 1) Residuals: Min 1Q Median 3Q Max -4.7722 -1.3619 -0.3593 1.5944 3.4297 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 4.9685 0.3347 14.85 4.39e-15 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 1.833 on 29 degrees of freedom

Fazendo um modelo linear com a distribuição normal, onde so temos média, por isso é o ~1. Veja que então temos o intercept, que é a média e temos um tal de Std. Error, que é exatamente o erro padrão que calculamos, o erro da média, e isso quer dizer que podemos, a partir desse número calcular o intervalo de confiança desse parâmetro para ver o que ele deve ser.

Na maioria dos modelos temos muitos e muitos parâmetros, e isso é importante, porque muito mais que apenas ficar valor que estimamos pro parâmetro, podemos querer saber o quanto no máximo ele vai ser e o quanto no mínimo, porque por exemplo, em evolução, talvez para alguma teoria, tenhamos mais interesse no máximo, ou no mínimo que esse vai vai ser na natureza, do que o valor médio, para saber se algo é possível ou não.

E de modo geral, todo santo modelo que você calcula, sempre que estima parâmetros, esse valor vai estar la, então bom ter uma intuição do que ele significa, e com essa intuição, quem sabe, até melhorar a interpretação das analises estatísticas que utilizamos, sem acreditar em um número cegamente.

Bem é isso ai, acho que esse post é bom para refletir um pouco, ja que existe uma eterna discussão, que as vezes fica mais acirrada, as vezes menos sobre abandonar o uso de valores p e começar a reportar mais parâmetros estimados e intervalos de confiança, ciência sempre é difícil, ainda mais conseguir os dados certos numa quantidade suficiente, mas saber o que fazer com eles também, mas enquanto essa discussão perdura, o melhor que podemos fazer é saber o que é o valor p, o que são intervalos de confiança e erro padrão para poder defender melhor nossos trabalhos e conclusões.

O script vai estar la no repositório recologia, e se eu escrevi alguma bobeira, algo errado, deixe um comentário corrigindo ou mande um e-mail e até mais.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
#http://recologia.com.br/2015/07/intervalo-de-confianca-e-erro-da-media
##################################
## Gerando dados
##################################
set.seed(1)
universo<-rnorm(1000000,5,2)
 
mean(universo)
sd(universo)
 
amostra<-sample(universo,30)
 
mean(amostra)
sd(amostra)
 
##################################
## Intervalo de confiança
##################################
#jpeg("01.jpg")
plot(0,mean(amostra),pch=19,xlim=c(0,2.1),ylim=c(-2,12),frame=F,
     ylab="Universo",xlab="")
arrows(x0=0,y0=mean(amostra),y1=mean(amostra)-1.96*sd(amostra),
       length = 0.05, angle = 90)
arrows(x0=0,y0=mean(amostra),y1=mean(amostra)+1.96*sd(amostra),
       length = 0.05, angle = 90)
#dev.off()
 
##################################
## Amostras do universo
##################################
#jpeg("02.jpg")
plot(0,mean(amostra),pch=19,xlim=c(0,2.1),ylim=c(-2,12),frame=F,
     ylab="Universo",xlab="")
arrows(x0=0,y0=mean(amostra),y1=mean(amostra)-1.96*sd(amostra),
       length = 0.05, angle = 90)
arrows(x0=0,y0=mean(amostra),y1=mean(amostra)+1.96*sd(amostra),
       length = 0.05, angle = 90)
 
for(i in seq(0.05,0.95,length=100)){
    unidade<-sample(universo,1)
    if((unidade>mean(amostra)-1.96*sd(amostra)) & (unidade<mean(amostra)+1.96*sd(amostra))){
        points(i,unidade,pch=19,col="blue",cex=0.5)
    } else {
        points(i,unidade,pch=19,col="red",cex=0.5)
    }
}
#dev.off()
 
##################################
## Muitos intervalos de confiança
##################################
#jpeg("03.jpg")
plot(0,mean(amostra),pch=19,xlim=c(0,2.1),ylim=c(-2,12),frame=F,
     ylab="Universo",xlab="")
arrows(x0=0,y0=mean(amostra),y1=mean(amostra)-1.96*sd(amostra),
       length = 0.05, angle = 90)
arrows(x0=0,y0=mean(amostra),y1=mean(amostra)+1.96*sd(amostra),
       length = 0.05, angle = 90)
 
for(i in seq(0.05,0.95,length=100)){
    unidade<-sample(universo,1)
    if((unidade>mean(amostra)-1.96*sd(amostra)) & (unidade<mean(amostra)+1.96*sd(amostra))){
        points(i,unidade,pch=19,col="blue",cex=0.5)
    } else {
        points(i,unidade,pch=19,col="red",cex=0.5)
    }
}
 
media<-mean(amostra)
amostra_inicial<-amostra
j=2
for(i in seq(1.05,1.95,length=30)) {
    amostra<-sample(universo,30)
    media[j]<-mean(amostra)
    j=j+1
    points(i,mean(amostra),pch=19)
    arrows(x0=i,y0=mean(amostra),y1=mean(amostra)-1.96*sd(amostra),
           length = 0.05, angle = 90)
    arrows(x0=i,y0=mean(amostra),y1=mean(amostra)+1.96*sd(amostra),
           length = 0.05, angle = 90)
}
#dev.off()
 
 
##################################
## Intervalo de confiança real
##################################
#jpeg("04.jpg")
plot(0,mean(amostra),pch=19,xlim=c(0,2.1),ylim=c(-2,12),frame=F,
     ylab="Universo",xlab="")
arrows(x0=0,y0=mean(amostra),y1=mean(amostra)-1.96*sd(amostra),
       length = 0.05, angle = 90)
arrows(x0=0,y0=mean(amostra),y1=mean(amostra)+1.96*sd(amostra),
       length = 0.05, angle = 90)
 
for(i in seq(0.05,0.95,length=100)){
    unidade<-sample(universo,1)
    if((unidade>mean(amostra)-1.96*sd(amostra)) & (unidade<mean(amostra)+1.96*sd(amostra))){
        points(i,unidade,pch=19,col="blue",cex=0.5)
    } else {
        points(i,unidade,pch=19,col="red",cex=0.5)
    }
}
 
media<-mean(amostra)
amostra_inicial<-amostra
j=2
for(i in seq(1.05,1.95,length=30)) {
    amostra<-sample(universo,30)
    media[j]<-mean(amostra)
    j=j+1
    points(i,mean(amostra),pch=19)
    arrows(x0=i,y0=mean(amostra),y1=mean(amostra)-1.96*sd(amostra),
           length = 0.05, angle = 90)
    arrows(x0=i,y0=mean(amostra),y1=mean(amostra)+1.96*sd(amostra),
           length = 0.05, angle = 90)
}
 
points(2.05,mean(universo),pch=19,col="green")
arrows(x0=2.05,y0=mean(universo),y1=mean(universo)-1.96*sd(universo),
       length = 0.02, angle = 90,col="green")
arrows(x0=2.05,y0=mean(universo),y1=mean(universo)+1.96*sd(universo),
       length = 0.02, angle = 90,col="green")
#dev.off()
 
 
##################################
## Erro da média
##################################
#jpeg("05.jpg")
plot(0,mean(amostra),pch=19,xlim=c(0,2.1),ylim=c(-2,12),frame=F,
     ylab="Universo",xlab="")
arrows(x0=0,y0=mean(amostra),y1=mean(amostra)-1.96*sd(amostra),
       length = 0.05, angle = 90)
arrows(x0=0,y0=mean(amostra),y1=mean(amostra)+1.96*sd(amostra),
       length = 0.05, angle = 90)
 
for(i in seq(0.05,0.95,length=100)){
    unidade<-sample(universo,1)
    if((unidade>mean(amostra)-1.96*sd(amostra)) & (unidade<mean(amostra)+1.96*sd(amostra))){
        points(i,unidade,pch=19,col="blue",cex=0.5)
    } else {
        points(i,unidade,pch=19,col="red",cex=0.5)
    }
}
 
media<-mean(amostra)
amostra_inicial<-amostra
j=2
for(i in seq(1.05,1.95,length=30)) {
    amostra<-sample(universo,30)
    media[j]<-mean(amostra)
    j=j+1
    points(i,mean(amostra),pch=19)
    arrows(x0=i,y0=mean(amostra),y1=mean(amostra)-1.96*sd(amostra),
           length = 0.05, angle = 90)
    arrows(x0=i,y0=mean(amostra),y1=mean(amostra)+1.96*sd(amostra),
           length = 0.05, angle = 90)
}
 
points(2.05,mean(universo),pch=19,col="green")
arrows(x0=2.05,y0=mean(universo),y1=mean(universo)-1.96*sd(universo),
       length = 0.02, angle = 90,col="green")
arrows(x0=2.05,y0=mean(universo),y1=mean(universo)+1.96*sd(universo),
       length = 0.02, angle = 90,col="green")
 
points(1.98,mean(media),pch=19,col="darkgray",cex=0.8)
arrows(x0=1.98,y0=mean(media),y1=mean(media)-1.96*sd(media),
       length = 0.02, angle = 90,col="darkgray")
arrows(x0=1.98,y0=mean(media),y1=mean(media)+1.96*sd(media),
       length = 0.02, angle = 90,col="darkgray")
#dev.off()
 
 
##################################
## Erro da média usando apenas a amostra
##################################
#jpeg("06.jpg")
plot(0,mean(amostra),pch=19,xlim=c(0,2.1),ylim=c(-2,12),frame=F,
     ylab="Universo",xlab="")
arrows(x0=0,y0=mean(amostra),y1=mean(amostra)-1.96*sd(amostra),
       length = 0.05, angle = 90)
arrows(x0=0,y0=mean(amostra),y1=mean(amostra)+1.96*sd(amostra),
       length = 0.05, angle = 90)
 
for(i in seq(0.05,0.95,length=100)){
    unidade<-sample(universo,1)
    if((unidade>mean(amostra)-1.96*sd(amostra)) & (unidade<mean(amostra)+1.96*sd(amostra))){
        points(i,unidade,pch=19,col="blue",cex=0.5)
    } else {
        points(i,unidade,pch=19,col="red",cex=0.5)
    }
}
 
media<-mean(amostra)
amostra_inicial<-amostra
j=2
for(i in seq(1.05,1.95,length=30)) {
    amostra<-sample(universo,30)
    media[j]<-mean(amostra)
    j=j+1
    points(i,mean(amostra),pch=19)
    arrows(x0=i,y0=mean(amostra),y1=mean(amostra)-1.96*sd(amostra),
           length = 0.05, angle = 90)
    arrows(x0=i,y0=mean(amostra),y1=mean(amostra)+1.96*sd(amostra),
           length = 0.05, angle = 90)
}
 
points(2.05,mean(universo),pch=19,col="green")
arrows(x0=2.05,y0=mean(universo),y1=mean(universo)-1.96*sd(universo),
       length = 0.02, angle = 90,col="green")
arrows(x0=2.05,y0=mean(universo),y1=mean(universo)+1.96*sd(universo),
       length = 0.02, angle = 90,col="green")
 
points(1.98,mean(media),pch=19,col="darkgray",cex=0.8)
arrows(x0=1.98,y0=mean(media),y1=mean(media)-1.96*sd(media),
       length = 0.02, angle = 90,col="darkgray")
arrows(x0=1.98,y0=mean(media),y1=mean(media)+1.96*sd(media),
       length = 0.02, angle = 90,col="darkgray")
 
erro_media<-sd(amostra_inicial)/sqrt(30)
points(2.02,mean(amostra_inicial),pch=19,col="lightgray",cex=0.8)
arrows(x0=2.02,y0=mean(amostra_inicial),y1=mean(amostra_inicial)-1.96*erro_media,
       length = 0.02, angle = 90,col="lightgray")
arrows(x0=2.02,y0=mean(amostra_inicial),y1=mean(amostra_inicial)+1.96*erro_media,
       length = 0.02, angle = 90,col="lightgray")
#dev.off()
 
##################################
## Exemplo do summary
##################################
erro_media
summary(lm(amostra_inicial~1))

Regressão linear ajustada com lm, bbmle, OpenBugs, Jags e Stan no R

No último post a gente falou do pacote sads do R. Muito legal, e ele faz vasto uso do pacote bbmle.

O pacote bbmle é bem legal, ele permite você escrever uma função de verosimilhança qualquer e então ele ajusta ela pra você. A gente tem um monte de post, por exemplo aqui, aqui e aqui onde o que a gente faz no fim das contas é escrever a função de verosimilhança e ajustar. Mas a gente sempre tinha ajustado usando estatística bayesiana, usando o OpenBugs pra ser mais preciso, e para conferir usando estatística frequentista usando o lme4, ou algum outro pacote que fazia o modelo pra gente.

Mas o bblme vai funcionar mais ou menos na mesma estrategia da linguagem bugs, escreve seu modelo e o bbmle ajusta pra você.

Vamos pegar um modelo simples aqui, a boa e velha regressão linear, e criar alguns dados.

01

Beleza, o mais simples de fazer é usar a função lm que ajusta modelos lineares no R a ajustar o modelo e olhar o resultado e pronto, simples assim.

> fit.lm <- lm(resposta~preditor) > summary(fit.lm) Call: lm(formula = resposta ~ preditor) Residuals: Min 1Q Median 3Q Max -1.0381 -0.4135 0.1270 0.3737 0.7836 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 2.41354 0.32834 7.351 5.28e-08 *** preditor 0.47595 0.01535 31.004 < 2e-16 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 0.491 on 28 degrees of freedom Multiple R-squared: 0.9717, Adjusted R-squared: 0.9707 F-statistic: 961.3 on 1 and 28 DF, p-value: < 2.2e-16

Certo, agora vamos lembrar como é a formula da reta da regressão.
O que a gente fez é:

resposta = \alpha + \beta \cdot preditor + erro

E esse erro é a parte onde entra a distribuição normal, gaussiana ou sei la, mas o erro então é o seguinte.

erro = Normal(0,\sigma^2)

A distribuição com média zero e desvio sigma.

Daria na mesma o que a gente escreveu acima estar tudo junto:

resposta = \alpha + \beta \cdot preditor + Normal(0,\sigma^2)

Ou ainda:

resposta = Normal(\alpha + \beta \cdot preditor,\sigma^2)

Ou as vezes ao invés daquele escrito normal acho que a gente vê escrito p ou P, mas as pessoas usam p pra tanta coisa em estatística que não é a toa que tudo fica confuso, e não sou só eu que fica confuso, veja so aqui.

Mas a gente consegue perceber, pelo gráfico, que a resposta e o preditor são nossos dados, que temos valores, e o que sobra são os chamados parâmetros nessa formula toda. Na regressão, a gente está lidando com três deles, dois que determinam a reta, o alpha (\alpha) que é onde ela cruza o eixo de resposta e o beta (\beta) que é a inclinação. Além desses dois caras temos o desvio padrão, que é o sigma (\sigma), esse cara diz o quanto afastado daquela reta as nossas amostras podem estar, não exatamente isso mas mais ou menos isso, a dispersão dos pontos em relação a reta.

Ok, por mais confuso que isso possa parecer, se esse é o modelo, e resposta e preditor nos temos, sempre esses três caras tem que sair do modelo, que no caso a gente otimiza para conseguir o melhor ajuste.

Então o bbmle, assim como o OpenBugs e companhia vai precisar de um modelo construído, que diga pra ele como esta a qualidade do ajuste, verosimilhança do modelo, e ele te diz quais os melhores valores para esses três parâmetros.

Na saida da regressão ali em cima, o alpha é o intercept (2.41354), o beta é o que ele chamou de preditor(0.47595) e por último, o sigma é o Residual standard error (0.491). Os parâmetros alpha e beta não se tem como ter certeza absoluta, so quem tem certeza é quem gerou os dados, como na natureza é ela que gera os dados, então a gente nunca tem certeza, mas a gente acha que esses são os melhores valores (lembrando que cada valor tem seu erro aceitável também, ali do lado).

Mas então, o lm faz todo esse processo pra gente, e não tem porque não usar ele para modelos assim, mas só para treinar, ou curiosidade, vamos ver como ficaria esse modelo no bbmle.

No bbmle a gente tem que fazer nossa própria função de verosimilhança.

funcao.likelihood <- function(alpha=1,beta=1,desvio=1) {
    -sum(dnorm(resposta,mean=alpha+beta*preditor,sd=desvio,log=TRUE))
}

Antes de pensar que essa função faz muita coisa, vamos só usar ela, pra ver como é mais simples que parece, ela simplesmente calcula um número, que é a verosimilhança do modelo, na verdade o negative log-likelihood, veja que tem um menos antes do sum ali.

Então você fala um valor pra cada parâmetro e a função fala se o quão bom é o ajuste, em termos de verossimilhança.

> funcao.likelihood(alpha=1,beta=1,desvio=1) [1] 1488.255

Se a gente muda esses valores, o ajuste muda

> funcao.likelihood(alpha=1.5,beta=1,desvio=1) [1] 1632.549

Quanto menor melhor, quanto mais a gente aproxima todos os valores dos valores que usamos para gerar os dados, menor ele vai ficar.

> funcao.likelihood(alpha=1.5,beta=0.9,desvio=1) [1] 1038.321

E se a gente usar parâmetros muito nada a ver, esse número fica gigante.

> funcao.likelihood(alpha=8,beta=0.9,desvio=1) [1] 3195.411

E isso é tudo que o bbmle precisa. Agora para ajustar os parâmetros é so usar essa função de verosimilhança na função mle do pacote bbmle.

> fit.mle <- mle2(funcao.likelihood) Houve 11 avisos (use warnings() para vê-los) > summary(fit.mle) Maximum likelihood estimation Call: mle2(minuslogl = funcao.likelihood) Coefficients: Estimate Std. Error z value Pr(z) alpha 2.413540 0.317207 7.6087 2.768e-14 *** beta 0.475945 0.014830 32.0925 < 2.2e-16 *** desvio 0.474314 0.061234 7.7459 9.487e-15 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 -2 log L: 40.38304

E pimba, estimamos tudo do modelo. O desvio, ou sigma aqui fica menor que na função lm, porque o método usado é diferente, mas basicamente temos o mesmo resultado. Recebemos um aviso também, que acho que porque pode ser tentando valores negativos de desvio na função, mas não existe desvio negativo, dai a função produz um NaN, tipo assim.

> funcao.likelihood(alpha=8,beta=0.9,desvio=-1) [1] NaN Mensagens de aviso perdidas: In dnorm(resposta, mean = alpha + beta * preditor, sd = desvio, : NaNs produzidos

Mas o mle acha o caminho certo e ajusta. Ele usa o método do Nelder and Mead como default, a gente falou do método do Newton aqui, e vimos o gradient descent aqui, esse é outro método que serve para achar os melhores valores para os parâmetros.

Mas o bbmle é legal que da para usar outros métodos de otimização, só citar qual método você quer usar.

> fit.mle2<-mle2(funcao.likelihood,method="L-BFGS-B",lower=c(alpha=0.0001,beta=0.0001,desvio=0.0001)) > summary(fit.mle2) Maximum likelihood estimation Call: mle2(minuslogl = funcao.likelihood, method = "L-BFGS-B", lower = c(alpha = 1e-04, beta = 1e-04, desvio = 1e-04)) Coefficients: Estimate Std. Error z value Pr(z) alpha 2.413540 0.317207 7.6087 2.768e-14 *** beta 0.475945 0.014830 32.0924 < 2.2e-16 *** desvio 0.474314 0.061234 7.7459 9.488e-15 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 -2 log L: 40.38304

Ou ainda usar outra função para otimizar, que não o optmin que é o otimizador default do R, eu acho.

> fit.mle3<-mle2(funcao.likelihood,optimizer="nlminb",lower=c(alpha=0.0001,beta=0.0001,desvio=0.0001)) > summary(fit.mle3) Maximum likelihood estimation Call: mle2(minuslogl = funcao.likelihood, optimizer = "nlminb", lower = c(alpha = 1e-04, beta = 1e-04, desvio = 1e-04)) Coefficients: Estimate Std. Error z value Pr(z) alpha 2.413539 0.317206 7.6087 2.768e-14 *** beta 0.475945 0.014830 32.0925 < 2.2e-16 *** desvio 0.474313 0.061234 7.7460 9.486e-15 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 -2 log L: 40.38304

Veja que esses outros métodos nos permitiram fazer mais coisas, uma delas foi delimitar limites para os parâmetros. No caso colocamos que tudo tem que ser maior que 0.0001, isso impossibilita um valor negativo pro desvio e veja que depois disso a gente não recebe mais warnings.
Eu acho que a parte boa do bbmle, é que se amanha a gente lê um artigo, com um modelo muito louco que alguém acabou de inventar, se a gente conseguir passar ele para uma função de verossimilhança, a gente consegue ajustar valores pra ele.

Alias o cara que escreveu esse pacote, o Ben Bolker, que é um cara que parece ser, além de inteligente, bem legal, já que ele responde todo mundo nos fóruns de R, e em tudo que é fórum da internet, bem ele escreveu um livro chamado Ecological Models and Data in R, que é muito legal, pois ao invés de ensinar sobre testes e mais testes com seus nomes complicados, ele ensina que existe funções determinísticas, distribuições estocásticas,e você juntas essas coisas e descreve o mundo, ai depois de muitos capítulos ele mostra que você continua fazendo a mesma coisa que fazia nas anovas da vida, mas usando o modelo que você pensa e não pensando o que pode testar com os dados que tem.

Mas já que estamos aqui, veja que na estatística bayesiana, apesar de partir de pressupostos diferentes, paradigmas diferente e tal, no final a gente quer descrever a mesma coisa, e vai ser a mesma coisa na pratica, escrever um modelo, ajustar os dados, essa mesma regressão é feita da seguinte forma no OpenBugs

Escrevemos o modelo

sink("linreg.txt")
cat("
model {
# Prior
 alpha ~ dnorm(0,0.001)
 beta ~ dnorm(0,0.001)
 sigma ~ dunif(0, 100)
 tau <- 1/ (sigma * sigma)
 
# Verossimilhança
 for (i in 1:n) {
    y[i] ~ dnorm(mu[i], tau)
    mu[i] <- alpha + beta*x[i]
 }
}
",fill=TRUE)
sink()

Depois temos que preparar o terreno, preparar os dados para exportar, criar uma função para inicializar a cadeia, definir os parâmetros da cadeia e então otimizar.

> bugs.data <- list(x=preditor,y=resposta,n=length(resposta)) > inits <- function(){ list(alpha=rnorm(1), beta=rnorm(1), sigma = rlnorm(1))} > params <- c("alpha","beta","sigma") > nc = 3 #Numero de cadeias > ni=1200 #Tamanho da cadeira > nb=200 #Numero de simulação que serão descartadas > nt=1 #Thinning rate > modelo1.bugs <-bugs(data = bugs.data, inits = inits, parameters = params,model = "linreg.txt",n.thin = nt,n.chains = nc,n.burnin = nb, n.iter = ni) > print(modelo1.bugs, dig = 3) Inference for Bugs model at "linreg.txt", Current: 3 chains, each with 1200 iterations (first 200 discarded) Cumulative: n.sims = 3000 iterations saved mean sd 2.5% 25% 50% 75% 97.5% Rhat n.eff alpha 2.460 0.354 1.812 2.223 2.438 2.680 3.230 1.013 170 beta 0.474 0.017 0.436 0.464 0.475 0.485 0.504 1.017 140 sigma 0.514 0.071 0.398 0.462 0.506 0.559 0.674 1.002 1600 deviance 43.640 2.655 40.630 41.700 42.960 44.780 50.750 1.003 1600 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 = 2.9 and DIC = 46.5 DIC is an estimate of expected predictive error (lower deviance is better).

Aqui o sigma fica um pouquinho maior, mas ainda temos basicamente o mesmo resultado. O bbmle é um pacote do R, aqui é meio diferente, porque o OpenBugs é um programa, de vida própria, que do R a gente um pacote chamado R R2OpenBUGS que tudo que faz é mandar todas essas informações que separamos no R para o OpenBugs, o OpenBugs processa elas e depois devolve bonitinho pro R, tudo organizadinho para então a gente interpretar, e fazer o que mais quiser como gráficos.

O OpenBugs, que na verdade vem da linguagem Bugs, de fazer inferência bayesiana, no caso Bugs significa Bayesian inference Using Gibbs Sampling, que usa um tipo de amostrador para tentar descobrir os valores reais dos parâmetros, e chegar ao melhor chute pelo menos. No final a gente queria saber o valor real do parâmetro, mas nunca saberemos, tem um post sobre Gibbs Sampler aqui, se a curiosidade bater.

Mas o OpenBugs não é o único a usar a linguagem Bugs para inferência, o jeitão de escrever o modelo etc, temos o WinBugs, que não esta mais em desenvolvimento, e o Jags, Just another Gibbs Sampler, que está a todo vapor eu acho. Na verdade ajustar modelos com o Jags, que é um programa externo também, é bem parecido, a gente pode ajustar o mesmo modelo denovo, da seguinte forma.

> jags <- jags.model(file='linreg.txt',data=bugs.data,inits=inits,n.chains = 3,n.adapt=100) Compiling model graph Resolving undeclared variables Allocating nodes Graph Size: 130 Initializing model |++++++++++++++++++++++++++++++++++++++++++++++++++| 100% > update(jags, 1000) |**************************************************| 100% > jags.samples(jags,c('alpha', 'beta','sigma'),1000) |**************************************************| 100% $alpha mcarray: [1] 2.380524 Marginalizing over: iteration(1000),chain(3) $beta mcarray: [1] 0.4774447 Marginalizing over: iteration(1000),chain(3) $sigma mcarray: [1] 0.5097927 Marginalizing over: iteration(1000),chain(3)

E temos nossos três parâmetros denovo, veja que usamos tudo igual ao que usamos no OpenBugs. O modelo, o jeito de mandar os dados, numa lista e tal. A gente ouve muito falar do Jags, eu acho que uso mais o OpenBugs porque os livros que olho sempre o pessoal ta usando ele, mas migrar para os modelos pro Jags é relativamente fácil, como acabamos de ver.

Agora uma outra nova opção que surgiu foi o Stan. Stan é um programa para ajustar modelos bayesianos, e mais novo, acho que tem mais suporte a computação paralela, mas ao contrario do OpenBugs e do Jags que usam o Gibbs Sampler, ele usa um negocio chamado Hamiltonian Monte Carlo, esse eu não tenho nem ideia de como funciona, mas a gente pode ajustar a mesma regressão com ele, o que não é tão difícil, e ver que funciona igual.

O jeito de construir os modelos nele é bem diferente, essa nossa regressão fica algo assim.

linreg_stan <- '
data{
  int n;
  real y[n];
  real x[n];
}
 
parameters{
  real alpha;
  real beta;
  real<lower=0, upper=100> sigma;
}
 
transformed parameters{
    real mu[n];
    for(i in 1:n){
    mu[i] <- alpha + beta*x[i];
  }
}
 
model{
  //Priors
  alpha ~ normal(0, 100);
  beta ~ normal(0, 100);
  sigma ~ uniform(0, 100);
 
  //Likelihood
  for(i in 1:n){
    y[i] ~ normal(mu[i], sigma);
}
}
'

O que no final das contas nem é tão diferente assim do Bugs, é definir a função de verosimilhança, como é o framework bayesiano, a gente tem os prior, só que além disso la no inicio, a gente tem que inicializar todos os valores que vamos usar, como no c++, linguagem que ele é feito. Dai é so preparar os dados para mandar, numa lista igual ao bugs e mandar ajustar.

> stan.data <- list(n = length(resposta), x = preditor, y = resposta) > fit.stan<- stan(model_code=linreg_stan,data=stan.data,iter = 1000,chains = 3,thin = 1) TRANSLATING MODEL 'linreg_stan' FROM Stan CODE TO C++ CODE NOW. COMPILING THE C++ CODE FOR MODEL 'linreg_stan' NOW. SAMPLING FOR MODEL 'linreg_stan' NOW (CHAIN 1). Iteration: 1 / 1000 [ 0%] (Warmup) Iteration: 100 / 1000 [ 10%] (Warmup) Iteration: 200 / 1000 [ 20%] (Warmup) . . . Um monte de mensagens . . . Iteration: 900 / 1000 [ 90%] (Sampling) Iteration: 1000 / 1000 [100%] (Sampling) # Elapsed Time: 0.166985 seconds (Warm-up) # 0.131985 seconds (Sampling) # 0.29897 seconds (Total)

E podemos olhar nossos resultado de forma parecida ao output do OpenBugs

> print(fit.stan,digits=3) Inference for Stan model: linreg_stan. 3 chains, each with iter=1000; warmup=500; thin=1; post-warmup draws per chain=500, total post-warmup draws=1500. mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat alpha 2.420 0.020 0.369 1.668 2.184 2.419 2.656 3.127 324 1.005 beta 0.476 0.001 0.017 0.442 0.465 0.475 0.487 0.511 331 1.003 sigma 0.520 0.004 0.077 0.395 0.464 0.513 0.567 0.693 454 1.003 lp__ 4.950 0.078 1.449 1.411 4.339 5.278 5.954 6.517 347 1.002 Samples were drawn using NUTS(diag_e) at Thu Jul 17 23:01:30 2014. For each parameter, n_eff is a crude measure of effective sample size, and Rhat is the potential scale reduction factor on split chains (at convergence, Rhat=1).

E denovo os mesmo valores de alpha, beta e sigma.
Bem a gente começou olhando o bbmle, usando o lm pra ver se estávamos chegando a valores razoáveis para os parâmetros da regressão e acabamos olhando também como ajusta usando o OpenBugs, Jags e Stan.

A primeira pergunta agora seria, qual caminho seguir?

A minha melhor resposta é:
Eu não faço a menor ideia, mas é legal pra caralho chegar quase no mesmo valor de alpha, beta e sigma de todos os jeitos.
Na pior da hipóteses, se o modelo que estivermos tentando ajustar for certo, vai dar certo. Claro que para uma regressão linear a gente vai usar o lm na verdade, acho que não faz sentido adicionar toda a complexidade desses outros métodos para algo simples. Mas para muitos outros modelos de nosso interesse, que não da para ajustar com lm, ou mesmo glm, o jeito é começar a embrenhar nesses métodos.

Claro que se da para encontrar as coisas prontas, como os modelos prontinhos no pacote sads, vamos usar la e não tentar reinventar a roda. Mas as vezes a gente quer mudar algo um pouquinho, e ai saber como tudo funciona é legal. Além disso, como no livro do Bolker, um comentário que acho muito interessante, não interessa o método que você usa, é sempre bom saber um pouquinho de cada método para conseguir ler melhor a literatura, já que sempre cada pessoa vai fazer as coisas de um jeito, e não vamos deixar de ler um artigo porque ele usa estatística bayesiana ou frequentista.

Bem e esse post é legal, porque eu sempre quis ver alguém ajustando algo que eu entenda, como uma regressão linear com todos esses programas e pacotes, como ninguém fez isso pra mim ler, eu mesmo fiz ^^. Agora eu tenho uma referencia quando quiser usar algo diferente do OpenBugs.

E é isso. o script vai estar la no repositório recologia, se eu escrevi alguma bobeira, algo errado, deixe um comentário corrigindo, ou mande um e-mail e ficamos por aqui.

set.seed(171)
preditor<-runif(30,10,30)
resposta<-rnorm(30,2+0.5*preditor,0.5)
 
plot(resposta~preditor)
 
############################
#lm do stats
############################
fit.lm <- lm(resposta~preditor)
summary(fit.lm)
 
############################
#bbmle
############################
library(bbmle)
 
funcao.likelihood <- function(alpha=1,beta=1,desvio=1) {
    -sum(dnorm(resposta,mean=alpha+beta*preditor,sd=desvio,log=TRUE))
}
 
fit.mle <- mle2(funcao.likelihood)
summary(fit.mle)
warnings()
 
fit.mle2<-mle2(funcao.likelihood,method="L-BFGS-B",lower=c(alpha=0.0001,beta=0.0001,desvio=0.0001))
summary(fit.mle2)
 
fit.mle3<-mle2(funcao.likelihood,optimizer="nlminb",lower=c(alpha=0.0001,beta=0.0001,desvio=0.0001))
summary(fit.mle3)
 
############################
#OpenBugs
############################
library(R2OpenBUGS)
 
sink("linreg.txt")
cat("
model {
# Prior
 alpha ~ dnorm(0,0.001)
 beta ~ dnorm(0,0.001)
 sigma ~ dunif(0, 100)
 tau <- 1/ (sigma * sigma)
 
# Verossimilhança
 for (i in 1:n) {
    y[i] ~ dnorm(mu[i], tau)
    mu[i] <- alpha + beta*x[i]
 }
}
",fill=TRUE)
sink()
 
bugs.data <- list(x=preditor,y=resposta,n=length(resposta))
 
# Função de inicialização
inits <- function(){ list(alpha=rnorm(1), beta=rnorm(1), sigma = rlnorm(1))}
 
# Escolha os parametros que quer estimar
params <- c("alpha","beta","sigma")
 
# Caracteristicas do MCMC
nc = 3 #Numero de cadeias
ni=1200 #Tamanho da cadeira
nb=200 #Numero de simulação que serão descartadas
nt=1 #Thinning rate
 
# Inicie o Amostrador
 
modelo1.bugs <-bugs(data = bugs.data, inits = inits, parameters = params,model = "linreg.txt",n.thin = nt,
                     n.chains = nc,n.burnin = nb, n.iter = ni)
 
print(modelo1.bugs, dig = 3)
 
############################
#Jags
############################
library(rjags)
 
jags <- jags.model(file='linreg.txt',data=bugs.data,inits=inits,n.chains = 3,n.adapt=100)
update(jags, 1000)
jags.samples(jags,c('alpha', 'beta','sigma'),1000)
 
 
############################
#Stan
############################
library(rstan)
 
linreg_stan <- '
data{
  int n;
  real y[n];
  real x[n];
}
 
parameters{
  real alpha;
  real beta;
  real<lower=0, upper=100> sigma;
}
 
transformed parameters{
    real mu[n];
    for(i in 1:n){
    mu[i] <- alpha + beta*x[i];
  }
}
 
model{
  //Priors
  alpha ~ normal(0, 100);
  beta ~ normal(0, 100);
  sigma ~ uniform(0, 100);
 
  //Likelihood
  for(i in 1:n){
    y[i] ~ normal(mu[i], sigma);
}
}
'
 
stan.data <- list(n = length(resposta), x = preditor, y = resposta)
 
fit.stan<- stan(model_code=linreg_stan,data=stan.data,iter = 1000,chains = 3,thin = 1)
print(fit.stan,digits=3)

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")

Usando Offset em modelos lineares

A gente falou aqui sobre modelos com Overdispersion e aqui sobre modelos com Zeros-Inflados.

Usando a linguagem Bugs, por um lado a gente pode sofre um pouco por ser obrigado a definir o modelo, mas por outro lado, a gente ganha no melhor entendimento do que está acontecendo. Overdispersion é um exemplo, que somente depois que eu vi o modelo no Bugs que eu fui entender o que diabos era aquele número que saio no glm quando a gente usava a family=quasipoisson.

Agora ta na hora de entender o que significa offset, aquela opção que você pode usar nos modelos lineares da função lm, nos modelos gerais lineares na função glm, afinal, em quase todo lugar.

Se a gente olhar a documentação da função com ?lm, a gente vai ver a descrição o argumento offset:

“this can be used to specify an a priori known component to be included in the linear predictor during fitting. This should be NULL or a numeric vector of length equal to the number of cases. One or more offset terms can be included in the formula instead or as well, and if more than one are specified their sum is used. See model.offset.”

No meu inglês nórdico, podemos ler:

“Este pode ser usado para especificar um componente conhecido a priori para ser incluído como preditor linear no ajuste. Este deve ser NULL ou um vetor numérico de tamanho igual ao número de casos. Um ou mais termos de offset podem ser incluídos na formula ao invés ou tal como, e se mais de um for especificado, sua soma será usada, olhe model.offset.”

Se você olhar a função model.offset, vai cair nas funções que transformam formula em matrizes, álgebra linear e não da mais para entender nada.

Mas vamos pensar de maneira mais simples. No Glm de Poisson, como estamos olhando, nos assumimos que as contagens esperadas são adequadamente descritas pelos efeitos das covariáveis. Porém, frequentemente, nos temos algo como uma “janela de contagem” que não é constante. Por exemplo, áreas de estudos que não tem o mesmo tamanho ou amostras que tem duração diferentes, a gente pode não determinar antes o tempo que vai olhar uma lagoa, simplesmente vamos olhando e contando quantos sapinhos vemos la. Para adicionar essa variação no modelo de Poisson, por exemplo se temos contagens de coelhos em vários campos de tamanho diferentes, nos podemos definir a área como offset, mais especificamente, como o modelo de Poisson a gente usa função de ligação logarítmica, a gente define o log da área de cada amostras. Simplificando, a gente modela a densidade como a resposta, e não a contagem per se.

Vamos fazer nosso exemplo aqui, vamos pensar em coelhinhos. Estamos contando coelhinhos em pastos e plantações, vamos fazer varias amostras, mais ou menos assim:

01

Cada pontinho é uma amostra. Mas cada uma dessas amostras são de areas cercados, com curva de nível, varias coisas no campo que dividem elas, então não são contínuos iguais, temos tamanhos diferentes para cada amostras. Podemos representar esses tamanhos aqui com o tamanho da bolinhas sendo equivalente ao tamanho da amostra.

02

No modelo de Poisson a gente tem:

C_i \sim Poisson(\lambda_i)

Cada amostra vem de uma distribuição de Poisson com média \lambda.
Agora como cada amostra foi feita num pequeno pasto de área A, nosso modelo então é na verdade:

C_i \sim Poisson(A_i \cdot \lambda_i)

Dessa forma, o preditor linear fica assim:

log(A_i \cdot \lambda_i)

E pelas regras de logaritmo, a gente pode separar essa multiplicação dentro do log em dois logs.

log(A_i) + log(\lambda_i)

Agora esse \lambda pode vir de um preditor linear, por exemplo com a variável x sendo o tipo de pasto que coletamos. Então

log(A_i \cdot \lambda_i) = log(A_i) + \alpha + \beta \cdot x_i

Veja que ficou assim porque podemos fazer aquela separação da área e do \lambda, por causa da multiplicação no log.

De forma mais simples, estamos pensando nesse log da área com um coeficiente 1 multiplicando ele, ou seja, a área fazendo o efeito dela proporcionalmente, mas isso pode ser mudado, usando um coeficiente para a área.

log(A_i \cdot \lambda_i) = \beta_0 \cdot log(A_i) + \alpha + \beta \cdot x_i

Agora é so passar isso para a linguagem Bugs da seguinte forma.

 for (i in 1:n) {
    C[i] ~ dpois(lambda[i])
    log(lambda[i]) <- 1 * logA[i] + alpha + beta *x[i]
 }

São duas linhas, a linha que tem a ~ é a parte estocástica, e a linha com <- é a parte determinística, que está sendo somado, além do preditor linear, o offset da área, ja logaritimizado, a gente não ta mandando pro bugs a área real, a gente já calcula o log no R e só mando o valor do log. Então é juntar os dados, configurar o MCMC, fazer a função de inicialização das cadeias e rodar o modelo. E para frisar, veja que a gente ta mandando o termo logA que é o log(A).

# Juntando os dados numa lista para o OpenBugs
bugs.data <- list(C = C, x = as.numeric(x)-1, logA = log(A), n = length(x))

Rodamos o modelo, mesmo com cadeias pequenas ele converge rápido, é um modelo relativamente simples e vemos o seguinte.

> print(out, dig = 3) Inference for Bugs model at "Offset.txt", Current: 3 chains, each with 1100 iterations (first 100 discarded), n.thin = 2 Cumulative: n.sims = 3000 iterations saved mean sd 2.5% 25% 50% 75% 97.5% Rhat n.eff alpha 0.579 0.126 0.323 0.493 0.583 0.669 0.809 1.003 3000 beta 0.966 0.148 0.685 0.861 0.966 1.064 1.255 1.001 3000 deviance 97.335 1.974 95.430 95.960 96.745 98.080 102.902 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 = 2.0 and DIC = 99.3 DIC is an estimate of expected predictive error (lower deviance is better).

Todos os parâmetros com um Rhat perto de 1, e vemos a diferença entre os dois tipos de locais que observamos.
Lembrando que beta aqui é a diferença entre os dois locais, isso pela forma como escrevemos o preditor linear, de uma olhada aqui se não entendeu essa parte. Mas vemos que essa diferença tem média 0.966, mas mais importante que isso o intervalo de confiança de 95% dessa distribuição vai de 0.685 a 1.255, ou seja , não engloba 0 (nenhuma) diferença, alias, com 0.148 de desvio padrão para esse parâmetro, o valor de zero diferença está muitos desvios padrão de distancia da média, temos muita certeza da diferença.

Agora a titulo de curiosidade vamos explorar aqui o que acontece ao ignorarmos o offset, ou seja construir o modelo de forma errada.

> glm.fit.sem.offset <- glm(C ~ x, family = poisson) > summary(glm.fit.sem.offset) Call: glm(formula = C ~ x, family = poisson) Deviance Residuals: Min 1Q Median 3Q Max -1.4298 -0.7816 -0.1419 0.8385 1.7360 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) 1.8245 0.1270 14.367 < 2e-16 *** xPlantação 0.9355 0.1499 6.242 4.31e-10 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 (Dispersion parameter for poisson family taken to be 1) Null deviance: 60.948 on 19 degrees of freedom Residual deviance: 17.615 on 18 degrees of freedom AIC: 103.45 Number of Fisher Scoring iterations: 4

Veja que no modelo bayesiano com offset, o alpha ou intercepto, como quiser chamar, ficou 0.579, sem o offset aqui ele ficou como 1.8245, ou seja muito maior, super-estimamos muito a contagem de coelhos, esse valor vai dar quantos coelhos esperamos achar no pasto e nas plantações, lembrando que nas plantações o lambda vai ser alpha mais beta, ou aqui intercept + xPlantação, então por mais que a diferença até que fica parecida, temos os valores de lambda bem super-estimado.

Mas não confunda com o problema ser modelo bayesiano ou não, se ajustarmos o modelo frequentista temos praticamente os mesmo parametros do modelo bayesiana, com o mesmo erro.

> glm.fit.com.offset <- glm(C ~ x, family = poisson, offset = log(A)) > summary(glm.fit.com.offset) Call: glm(formula = C ~ x, family = poisson, offset = log(A)) Deviance Residuals: Min 1Q Median 3Q Max -1.68068 -0.44061 -0.06695 0.31836 1.69544 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) 0.5811 0.1270 4.575 4.75e-06 *** xPlantação 0.9658 0.1499 6.445 1.16e-10 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 (Dispersion parameter for poisson family taken to be 1) Null deviance: 59.836 on 19 degrees of freedom Residual deviance: 13.535 on 18 degrees of freedom AIC: 99.376 Number of Fisher Scoring iterations: 4

Apesar da grande diferença nas estimativas dos parâmetros, o uso ou não de offset, não é facil detectar que ele esta faltando, temos que ter pensado isso a priori. Se a gente olhar os resíduos de ambos os modelos não vemos nenhuma tendencia, pelo menos eu não vejo.

03

Temos um desvio residual menor no modelo empregando o offset, 13.535 com offset contra 17.615 sem offset, mas a gente nunca saberia disso se não fizer o ajuste com o offset, a novamente, não pegaria a necessidade disso olhando somente um plot dos resíduos sem offset.

Bem, é isso, mais um conceito legal para se entender, e facilitar a vida ao ler artigos e quem saber resolver os problemas de alguém. O repositório recologia esta la. E é isso, até o próximo post.

Referência:
Marc Kéry (2010) Introduction to WinBUGS for Ecologists. Academic Press, Burlington.

### Gerando dados
set.seed(1234)
n.locais <- 10
A <- runif(n = 2*n.locais, 2,5)		# Area de 2 a 5 km2 
x <- gl(n = 2, k = n.locais, labels = c("Pasto", "Plantação"))
preditor.linear <- log(A) + 0.69 +(0.92*(as.numeric(x)-1))
lambda <- exp(preditor.linear)
C <- rpois(n = 2*n.locais, lambda = lambda)
 
#Figura 1
stripchart(C~x,vertical=T,pch=19,method="stack",offset=0.8,at=c(1.2,2.2),frame=F,xlab="Contagem",ylab="Local")
 
#Figura 2
stripchart(C~x,vertical=T,pch=19,method="jitter",offset=0.8,at=c(1,2),frame=F,xlab="Contagem",ylab="Local",cex=A/2)
 
### Análise com OpenBugs
# Definindo o modelo
sink("Offset.txt")
cat("
model {
# Priors
 alpha ~ dnorm(0,0.001)
 beta ~ dnorm(0,0.001)
 
# Likelihood
 for (i in 1:n) {
    C[i] ~ dpois(lambda[i])
    log(lambda[i]) <- 1 * logA[i] + alpha + beta *x[i]
 }
}
",fill=TRUE)
sink()
 
# Juntando os dados numa lista para o OpenBugs
bugs.data <- list(C = C, x = as.numeric(x)-1, logA = log(A), n = length(x))
 
# Função para iniciar as cadeias com valores ao acaso
inits <- function(){ list(alpha=rlnorm(1), beta=rlnorm(1))}
 
# Parametros a estimar
params <- c("alpha","beta")
 
# MCMC settings
nc <- 3    # Número de cadeias
ni <- 1100 # Tamanho de cada cadeia
nb <- 100  # Número de amostras a descartar como burn-in
nt <- 2    # Thinning rate
 
# Inicie o Gibbs-sampler
library(R2OpenBUGS)
 
out <- bugs(data=bugs.data, inits=inits, parameters.to.save=params, model.file="Offset.txt",
            n.thin=nt,n.chains=nc,n.burnin=nb, n.iter=ni)
 
#Resultado
print(out, dig = 3)
 
### Analise usando R
 
glm.fit.sem.offset <- glm(C ~ x, family = poisson)
glm.fit.com.offset <- glm(C ~ x, family = poisson, offset = log(A))
summary(glm.fit.sem.offset)
summary(glm.fit.com.offset)
 
#figura 3
par(mfrow=c(2,1))
plot(resid(glm.fit.sem.offset),ylim=c(-2.5,2),main="Sem Offset",pch=19,ylab="Resíduos",xlab="Amostras",cex=1.3)
plot(resid(glm.fit.com.offset),ylim=c(-2.5,2),main="Com Offset",pch=19,ylab="Resíduos",xlab="Amostras",cex=1.3)

Teste t de Poisson

O titulo de post é um nome bem esquisito, e errado também, teste t usa a estatística t, então não da para ter um teste t de Poisson exatamente. Mas a gente acaba associando alguns nomes a situações.
Teste t para avaliar um fator de dois níveis, anova para um fator de vários níveis, ancova para variáveis continuas e variáveis discretas e assim vai.
A questão é que todos esses métodos estão na verdade dentro de uma categoria, existe uma teoria que unifica todos esses métodos, esses e outros, como aqueles que já comentamos usando distribuição binomial, ou mesmo com distribuição de poisson como vemos aqui são todos casos da mesma “estratégia” em estatística aplicada, se meus termos aqui não são ruins. Mas em 1972, John Nelder e Robert Wedderburn unificaram mais ou menos tudo que falamos ali em cima na ideia geral de Modelos gerais linearizados.. Eles mostraram que todas essas técnicas, representadas como tipos diferentes de analises, incluindo chi-quadrados até que a gente ve muito em génetica, podiam todos serem representados como casos especiais de modelos lineares. Quando eu fiz estatística na faculdade, na graduação, eu nunca fui apresentado a essa ideia que tudo é “modelo linear”, que tudo esta interligado e por isso muitas ideias são comuns, teste t e anova eram coisas distintas, técnicas que se você tem a situação a você aplica isso, se tem a situação b, aplica aquilo e se por ventura tem a situação c você se ferrou, porque não tem um teste pronto pra você. E isso gera o problema, que as vezes as pessoas abandonam a ideia, pressupostos da biologia por exemplo, para forçar os dados serem aplicáveis a um teste e temos que lidar com um monte de problemas decorrentes disso. Além é claro de isso muitas vezes diminuir o avanço cientifico, se é que alguém ainda se importa com isso.

Mas então, um modelo geral linear vai ter mais ou menos três componentes.

1. Uma distribuição estatística, usada para descrever a parte aleatória da resposta do modelo, o eixo y dos gráficos que a gente normalmente faz. Essa é a parte estocástica do sistema que estamos tentado descrever.
2. Uma função de ligação, que é aplicada ao valor esperado, um E(y).
3. E um preditor linear, os a+bx, intercepto e inclinação, as variáveis preditoras que na maioria das vezes estamos interessados. se a função de ligação é a função g, então um g(E(y)). Essa é a parte determinística do modelo, veja que aqui entra qualquer função matemática na verdade, mas essa ai é bem comum.

Nesse sentido a regressão linear que a gente faz é o seguinte.

1. Distribuição: y_i \sim Normal(\mu_i , \sigma^2)
2. Função de ligação: Identidade, \mu_i = E(y_i) que é o preditor linear.
3. Preditor linear: \alpha + \beta \cdot x_i

Se você pensa que nunca fez isso, veja que todas as simulações que fazemos aqui a gente meramente vai do item 3 pro item 1, a gente gera os valores determinísticos, transforma pela função de ligação e depois gera a parte aleatória dos dados, e depois roda uma analise para ver se recuperamos os valores do preditor linear.

Mas e um modelo usando Poisson?

Ai teremos o seguinte:

1. Distribuição: C_i \sim Poisson(\lambda_i)
2. Função de ligação: \log{(\lambda_i)} ou \log({E(C_i)})
3. Preditor linear: \alpha + \beta \cdot x_i

Viu, a mesma ideia, só mudamos a função de ligação, que antes era a identidade, que nada mais é que o mesmo valor do número, agora usamos o log, log na base natural, que é o número do euler la, que a gente ja esbarrou por ele aqui.

Bem agora a analise, primeiro vamos entender os dados.

Vamos supor que queremos estudar o número de aranhas que encontramos em arbustos. Então a gente pega um arbusto inteiro, poe ele num sacão e leva para o laboratório, la vasculhamos ele e contamos quantas aranhas conseguimos achar. Particularmente os arbustos são de uma espécie de planta que ocorre tanto na zona rural como na zona urbana, em terrenos baldios. Mas nos achamos que aranhas não conseguem dispersar para dentro da zona urbana, logo esperamos que elas sejam restritas a áreas rurais, certo? Claro que uma ou outra aranha vai acabar chegando na cidade, ainda mais que estamos falando “aranhas” sem considerar as espécies delas. Mas de modo geral esperamos menos aranhas na cidade, ou mais claramente falando, esperamos que em média, arbustos na cidade tenham menos aranhas, ou seja em média as contagens de aranhas serão menor na cidade quando comparado a zona rural.

Então temos a pergunta, existe diferença no número de aranhas entre a área urbana e a área rural?

Bem sabendo o que queremos saber, fica fácil, é so pegar arbustos na cidade e na área rural e contar a aranhas e ver se existe diferença nas contagens de aranhas em cada área. Mas porque não usar um teste t normal? Porque primeiro aranhas so vem inteiras, não existe meia aranha nem um quarto de aranha, como a distribuição normal é continua, ela vai prever uma variável continua, o que é ruim. Mas além disso, a distribuição normal, se eu ver um número muito baixo de aranhas, algo que deixe a média perto de zero, eu vou prever números negativos de aranhas em arbusto, o que não pode existir, se num local existem menos aranhas, eu tenho que ver mais zeros e não números negativos. Esses dois problemas são solucionados usando a distribuição de Poisson, existem mais coisas a se considerar, mas por agora vamos em frente, outras considerações ficam como cenas dos próximos capítulos.

Então fomos la no mato, coletamos os arbustos e contamos a aranhas, e vamos ver o que conseguimos:

02

Olha ai, pegamos 20 arbustos no total, 10 na zona rural e 10 na zona urbana. Como ja suspeitávamos, na cidade não vimos mais que 4 aranhas por arbusto, enquanto no campo olha la, a gente tinha muito mais, até 8 aranhas. Podemos fazer um box-plot.

01

Certo, são os mesmo dados, mas agora acho que visualmente da para ver outra questão, aqui temos uma variação muito maior dos dados na área rural e isso quebra uma das premissas de modelos lineares em geral, como o teste t, que é homogeniedade de variâncias. Mas não é difícil pensar que se a gente sempre vê em média uma aranha, não tem como ter uma variação igual quando a gente tem muitas, pense que não tem como contar algo para menos que zero então um número pequeno implica de certo modo que só podemos ter uma variação pequena, agora se a gente contar em média 100 aranhas nos arbustos, seria fácil pensar por exemplo que erramos e contamos 90, ou 110, e isso faz um desvio muito grande em relação a uma contagem pequena, o ponto é que o desvio tem uma relação com o quanto a gente conta, ou seja com o \lambda, mas na distribuição de poisson é exatamente assim que funciona, o \lambda determina a média e a variação dos dados, no caso quanto maior a contagem média, maior a variação. Temos como mexer um pouco nessa relação, mas isso é outro post, por agora já da para ter uma noção de outra característica do modelo que estamos fazendo, que é em modelos usando distribuição de poisson, a média e o desvio estão relacionados entre si.

Mas podemos ver as médias das áreas também, para ver como ficaram.

> aggregate(C, by = list(x), FUN = mean) Group.1 x 1 Urbana 2.4 2 Rural 5.3

E a média da área rural é quase o dobro da média na área urbana. Agora podemos tentar definir isso como parâmetros de um modelo linear generalizado usando a distribuição de poisson.

Primeiro temos que definir o modelo na linguagem bugs.

# Definindo o modelo na linguagem Bugs
sink("Poisson.t.txt")
cat("
model {
# Priors
 alpha ~ dnorm(0,0.001)
 beta ~ dnorm(0,0.001)
# Likelihood
 for (i in 1:n) {
    C[i] ~ dpois(lambda[i])
    log(lambda[i]) <- alpha + beta *x[i]
# Ajuste
    residuo[i] <- (C[i] - lambda[i]) / sqrt(lambda[i])
}
}
",fill=TRUE)
sink()

Além do modelo em si, vamos salvar os resíduos.

Agora se pararmos um momento para olhar o likelihood, podemos ver os itens que definimos la em cima.
O modelo é:

# Likelihood
 for (i in 1:n) {
    C[i] ~ dpois(lambda[i])
    log(lambda[i]) <- alpha + beta *x[i]
}

E olha aqui, temos um loop para processar todos as amostras, e dentro dele temos:

1. A distribuição de poisson

C[i] ~ dpois(lambda[i])

2. A função de ligação

log(lambda[i])

3. A parte determinística do modelo.

alpha + beta *x[i]

Esta numa linguagem de computador, mas é exatamente o que definimos la em cima.
Mas agora com o modelo pronto, definimos as configurações do MCMC e rodamos ele.
Com os resultados em mão, antes de fazer qualquer inferência, primeiro de tudo a gente olha se as cadeias convergiram. Existem muitas formas de ver isso, mas primeiro podemos olhar para os valores de Rhat, eles são simples de interpretar, porque é so buscar se estes estão inferiores a 1.1, o valor limite, acima disso temos que pensar bem como prosseguir, mas nesse caso…

03

Todos menores que 1.1.

Podemos fazer um teste rapido, e simplesmente testar se alguém ficou acima de 1.1

which(out$summary[,8] > 1.1)

Certo, outra coisa, podemos olhar agora os resíduos, estes tem que ter uma distribuição normal, se eles não tiverem uma distribuição normal, significaria que alguma coisa muito importante que não estamos levando em conta esta acontecendo, que tem um efeito claramente mudando muito tudo. Se os resíduos estiveram normal, com média 0, quer dizer que estamos errando igualmente tanto para mais quanto para menos.

04

Olha ai, passamos uma linha no zero, e temos mais ou menos a mesma quantidade de pontos pra cima e para baixo, e nenhuma tendência aparente. Podemos olhar também esses resíduos na forma de um histograma.

05

E apesar de umas barras talvez esquisitas, podemos ver também a densidade em vermelho, que é um “histograma” infinito, mais suave, e vemos bem uma distribuição normal, com média perto do zero, marcado em azul na figura.

Existem também testes para testar a normalidade, como o teste de shapiro-wilks.

> shapiro.test(out$mean$residuo) Shapiro-Wilk normality test data: out$mean$residuo W = 0.9374, p-value = 0.2144

Mas a gente tem que tomar cuidado para não entrar num paradoxo, se a gente começar a testar o teste do teste do teste, a gente não vai parar nunca, e nunca voltaremos para a questão que realmente nos interessa.
Mas pelo menos eu ja estou mais que convencido que o modelo teve um ajuste bom, não tem problemas, então poderíamos olhar as conclusões que ele nos traz, quais inferências podemos fazer.

06

E exatamente o \beta é o que nos interessa, ele é a diferença entre a área rural e a área urbana.
Voltando a nossa pergunta inicial, que era se existia diferença no número de aranhas entre a área urbana e a área rural, um valor de \beta=0 significaria que não há diferença. Como em todas as amostras do MCMC, não temos nenhum 0 nesse caso, nem nada menor que zero estamos fortemente propensos a acreditar que existe uma diferença. Podemos calcular o intervalo de confiança de 95%, ou seja, qual intervalo mais comum de valores e falar que este é maior que zero. Poderíamos sei la, calcular um intervalo de 99%, mas aqui o que interessa é que podemos fazer uma afirmação mais ou menos assim.

“Temos 95% de certeza que deve haver uma diferença no número de aranhas presente em arbustos na zona rural e zona urbana.”

Podemos apresentar então quantas aranhas esperamos encontrar em cada local, e a partir daqui é tentar entender melhor porque existe essa diferença.

07

E assim concluímos nosso primeiro glm usando distribuição de Poisson, já que esse post está ficando grande mais. Mas da para ver como os glms abrem uma grande gama de perguntas que podemos abordar. Espero que não tenha falado muita bobeira no caminho, mas modelos lineares usando distribuição normal, binomial e de poisson cobrem a grande maioria das coisas que lemos comumente na literatura científica relacionada a biologia.

Bem o script está disponível no repositório recologia no github, além de aqui em baixo, e acho que é isso, até o próximo post.

Referencias:
Marc Kéry (2010) Introduction to WinBUGS for Ecologists. Academic Press, Burlington.

set.seed(123)
### Gerando os dados
#Gerando o preditor linear, o item 3 do post.
n.arbustos <- 10
x <- gl(n = 2, k = n.arbustos, labels = c("Urbana", "Rural"))
n <- 2*n.arbustos
 
#Aplicando a função de ligação, exponencial é o contrario de logaritimo,
#estamos no caminho contrario aqui, esse é o item 2 do post
lambda <- exp(0.69 + 0.92*(as.numeric(x)-1))
 
#Adicionando a parte estocastica, o item 1 do post
C <- rpois(n = n, lambda = lambda)
 
#Vamos olhar o que criamos
aggregate(C, by = list(x), FUN = mean)
 
boxplot(C ~ x, col = "grey", xlab = "Tipo de ambiente", ylab = "Contagem de aranhas", las = 1,frame=F)
 
stripchart(C ~ x,xlab = "Tipo de ambiente", ylab = "Contagem de aranhas", las = 1,frame=F,
           method="stack",vertical=T,pch=19,col=1,at=c(1.4,2.1),offset = 0.8)
legend("topleft",pch=19,legend="Um arbusto")
 
 
### Analise Bayesiana
# Definindo o modelo na linguagem Bugs
sink("Poisson.t.txt")
cat("
model {
 
# Priors
 alpha ~ dnorm(0,0.001)
 beta ~ dnorm(0,0.001)
 
# Likelihood
 for (i in 1:n) {
    C[i] ~ dpois(lambda[i])
    log(lambda[i]) <- alpha + beta *x[i]
 
# Ajuste
    residuo[i] <- (C[i] - lambda[i]) / sqrt(lambda[i])
 
}
 
}
",fill=TRUE)
sink()
 
# Juntando os dados para mandar para o Openbugs
dados.bugs<- list(C = C, x = as.numeric(x)-1, n = length(x))
 
# Função geradora de valores iniciais
inits <- function(){ list(alpha=rlnorm(1), beta=rlnorm(1))}
inits()
 
# Parametros a estimar
params <- c("alpha", "beta","residuo")
 
# Configurações do MCMC
nc <- 3
ni <- 3000
nb <- 1000
nt <- 2
 
# Iniciando o Gibss Sampler
library(R2OpenBUGS)
out <- bugs(data=dados.bugs, inits=inits, parameters.to.save=params,model.file="Poisson.t.txt",
            n.thin=nt, n.chains=nc, n.burnin=nb,n.iter=ni)
 
 
### Checando convergencia e ajuste do modelo
 
hist(out$summary[,8],col="grey",main="Valores de Rhat")
which(out$summary[,8] > 1.1)
 
plot(out$mean$residuo, las = 1,frame=F,ylab="Resíduos",xlab="Amostras",xlim=c(0,20))
abline(h = 0)
 
hist(out$mean$residuo,xlab="Resíduos",freq = F)
abline(v = 0,lwd=2,col="blue",lty=3)
points(density(out$mean$residuo),type="l",lty=2,lwd=2,col="red")
 
shapiro.test(out$mean$residuo)
 
### Inferencia baseada no modelo que criamos
print(out, dig = 3)
 
#Função do Livro Doing Bayesian Analysis
HDIofMCMC <- function( sampleVec , credMass=0.95 ) {
    sortedPts = sort( sampleVec )
    ciIdxInc = floor( credMass * length( sortedPts ) )
    nCIs = length( sortedPts ) - ciIdxInc
    ciWidth = rep( 0 , nCIs )
    for ( i in 1:nCIs ) {
        ciWidth[ i ] = sortedPts[ i + ciIdxInc ] - sortedPts[ i ]
    }
    HDImin = sortedPts[ which.min( ciWidth ) ]
    HDImax = sortedPts[ which.min( ciWidth ) + ciIdxInc ]
    HDIlim = c( HDImin , HDImax )
    return( HDIlim )
}
 
int95<-HDIofMCMC(out$sims.list$beta)
hist(out$sims.list$beta, col = "grey", las = 1, xlab = "Coeficiente para area rural", main = "",freq=F)
points(density(out$sims.list$beta),type="l",lty=2,lwd=3,col="red")
arrows(int95[1],0.25,int95[2],0.25,col="blue",angle =90,lwd=3,lty=3,code=3)
legend("topleft",legend="Intervalo de confiança de 95%",lwd=3,lty=3,col="blue",bty="n")
 
par(mfrow = c(2,1))
int95<-HDIofMCMC(exp(out$sims.list$alpha))
hist(exp(out$sims.list$alpha), main = "Área urbana",col = "grey", xlab = "", xlim = c(0,10), breaks = 20,freq=F)
points(density(exp(out$sims.list$alpha)),type="l",lty=2,lwd=3,col="red")
arrows(int95[1],0.1,int95[2],0.1,col="blue",angle =90,lwd=3,lty=3,code=3)
legend("topright",legend="Intervalo de confiança de 95%",lwd=3,lty=3,col="blue",bty="n")
 
int95<-HDIofMCMC(exp(out$sims.list$alpha + out$sims.list$beta))
hist(exp(out$sims.list$alpha + out$sims.list$beta), main = "Área rural",freq=F,
     col = "grey", xlab = "Número de aranhas esperado", xlim = c(0,10), breaks = 20)
points(density(exp(out$sims.list$alpha + out$sims.list$beta)),type="l",lty=2,lwd=3,col="red")
arrows(int95[1],0.1,int95[2],0.1,col="blue",angle =90,lwd=3,lty=3,code=3)
 
### Analise Frequentista
poisson.t <- glm(C ~ x, family = poisson)
summary(poisson.t)
anova(poisson.t, test = "Chisq")