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

Noções de notação matemática para o teste t

Opa, após alguma semanas parado apanhando pra conseguir voltar o blog, surge um novo post.

Vamos ver aqui um pouquinho de como é a notação desses modelos que a gente faz, como escrever eles com aquelas letrinhas bonitinhas da matemática que normalmente confunde as pessoas.
A parte boa de tentar entender essa notação, é que ela pode ajudar nosso entendimento dos modelos, assim como melhorar nossa intuição do que deve estar acontecendo.

Bem teste-t é uma análise bem básica, que a gente ja viu milhares de vezes, como aqui e aqui, para ver diferença de médias.
Mas vemos ver de novo ele para os seguintes dados.

massa<-c(6,8,5,7,9,11)
regiao<-factor(c("Campo Grande","Campo Grande","Campo Grande","Campo Grande","Dourados","Dourados"))

Podemos fazer um gráfico rápido para facilitar o entendimento.

01

Então temos um fator, de dois níveis, que são as cidades, campo grande e dourados, e temos a massa de indivíduos de alguma espécie qualquer.

No R, fazer o teste t significa digitar

t.test(massa~regiao,var.equal=T)

E temos como resposta…

Two Sample t-test data: massa by regiao t = -3.0551, df = 4, p-value = 0.03784 alternative hypothesis: true difference in means is not equal to 0 95 percent confidence interval: -6.6808175 -0.3191825 sample estimates: mean in group Campo Grande mean in group Dourados 6.5 10.0

Certo, mas isso na verdade a gente está fazendo um modelo linear, é algo como um caso especial, não um caso especial, simplesmente um caso de modelo linear com 1 fator de 2 níveis, falar um caso significa que queremos ver o teste-t dentro da teoria que usamos para fazer todos os modelo lineares ou modelos gerais linearizados.

Bem um modelo linear no R é feito simplesmente usando o comando lm.

lm(massa~regiao,var.equal=T)
Call: lm(formula = massa ~ regiao) Coefficients: (Intercept) regiaoDourados 6.5 3.5

Bem se a gente calcular a média das massas para as duas cidades a gente vai ver que:

> aggregate(massa,list(regiao),mean) Group.1 x 1 Campo Grande 6.5 2 Dourados 10.0

Agora esse 6.5, que está indicado como intercept da para saber o que é, é a media de campo grande, agora e esse 3.5 de dourados, e aonde estão os valores p?

Primeiro, temos que perceber que antes de tudo temos que estimar parâmetros, para construir um modelo, mas o que diabo é o modelo? Agora quando a gente tentar escrever a notação matemática, vamos conseguir talvez entender melhor, visualizar o quadro maior.

Quando a gente fala modelo, a gente ta falando em fazer uma função matemática, que representa a situação que a gente observa, qual a melhor função matemática para descrever o que a gente ta vendo?
Uma função exata não da para fazer quase sempre, um função probabilística sempre da
O porque não conseguimos fazer uma função matemática perfeita para os dados? Porque existe quase sempre um componente estocástico, um componente chamado acaso, que a gente pode esperar algo, mas não temos como saber exatamente o que. Mesmo que você tenha todos os dados dos mundo, a posição exata de cada atomo, sua carga, toda a informação, ainda sim não da para prever o tempo perfeitamente. Perfeitamente não, mas razoavelmente bem sim.

Ta ficando confuso, mas muita calma nessa hora.

Matematicamente, o que estamos fazendo aqui é uma função assim:

massa_i = \alpha + \beta \cdot regiao_i + \epsilon_i

Poxa, a resposta então é a massa, ali a gente tem a equação de reta, a+bx e um erro, esse sendo representado por esse simbolo epsilon (\epsilon_i).
Os i nas são referentes aos índices, índices do que? Dos casos, das amostras que temos, vamos olhar a massa, no R a gente tem:

> massa [1] 6 8 5 7 9 11 > massa[1] [1] 6 > massa[3] [1] 5

Então a gente tem 6 amostras, o i é o índice, quando o índice = 1, o valor de massa que observamos é 6, quando o índice=3, o valor observado é 5, isso é o i, índice de posição, o mesmo vale para a região.

Beleza, mas a notação completa do modelo é assim;

massa_i = \alpha + \beta \cdot regiao_i + \epsilon_i
\epsilon_i \sim Normal(0,\sigma^2)

Ou as vezes assim:

massa_i = Normal(\alpha + \beta \cdot regiao_i,\sigma^2)

Então a primeira linha é a equação determinística, mas nela ta somado um número que vem de uma distribuição normal com média zero, a segunda linha, no segundo caso a gente so juntou tudo isso numa linha, mas o primeiro jeito eu acho que é mais simples de entender.

Mas como a gente ve isso em relação aos nossos dados? Assim.

  6  = \alpha \cdot 1 + \beta \cdot 0 + \epsilon_1 \\  8  = \alpha \cdot 1 + \beta \cdot 0 + \epsilon_2 \\  5  = \alpha \cdot 1 + \beta \cdot 0 + \epsilon_3 \\  7  = \alpha \cdot 1 + \beta \cdot 0 + \epsilon_4 \\  9  = \alpha \cdot 1 + \beta \cdot 1 + \epsilon_5 \\  11 = \alpha \cdot 1 + \beta \cdot 1 + \epsilon_6 \\

Essencialmente esse é o meu modelo, eu to tentando ver ele em relação a cada amostra de dados que tenho.
Quando a gente estima parâmetros, o que eu quero saber é qual valor eu posso colocar para alpha(\alpha) e para beta(\beta) para solucionar esse monte de equação ai. So que lembre-se que o epsilon (\epsilon) é um valor que vem de uma distribuição normal, uma distribuição de eventos estocásticos, dados ao acaso, então ele pode ser qualquer valor, eu posso colocar qualquer número pra ele, qualquer número mesmo, ou seja, eu sempre vou conseguir solucionar todas essa equações ae em cima, para qualquer valores de alpha e beta. Por exemplo , eu posso colocar o valor de alpha como 100, beta como um milhão e epsilon como -94 e ta resolvido.

 6  = 100 \cdot 1 + 1000000 \cdot 0 + (-94)

Espero não ter errado, mas 100 vezes 1 da 100, dai um milhão vezes zero da zero, e tirando 94 da 6. So que essa solução é ruim, porque esses valores de alpha e beta fazem a gente ter que usar valores muito grandes para epsilon, mais precisamente, muito diferentes para cada equação. Imagina la em baixo quando o um milhão é multiplicado por 1, a gente vai precisar de um epsilon de mais de um milhão la, enquanto usamos 94 na primeira equação, ou seja para conseguir valores assim precisamos de um desvio muitooooo grande para o epsilon, e a pergunta é, quais valores eu ponho em alpha e beta para usar o menor desvio possível no epsilon. Pra variar o mínimo possível os valores que vou ter que colocar no epsilon para satisfazer todas essas equações ai de cima.

Certo, até aqui tudo bem, mas como isso está no R?

Bem o calculo não é feito com aquele monte equação, tem um jeito mais simples de fazer isso no R, que é usando álgebra linear(ou não, dependendo do tamanho do seu ódio para com álgebra linear). Mas a gente usa operações com matrizes.

Fica assim a conta para o nosso exemplo.

 \left(\begin{array}{c} 6 \\ 8 \\ 5 \\ 7 \\ 9 \\ 11 \end{array} \right) = \left(\begin{array}{cc} 1&0 \\ 1&0 \\ 1&0 \\ 1&0 \\ 1&1 \\ 1&1 \end{array} \right) \ast  \left(\begin{array}{c} \alpha \\ \beta \end{array} \right) +  \left(\begin{array}{c} \epsilon_1 \\ \epsilon_2 \\ \epsilon_3 \\ \epsilon_4 \\ \epsilon_5 \\ \epsilon_6 \end{array} \right)

Se a gente resolver essa multiplicação e a soma de matrizes e vetores, a gente vai resolver exatamente aquele conjunto de equaçãozinha la em cima.

Agora é so pensar, antes do igual, é nosso vetor chamado massa:

> massa [1] 6 8 5 7 9 11

Essa é a primeira matriz, antes do sinal de igual, a matriz de 1 coluna e 6 linhas(o número de linhas é o número de amostras) de resposta.

A segunda matriz, é a matriz do modelo. Pra ver ela é so ver como é guardado os dados da região.
Ele é um fator, que tem 2 componentes, os níveis e como eles estão distribuídos.

Com o comando levels a gente ve os níveis no R

> levels(regiao) [1] "Campo Grande" "Dourados" > as.integer(regiao) [1] 1 1 1 1 2 2

E temos números que representam onde eles estão, é assim porque a gente faz conta com número, para nos os nomes tem significado, mas para a álgebra das matrizes não.

Mas é so unir essas duas informações que temos nossa informação devolta

> levels(regiao)[as.integer(regiao)] [1] "Campo Grande" "Campo Grande" "Campo Grande" "Campo Grande" "Dourados" "Dourados"

Bem com essa informação, o R constrói aquela matriz que de 0 e 1 do modelo

> model.matrix(~regiao) (Intercept) regiaoDourados 1 1 0 2 1 0 3 1 0 4 1 0 5 1 1 6 1 1 attr(,"assign") [1] 0 1 attr(,"contrasts") attr(,"contrasts")$regiao [1] "contr.treatment"

E assim ja temos a nossa função, agora é só minimizar os quadrados para obter a resposta. Isso é os outros dois componente la de cima.

O comando lm ta achando aquele vetor que tem o alpha e o beta:

> lm(massa~regiao) Call: lm(formula = massa ~ regiao) Coefficients: (Intercept) regiaoDourados 6.5 3.5

Se você não se ligou ainda, essas são as médias, do que todas são dadas a partir do intercepto, ou seja a média de campo grande é 6.5, e a média de dourados é 10, que é o 6.5+3.5, aqui é quando o beta é multiplicado por 1 na matriz la em cima.

E com os menores resíduos possíveis, resido é os números que temos que atribuir a epsilon, esses aqui:

> resid(lm(massa~regiao)) 1 2 3 4 5 6 -0.5 1.5 -1.5 0.5 -1.0 1.0

Se a gente substituir esses valores la, temos a solução de todas as equações, usando o menor desvio possível.

Bem, então o nosso exemplo é assim:

 \left(\begin{array}{c} 6 \\ 8 \\ 5 \\ 7 \\ 9 \\ 11 \end{array} \right) = \left(\begin{array}{cc} 1&0 \\ 1&0 \\ 1&0 \\ 1&0 \\ 1&1 \\ 1&1 \end{array} \right) \ast  \left(\begin{array}{c} 6.5 \\ 3.5 \end{array} \right) +  \left(\begin{array}{c} -0.5 \\ 1.5 \\ -1.5 \\ 0.5 \\ -1 \\ 1 \end{array} \right)

Essa operação ai, é exatamente esse conjunto de equações.

  6  = 6.5 \cdot 1 + 3.5 \cdot 0 + -0.5 \\  8  = 6.5 \cdot 1 + 3.5 \cdot 0 + 1.5 \\  5  = 6.5 \cdot 1 + 3.5 \cdot 0 + -1.5 \\  7  = 6.5 \cdot 1 + 3.5 \cdot 0 + 0.5 \\  9  = 6.5 \cdot 1 + 3.5 \cdot 1 + -1 \\  11 = 6.5 \cdot 1 + 3.5 \cdot 1 + 1 \\

Por exemplo, o primeiro caso, 6.5 vezes 1 da 6.5, mais 3.5 vezes zero da zero, e continuamos com 6.5, tirando -0.5 da 6, que é o nosso valor observado.

E qual diabos é a função ai? Podemos escrever ela assim.

   f(x) = \left\{       \begin{array}{lr}         6.5 + \epsilon & : x = 1\\         6.5 + 3.5 + \epsilon & : x = 2       \end{array}     \right.

Onde x é nada mais que a região. E lembre-se que a gente tem a região representada assim

> regiao [1] Campo Grande Campo Grande Campo Grande Campo Grande Dourados Dourados Levels: Campo Grande Dourados

Mas o computador le assim:

> levels(regiao) [1] "Campo Grande" "Dourados" > as.integer(regiao) [1] 1 1 1 1 2 2

E dai vem o 1 e 2 da equação la em cima, x so pode assumir esses dois valores para essa equação. Agora imagine que para generalizar aquela equação para 3 níveis, como numa anova, é so adicionar uma linha, simples certo?

Agora da para visualizar como o teste-t é um caso especial de anova com apenas 2 níveis.

Bem vamos olhar mais uma vez o resultado de um modelo linear no R

> summary(lm(massa~regiao)) Call: lm(formula = massa ~ regiao) Residuals: 1 2 3 4 5 6 -0.5 1.5 -1.5 0.5 -1.0 1.0 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 6.5000 0.6614 9.827 0.000601 *** regiaoDourados 3.5000 1.1456 3.055 0.037841 * --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 1.323 on 4 degrees of freedom Multiple R-squared: 0.7, Adjusted R-squared: 0.625 F-statistic: 9.333 on 1 and 4 DF, p-value: 0.03784

Bem agora a gente sabe o que são aqueles resíduos, os coeficientes ou parâmetros estimados, e aquele itemzinho chamado “Residual standard error” é o desvio usado no Epsilon.

> sqrt(deviance(lm(massa~regiao))/df.residual(lm(massa~regiao))) [1] 1.322876

A menor quantidade de erro para satisfazer o modelo com esses parâmetros ae.

Assim se não houve-se diferença na massa entre dourados e campo grande, esperaríamos um valor de beta próximo de zero. A nossa equação ficou.

   f(x) = \left\{       \begin{array}{lr}         6.5 + \epsilon & : x = 1\\         6.5 + 3.5 + \epsilon & : x = 2       \end{array}     \right.

Deveria ter a mesma média para x=1 e x=2, o termo de diferença, que é o beta, deveria ser zero.
Como não é, é porque tem essa diferença de 3.5, lembrando que é uma diferença de mais ou menos 3.5, não exatamente, porque o epsilon pode assumir qualquer valor.

Bem, ficamos por aqui, por enquanto. Mas acho que agora talvez esses modelos façam mais sentido para todo mundo. Eu me pergunto porque eu nunca vi uma aula as pessoas mostrando esses modelos desse jeito, acho que eu entenderia melhor, mas ta bom, ta ruim mas ta bom.

Bem apesar do pau no site, o repositório recologia sempre esteve la. E é isso, até o próximo post.

massa<-c(6,8,5,7,9,11)
regiao<-factor(c("Campo Grande","Campo Grande","Campo Grande","Campo Grande","Dourados","Dourados"))
 
plot(massa~regiao,frame=F)
 
regiao
 
levels(regiao)
 
as.integer(regiao)
 
levels(regiao)[as.integer(regiao)]
 
t.test(massa~regiao,var.equal=T)
 
lm(massa~regiao)
 
resid(lm(massa~regiao))
 
model.matrix(~regiao)
 
summary(lm(massa~regiao))
 
model.matrix(~regiao-1)
 
summary(lm(massa~regiao-1))
 
aggregate(massa,list(regiao),mean)
 
sqrt(deviance(lm(massa~regiao))/df.residual(lm(massa~regiao)))

Existe mais de uma maneira de se fazer a mesma coisa, um exemplo com o modelo de Cinética enzimática de Michaelis-Menten

Bem, tem um monte de gente que diz frases legais sobre o R. Essa é uma:

This is R. There is no if. Only how.

O autor é o “Simon ‘Yoda’ Blomberg, R-help (April 2005)”

Na minha tradução usando meu inglês nórdico fica…

No R não existe se é possível, apenas como é possível.

Certo, mas vamos olhar uma coisa eu acho bem legal. O modelo do Michaelis-Menten, eu particularmente nunca fui bom de química, além de ser ruim em muitas outras disciplinas, mas eu acho esse modelo legal porque foi um dos primeiros modelos não lineares que eu vi, que fez sentido, que eu entendi de verdade, eu acho.
Basicamente ele relaciona como a velocidade que uma reação acontece com a concentração de uma enzima. Se você tiver um monte de reagente numa solução, ela ocorre devagarinho, mas a gente pode adicionar uma enzina para aumentar a velocidade, lembrando que uma enzina não é consumida na reação. Mas se a gente adicionar so um pouquinho de enzina, a gente aumenta um pouquinho, se a gente adiciona bastante, a gente aumenta mais a velocidade da reação, mas até um limite. So que no começo adicionar enzimas, aumentar a concentração desta na solução, aumenta bem a velocidade da reação, mas conforme a gente coloca mais e mais, a solução satura, o efeito de aumentar a velocidade da reação satura, ou seja tem um limite, um teto de aumento de velocidade, acima disso não adianta mais adicionar enzina que não da em nada.
Essa idéia pode ser descrita na forma de uma equação, o modelo de Michaelis-Menten, não tenho certeza se meu exemplo é muito bom, mas a idéia é por ai.

A modelo é o seguinte:

 {Vel} = {{Vmax\cdot Con}\over{Km+Con}}

Certo, então a gente tem duas constantes, Vmax e Km, Vmax é a maior velocidade de reação possível, enquanto Km é a constante que diz o quanto o aumento da concentração de algo, uma enzima por exemplo, aumenta a velocidade da reação, até o limite de Vmax.

Vamos aquele velho esquema de de gerar alguns dados, onde Vmax=3.4 e Km=0.4, os mesmo valores do artigo do wikipedia sobre o modelo.
Com os dados em mão, vamos fazer um gráfico para dar uma olhada.

01

Certo, primeiro vamos olhar aquelas características que falamos, Vmax é o teto do modelo, a tendencia máxima da reação. O segundo parâmetro, Km, é o momento também onde a velocidade da reação é metade do Vmax, algumas características legais.
Outra coisa é que esse modelo não precisa só ser pensado em nível de reações químicas com enzimas, por exemplo, uma vez eu tentei usar o modelo de Michaelis-Menten para relacionar o número de interações de parasitas com o número de estudos deles, ou seja quanto mais estudos, em mais hospedeiros esses os parasitas eram descobertos, só que isso tinha um limite, chega uma hora que não interessa quantos estudos mais sejam feitos, a gente ja sabe quem uma determinada espécie de parasita ataca, logo há uma saturação, e mais estudos dificilmente vão descobrir hospedeiros novos, mas se ha poucos estudos para uma especie, é bem provável que novas espécies de hospedeiros sejam descobertas. Pensando desse jeito, o Vmax é o total de espécies de hospedeiros que uma espécie de parasita pode parasitar e Km seria o quanto um estudo contribui para a descoberta do total de hospedeiros possíveis de um parasita, sendo que essa contribuição satura, diminui conforme a gente acumula mais estudos.

Certo mas chega de justificar utilidades, temos dados, temos um modelo e queremos ajustar valores para os parâmetros, para as constantes dele.
De cara a gente vê que o modelo não é linear, não é uma linha, logo a gente precisa ajustar um modelo não linear. Vamos fazer isso usando um modelo não linear, no R usando o comando nls, que serve para ajustar modelos não lineares.
Bem ele tem várias coisas a se estudar, mas a principio vamos simplesmente ajustar o modelo, ele não é muito diferente de usar o comando lm. A gente diz a formula, fornece os dados e a função faz o resto. Mas como a função otimiza os valores das constantes, a gente tem que partir de algum ponto, ou seja tem que fornecer um ponto de inicio para começar a otimizar. Vamos fornecer valores quaisquer sem pensar muito, mas existem formas elegantes de conseguir valores iniciais.

nls.menten<-nls(Vel~(Vmax * Con)/(Km + Con),
                data=dados,start=list(Vmax=1,Km=1),algorithm="default")
summary(nls.menten)

Notem que o modelo que a gente usou na formula:

Vel~(Vmax * Con)/(Km + Con)

É exatamente o modelo la de cima, so que escrito do jeito do R, e usamos Con como nome da concentração, além disso fornecemos os dados, o start é uma lista com valores iniciais para todas as constantes, temos duas aqui Vmax e Km, e o algorithm é a estrategia de otimizar os valores, a gente ta usando o default, mas existe varias maneiras de otimizar valores, no nls por padrão vem com três, mas existe outras ainda, mas isso é outra historia.

Mas continuando, temos como resposta:

Formula: Vel ~ (Vmax * Con)/(Km + Con) Parameters: Estimate Std. Error t value Pr(>|t|) Vmax 3.48185 0.03615 96.32 <2e-16 *** Km 0.44281 0.02307 19.20 <2e-16 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 0.07808 on 28 degrees of freedom Number of iterations to convergence: 7 Achieved convergence tolerance: 3.189e-07

Olha ae, as constantes estão até com o nome original da formula, bem parecido com o que usamos para gerar os dados, que foi Vmax=3.4 e Km=0.04. O parâmetros convergiram legal para esses valores, e não deveria ser diferente, ja que geramos esses dados a partir exatamente desse modelo e esses valores de constantes. Se considerar os desvios, os valores originais estão dentro do intervalo de confiança de 95%.
Podemos colocar esse modelo no gráfico, junto aos dados e ver como ele fica.

02

Ta mas agora vamos olhar uma coisa interessante. Podemos transformar esse modelo, os dados para um modelo linear, uma reta.
Difícil pensar né, como uma curva dessas pode virar uma reta. Mas tudo o que a gente tem que fazer é inverter o valor de concentração e de Velocidade da reação.

Primeiro vamos fazer um gráfico com 1/Vel e 1/Con e ver como fica, para avaliar isso visualmente.

03

Olha ae que louco, agora é so fazer uma regressão ai para estimar a equação da reta.
Mas vamos fazer o seguinte, vamos usar um glm, que é modelo linear genearilizado em português, eu acho, ou generalized linear model em inglês.

Bem, ajustar isso no r também é bem simples, basta usar a função glm(), da seguinte forma.

glm.menten<-glm(Vel ~c(Con^-1), data = dados,family = gaussian(link = "inverse"))
summary(glm.menten)

Agora você deveria perguntar, ue mas vc não vai inverter as coisas? Eu estou invertendo a concentração aqui “c(Con^-1)”.
Primeira coisa, eles esta dentro do concatenar “c()” para nao enviar ele (^1) para a formula, o valor para ser invertido, e a formula continuar igual, a+bx.
Outra coisa, olhem que eu estou elevando ela a -1, se você não lembra de exponenciação, de uma revisada rápida no wikipedia, mas essa regra é a seguinte…

{a^{-1}} = {1\over{a}}

A ultima questão é, porque isso só foi feito para a concentração (con) e não para a resposta Vel?

A resposta esta aqui “family = gaussian(link = “inverse”)”, todo glm usa uma distribuição e uma função de ligação, uma das muitas possibilidades da função de ligação é a inversa, ou seja 1/valor, ou seja aqui eu defini que usaremos a distribuição Gaussiana, que é a mesma coisa que distribuição normal que até hoje eu acho meio confuso os muitos nomes dela, mas então, aqui falamos que queremos que a resposta seja invertida, a função de ligação é a inversa, ou seja a função vai inverter os dados de Vel, o glm ta fazendo isso para a gente.

Mas sabendo que ta tudo invertido, o resultado fica.

Call: glm(formula = Vel ~ c(Con^-1), family = gaussian(link = "inverse"), data = dados) Deviance Residuals: Min 1Q Median 3Q Max -0.14465 -0.05057 -0.02119 0.05772 0.13804 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 0.287204 0.002982 96.32 <2e-16 *** c(Con^-1) 0.127178 0.005528 23.01 <2e-16 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 (Dispersion parameter for gaussian family taken to be 0.006095741) Null deviance: 15.60591 on 29 degrees of freedom Residual deviance: 0.17068 on 28 degrees of freedom AIC: -63.938 Number of Fisher Scoring iterations: 4

Agora temos uma equação de reta, a+bx, intercepto e inclinação, ou como queira chamar. Mas com uma rápida batida de olho, a gente vê que esses valores não se parecem nem um pouco com Vmax ou Km que usamos.
Mas se a gente olhar o ajuste deles aos dados…

04

Fica muito bom. Mas o que diabos significa esses valores de intercepto e inclinação aqui?

Aqui a jogada é a seguinte, os valores foram invertidos certo? Então vamos relembrar o modelo original.

 {Vel} = {{Vmax \cdot Con} \over {Km + Con}}

So que Vel ta escrito na forma de  Vel^{-1} que é o mesmo que  1\over{Vel} .

Como estamos olhando para uma igualdade, temos que fazer a mesma coisa dos dois lados certo? Logo…

 {Vel^{-1}} = ({{Vmax \cdot Con} \over {Km + Con}})^{-1}

Que é o mesmo que:

 {1\over{Vel}} =  {{Km + Con} \over {Vmax \cdot Con}}

Certo, nenhuma algebra muito intensa aqui, a gente so inverteu os valores dos dois lados, so que do outro lado como ja era um fração, ela inverteu.

Mas agora é uma soma em cima, ou seja:

 {1\over{Vel}} =  {{Km} \over {Vmax \cdot Con}} + {{Con} \over {Vmax \cdot Con}}

Para facilitar vamos somente mudar as coisas de lugar:

 {1\over{Vel}} =  {{Con} \over {Vmax \cdot Con}} + {{Km} \over {Vmax \cdot Con}}

E vamos cortar aquele Con dividido por Con

 {1\over{Vel}} =  {1 \over {Vmax}} + {{Km} \over {Vmax \cdot Con}}

E do outro lado o que a gente pode ver? Bem é a razão de Km por Vmax que multiplica o inverso da Concentração.

 {1\over{Vel}} =  {1 \over {Vmax}} + {{Km} \over {Vmax}}\cdot{1\over{Con}}

Agora com o que isso te parace, imagine que 1/con é um número x qualquer, num é uma coisa + outra coisa que multiplica x?
Isso não é a+bx, ou intercepto + inclinação?
Vamos substituir ali o intercepto e inclinação por alpha e beta para a gente ver melhor.

 {1\over{Vel}} =  \alpha + \beta\cdot{1\over{Con}}

Fechou, é uma equação da reta, simples assim.
Agora a gente sabe o que fazer com os coeficientes la em cima.
So pensar pelo que substituimos aqui o alpha(\alpha) e o beta(\beta).

O alpha é:
\alpha={1 \over {Vmax}}

E o beta é:
\beta= {{Km} \over {Vmax}}

Então ali no modelo glm a gente tem que fazer o que para ver o Km e o Vmax?

Vmax vai ser:

Vmax={1 \over {\alpha}}

E Km vai ser:

Km = \beta \cdot {Vmax}

Que ja temos, mas mas podemos pensar também como

Km = \beta \cdot {\alpha^{-1}}

Agora no R é so a gente pegar os valores de coeficientes do modelo glm, com a função coef e fazer essas continhas, que para Vmax e Km fica:

coef(glm.menten)[1]^-1
coef(glm.menten)[2]*(coef(glm.menten)[1]^-1)

Isso da:

> coef(glm.menten)[1]^-1 (Intercept) 3.481847 > coef(glm.menten)[2]*(coef(glm.menten)[1]^-1) c(Con^-1) 0.4428143

Voala, quase como mágica, os mesmos valores estimados com o modelo não linear, usando nls.
a gente pode ainda agora ver os valores, usando o inverso, que entendemos o porque disso agora, e ver como fica os dois modelos ajustados com nls e glm aos dados.

05

E são a mesma coisa.
Muito legal né. Mas o que eu queria mostrar de verdade é como as vezes vemos uma pessoa trabalhando de um jeito, outra trabalhando de outro, mas as vezes as pessoas estão fazendo a mesma coisa. Nesse caso a gente ajustou o mesmo modelo, com uma função não linear, e depois com uma função linear que a gente conhece bem, o a+bx (conhecida tambem como \alpha +\beta \cdot x), só que é preciso pensar um pouco para visualizar onde foram parar os parâmetros que temos interesse no caso da função linear. Mas de qualquer forma ajustamos exatamente a mesma coisa. Fizemos a mesma coisa de duas formas completamente diferentes, para chegar ao mesmo resultado..

Eu acho um pouco legal pensar também, como algumas vezes, a gente pode estar advogando um modelo contra o outro, e eles não podem ser algo assim, simplesmente a mesma coisa escrita de outra forma.

Foi mais uma curiosidade, mas serviu para ver também como uma função de ligação, como simplesmente inverter uma função faz o milagre do enrretamento das coisas em alguns casos.

Eu vi essa curiosidade numa apostilinha de glm, já que nunca teria essa sacada na vida hehehe, no fim do texto eu deixo o link e a referência, mas acho que é apenas um curso de R sobre GLM, nada que de para citar, mas é legal a apostila.

Bem o script vai estar ainda no final aqui, mas eu vou deixar ele la no repositório Recologia do Github, junto com os outros, fica um lugar fácil para achar somente o script, como as vezes eu mesmo procuro coisas minhas para copiar ou lembrar comandos, fica mais fácil de achar tudo.

Bem é isso, até o próximo post.

Referência:

Ulrich Halekoh, Søren Højsgaard 2008 Generalized Linear Models (GLM) 32 pp.

#Gerando dados
set.seed(21)
m_menten <- function(C,Vmax,Km) { return( (Vmax * C)/( Km + C) ) }
con<-runif(30,0,4)
v<-m_menten(con,Vmax=3.4,Km=0.4)
v<-v+rnorm(30,0,0.08)
dados<-data.frame(Con=con,Vel=v)
 
#Figura 1, Caracteristicas do modelo
plot(Vel~Con,frame=F,ylim=c(0,4),xlim=c(0,4.5),ylab="Velocidade",xlab="Concentração",pch=19,data=dados)
abline(h=3.4,lty=2,col="red")
text(4,3.5,"Vmax")
abline(h=0.5*(3.4),lty=2,col="red")
text(4,(0.5*(3.4))-0.1,expression("1/2 Vmax"))
lines(x=c(0.4,0.4),y=c(0,0.5*(3.4)),lty=2,col="red")
text(0.55,1,"Km")
 
#Ajustando um modelo não linear para os dados
nls.menten<-nls(Vel~(Vmax * Con)/(Km + Con),data=dados,start=list(Vmax=1,Km=1),algorithm="default")
summary(nls.menten)
 
#Figura 2 - ajuste modelo não linear
plot(Vel~Con,frame=F,ylim=c(0,4),xlim=c(0,6),ylab="Velocidade",xlab="Concentração",pch=19,data=dados)
points(seq(0,6,0.01),predict(nls.menten,newdata=data.frame(Con=seq(0,6,0.01))),type="l",lty=2,col="blue")
legend("topleft",lty=2,col="blue",legend="Ajuste com nls()",bty="n")
 
#Figura 3.
plot(c(1/Vel)~c(1/Con),frame=F,ylim=c(0,1),xlim=c(0,6),ylab="1/Velocidade",xlab="1/Concentração",pch=19,data=dados)
 
#Ajustando um modelo Geral Linear
glm.menten<-glm(Vel ~c(Con^-1), data = dados,family = gaussian(link = "inverse"))
summary(glm.menten)
 
#Retornando os valores de Vmax e Km
coef(glm.menten)[1]^-1
coef(glm.menten)[2]*(coef(glm.menten)[1]^-1)
 
#Figura 4 - Ajuste com modelo linear
plot(c(1/Vel)~c(1/Con),frame=F,ylim=c(0,1),xlim=c(0,6),ylab="1/Velocidade",xlab="1/Concentração",pch=19,data=dados)
points(c(seq(0.1,6,0.01)),predict(glm.menten,newdata=data.frame(Con=c(1/seq(0.1,6,0.01)))),type="l",lty=2,col="green")
legend("topleft",lty=2,col="green",legend="Ajuste com glm()",bty="n")
 
 
#Figura 5 - Comparando os dois ajustes
plot(Vel~Con,frame=F,ylim=c(0,4),xlim=c(0,6),ylab="Velocidade",xlab="Concentração",pch=19,data=dados)
points(seq(0,6,0.01),predict(nls.menten,newdata=data.frame(Con=seq(0,6,0.01))),type="l",lty=2,col="blue")
points(c(seq(0.1,6,0.01))^-1,predict(glm.menten,newdata=data.frame(Con=c(1/seq(0.1,6,0.01))))^-1,type="l",lty=2,col="green")
legend("topleft",lty=2,col=c("blue","green"),legend=c("Ajuste com nls","Ajuste com glm"),bty="n")

Crescimento populacional de aves e algumas simulações.

Bem já tivemos alguns posts sobre crescimento populacional por aqui. Muita teoria, mas sempre eu ficava com dúvida quando e como esses modelos seriam confrontados com dados reais? O que adianta saber esse monte de formula se são inúteis na vida real?
Bem vamos ver um pequeno exemplo que peguei do livro “Primer of Ecology with R” do M. Henry H. Stevens. Ele tem cara de ser um cara legal hehehe.

Primeiro temos que instalar o pacote primer, que ele fez junto do livro para servir de apoio para o livro, que vem com várias funções legais, dados e outras coisas para explorarmos. Mas por agora vamos olhar o dataset chamado sparrows que vem no pacote.

01

O dados consistem da contagem de passarinhos ao longo de alguns anos. Bem como vimos o crescimento da população pode ser medido como o tamanho da população de um ano dividido pela população do ano anterior, ja que o lambda vai ser a divisão.

Vamos lembrar que podemos descrever o crescimento discreto como:

 N_{t+1}= N_t\cdot(1+r)

E esse 1+r é exatamente o

 N_{t+1}= N_t\cdot\lambda

Assim se ele é maior que 1, a população ta crescendo, se ele é menor que 1, a população esta diminuindo e caso ele seja exatamente 1, a população não mudaria de um ano para o outro.

No parte inferior do gráfico anterior estamos vendo exatamente isso, o \lambda da mudança entre os anos, mas ele parece variar bastante de um ano para o outro.

Mas se esses são os valores de \lambda, e eles se repetem, poderiamos pensar que o próximo ano, seria a população no último ano vezes um valor parecido com algum desses ae. Então podemos tentar ter uma idéia de como será a população no próximo ano se sorteamos um desses valores de crescimento e multiplicar última contagem da população.

Mas não precisamos parar por ai, podemos pegar o resultado dessa simulação para o ano seguinte e fazer denovo, sortear outro valor de crescimento e ver o que acontece, e denovo e denovo, umas 50 vezes por exemplo.

Então teríamos isso:

02

Uma previsão de como estará a população daqui a 50 anos, mas essa é apenas uma simulação dependendo da nossa sorte de ter pego somente número altos ou azar de pego muitos números baixos, veríamos a população subir ou descer, então deveríamos fazer isso denovo para ver o que vai acontecer, repetir o mesmo processo outra vez.
Alias, como as contas residem nas costas do computador, vamos fazer isso 10 vezes, e olhar um gráfico como o acima para ver o que esperar.

03

Agora vejam como aconteceu um caso, onde a população cresceu muito, muito mesmo, a linha la em cima, mas notem que ela esta la isolada, a maioria esta dentro de um mesmo patamar, mas temos também outra simulação que parece que decaiu muito, culminando na extinção da população. Então temos vários cenários possíveis, a população ficando gigante, a população extinguindo, mas em media não é isso que vemos, vemos algo entre 10 e mil indivíduos como expectativa.
Mas olhômetro tem seu limites, então o que podemos fazer é realmente muitas, mas muitas simulações mesmo e ver em média, qual o tamanho da população devera ter daqui 50 anos.

04

Então utilizamos uma função para fazer esse serviço e geramos 1000 simulações, agora é so pegar o último valor da simulação, o tamanho da população que esperamos para daqui 50 anos e olhar. Vamos fazer um histograma para entender melhor o resultado.
Bem primeiro vemos que olhar diretamente o histograma é ruim, ja que existem alguns valores gigantes, lembram do primeiro gráfico que tinha um valor de lambda de mais ou menos 3? Imagina que esse cara sai várias vezes na simulação, a população ficaria enorme, e esse são os extremos que estamos vendo. Isso faz o histograma ficar ruim de ver daquele jeito.

Mas podemos usar a escala logaritima. So uma coisa é que aqui usamos log do tamanho da população + 1, isso porque logo de 0 não existe, mas log de 1 é 0, então não temos problemas com os zeros dos dados e conseguimos ver como esta a população numa escala logaritima.

Agora vemos algo como uma distribuição normal, mas alguns valores gigantes, so que esses são poucos. Aqui podemos calcular o intervalo de confiança de 95% dos dados e colocar uma linhas ali no histograma identificando ele, o intervalo de confiança de 95%. Então esperamos que daqui 50 anos a população esteja por ai. Vemos muitos casos onde a população foi para zero, isso pode ser perigoso, isso quer dizer que a qualquer momento que a espécie baixar muito a população deveríamos preocupar, poderíamos nessa simulação ainda ver por exemplo quantos anos para frente começam a aparecer muitos casos de extinção da espécie, para ter uma idéia o intervalo mínimo que deveríamos ir olhar como a espécie esta.

Agora fizemos uma simulação, considerando o crescimento como um efeito estocástico, mas qual a confiança nessas simulações?

Bem podemos tentar olhar como é a distribuição dos crescimentos que usamos na simulação, e melhor que isso, podemos usar uma distribuição t, que leva em consideração o número de valores que temos, para ver se uma distribuição t pode descrever bem esses dados.
Então calculamos o quantiles e determinamos, usando a distribuição t pelos motivos que colocamos nesse post, o intervalo de confiança do que deve ser o lambda.

05

Mas quando fazemos um Q_Qplot vemos que temos um valor muito grande e dois valores muitos pequenos, que parecem não se enquadrar na distribuição, talvez outliers, e esses podem ter grande impacto na simulação, tanto gerando os valores extremos como aumentando em demasia os casos em que a população foi para zero, lembrando que uma vez que a população chega a zero, zero vezes qualquer número é sempre zero, então ela fica zero para o resto da simulação. Então temos que olhar com cuidado esses resultados.

Mas ainda sim temos aqui uma previsão nua do que deve acontecer nos próximos anos dessa espécie. Por um lado podemos encarar essa simulação como irreal, já que não consideramos um limite para o crescimento, como no caso do crescimento logístico, ou competição ou qualquer outra coisa. Mas essa simulação feita aqui pode ser legal exatamente por isso, estamos pensando o vai acontecer sem fazer muitas premissas, assumindo que o mundo é extremamente simples, e temos um primeiro resultado, uma primeira perspectiva para o futuro. E essa medida aqui adquirida pode ainda ser confrontada com o que se conhece da espécie, o que achamos que pode acontecer quanto a interações bióticas e como a espécie deveria estar, Algo como quando usamos o teorema de Hardy-Weinberg, onde o interessante não é o resultado em sim, mas porque o que observamos é diferente desse resultado e como é essa diferença, e isso nos traz muita informação.

Referencia: Stevens, M. H. H. A Primer of Ecology with R. UseR! Series, Springer (2009)

#####################################################################
 
set.seed(3)
library(primer)
data(sparrows)
names(sparrows)
 
#Figura 1
jpeg("01.jpg")
layout(matrix(1:2, nc=1))
plot(Count ~ Year, type="b",data=sparrows,ylab="Contagem",xlab="Ano"
,ylim=c(0,125),frame.plot=F)
obs.R <- sparrows$Count[-1]/sparrows$Count[-length(sparrows$Count)]
plot(obs.R ~ sparrows$Year[-length(sparrows$Count)],ylab="Crescimento",
xlab="Ano",ylim=c(0,3),frame.plot=F)
abline(h=1, lty=3)
dev.off()
 
#Primeira Simulação
years<-50
sim.Rs<-sample(x=obs.R,size=years,replace=TRUE)
output<-numeric(years+1)
output[1]<-sparrows$Count[sparrows$Year==max(sparrows$Year)]
 
for(t in 1:years ) {
output[t+1]<-output[t] * sim.Rs[t]
}
 
plot(0:years, output, type="l",xlab="Anos",ylab="Contagem",main="Simulação",
frame.plot=F)
 
#10 Simulações
sims<-10
sim.RM<-matrix(sample(obs.R, sims * years, replace=TRUE),
nrow=years, ncol=sims)
output[1]<-sparrows$Count[sparrows$Year==max(sparrows$Year)]
outmat<-sapply(1:sims,function(i) {
for( t in 1:years ) {
output[t+1] <- output[t] * sim.RM[t,i]
}
output
}
)
 
options(scipen=10)
matplot(0:years,outmat, type="l",xlab="Anos",ylab="Contagem",
frame.plot=F,log="y")
 
#Função de fazer Simulações
PopSim <- function(Rs, N0, years=50, sims=10) {
sim.RM<-matrix( sample(Rs,size=sims*years,replace=TRUE),nrow=years,
ncol=sims)
output<-numeric(years+1)
output[1]<-N0
outmat<-sapply(1:sims, function(i) {
for( t in 1:years ) {
output[t+1]<-round(output[t]*sim.RM[t,i],0)
}
output
})
return(outmat)
}
 
#Agora Gerando 1000 simulações
output<-PopSim(Rs=obs.R, N0=43, sims=1000)
dim(output)
 
#Expectativa da população daqui 50 anos
N.2053<-output[51,]
summary(N.2053, digits=6)
quantile(N.2053, prob=c(0.0275, .975) )
 
#Figura 4, histogramas
par(mfrow=c(2,1))
hist(N.2053, main="N",ylab="Frequencia")
hist(log10(N.2053+1), main="log(N+1)",ylab="Frequencia" )
abline(v=log10(quantile(N.2053, prob=c(0.0275, .975) )+1), lty=3)
 
#Calculando o intervalo de confiança com distribuição t
logOR<-log(obs.R)
n<-length(logOR)
t.quantiles<-qt(c(0.025, 0.975),df=n-1)
 
se<-sqrt(var(logOR)/n)
CLs95<-mean(logOR)+t.quantiles*se
 
R.limits<-exp(CLs95)
R.limits
 
N.Final.95<-sparrows$Count[sparrows$Year==max(sparrows$Year)]*R.limits^50
round(N.Final.95)
 
#Figura 5 QQplot
qqplot(qt(ppoints(n), df=n-1), scale(logOR))
qqline(scale(logOR))
 
#########################################################################

Distribuição Normal, distribuição t e uma interessante propriedade.

Bem, a gente vive ajustando modelos aqui e sempre há uma preocupação com número de amostras, número de pontos, graus de liberdade e esse tipo de paranóia. Mas recentemente eu vi algo muito interessante relacionado a distribuição t e decidi olhar uma coisa que mostra como essas são preocupações importantes.

Vamos simular alguns dados normais com a função rnorm do R, para quem não sabe, da para a gente fornecer um valor qualquer de média e desvio padrão que o r “inventa” vários valores que respeitam essas características que desejamos em uma distribuição, no caso a normal aqui.
Então eu vou simular aqui 30 valores com média 15 e desvio padrão de 3.

> medida<-rnorm(30,15,3) > medida [1] 10.558297 19.731508 12.129767 12.239984 9.007074 14.183112 14.053954 13.115234 14.680608 16.284044 [11] 12.666841 11.118353 12.661300 15.035855 14.542751 12.889607 18.566637 16.021537 16.520905 14.120085 [21] 15.670924 21.021604 18.035937 14.092622 11.924265 14.197846 14.402683 15.393368 15.437400 16.086194

Tem 30 valores ae em cima, podem contar, e para facilitar olhar a distribuição deles, vamos olhar um histograma desses valores.

01

Uma propriedade interessante é a gente poder centralizar esses valores onde a gente quiser, nesse caso o centro deles é a média, mas se a gente diminuir eles pela média e dividir o resultado disso pelo desvio padrão, a gente vai ter que o centro vai ser a 0 e 95% dos valores estarão entre 2 desvios padrão. Então fazendo essa conta temos o seguinte.

> media<-mean(medida) > desvio<-sd(medida) > medida.ajus<-(medida-media)/desvio > medida.ajus [1] -1.515533285 1.970461267 -0.918345074 -0.876460256 -2.105027691 -0.138034844 -0.187117338 -0.543848613 [9] 0.051023271 0.660359440 -0.714246487 -1.302701405 -0.716352088 0.186023854 -0.001365031 -0.629591181 [17] 1.527788175 0.560601567 0.750370783 -0.161986421 0.427362130 2.460722243 1.326112048 -0.172422612 [25] -0.996439389 -0.132435779 -0.054593635 0.321885401 0.338618360 0.585172589

Mas para facilitar vamos olhar denovo no gráfico, vejam como ele não muda nada em relação ao anterior, a única coisa que o centro agora é zero, ao invés de 15 como era antes. Vamos ainda ajustar uma distribuição normal com média 0 e desvio 1, que notem que fica parecido com as barras, e notem como as “perninhas” delas englobam todas as barrinhas. Ou seja antes de ficar bem retinho os “pezinhos” da distribuição, em baixo dela esta todas as barrinhas.

02

Outra coisa é que essa linha é idealmente o mesmo que a gente aumentar o número de barras, a gente poderia quebrar as barras ao invés de 1 em 1 unidade, de 0.5 em 0.5 e dobrar o número de barras, ou triplicar, se aumentarmos isso infinitamente teremos mais ou menos essa linha.

Mas como queremos olhar a distribuição t, vamos primeiro olhar a distribuição normal.

03

E agora vamos fazer o mesmo gráfico com a distribuição t com 5 graus de liberdade.

04

Hummm parece que são a mesma coisa, pelo menos para mim parece.
Mas para comparar melhor vamos olhar as 2 juntas.

05

Da para ver que a distribuição t nesse caso (com 5 graus de liberdade) tem as perninhas um pouco mais abertas que a distribuição normal, alias se olhar é isso que os graus de liberdade fazem…

06

Olha so, quanto mais a gente aumenta os graus de liberdade, mais parecida ela fica com a distribuição normal, mais ela fecha as perninhas. Agora o que acontece é que essa abertura de perninha que eu estou citando é exatamente o intervalo de confiança, então o que acontece com poucos graus de liberdade? A gente tem um intervalo de confiança mais largo comparando com um com menos, so olhar o gráfico. Isso quer dizer que com poucos graus de liberdade eu até aceito valores mais extremos, mas com mais graus de liberdade eu não aceitaria.

Vamos tentar fazer um exemplo aqui.

Vamos gerar 10 valores apenas, com a média 15 e o desvio padrão de 3, mas agora vamos adicionar um outlier, vamos adicionar um valor de 50. Esse valor faz com que a média suba um monte a mais que 15, tente calcular.

> x<-c(rnorm(10,15,3),50) > x [1] 17.02194 21.21611 13.37691 11.78852 13.88263 13.54458 15.82435 13.56146 17.39432 11.98665 50.00000

07

Agora vamos então tentar ajustar uma distribuição normal nesses dados, da para fazer isso com a função fitdistr, que você fala a distribuição e ele ajusta para a gente.

> fitdistr(x,"normal") mean sd 18.145225 10.413227 ( 3.139706) ( 2.220107)

Veja como a média fica longe do valor “real” que é 15, e o desvio então fica 10, sendo que era para ser algo em torno de 3. Mas olhar uma coisa só não da parâmetro para comparação, então vamos ajustar a mesma coisa com a distribuição t. Nesse caso a gente constrói nossa própria função t, que alias é um dos exemplos dessa função se você digitar ?fitdistr. Então não fui eu que tive a sacada de escrever isso, só estou usando :), mas vamos la.

> mydt <- function(x, m, s, df) dt((x-m)/s, df)/s > fitdistr(x, mydt, list(m = 0, s = 1), df = length(x)-1, lower = c(-Inf, 0)) m s 15.737086 5.881490 ( 1.903449) ( 2.061368)

E pimba, ambos média que fica quase 15 e o desvio que fica 5.88, quase 6 são valores muito melhor para estimar essa distribuição da qual os dados partiram que o que conseguimos com a distribuição normal. Podemos colocar essas duas distribuições no histograma e ver melhor que com números o que esta acontecendo.

08

Bem as perninhas abertas da distribuição t encaixaram muito melhor o outlier, o valor de 50 que comparando com a distribuição normal, ela teve que mover o centro e abrir um monte as perninhas para caber esse valor, enquanto com a distribuição t, ficou bem mais parecido com o que deveria ser, o centro bem perto do valor de 15 e um desvio não tão gigante. Por isso a preocupação com graus de liberdade corretamente descrito, número de amostras e pseudo-replicação são interessantes, aqui é só um pequeno exemplo como um ajuste muito melhor pode ser alcançado quando esse tipo de preocupação é levado em conta. E a distribuição t, além de ser o modelo nulo do famoso teste t é muito mais interessante do que esperamos.

#Pacote do fitdistr e seed para gerar valores iguais
library(MASS)
set.seed(12)
 
#dados exemplo
medida<-rnorm(30,15,3)
medida
 
hist(medida,xlim=c(5,25),main="Histograma da Medida")
 
#normalizando os dados
media<-mean(medida)
desvio<-sd(medida)
medida.ajus<-(medida-media)/desvio
medida.ajus
 
#histograma com distribuição normal
hist(medida.ajus,xlim=c(-5,5),prob=T,main="Histograma da Medida")
curve(dnorm(x,mean=0,sd=1),-4,4,lwd=2,col="red",add=T)
 
#distribuição normal
curve(dnorm(x,mean=0,sd=1),-5,5,lwd=2,main="Distribuição Normal",
frame.plot=F,col="red")
 
#distribuição t
curve(dt(x,df=5),-5,5,lwd=2,main="Distribuição t com 5 graus de liberdade",
frame.plot=F,lty=2,col="blue")
 
#Distribuição normal e t juntas
curve(dnorm(x,mean=0,sd=1),-5,5,lwd=2,main="Distribuição Normal",
frame.plot=F,col="red")
curve(dt(x,df=1),-5,5,lwd=2,main="Distribuição Normal",
lty=2,col="blue",add=T)
 
#distribuição normal e t com varios df
curve(dnorm(x,mean=0,sd=1),-5,5,lwd=2,main="Distribuição Normal",
frame.plot=F,col="red")
curve(dt(x,df=1),-5,5,lwd=1,lty=2,col="blue",add=T)
curve(dt(x,df=3),-5,5,lwd=1,lty=2,col="green",add=T)
curve(dt(x,df=10),-5,5,lwd=1,lty=2,col="black",add=T)
curve(dt(x,df=30),-5,5,lwd=1,lty=2,col="yellow",add=T)
legend("topright",col=c("blue","green","black","yellow"),lty=2,lwd=2,
legend=c(1,3,10,30),title="Graus de Liberdade")
 
#segundo exemplo com um outlier
x<-c(rnorm(10,15,3),50)
x
 
x.s<-scale(x)
 
hist(x.s,breaks=seq(-5,5,0.1),prob=F)
 
#ajustando distribuição normal
fitdistr(x,"normal")
 
#ajustando distribuição t
mydt <- function(x, m, s, df) dt((x-m)/s, df)/s
fitdistr(x, mydt, list(m = 0, s = 1), df = length(x)-1, lower = c(-Inf, 0))
 
#histograma comparando os ajustes
hist(x,breaks=seq(0,50,0.5),prob=T,ylim=c(0,0.4))
curve(dnorm((x-18.145225)/10.413227),0,50,lwd=2,col="red",add=T)
curve(dt((x-15.737086)/5.881490,df=10),0,50,lwd=2,lty=2,col="blue",add=T)
legend("topright",col=c("red","blue"),lty=c(1,2),lwd=2,
legend=c("normal","t com df 10"),title="Distribuição")