CHP7266 - Estatástica e Experimentação (em R)
Aula 8
Prof. Lorenzo Zanette - l.zanette@ufc.br
Como analisamos relações entre duas variáveis contínuas?
A regressão para a média (de K. Pearson)
Regression towards mediocrity in hereditary stature
Problema 1 : definir relação causa e efeito
Problema 0 : Definir variáveis e relação entre variáveis
No mundo real, em geral…
Ou ainda multivariadamente pior:
Um exemplo: herbívoros x predadores
## [1] 13.34 25.22 31.68 4.23 26.43 27.05 18.40 18.63 18.48 15.88## [1] 59.36 89.77 111.36 37.45 112.99 104.14 74.20 70.87 71.18 95.90##
## Pearson's product-moment correlation
##
## data: herbivoros and predadores
## t = 5.9282, df = 8, p-value = 0.0003504
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.6321222 0.9769825
## sample estimates:
## cor
## 0.9025378No entanto, não inferimos sobre causalidade…
Mas, se construíssemos um modelo para testar um relação linear possível?
variável resposta ~ variável explicativa
predadores ~ herbívoros
“Traduzindo” se predadores são y e herbívoros são x , a relação entre eles pode ser aproximada por algo parecido com:
y = a + bx + erro
Lembram!!!
Precisamos estimar valores de b e a
plot(predadores~herbivoros,pch=16)
abline(h=mean(predadores),col="red")
abline(v=mean(herbivoros),col="green")
for(i in 1:10) lines(c(herbivoros[i],mean(herbivoros)),c(predadores[i],predadores[i]),col="green",lty=2)
points(y=mean(predadores),x=mean(herbivoros),pch=16,col="red")
for(i in 1:10) lines(c(herbivoros[i],herbivoros[i]),c(mean(predadores),predadores[i]),col="red",lty=3)
Mas e o b, o b…o grande b
#soma dos produtos (dos desvios)
D_herb<-(herbivoros-mean(herbivoros))
D_pred<-(predadores-mean(predadores))
SSXY<-sum(D_herb*D_pred)
SSXY## [1] 1579.556## [1] 571.4448## [1] 2.764144Temos b, falta a
Agora faltam os outros valores para formar a reta:
## [1] 64.49523 97.33327 115.18964 39.31388 100.67788 102.39165 78.48180
## [8] 79.11756 78.70293 71.51616plot(predadores~herbivoros,pch=1)
points(x=herbivoros,y=pred_exp,col="red",pch=16)
points(x=herbivoros,y=pred_exp,col="red",pch=16,type='l')
# Os resíduos da regressão...
for(i in 1:10) lines(c(herbivoros[i],herbivoros[i]),c(pred_exp[i],predadores[i]),col="red",lty=3)
Agora tomemos o caminho mais curto…
Um modelo linear no R:
Testamos a normalidade para ficar em uma regressão simples:
E finalmente…
lm(predadores~herbivoros)->regre1
#coeficientes estimados para regressão (**a** e **b**)
coef(regre1)## (Intercept) herbivoros
## 27.621554 2.764144
##
## Call:
## lm(formula = predadores ~ herbivoros)
##
## Residuals:
## Min 1Q Median 3Q Max
## -8.2476 -6.9260 -4.0557 0.8453 24.3838
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 27.6216 9.9405 2.779 0.02397 *
## herbivoros 2.7641 0.4663 5.928 0.00035 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 11.15 on 8 degrees of freedom
## Multiple R-squared: 0.8146, Adjusted R-squared: 0.7914
## F-statistic: 35.14 on 1 and 8 DF, p-value: 0.0003504E jamais, esqueçam de olhar para os resíduos…eles é que precisam ser normais e sem padrões
Para ver mesmo…simulemos um amostra enorme!!!
set.seed(1234)
herbi<-rnorm(300,23,8)
herbivoros<-round(herbi,2)
erro<-rnorm(300,0.1,10)
alpha<-24
beta1<-3
predadores<-alpha+(beta1*herbi)+erro
predadores<-round(predadores,2)
lm(predadores~herbivoros)->regre2
coef(regre2)## (Intercept) herbivoros
## 22.946362 3.025908
Cuidado com a joaninhas
