CHP7266 - Estatástica e Experimentação (em R)

---

---

Aula 9

Prof. Lorenzo Zanette - l.zanette@ufc.br

---

---

Ainda nosso exemplo: herbívoros x predadores

# dano causado po herbivoros
herbivoros

##  [1] 13.34 25.22 31.68  4.23 26.43 27.05 18.40 18.63 18.48 15.88

# densidade predadores po planta
predadores

##  [1]  59.26  89.67 111.26  37.35 112.89 104.04  74.10  70.77  71.08  95.80

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

No “olhometro”

(max(predadores)-min(predadores))/(max(herbivoros)-min(herbivoros))

## [1] 2.751913

Como minimizar a soma de desvios

Encontrar valor mínimo na curva de “beta”

Resultado:

lm(predadores~herbivoros)->regre1
coef(regre1)

## (Intercept)  herbivoros 
##   27.521554    2.764144

Comparemos os resíduos…

summary(regre1)

## 
## 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.5216     9.9405   2.769  0.02435 *  
## 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.0003504

Mas, lembrando…

y = a + bx + erro

Falta estimar o erro. Assumimos que ele é “normal” com média = 0 e var = ?

Podemos “desconstruir” a variância em seus componentes:

• SSY -> soma dos desvios da variável resposta (ex.predadores) = variação total que queremos explicar

• RSS -> soma dos resíduos do modelo linear = variação que o modelo não explica

• SS_mod -> soma dos desvios que atribuída a relação y = a + bx

# Como obter a soma dos desvios do modelo 

sum((predict(regre1)-predadores)^2)->RSS
RSS

## [1] 993.881

sum((mean(predadores)-predadores)^2)->SSY
SSY

## [1] 5360.001

SS_mod<-SSY-RSS
SS_mod  

## [1] 4366.12

#ou 
#soma de desvios do resíduos
RSS<-deviance(lm(predadores~herbivoros))
RSS

## [1] 993.881

#soma de desvios de predadores (total)
SSY<-deviance(lm(predadores~1))
SSY

## [1] 5360.001

#soma de desvios do modelo
SS_mod<-SSY-RSS
SS_mod

## [1] 4366.12

#ou
D_herb<-(herbivoros-mean(herbivoros))
D_pred<-(predadores-mean(predadores))
SSXY<-sum(D_herb*D_pred)
SSXY

## [1] 1579.556

coef(regre1)

## (Intercept)  herbivoros 
##   27.521554    2.764144

SS_mod<-SSXY*unname(coef(regre1)[2])
SS_mod

## [1] 4366.12

Lembrando que com essas medidas podemos estimar “quão bom” é nosso modelo linear

Coeficiente de determinação:

SSY

## [1] 5360.001

SS_mod

## [1] 4366.12

RSS

## [1] 993.881

# Quanto a var atribuida ao modelo representa da var total?
CoefD<-SS_mod/SSY
CoefD

## [1] 0.8145745

#ou
summary(regre1)

## 
## 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.5216     9.9405   2.769  0.02435 *  
## 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.0003504

Mas, queremos uma estimativa da variação do y=a+bx+erro

summary.aov(regre1)

##             Df Sum Sq Mean Sq F value  Pr(>F)    
## herbivoros   1   4366    4366   35.14 0.00035 ***
## Residuals    8    994     124                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# Médias de vars. = variâncias...
summary.aov(regre1)[[1]][3]

##             Mean Sq
## herbivoros   4366.1
## Residuals     124.2

sqrt(summary.aov(regre1)[[1]][3])

##             Mean Sq
## herbivoros   66.077
## Residuals    11.146

Aula 10

Analisando variâncias

Variável(eis) categórica como explicativa = fator , com 2 ou mais níveis

Ex. (Crawley, 2013)

One-way ANOVA

Sucesso na germinação de sementes em diferentes tipos de solo

Variável resposta ~ fator (com diferentes níveis)

Germinação de sementes ~ solo

areia<-c(6,10,8,6,14,17,9,11,7,11)
terra<-c(13,16,9,12,15,16,17,13,18,14)
argila<-c(17,15,3,11,14,12,12,8,10,13)

data.frame(areia,argila,terra)->solos

Olhemos para as médias

sapply(list(areia,argila,terra),mean)

## [1]  9.9 11.5 14.3

(tab_solo<-stack(solos))

##    values    ind
## 1       6  areia
## 2      10  areia
## 3       8  areia
## 4       6  areia
## 5      14  areia
## 6      17  areia
## 7       9  areia
## 8      11  areia
## 9       7  areia
## 10     11  areia
## 11     17 argila
## 12     15 argila
## 13      3 argila
## 14     11 argila
## 15     14 argila
## 16     12 argila
## 17     12 argila
## 18      8 argila
## 19     10 argila
## 20     13 argila
## 21     13  terra
## 22     16  terra
## 23      9  terra
## 24     12  terra
## 25     15  terra
## 26     16  terra
## 27     17  terra
## 28     13  terra
## 29     18  terra
## 30     14  terra

names(tab_solo)<-c("germina","solo")

attach(tab_solo)

Como comparar as médias comparando a variação ???

Vejamos como seriam médias iguais para um fator de dois níveis:

Vejamos como seriam níveis com médias diferentes, mas com a mesma variabilidade

Considerando nosso exemplo:

summary(solos)

##      areia           argila          terra     
##  Min.   : 6.00   Min.   : 3.00   Min.   : 9.0  
##  1st Qu.: 7.25   1st Qu.:10.25   1st Qu.:13.0  
##  Median : 9.50   Median :12.00   Median :14.5  
##  Mean   : 9.90   Mean   :11.50   Mean   :14.3  
##  3rd Qu.:11.00   3rd Qu.:13.75   3rd Qu.:16.0  
##  Max.   :17.00   Max.   :17.00   Max.   :18.0

Quantificando a variação com as somas dos quadrados…

SST<-sum((germina-mean(germina))^2)
SSE<-sum((areia-mean(areia))^2)+sum((argila-mean(argila))^2)+sum((terra-mean(terra))^2)
SST

## [1] 414.7

SSE

## [1] 315.5

SSA<-SST-SSE
SSA

## [1] 99.2

Para comparar as vars no nosso experimento com 3 tipos de solo, precisamos que a amplitude de variação em cada nível do fator não sejam muito diferentes (i.e. significativamente dif.)

tapply(germina,solo,var)

##     areia    argila     terra 
## 12.544444 15.388889  7.122222

fligner.test(germina~solo)

## 
##  Fligner-Killeen test of homogeneity of variances
## 
## data:  germina by solo
## Fligner-Killeen:med chi-squared = 0.36507, df = 2, p-value = 0.8332

Finalmente, rodando uma anova one-way no R

lm(germina~solo)->mod_aov

summary.aov(mod_aov)

##             Df Sum Sq Mean Sq F value Pr(>F)  
## solo         2   99.2   49.60   4.245  0.025 *
## Residuals   27  315.5   11.69                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

SSA

## [1] 99.2

SSE

## [1] 315.5

SST

## [1] 414.7

summary(mod_aov)

## 
## Call:
## lm(formula = germina ~ solo)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
##   -8.5   -1.8    0.3    1.7    7.1 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    9.900      1.081   9.158 9.04e-10 ***
## soloargila     1.600      1.529   1.047  0.30456    
## soloterra      4.400      1.529   2.878  0.00773 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.418 on 27 degrees of freedom
## Multiple R-squared:  0.2392, Adjusted R-squared:  0.1829 
## F-statistic: 4.245 on 2 and 27 DF,  p-value: 0.02495

tapply(germina,solo,mean)

##  areia argila  terra 
##    9.9   11.5   14.3

plot(germina~solo,col="orange")