Inferência estatística

Inferência estatística nada mais é do que tentar fazer afirmações sobre as características de uma população, com base em informações dadas por amostras. Para isso, nós vamos estimar os parâmetros reais da população através de nossa amostra. Nesse caso o termo população é qualquer conjunto de elementos que você está observando, ou seja, estamos nos referindo a uma população de números que podem ou não estar associados a uma população biológica. Amostra é qualquer subconjunto desta população.

A altura média das árvores da floresta é um parâmetro que não sabemos e então usaremos a amostra para estimar este parâmetro. Podemos estimar a distribuição das alturas na amostra e esperamos que esta reflita a distribuição de altura de todas as árvores. Temos então como calcular os descritores da amostra e como estimativa dos descritores reais da população. Porém, afirmar algo sobre uma população baseado em uma amostra pode ser perigoso e nos fazer cometer erros.

O que nós fazemos então é ajustar um modelo teórico de distribuição de probabilidades nos nossos dados para poder inferir sobre os parâmetros da população. Por exemplo, observando o histograma da distribuição das alturas, os valores estimados da média e variância dos dados, modelamos a ‘altura das árvores’ como uma variável aleatória que tem uma distribuição de probabilidades normal com média (μ) e desvio padrão (σ).

Baixe o arquivo de dados arvores.csv

# dados da altura das árvores coletadas pelo pesquisador
arv <-read.table("arvores.csv", sep=";")

hist(arv$alt)

mean(arv$alt)
var(arv$alt)
sd(arv$alt)

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). A palavra aleatória indica que em geral só conhecemos aquele valor depois do experimento ser realizado.

Exemplos de variáveis aleatórias:

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.

Variável aleatória discreta: número ou a quantidade observada na unidade experimental ou tentativa.
- Representada por números inteiros (0, 1, 2, 3, 4…);
- Não pode conter números negativos;
- Número finito de possibilidades;
- Podemos achar a probabilidade de cada evento.

Variável aleatória contínua: usualmente medidas contínuas como peso, altura, distância, pH, biomassa, etc.
- Representada por números não inteiros (1.3, - 1.54, - 1.7);
- Pode conter números negativos;
- Número infinito de possibilidades;
- Probabilidade de cada evento é zero.

Estas variáveis aleatórias terão sempre uma função de probabilidade associando cada possível valor da variável (X) à sua probabilidade de ocorrência P(X). Quando conhecemos todos os valores de uma variável aleatória, juntamente com suas respectivas probabilidades, temos uma distribuição de probabilidades. As distribuições de probabilidades discretas são conhecidas como função massa de probabilidade, enquanto que distribuições de probabilidade contínuas são conhecidas como função de densidade de probabilidade. A diferença está no fato de que nas distribuições discretas temos a probabilidade para cada valor de X, enquanto que nas distribuições contínuas temos a probabilidade para um intervalo.

Hoje em dia você pode usar diversas distribuições de probabilidades que melhor se ajustam aos seus dados (e até mesmo fazer testes que não levam em conta estas distribuições teóricas, mas não veremos isso aqui). Falaremos brevemente, sem mencionar os aspectos matemáticos de algumas distribuições mais usadas por biólogos. Para saber mais sobre estas distribuições recomendo ler a apostila “Estatística aplicada à ecologia usando o R”, de onde vieram algumas informações deste roteiro.

Bernouli

Quando jogamos uma moeda, os resultados possíveis são cara ou coroa. Se a moeda não é viciada, a probabilidade de dar é cara a mesma probabilidade de dar coroa, ou seja 50%. Neste caso estamos interessados na ocorrência de sucesso ou fracasso. Assim, podemos definir uma variável aleatória X que assume apenas dois valores: 1 se ocorrer sucesso, 0 se ocorrer fracasso. Seu único parâmetro será p que indica a probabilidade de sucesso.

Seu valor esperado na distribuição será o próprio p

E[X] = p

e a sua variância

var[X] = p(1-p)

Usamos a notação X ~ Bernouli(p) para indicar esta variável.

plot(dbinom(0:1, size=1, prob=0.3), 
    xlim=c(0,3), ylim=c(0,1), 
    xlab="resultados previstos",
    ylab="Probabilidade", 
    main="Distribuição Bernoulli p=0.3",
    type="h", xaxt="n")
axis(1, at=c(1,2), labels = c(0,1))

Binomial

Agora imagine que repetimos um ensaio de bernoulli n vezes. Esta amostra será constituída de uma sequência de sucessos e fracassos. A distribuição binomial é a distribuição de probabilidade discreta do número de sucessos em uma sequência de n tentativas tal que:

1. as tentativas são independentes;

2. cada tentativa resulta apenas em duas possibilidades, sucesso ou fracasso; e

3. a probabilidade de cada tentativa, p, permanece constante.

Se X ~ Bin(n, p), isto é, X é uma variável aleatória distribuída binomialmente, então o valor esperado de X é:

E[X] = np

e a variância é

var[X] = np(1-p)

plot(y=dbinom(0:6, size=6, prob=0.3), 
    x=c(0:6),
    xlab="Número de girinos predados",
    ylab="Probabilidade", 
    main="Distribuição binomial:  n=6, p=0.3",
    type="h")

Poisson

Na teoria da probabilidade e na estatística, a distribuição de Poisson é uma distribuição de probabilidade discreta. Expressa a probabilidade de uma série de eventos ocorrem em um período fixo de tempo, área, volume, quadrante, etc. Esta distribuição segue as mesmas premissas da distribuição binomial:

1. as tentativas são independentes;

2. a variável aleatória é o número de eventos em cada amostra; e

3. a probabilidade é constante em cada intervalo.

Se X ~ Pois(λ), isto é, X é uma variável aleatória com distribuição Poisson, então o valor esperado de X é

E[X] = (λ)

e a variância é

var[X] = (λ)

#probabilidade de que 5 borboletas visitem a flor:
dpois(5,10)

# histograma da distribuição
plot(dpois(seq(1,20, by =1), lambda = 10), type ="h",
     xlab = "Número de visitas", 
     ylab = "Probabilidade", 
     main = "Função massa de probabilidade")

Normal

A distribuição normal é uma das mais importantes distribuições com probabilidades contínuas. Esta distribuição é descrita pelos parâmetros média (μ) e desvio padrão (σ), ou seja, conhecendo-se estes parâmetros consegue-se determinar qualquer probabilidade em uma distribuição Normal.

A importância da distribuição normal como um modelo de fenômenos quantitativos é devido em parte ao Teorema Central do limite . O teorema afirma que “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.

Variáveis aleatórias com distribuição aproximadamente normal apresentam as seguintes propriedades:

– Metade (50%) está acima (e abaixo) da média

– Aproximadamente 68% está dentro de 1 desvio padrão da média

– Aproximadamente 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

# calculando a probabilidade de ter um peixe maior ou igual a 20.15
prob <- pnorm(q = 20.15, mean = 17.1, sd = 1.21, lower.tail = F)
prob

#histograma
x<-seq(13, 22, 0.1)
y=dnorm(x, mean = 17.1, sd=1.21)
plot(x, y, type="l", lwd=2, col="red", 
     ylab = "Probabilidade",
     xlab= "Comprimento dos peixes",
     main ="Função densidade de probabilidade")
abline(v = 17.1, lty=2)
abline(v = 20.15, lty=3)

Graus de liberdade

Antes de partirmos para o teste t em si, precisamos entender o conceito de graus de liberdade. Para encontrar a probabilidade associda com qualquer valor em partircular de uma estatística (no nosso exemplo, a estatística t), precisamos conhecer quantas observações independentes foram usadas para calculá-las. Usualmente se diz que os graus de iberade dos testes estatísticos são o número de observações independentes menos o número de parâmetros estimados. Em nosso exemplo, tínhamos 100 árvores amostradas e estimamos a média e o desvio padrão para construir a distribuição de resultados esperados quando a hipótese nula era verdadeira. Portanto, o grau de liberdade para a estatística t em nosso exemplo é 100-2 = 98.

Teste t

Não entraremos em detalhes de como se calcula a estatística t e como é a sua distribuição de probabilidades (recomendo ver na wikipedia). O importante a saber aqui é que tendo as duas amostras modeladas como uma distribuição normal, calculamos a estatística t e com um certo valor de graus de liberdade, nós podemos recorrer às tabelas da distribuição t para ver qual a probabilidade de se obter aquele valor que obtivemos dada que a hipótese nula é verdadeira. Se esta probabilidade for muito pequena, a chance de estarmos caindo no erro do tipo 1 é pequena.

No R, a função usada para fazer um teste t é t.test.

t.test(arv$alt~arv$sp, var.equal=T)
#o argumento var.equa=T diz que as variâncias das duas espécies são as mesmas
# podemos fazer o teste com a premissa de que as variâncias não sejam iguais,
# usando var.equal=F, isso vai "consumir" graus de liberdade.

Vamos interpretar este resultado: temos o valor de t calculado e os graus de liberade (já calculado pela gente antes). Com estes 2 valores, o que o R faz é calcular a probabilidade de que a hipótese nula seja verdadeira. Esse é o valor-p, que neste exemplo deu um valor beem pequeno. Vemos então que podemos rejeitar H0, com uma probabilidade muito baixa de estarmos incorrendo no erro de tipo 1, e desta forma acreditamos que as espécies de árvores tenham distribuições de alturas diferentes (médias diferentes).

O teste t possui três variantes básicas: a mais usada é normalmente chamada de teste t para amostras independentes, e permite a comparação de duas médias (a que acabamos de fazer) (mais roteiros aqui e aqui ); o teste t para uma amostra compara uma média com um valor fixo (roteiro aqui ); e o teste t pareado (ou teste t para amostras dependentes) compara médias de unidades amostrais dependentes aos pares (roteiro aqui ). Todos estas variações do teste t podem ser executadas no R pela função t.test, e o uso do teste t pareado e do teste t para uma amostra pode ser ativado com o uso dos argumentos paired e mu, respectivamente. Para saber mais como realizar estes testes no R, veja o capítulo 3 da apostila Introdução ao uso do software R para as Ciências Biológicas.

Recomendo também este roteiro unicaulda de teste t, que explica também a diferença entre um teste unicaudal e bicaudal. No teste unicaudal sua hipótese tem uma direção definida para a diferença entre as médias, por exemplo, esperamos que a sp1 seja maior do que a espécie 2. No teste bicaudal (o que acabamos de fazer), a hipótese diz respeito apenas se há diferença entre as médias e não qual é maior.

ANOVA

Vamos então incluir no nosso exemplo inicial uma terceira espécie amostrada (e vamos diminuir o número de amostras por espécie para ficar mais fácil de fazer os cálculos passo a passo):

alt.sp3 <- c(36.90076761,33.23131727,32.82767311,24.93410577,31.87626329,
         28.76824248,26.6436144,31.41238398,27.65383929,31.70425719,
         33.60752445,28.30263354,29.32210674,23.4790054,22.20793046,
         18.28797293,23.39044736,28.75820917,29.94344703,34.47430235,
         32.09820862,26.99845592,27.64858951,31.90557306,33.63236368,
         25.35966179,35.00885172,32.89125598,18.98802229,30.8516906,
         22.04191398,35.37232889,35.67150408,34.30257037,24.43038475,
         20.6274663,31.15717916,23.27023231,36.20568545,34.33068448,
         28.40769054,31.00578874,29.50989158,31.69915991,37.25975797,
         30.1279997,38.623371,35.36472829,29.88304259,31.00145491)

sp3 <- data.frame(alt=alt.sp3, sp=rep("sp3",50))
arv2 <- rbind(arv,sp3)

amostras <- seq(1,150,by=5)

arv3 <- arv2[ amostras, ]

Neste exemplo lidaremos com um modelo simples de anova chamado one way ANOVA, porque lida apenas com uma variável independente categórica (fator).

A ANOVA testará a hipótese nula de que as médias das alturas das 3 espécies não diferem.

Faremos agora a ANOVA “na unha”, para entendermos melhor cada passo da análise e a interpretação de seus resultados para depois usarmos a função do R que faz esta análise diretamente.

ANOVA na unha

Vamos calcular os valores da tabela de ANOVA. Começando com os desvios quadráticos, ou seja, quanto os dados desviam da média (idéia parecida com a variância). O ponto importante é que essa variação é aditiva e portanto, pode se decomposta. A variação total é decomposta em variação relacionada ao tratamento (entre grupos), no nosso caso às espécies, e uma variação interna dos grupos (chamada de erro). A estatística F é a razão entre essas variaões após dividir cada uma delas pelos seus respectivos graus de liberdade. Complicou? Vamos fazer os cálculos e ver os gráficos para ver se entendemos melhor:

A tabela de ANOVA que vamos construir é essa:

Primeiro vamos mudar a forma do nosso data frame para facilitar os comandos

arv4 <- data.frame(sp1=arv3$alt[arv3$sp=="sp1"],
                   sp2=arv3$alt[arv3$sp=="sp2"],
                   sp3=arv3$alt[arv3$sp=="sp3"])
arv4

boxplot(arv4)

Vamos calcular a média geral, as médias para cada espécie, e as diferenças entre a média geral e a média para cada espécie, para chegar ao valor da soma dos quadrados total:

media.geral <- mean(arv3$alt)
media.geral

medias.sps <- apply(arv4, 2, mean)
medias.sps

dif.geral <- arv4 - media.geral
ss.especies <- dif.geral^2
ss.especies

ss.total <- sum(ss.especies)
ss.total

Para calcular os desvios quadráticos totais, nós subtraímos cada altura das árvores pela média geral e elevamos ao quadrado. Veja o gráfico:

vetor.cor <- rep(1:3, each=10) #vetor de cores

plot(x = c(1:30), y = arv3$alt, ylim=c(10,40),
     pch=(rep(c(15,16,17), each=10)),
     col=vetor.cor,
     ylab="Variável Resposta", xlab="Observações",
     main="Variação total")
    for(i in 1:30)
    {
    lines(c(i,i),c(arv3$alt[i],mean(arv3$alt)),col=vetor.cor[i])
    }
    abline(h=media.geral)

Agora vamos fazer a somatória dos desvios quadráticos dentro de cada grupo (ss.intra).

O gráfico para entender esse cálculo:

vetor.medias<-rep(medias.sps, each=10)

plot(c(1:30), arv3$alt, ylim=c(10,40),
     pch=(rep(c(15,16,17),each=10)),
     col=vetor.cor,
     main="Variação Intra Grupos",
     ylab="Variável Resposta", xlab="Observações")
    for(i in 1:30)
    {
    lines(c(i,i),c(vetor.medias[i],arv3$alt[i]),col=vetor.cor[i])
    }
    lines(c(1,10),c(medias.sps[1],medias.sps[1]),col=1)
    lines(c(11,20),c(medias.sps[2],medias.sps[2]),col=2)
    lines(c(21,30),c(medias.sps[3],medias.sps[3]),col=3)

Cálculo dos valores:

medias.sps
ss.sp1=sum((arv4$sp1-medias.sps["sp1"])^2)
ss.sp1
ss.sp2=sum((arv4$sp2-medias.sps["sp2"])^2)
ss.sp2
ss.sp3=sum((arv4$sp3-medias.sps["sp3"])^2)
ss.sp3
ss.intra=ss.sp1+ss.sp2+ss.sp3
ss.intra

A soma dos quadrados entre grupos:

plot(c(1:30), vetor.medias, ylim=c(10,40), 
  pch=(rep(c(15,16,17),each=10)),
  col=vetor.cor, 
  main="Variação Entre Grupos", 
  ylab="Variável Resposta", xlab="Observações")
 for(i in 1:30)
    {
    lines(c(i,i),c(vetor.medias[i],mean(vetor.medias)),col=vetor.cor[i])
    }
    abline(h=media.geral)
    points(c(1:30),arv3$alt, ylim=c(10,50), 
           pch=(rep(c(0,1,2),each=10)), col=vetor.cor, cex=0.5)

medias.sps
media.geral
ss.entre=10*sum((medias.sps-media.geral)^2)
ss.entre

#conferindo os cálculos
ss.intra+ss.entre
ss.total

Cálculo do F

# Desvios médios
ms.entre=ss.entre/2
ms.intra=ss.intra/27
ms.entre
ms.intra

#F - razão das variâncias
F.sps=ms.entre/ms.intra
F.sps

Distribuição F

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

curve(expr=df(x, 2,27),
      main="Distribuição F de Fisher (df=2,27)", 
      xlab="Valor F",
      ylab="Densidade Probabilística (df)",
      xlim=c(2,12))
abline(v=F.sps, col="red")
abline(h=0, lty=2)

xf=seq(F.sps, 12, 0.01)
ydf=df(xf, 2, 27)
polygon(c(F.sps,xf),c(0,ydf),col="red")

text(x= 7,y=0.08,paste("pf(x) =",
     round(pf(F.sps,2,27,lower.tail=F),4)), 
     cex=0.8, col="red")

Cálculo do P

p.sps=pf(F.sps, 2, 27, lower.tail=FALSE)
p.sps

Para mais informações sobre o F e as comparações com o teste t, veja esse roteiro.

A tabela final

Colocando os dados calculados em nossa tabela de anova

anova.sps <- aov(alt~sp, data=arv3)
summary(anova.sps)

TukeyHSD(anova.sps)

#revendo gráfico dos dados
stripchart( arv3$alt~arv3$sp,
           vertical = TRUE, pch = 16, 
           col = c("black", "red", "green"))

#boxplot
boxplot(arv4)

lm.anova.sp <- lm(alt~sp, data=arv3)
summary.aov(lm.anova.sp) #mesmo que o summary do aov

# mostra o resultado de forma um pouco diferente, consegue entender?
summary(lm.anova.sp) 

Regressão linear simples

Vamos pensar agora se nossa variável independente (explicativa) fosse uma variável contínua, por exemplo a quantidade de nutrientes no solo. Nós criamos uma variável de quantidade de nutrientes (completamente fictícia como todos os dados), e queremos agora saber se a quantidade de nutrientes no solo afeta a altura das árvores. Nossa hipótese é de que quanto mais nutrientes, mais alta a árvore pode ficar. Nesse caso, chamamos esse modelo linear classicamente como regressão.

Inserindo os dados de quantidade de nutrientes:

# criando a variável nutrientes
nutri <- c(1.09, 1.73, 1.49, 2.84, 3.08, 8.24, 8.67, 6.10, 8.37, 8.03,
           10.62, 14.26, 13.77, 14.81, 14.40, 17.27, 15.46, 18.66, 16.03,
           17.37, 26.80, 27.1, 29.09, 20.46, 24.51, 28.56, 24.87, 27.44,
           26.02, 22.94)

# ordenando arv3 por altura (para criar uma relação com nutrientes)
arv3 <-arv3[order(arv3$alt), ]

#colocando no data.frame de dados
arv3$nutri <- nutri

Agora vamos observar a relação entre a altura da ávore e a quantidade de nutrientes no solo:

plot(arv3$alt ~ nutri)

Vamos então modelar a relação e ver os resultados:

alt.nutr <- lm(alt ~ nutri, data = arv3)
summary(alt.nutr)

Os resultados de uma regressão linear que rejeita a hipótese nula são comumente apresentados com uma representação gráfica que mostra os dados (em um gráfico de dispersão) e a equação (na forma de uma reta). Para adicionar a reta da equação ao gráfico criado pelo comando plot, basta, com o gráfico aberto, executar o comando abline, fazendo referência ao objeto que contém o resultado da regressão:

Vamos agora desenhar esta reta da regressão para entendermos melhor:

plot(arv3$alt ~ nutri)
abline(alt.nutr, col = "red") # reta da regressão

#reta da hipótese nula de ausência de efeito dos nutrientes
abline(h = mean(arv3$alt), col = "blue", lty = 2)  

Os resultados do summary são parecidos com os de uma ANOVA, mas conta com alguns detalhes a mais. No exemplo acima, a linha do (Intercept) refere-se ao intercepto da equação linear (o a, na equação y = a + bx), e verifica a hipótese nula de que ele seria igual a zero; de uma forma geral, nós raramente nos preocupamos com a esta parte do resultado, pois o intercepto normalmente não faz parte da nossa pergunta. Já a linha debaixo, com o nome da variável explicativa (nutri), é sempre do nosso interesse: ela verifica a hipótese nula de que o parâmetro b da equação (inclinação da reta) seja igual a zero; ou seja, ao rejeitarmos a hipótese nula nós concluímos que de fato há uma relação estatisticamente significativa entre as duas variáveis, e que o efeito desta relação é medido por este parâmetro. Por fim, a interpretação da correlação, que está representada no parâmetro Adjusted R-squared também é importante, e mede o quanto da variação na variável resposta foi de fato explicada pela variável explicativa.

Confuso? Então vejamos os resultados do exemplo acima: o teste rejeita a hipótese nula de que a altura das árvores não varia em função da quantidade de nutrientes (o valor de p foi menor do que 0.05). O resultado da regressão permite concluirmos que, em média, cada incremento de uma unidade de nutriente aumenta a altura em 0,58 (este é o valor de b da equação). Por fim, o valor de R² ajustado nos diz que 76,3% da altura da árvore é explicada pela quantidade de nutrientes do local.

ANCOVA

Com os dados das árvores nós testamos independentemente duas hipóteses: se a altura média das árvores era diferente para cada espécie e se a altura das árvores dependia da quantidade de nutrientes no solo. Mas porque fazer estas análises separadamente? Ambas as variáveis independentes podem influenciar na altura das árvores e, para isso, precisamos modelá-las em conjunto.

O modelo linear proposto abaixo é muitas vezes chamado de ANCOVA, ou seja, Análise de Covariância, que tem duas variáveis explicativas, uma categórica (no nosso caso a espécie de planta) e uma contínua (os nutrientes).

alt.sp.nutr <- lm(alt ~ nutri + sp, data = arv3)
summary(alt.sp.nutr)

O modelo proposta acima é o que chamamos de modelo aditivo, ou seja, as médias das alturas das árvores podem diferir, mas o efeito do nutriente é o mesmo para todas as espécies.

Vejamos o gráfico de resultados do modelo:

# tirando os valores de intercepto do modelo para cada especie
coef(alt.sp.nutr)
int.sp1 <- coef(alt.sp.nutr)[1]
int.sp2 <- coef(alt.sp.nutr)[1]+coef(alt.sp.nutr)[3]
int.sp3 <- coef(alt.sp.nutr)[1]+coef(alt.sp.nutr)[4]

plot(alt~nutri, data=arv3)

# reta de regressão para espécie 1
abline(a=int.sp1, b=coef(alt.sp.nutr)[2], col="black")
# reta de regressão para sp 2
abline(a=int.sp2, b=coef(alt.sp.nutr)[2], col="red")
# reta de regressão para ps 3
abline(a=int.sp3, b=coef(alt.sp.nutr)[2], col="green")

Mas e se nós acreditamos que as espécies respondem de maneira diferente à quantidade de nutrientes no solo? Ou seja, uma espécie consegue crescer muito mais com a mesma quantidade de nutrientes do que as outras.

Para testar esse segundo modelo, inserimos então a interação entre a variável qualitativa espécie e a quantitativa nutrientes. Essa interação nada mais é do que um novo fator resultado da multiplicação entre as variáveis. Esse novo modelo é chamado de modelo com interação.

Existem duas notações de fórmula diferentes para este modelo mas que significam a mesma coisa:

alt.sp.nutr2 <- lm(alt ~ nutri + sp + nutri:sp, data=arv3)
summary(alt.sp.nutr2)

alt.sp.nutr3 <- lm(alt ~ nutri*sp, data=arv3)
summary(alt.sp.nutr3)

Vejamos o gráfico de resultado

# interceptos
int.sp1 <- coef(alt.sp.nutr2)[1]
int.sp2 <- coef(alt.sp.nutr2)[1] + coef(alt.sp.nutr2)[3]
int.sp3 <- coef(alt.sp.nutr2)[1] + coef(alt.sp.nutr2)[4]

#coeficientes
coef.sp1 <- coef(alt.sp.nutr2)[2]
coef.sp2 <- coef(alt.sp.nutr2)[2] + coef(alt.sp.nutr2)[5]
coef.sp3 <- coef(alt.sp.nutr2)[2] + coef(alt.sp.nutr2)[6]

plot(arv3$alt~nutri, col=rep(c("black", "red", "green"),each=10))

abline(a=int.sp1, b=coef.sp1, col="black", xlim=c(0.5))
abline(a=int.sp2, b=coef.sp2, col="red")
abline(a=int.sp3, b=coef.sp3, col="green")

Olhando a tabela de ANOVA deste modelo vemos que a interação (nutri:sp) não foi significativa, e os nutrientes e as espécies são responsáves pela variação no na altura das árvores de maneira indepentende, ou seja, aditivamente:

summary.aov(alt.sp.nutr3)

Se comparármos o modelo aditivo com o interativo através da anova, vemos que estes modelos nao são diferentes significativamente, e portanto, nós ficamos com o modelo mais simples (princípio da parcimônia):

anova(alt.sp.nutr3, alt.sp.nutr)

Portanto, o modelo mais parcimonioso para nossa ANCOVA é o modelo aditivo, em que vemos a influência dos nutrientes na altura das árvores assim como a diferença na altura sendo explicada pela espécie de planta.

GLM Poisson

Vamos testar se o número de atropelamentos de uma estrada está relacionado com a distância a um Parque Nacional. Nossa variável resposta será o número de atropelamentos de um local e, portanto, vamos modelá-la como pertencente a uma distribuição de Poisson.

Dados

atrop <- read.table("atropela.csv", header = T, sep = ";")
atrop

Vamos ver o gráfico de dispersão dos dados

plot( atrop$TOT.N ~ atrop$D.PARK)

Para modelar os dados vamos então fazer um modelo GLM, colocando como distribuição a Poisson, que tem por padrão a função de ligação “log”

mod <- glm(TOT.N ~ D.PARK, data = atrop, 
           family = poisson)

summary(mod)

Vamos testar se este modelos é diferente do modelo nulo

anova(mod, test = "Chisq")

Podemos ver que este modelo é significativamente diferente do modelo nulo, portanto, a inclinação da reta da relação entre atropelamentos e distância ao parque é diferente de zero.

Vamos plotar gráfico de resultado. Devemos lembrar que os coeficientes que aparecem no summary do modelos estão na escala do link function (nesse caso em log), por isso, antes de plotar a reta precisamos des-logaritmizar

coefi <- coef(mod)
coefi

plot(atrop$TOT.N~atrop$D.PARK,
     xlim = c(0,25000), ylim = c(0, 110),
     xlab = "Distância ao parque", ylab = "Número de atropelamentos")

#essa função plota a curva do modelo, colocamos exp para voltar à escala original dos dados
curve(exp(coefi[1] + coefi[2]*x), 
      xlim = c(0,25000), ylim = c(0, 110),
      xlab = "", ylab = "", 
      add=T, #adicionando a curva no plot já desenhado
      col="red")

GLM Binomial

Vamos usar os dados de um trabalho que avaliou a quantidade de veados infectados com um parasita em cada local de amostragem. O pesquisador estava interessado em saber se a proporção de animais infectados era influenciada pela quantidade de um tipo de vegetação (OpenLand) nos locais de amostra.

veados <- read.table("veados_parasitas.csv", header=T, sep=";")
str(veados)
summary(veados)

Nestes dados temos o número de animais infectados, amostrados e a quantidade de habitat em cada local de amostragem.

A proporção de animais infectados é simplesmente o número de infectados do totoal de animais.

prop.veados <- veados$infectado/veados$amostrado
prop.veados

Vejamos o gráfico dessa relação:

plot(prop.veados ~ veados$openland,
     ylim=c(0,1))

De fato parece que a proporção de animais infectados diminui com o aumento da quantidade de ‘openland’ na paisagem.

Vamos fazer um modelo dos dados como pertencendo a uma distribuição binomial, pois diz respeito ao número de “sucessos” (infectados) em relação ao total de “tentativas” (animais amostrados). Usar os dados de proporção como a variável dependente estaria errado, pois em dados de proporção a distribuição não é normal e a variância não é homogênea (ela é muito maior em torno de 0.5).

Para modelar este dados a primeira coisa é criar a variáve Y como a combinação entre o número de infectados e não infectados usando a função cbind (junta colunas):

nao.infectado <- veados$amostrado- veados$infectado
y <- cbind(veados$infectado, nao.infectado)

Criando o modelo e vendo os resultados:

mod.veado <- glm(y ~  openland, data=veados, family="binomial")
summary(mod.veado)

Vamos ver se nosso modelo é significativamente diferente do modelo nulo:

anova(mod.veado, test = "Chisq")

Sim, podemos dizer que a quantidade de habitat ‘openland’ influencia na probabilidade dos animais serem infectados pelo parasita. Essa relação é negativa, ou seja, quanto maior a quantidade de habitat, menores os níveis de infecção na população de veados.

Plotando a curva do modelo

# salvando os coeficientes
coefs <- coef(mod.veado)
coefs

#lembre que estes coefs foram linearizados com a função de ligação
# para a binomial a link function é a logit
# então temos que transformar os y no inverso do logit, a função é essa:
invlogit = function(x) return(exp(x)/(1+exp(x)))

plot(prop.veados ~ veados$openland,
     ylim=c(0,1))
#usando o invlogit para desenhar a curva:
curve(invlogit(coefs[1] + coefs[2]*x), add=T, col="red")

Teste de normalidade

Um dos testes de normalidade muito utilizado é o teste de Shapiro-Wilks (shapiro.test). Neste teste a hipótese nula é de que os dados são normais, ou seja, se o resultado do teste der não significativo significa que os dados não se distanciam de uma distribuição normal.

shapiro.test(arv3$alt)

Portanto, podemos ficar um pouco mais seguros de ter modelado nossos dados de altura das árvores como tendo uma distribuição normal. Mas não nos esqueçamos de que quando fizemos o teste t, nós assumimos que a altura das árvores da sp1 e a altura da sp2 seguiam distribuição normal, então vamos fazer para cada uma separadamente:

shapiro.test(arv$alt[arv$sp == "sp1"])
shapiro.test(arv$alt[arv$sp == "sp2"])

Teste de homogeneidade de variâncias

Para saber se a variância dos dados é homogênea entre os grupos (no caso as variáveis categóricas X) podemos fazer o teste de Barltett ou Levene.

Teste de Bartlett:

bartlett.test(arv3$alt~arv3$sp)

Teste de Levene:

# carregando o pacote car que tem a função do teste
library(car)

leveneTest(arv3$alt~arv3$sp)

O resultado não significativo para ambos os testes diz que as variâncias das alturas das espécies podem ser consideradas homogêneas.

Ótimo! Vamos agora checar o ajuste de nossos modelos lineares através do diagnóstico dos resíduos.

Checagem do modelo - Resíduos

Depois de ajustar um modelo aos dados nós precisamos investigar quão bem o modelo descreve os dados. Em particular, nós deveríamos olhar se existe alguma tendência sistemática na qualidade do ajuste (goodness of fit) do modelo (lembre-se, nós ajustamos um modelo aos nossos dados e não nossos dados ao modelo!). Para isso nós trabalhamos com os resíduos do nosso modelo.

Os resíduos são os desvios que nosso dados observados tem do modelo ajustado. Lembre-se das figuras da “ANOVA na unha”, onde tinham os valores das médias para cada espécie e linhas ligando os pontos dos dados à este valor, essas linhas são os desvíos em relação à média, e portanto so resíduos do modelo.

Lembrando que nossa premissa de normalidade recai sobre os resíduos do modelo. Ou2 seja, são os resíduos que precisam ter distribuição normal. Nós checamos se a nossa variável resposta tem distribuição normal (por exemplo com o teste acima) porque se ela tem, provavelmente seu resíduos também o terão.

Podemos “pegar” os resíduos do nosso modelo usando a função residuals ou resid.

resid(anova.sps)
residuals(anova.sps) # o mesmo que a função de cima

Outra coisa importante para se ver durante a checagem do modelo é se existem dados (geralmente outliers) que estão influenciando mais o ajuste do modelo do que todos os outros. Numa regressão podemos pensar se existe alguns pontos que estão puxando a reta da regressão muito para baixo ou para cima, influenciando na estimação dos coeficientes da regressão.

Para analisar tudo isso nós usamos a função plot com o objeto do nosso modelo. Esse comando nos retorna 4 gráficos de diagnóstico. Vejamos nosso modelo da ANOVA:

# separando a área gráfica em 4 espaço para os 4 gráficos
par(mfrow=c(2,2))

plot(anova.sps)

O primeiro gráfico plota a relação entre os resíduos do modelo (resid(anova.sps)) e os valores ajustados (3 médias das alturas para as 3 espécies). O plot identifica valores que destoam dos valores médios. A linha vermelha que passa próxio ao zero dos resíduos indica se as variâncias dos resíduos estão homogêneas. Quanto mais “reta” essa linha, melhor o diagnóstico

O segundo gráfico já é conhecido (qqnorm) e mostra os ajustes dos quantis dos resíduos aos quantis de uma distribuição normal padrão (Z, média = 0, variância = 1). Quanto mais alinhados com a reta do plot, mais os dados se ajustam à normal.

O terceiro gráfico é uma variação do primeiro e compara os valores da raiz quadrada dos resíduos padronizados. Lembrando que os resíduos padronizados serão (residuo - media(residuo))/ sd(residuo).

O último gráfico nos mostra se nós temos dados “influentes”, ou seja que estão puxando as médias das espécies pra cima ou pra baixo e, portanto modificando os coeficiente dos modelos. Neste exemplo não parece haver tendência de dados que estejam fazendo isso.

Vamos observar os gráficos de diagnóstico para todos os modelos que criamos anteriormente

# regressão: alt ~ nutri
plot(alt.nutr)

# ANCOVA modelo aditivo: alt ~ nutri + sps
plot(alt.sp.nutr)

# ANCOVA modelo interativo: alt ~ nutri * sps
plot(alt.sp.nutr2)

# GLM poisson: atrop ~ dist.park
plot(mod)

# GLM binoimial: prop.infectado ~ openland
plot(mod.veado)

# para "desligar" o par inicial
par(mfrow=c(1,1))

Então, o que acharam dos ajustes dos outros modelos? Bom, ruim? Leva tempo para nos acostumarmos a olhar estes gráficos e tirarmos conclusões de nossos ajustes. Não se desespere, procure ajuda vendo exemplos em livros e roteiros das análises no R (“google it!”).

Parabéns por chegar até aqui!!!

Para entender melhor os temas abordados nesse roteiro, além dos links já passados, recomendo fortemente a leitura dos Capítulos 8, 9, 10 e 13 (pelo menos) do livro The R Book do Crawley (se vc ainda não baixou, aqui o link ).