Departamento de Ecologia IB-USP

Um problema

Qual é a alturas média das årvores em uma floresta?

InferĂȘncia estatĂ­stica

  • Fazer afirmaçÔes sobre as caracterĂ­sticas de uma população, com base em informaçÔes dadas por amostras.

  • População Ă© qualquer conjunto de elementos que vocĂȘ estĂĄ observando

  • Amostra Ă© qualquer subconjunto desta população.

Ajustando um modelo

  • Estimar parĂąmetros da população baseados na amostra, modelo teĂłrico de distribuição de probabilidades

Modelo de distribuição para altura das årvores

  • VariĂĄvel AleatĂłria de distribuição de probabilidades Normal com mĂ©dia (ÎŒ) e desvio padrĂŁo (σ).

## [1] 22.72789
## [1] 32.3065
## [1] 5.683881

Definindo variĂĄveis aleatĂłrias

A variĂĄvel aleatĂłria X Ă© uma variĂĄvel que tem um valor Ășnico (determinado aleatoriamente) para cada resultado de um experimento (lembre-se esse termo veio da teoria de probabilidade).

Exemplos:

  1. Altura das ĂĄrvores de uma floresta (!!)

  2. NĂșmero de presas capturadas por um predador em um determinado dia;

  3. Comprimento de um peixe adulto selecionado aleatoriamente.

As variĂĄveis aleatĂłrias podem ser discretas ou contĂ­nuas.

Função de probabilidade

Função que define as probabilidades P(X) para cada valor possível da variåvel X.

Distribuição de Probabilidades.

Algumas DistribuiçÔes de Probabilidade

Bernouli: X ~ Bernouli(p)

  • Discreta, valores possĂ­veis 0 e 1. p - probabilidade de sucesso
  • Jogar moedas: cara - 0, coroa - 1.
  • Presença/ausĂȘncia de uma espĂ©cie em locais amostrados

Binomial: X ~ Bin(n, p)

  • Discreta: nĂșmero de sucessos em n tentativas . p - prob sucesso.

  • NĂșmero sementes predadas em cada experimento contendo 20 semente.

  • NĂșmero de animais infectados em relação ao nĂșmero de capturados em cada local

Exemplo Binomial

  • Predação de girinos por odonatas, p = 0.30.
    Total girinos: 6 num lago.
  • Qual probabilidades de que 0, 1, 2, 3, 5 ou todos sejam predados?

Poisson: X ~ Pois(λ)

  • Discreta: Probabilidade de uma sĂ©rie de eventos ocorrem em um perĂ­odo fixo de tempo, ĂĄrea, volume, quadrante, etc.

  • Dados de contagem

  • AbundĂąncia de espĂ©cies em um local

  • NĂșmero de indivĂ­duos capturados por tempo

Exemplo Poisson

  • NĂșmero de visitas de borboletas em uma planta em 15 min: 10
  • Probabilidade de 5 visitas dpois(5,10)

Normal: X ~ N(ÎŒ, σ)

  • SimĂ©trica, contĂ­nua

  • Metade (50%) dos dados estĂĄ acima (e abaixo) da mĂ©dia

  • 68% dados estĂĄ dentro de 1 desvio padrĂŁo da mĂ©dia

  • 95% estĂĄ dentro de 2 desvios padrĂ”es da mĂ©dia

  • Virtualmente todos os valores estĂŁo dentro de 3 desvios padrĂ”es da mĂ©dia

Exemplo Normal

  • Qual Ă© a probabilidade de que um peixe capturado aleatoriamente tenha 20,15 cm ou mais? MĂ©dia da população Ă© 17,1 cm, desvio padrĂŁo 1,21 cm?

Teorema do Limite Central

"Toda soma de variĂĄveis aleatĂłrias independentes de mĂ©dia finita e variĂąncia limitada Ă© aproximadamente Normal, desde que o nĂșmero de termos da soma seja suficientemente grande."

Independentemente do tipo de distribuição da população, na medida em que o tamanho da amostra aumenta, a distribuição das médias amostrais tende a uma distribuição Normal.

Exemplo TLC

Voltando ao nosso problema

  • Dados de duas espĂ©cies
levels(arv$sp)
## [1] "sp1" "sp2"

Serå que a altura média de uma espécie nessa floresta é diferente da altura média da outra espécie?

Ou seja, serå que as espécies possuem uma mesma distribuição de probabilidades e as diferenças encontradas são devido ao acaso ou serå que cada espécie é uma variåvel aleatória com médias diferentes?

DistribuiçÔes de frequĂȘncias das alturas para cada espĂ©cie

Diferença entre as médias das alturas

mean(arv$alt[arv$sp=="sp2"]) - mean(arv$alt[arv$sp=="sp1"])
## [1] 5.109466

Posso afirmar que as médias são diferentes?

H0 - não hå diferença entre as alturas das duas espécies.
A média da altura da espécie 1 é igual à média da altura da espécie 2 e a diferença observada se deve ao acaso.

H1 - hå diferença entre as alturas das duas espécies.
A média da altura da espécie 1 é diferente da média da altura da espécie 2, elas vem de distribuiçÔes diferentes.

Erros associados a escolher cada hipĂłtese como verdadeira

AFIRMAÇÃO: nĂŁo hĂĄ diferença entre as alturas das espĂ©cies,
VERDADE: hå diferenças.

AFIRMAÇÃO: hĂĄ diferença entre as alturas das espĂ©cies,
VERDADE: não hå diferenças.

Definição formal dos erros:

ERRO DO TIPO I (α): rejeitar a hipótese nula dada que ela é verdadeira.

ERRO DO TIPO II (B): aceitar a hipĂłtese nula dado que ela Ă© falsa.

Testando as hipĂłteses

Contrastar a probabilidade de que a diferença que vocĂȘ encontrou nas mĂ©dias das alturas das duas espĂ©cies pode se considerada devido ao acaso ou nĂŁo.

Qual é a probabilidade de que as duas amostras pertencem à mesma população de medidas?

Podemos ou não rejeitar a hipótese nula (de que a diferença é ao acaso), baseado em alguma probabilidade de estarmos cometendo algum erro (α).

Graus de liberdade

"O nĂșmero de observaçÔes independentes menos o nĂșmero de parĂąmetros estimados."

Quantas observaçÔes independentes foram usadas para calcular a média e a variùncia dos dados sob a hipótese nula?

Nosso exemplo: 100 ĂĄrvores - 2 parĂąmetros = 98 graus de liberdade

EstatĂ­stica t

É uma medida que reflete a diferença esperada entre a a situação quando H0 Ă© correta e quando a H1 Ă© a correta.

Tendo o t e os graus de liberdade recorremos à distribuição t para saber a probabilidade de obter o valor observado dado que a H0 é verdade.

Se esta probabilidade for muito pequena, a chance de estarmos caindo no erro do tipo 1 Ă© pequena.

Teste t

t.test(arv$alt~arv$sp, var.equal=T)
## 
##  Two Sample t-test
## 
## data:  arv$alt by arv$sp
## t = -5.0125, df = 98, p-value = 2.387e-06
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -7.132311 -3.086621
## sample estimates:
## mean in group sp1 mean in group sp2 
##          20.17316          25.28262

Probabilidade na distribuição t

Modelos lineares

Modelos lineares

Grande família de modelos estatísticos que modelam uma relação linear entre a variåvel dependente Y e a(s) variåvel(is) independente(s) X.

  • teste t, ANOVA, ANCOVA, regressĂŁo, etc

  • LMs: y tem distribuição normal

  • GLMs: y tem outra distribuição (binomial, poisson) que Ă© linearizada atravĂ©s de uma funçÔa de ligação

ANOVA one-way

Comparando alturas de 3 espécies de plantas

H0: que as médias das alturas das 3 espécies não diferem.

ANOVA na unha

Soma dos quadrados total

Sutraímos cada altura das årvores pela média geral e elevamos ao quadrado.

Soma dos quadrados total

(media.geral <- mean(arv3$alt))
## [1] 23.28284
(medias.sps <- apply(arv4, 2, mean))
##      sp1      sp2      sp3 
## 20.44245 21.78334 27.62274
dif.geral <- arv4 - media.geral
(ss.especies <- dif.geral^2)
##                sp1       sp2        sp3
##  [1,]   4.00337202 17.121094 185.447876
##  [2,]  11.60801175 21.347420  30.089610
##  [3,]   0.07824755  3.323227 106.599051
##  [4,]  15.75374565 26.834280  24.948725
##  [5,]   9.15530094  3.247214  77.710675
##  [6,] 224.43179074 23.066347   4.313177
##  [7,]  15.28543337  7.866884   1.539904
##  [8,]  53.88994194 35.359990   7.051024
##  [9,]   7.21577099 15.631997  26.264064
## [10,]   1.52250306 25.688498  46.856173
(ss.total <- sum(ss.especies))
## [1] 1033.251

Soma dos quadrados dentro de cada grupo

Soma dos quadrados dentro de cada grupo

(ss.sp1=sum((arv4$sp1-medias.sps["sp1"])^2))
## [1] 262.2655
(ss.sp2=sum((arv4$sp2-medias.sps["sp2"])^2))
## [1] 157.0018
(ss.sp3=sum((arv4$sp3-medias.sps["sp3"])^2))
## [1] 322.4728
(ss.intra=ss.sp1+ss.sp2+ss.sp3)
## [1] 741.7401

Soma dos quadrados entre grupos

Soma dos quadrados entre grupos

medias.sps
##      sp1      sp2      sp3 
## 20.44245 21.78334 27.62274
media.geral
## [1] 23.28284
ss.entre=10*sum((medias.sps-media.geral)^2)
ss.entre
## [1] 291.5112
#conferindo os cĂĄlculos
ss.intra+ss.entre
## [1] 1033.251
ss.total
## [1] 1033.251

Desvios médios

Soma dos quadrados / graus de liberdade

graus de liberdade entre grupos = 3 - 1

graus de liberdade dentro dos brupos = 30 - 3

ms.entre=ss.entre/2
ms.intra=ss.intra/27
ms.entre
## [1] 145.7556
ms.intra
## [1] 27.47186

CĂĄlculo de F

F = desvio médio entre grupos / desvio médio intra grupos

#F - razĂŁo das variĂąncias
F.sps=ms.entre/ms.intra
F.sps
## [1] 5.305634

Distribuição F

Vamos ver na distribuição F qual a probabilidade de termos econtrado o valor F.sps ao acaso:

A tabela final

Colocando os dados calculados em nossa tabela de anova

ANOVA no R

A função que constrói esta tabela de ANOVA e faz o teste estatístico é a aov. Vamos comparar:

anova.sps <- aov(alt~sp, data=arv3)
summary(anova.sps)
##             Df Sum Sq Mean Sq F value Pr(>F)  
## sp           2  291.5  145.76   5.306 0.0114 *
## Residuals   27  741.7   27.47                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Agora nĂłs sabemos de ondem vem os nĂșmeros nesta tabela e como foi calculado o P, certo?

Testes à posteriori - onde estão as diferenças?

  • veja os grĂĄficos

  • faça teste Ă  posteriori - Tukey HSD

Teste Tukey

TukeyHSD(anova.sps)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = alt ~ sp, data = arv3)
## 
## $sp
##             diff        lwr       upr     p adj
## sp2-sp1 1.340893 -4.4708805  7.152667 0.8360276
## sp3-sp1 7.180300  1.3685259 12.992073 0.0132022
## sp3-sp2 5.839406  0.0276327 11.651180 0.0487468

Outra maneira de fazer a ANOVA no R

lm.anova.sp <- lm(alt~sp, data=arv3)
summary.aov(lm.anova.sp) #mesmo que o summary do aov
##             Df Sum Sq Mean Sq F value Pr(>F)  
## sp           2  291.5  145.76   5.306 0.0114 *
## Residuals   27  741.7   27.47                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Resultado ANOVA como LM

summary(lm.anova.sp) 
## 
## Call:
## lm(formula = alt ~ sp, data = arv3)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -12.1407  -3.0002  -0.0742   3.2719   9.2780 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   20.442      1.657  12.334 1.32e-12 ***
## spsp2          1.341      2.344   0.572  0.57202    
## spsp3          7.180      2.344   3.063  0.00492 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.241 on 27 degrees of freedom
## Multiple R-squared:  0.2821, Adjusted R-squared:  0.229 
## F-statistic: 5.306 on 2 and 27 DF,  p-value: 0.01139

RegressĂŁo linear simples

  • Y contĂ­nuo - X continuo

Modelo da RegressĂŁo

alt.nutr <- lm(alt ~ nutri, data = arv3);     summary(alt.nutr)
## 
## Call:
## lm(formula = alt ~ nutri, data = arv3)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -6.7176 -1.7171  0.0088  1.4124  9.2072 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 14.38711    1.05832  13.594 7.44e-14 ***
## nutri        0.58006    0.05971   9.714 1.82e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.906 on 28 degrees of freedom
## Multiple R-squared:  0.7712, Adjusted R-squared:  0.763 
## F-statistic: 94.37 on 1 and 28 DF,  p-value: 1.818e-10

Plot da RegressĂŁo

AnĂĄlise de CovariĂąncia

  • Y continuo

  • X1 contĂ­nuo - X2 categĂłrico

Modelo da ANCOVA: aditivo

alt.sp.nutr <- lm(alt ~ nutri + sp, data = arv3);     summary(alt.sp.nutr)
## 
## Call:
## lm(formula = alt ~ nutri + sp, data = arv3)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -5.6682 -0.9554 -0.1585  1.3581  7.6134 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  13.3971     1.1034  12.142 3.23e-12 ***
## nutri         0.5256     0.0561   9.369 8.09e-10 ***
## spsp2         1.6421     1.1423   1.437  0.16251    
## spsp3         3.8326     1.1965   3.203  0.00357 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.553 on 26 degrees of freedom
## Multiple R-squared:  0.836,  Adjusted R-squared:  0.817 
## F-statistic: 44.16 on 3 and 26 DF,  p-value: 2.401e-10

GrĂĄfico ANCOVA aditivo

Modelo ANCOVA com interação

alt.sp.nutr3 <- lm(alt ~ nutri*sp, data=arv3);        summary(alt.sp.nutr3)
## 
## Call:
## lm(formula = alt ~ nutri * sp, data = arv3)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -5.6724 -0.9444 -0.3313  1.0416  7.0702 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 13.40167    1.42727   9.390 1.66e-09 ***
## nutri        0.52527    0.08906   5.898 4.38e-06 ***
## spsp2        2.93717    1.97080   1.490    0.149    
## spsp3        0.43636    2.75666   0.158    0.876    
## nutri:spsp2 -0.10095    0.12423  -0.813    0.424    
## nutri:spsp3  0.17187    0.14350   1.198    0.243    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.474 on 24 degrees of freedom
## Multiple R-squared:  0.8578, Adjusted R-squared:  0.8282 
## F-statistic: 28.96 on 5 and 24 DF,  p-value: 2.003e-09

Gråfico ANCOVA com interação

Tabela de ANOVA do modelo

summary.aov(alt.sp.nutr3)
##             Df Sum Sq Mean Sq F value   Pr(>F)    
## nutri        1  796.8   796.8 130.179 3.53e-11 ***
## sp           2   66.9    33.5   5.467   0.0111 *  
## nutri:sp     2   22.6    11.3   1.846   0.1796    
## Residuals   24  146.9     6.1                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Modelos lineares generalizados

  • Outras distribuiçÔes de probabilidades alĂ©m da Normal

  • VariĂąncias nĂŁo homogĂȘneas

Mais comuns:

  • Poisson – Ășteis para dados de contagem

  • Binomial – Ășteis para dados com proporçÔes, ou presença/ausĂȘncia

Etapas do GLM

  1. Uma hipótese sobre a distribuição da variåvel resposta Yi.

  2. Especificação do modelo (variåveis X).

3.Função de ligação entre a média Yi e as variåveis X. Função que lineariza a regressão

FunçÔes de ligação

Distribuição link function
normal identity
poisson log
binomial logit
Gamma reciprocal

GLM Poisson

  • Y - NĂșmero de atropelamentos de anuros em uma estrada

  • X - distĂąncia a um Parque Nacional.

Modelo GLM Poisson

mod <- glm(TOT.N ~ D.PARK, data = atrop, family = poisson)
summary(mod)
## 
## Call:
## glm(formula = TOT.N ~ D.PARK, family = poisson, data = atrop)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -8.1100  -1.6950  -0.4708   1.4206   7.3337  
## 
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  4.316e+00  4.322e-02   99.87   <2e-16 ***
## D.PARK      -1.059e-04  4.387e-06  -24.13   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 1071.4  on 51  degrees of freedom
## Residual deviance:  390.9  on 50  degrees of freedom
## AIC: 634.29
## 
## Number of Fisher Scoring iterations: 4

Testando contra hipĂłtese nula

anova(mod, test = "Chisq")
## Analysis of Deviance Table
## 
## Model: poisson, link: log
## 
## Response: TOT.N
## 
## Terms added sequentially (first to last)
## 
## 
##        Df Deviance Resid. Df Resid. Dev  Pr(>Chi)    
## NULL                      51     1071.4              
## D.PARK  1   680.55        50      390.9 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

GrĂĄfico Modelo GLM Poisson

GLM Binomial

  • Y - Quantidade de veados infectados dentre os capturados

  • X - quantidade de vegetação aberta

Modelo GLM binomial

y <- cbind(veados$infectado, veados$amostrado - veados$infectado)
mod.veado <- glm(y ~  openland, data=veados, family="binomial")
summary(mod.veado)
## 
## Call:
## glm(formula = y ~ openland, family = "binomial", data = veados)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -3.8052  -1.1463   0.6707   1.3416   3.5483  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   2.4820     0.1736   14.29   <2e-16 ***
## openland     -6.0336     0.4432  -13.61   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 376.18  on 31  degrees of freedom
## Residual deviance: 121.30  on 30  degrees of freedom
## AIC: 204.65
## 
## Number of Fisher Scoring iterations: 5

Teste GLM binomial

anova(mod.veado, test = "Chisq")
## Analysis of Deviance Table
## 
## Model: binomial, link: logit
## 
## Response: y
## 
## Terms added sequentially (first to last)
## 
## 
##          Df Deviance Resid. Df Resid. Dev  Pr(>Chi)    
## NULL                        31     376.18              
## openland  1   254.88        30     121.30 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

GrĂĄfico GLM binomial

Checando as premissas dos modelos

Teste de normalidade

Teste de Shapiro-Wilks

  • H0: dados sĂŁo normais
shapiro.test(arv3$alt)
## 
##  Shapiro-Wilk normality test
## 
## data:  arv3$alt
## W = 0.96262, p-value = 0.3609
shapiro.test(arv$alt[arv$sp == "sp1"])
## 
##  Shapiro-Wilk normality test
## 
## data:  arv$alt[arv$sp == "sp1"]
## W = 0.9906, p-value = 0.9591
shapiro.test(arv$alt[arv$sp == "sp2"])
## 
##  Shapiro-Wilk normality test
## 
## data:  arv$alt[arv$sp == "sp2"]
## W = 0.97453, p-value = 0.3502

Teste de homogeneidade de variĂąncias

Teste de Bartlett:

bartlett.test(arv3$alt~arv3$sp)
## 
##  Bartlett test of homogeneity of variances
## 
## data:  arv3$alt by arv3$sp
## Bartlett's K-squared = 1.1109, df = 2, p-value = 0.5738

Teste homogeneidade de variĂąncias

Teste de Levene:

leveneTest(arv3$alt~arv3$sp)
## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  2   0.512  0.605
##       27

Checagem do ajuste do modelo - ResĂ­duos

RESÍDUO = Y(observado) - Y(dados ajustados)

A premissa de normalidade na verdade recai sobre os resĂ­duos do modelo.

resid(anova.sps)
##           1           6          11          16          21          26 
##   0.8395548  -0.5666556   3.1201253  -1.1287012   5.8661704 -12.1406501 
##          31          36          41          46          51          56 
##   6.7500566  -4.5005793   0.1541789   1.6065001  -2.6382600  -3.1208224 
##          61          66          71          76          81          86 
##   3.3224765  -3.6806771  -0.3024985   6.3022481  -1.3052922   7.4459310 
##          91          96         101         106         111         116 
##  -2.4542276  -3.5688778   9.2780228   1.1454976   5.9847796  -9.3347719 
##         121         126         131         136         141         146 
##   4.4754638  -2.2630830  -5.5808309  -6.9952785   0.7849457   2.5052549
#residuals(anova.sps) # o mesmo que a função de cima

ResĂ­duos ANOVA: alt ~ sp

ResĂ­duos regressĂŁo: alt ~ nutri

ResĂ­duos ANCOVA aditivo: alt ~ sp + nutri

ResĂ­duos ANCOVA interativo: alt ~ sp * nutri

ResĂ­duos GLM poisson: atrop ~ dist.park

ResĂ­duos LM com a poisson: atrop ~ dist.park

ResĂ­duos GLM binomial: prop.infec ~ openland

ResĂ­duos LM com a binomial: prop.infec ~ openland