###################################
## Q3 — ANCOVA                   ##
## Autora: Victória Souza        ##
###################################
##############################
#### Aula 1 Curso Bocaina ####
##############################

# Pacotes
library(RT4Bio)
## Carregando pacotes exigidos: MASS
## Carregando pacotes exigidos: survival
# Para essa aula vamos usar os arquivos da pasta "Aula 1 Bocaina"

# Tema da aula: Dados contínuos com distribuição Gamma ou Normal 

# Mudar diretório

#1- Criar tabela utilizando arquivo "a1.1.txt" 

#2- Faça todos os procedimentos iniciais para verificar se os dados da tabela estão ok com os dados

tab = read.table ("a1.1.txt",header=T) #header significa ler a primeira linha como titulos
tab
##    adubo   peso
## 1      1  2.822
## 2      2  2.359
## 3      3  3.091
## 4      4  2.530
## 5      5  3.475
## 6      6  3.649
## 7      7  6.700
## 8      8  8.583
## 9      9  5.321
## 10    10 10.130
## 11     1  5.236
## 12     2  5.448
## 13     3  6.123
## 14     4  4.405
## 15     5 10.823
## 16     6  9.084
## 17     7  8.100
## 18     8 11.029
## 19     9 14.502
## 20    10 12.414
attach(tab)
summary(tab)
##      adubo           peso       
##  Min.   : 1.0   Min.   : 2.359  
##  1st Qu.: 3.0   1st Qu.: 3.606  
##  Median : 5.5   Median : 5.785  
##  Mean   : 5.5   Mean   : 6.791  
##  3rd Qu.: 8.0   3rd Qu.: 9.345  
##  Max.   :10.0   Max.   :14.502
names(tab)
## [1] "adubo" "peso"
#3- Verifique graficamente se existe alguma relação entre as variaveis

plot(peso~adubo, las=1)

Interpretação: Dispersão sugere tendência linear positiva: à medida que a quantidade de adubo aumenta (variável explicativa), o peso da plântula tende a aumentar (variável resposta).

#4- Elabore uma hipótese

# Um aumento da concentração de adubo determina um aumento de peso da plântula

#5- Teste em uma função linear (y=a+bx)
# y = peso; x = adubo; "a" e o "b' parâmetros.  

m1=glm(peso~adubo) #peso em função do adubo
m1
## 
## Call:  glm(formula = peso ~ adubo)
## 
## Coefficients:
## (Intercept)        adubo  
##      1.9828       0.8743  
## 
## Degrees of Freedom: 19 Total (i.e. Null);  18 Residual
## Null Deviance:       247.3 
## Residual Deviance: 121.2     AIC: 98.78

Interpretação: Coeficiente Linear (y-intercept): a=1.98. Coeficiente Angular (inclinação da reta): b=0.87.

anova(m1,test="F")
## Analysis of Deviance Table
## 
## Model: gaussian, link: identity
## 
## Response: peso
## 
## Terms added sequentially (first to last)
## 
## 
##       Df Deviance Resid. Df Resid. Dev      F    Pr(>F)    
## NULL                     19     247.27                     
## adubo  1   126.11        18     121.15 18.737 0.0004046 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Interpretação (ANOVA): Significativo (p=0.0004046), indicando dependência do peso em relação à quantidade de adubo.

#6- construa agora a equação
summary(m1)
## 
## Call:
## glm(formula = peso ~ adubo)
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   1.9828     1.2532   1.582 0.131022    
## adubo         0.8743     0.2020   4.329 0.000405 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 6.730662)
## 
##     Null deviance: 247.27  on 19  degrees of freedom
## Residual deviance: 121.15  on 18  degrees of freedom
## AIC: 98.784
## 
## Number of Fisher Scoring iterations: 2

Interpretação (coeficientes): Reta ajustada: peso = 1,9828 + 0,8743 × adubo. Intercepto = 1,98 (peso quando adubo=0) e inclinação = 0,874 (a cada unidade de adubo, o peso aumenta em ~0,874).

# y = a +bx
# peso = 1.9828+0.8743adubo

#7- construir o grafico com a curva obtida no modelo
plot(peso~adubo,las=1,pch=16,ylim=c(0,15), xlim=c(0,12),bty="l")
curve(1.9828+0.8743*x,add=T,lty=2)
legend(0,15,legend=expression((Peso==1.9828+0.8743*Adubo)),bty="n")

Interpretação: A reta acompanha a tendência dos pontos, sugerindo a adequação do ajuste linear.

#8- Instale o pacote "RT4Bio"....

# Crítica ao modelo...

# A distribuição utilizada esta ok?!

#8.1 Dividir o valor de Residual deviance/degrees of freedom (Residual deviance) 
summary(m1)
## 
## Call:
## glm(formula = peso ~ adubo)
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   1.9828     1.2532   1.582 0.131022    
## adubo         0.8743     0.2020   4.329 0.000405 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 6.730662)
## 
##     Null deviance: 247.27  on 19  degrees of freedom
## Residual deviance: 121.15  on 18  degrees of freedom
## AIC: 98.784
## 
## Number of Fisher Scoring iterations: 2
121.15/18 
## [1] 6.730556
# 6.730556 é muito próximo à 6.730662, indica que o modelo está utilizando distribuição adequada

#8.2
hist(peso)

#8.3
library(RT4Bio)
rdiagnostic(m1)

Interpretação (diagnóstico): Não há evidências de que os dados não atendam aos parâmetros de normalidade/homocedasticidade; o modelo linear é aceitável.

# Distribuição ok!!


# Abrir o script salvo, e incluir no objeto trabalhado no modelo anterior uma nova variavel, "variedade da semente", que apresenta duas cores (vermelho e azul)
varied=rep(c("vermelha","azul"),c(10,10))
varied
##  [1] "vermelha" "vermelha" "vermelha" "vermelha" "vermelha" "vermelha"
##  [7] "vermelha" "vermelha" "vermelha" "vermelha" "azul"     "azul"    
## [13] "azul"     "azul"     "azul"     "azul"     "azul"     "azul"    
## [19] "azul"     "azul"
summary(varied)
##    Length     Class      Mode 
##        20 character character
tab2 = data.frame (tab,varied) 
summary(tab2)
##      adubo           peso           varied         
##  Min.   : 1.0   Min.   : 2.359   Length:20         
##  1st Qu.: 3.0   1st Qu.: 3.606   Class :character  
##  Median : 5.5   Median : 5.785   Mode  :character  
##  Mean   : 5.5   Mean   : 6.791                     
##  3rd Qu.: 8.0   3rd Qu.: 9.345                     
##  Max.   :10.0   Max.   :14.502
# Agora temos uma variavel contínua e outra categórica, seiria considerado uma ANCOVA, mas qui continua sendo um GLM

# incluindo outra hipótese e a interação
# adubo, varied, adubo:varied


m2=glm(peso~adubo*varied)
mn=glm(peso~1) #modelo nulo
anova(m2,mn,test="F")
## Analysis of Deviance Table
## 
## Model 1: peso ~ adubo * varied
## Model 2: peso ~ 1
##   Resid. Df Resid. Dev Df Deviance      F    Pr(>F)    
## 1        16     45.146                                 
## 2        19    247.267 -3  -202.12 23.878 3.773e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Interpretação: O modelo com adubo + cor + interação ajusta melhor que o nulo (p=3.773e-06; altamente significativo); os preditores explicam variação do peso.

# 0.000003773

# descartar o modelo nulo

m2=glm(peso~adubo*varied)
anova(m2,test="F")
## Analysis of Deviance Table
## 
## Model: gaussian, link: identity
## 
## Response: peso
## 
## Terms added sequentially (first to last)
## 
## 
##              Df Deviance Resid. Df Resid. Dev       F    Pr(>F)    
## NULL                            19    247.267                      
## adubo         1  126.115        18    121.152 44.6957 5.235e-06 ***
## varied        1   74.128        17     47.024 26.2713 0.0001017 ***
## adubo:varied  1    1.878        16     45.146  0.6656 0.4265816    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Interpretação: A interação adubo:varied não é significativa (p=0.4265816; p>0,05); ou seja, adubo tem efeito similar nas duas variedades (azul e vermelha). A interação pode ser retirada do modelo. As inclinações podem ser consideradas iguais entre as duas cores (retas paralelas).

#simplificando o modelo
m2=glm(peso~adubo*varied)
m2.1=glm(peso~adubo+varied) #sem interação
anova(m2,m2.1,test="F")
## Analysis of Deviance Table
## 
## Model 1: peso ~ adubo * varied
## Model 2: peso ~ adubo + varied
##   Resid. Df Resid. Dev Df Deviance      F Pr(>F)
## 1        16     45.146                          
## 2        17     47.024 -1   -1.878 0.6656 0.4266

Interpretação: Diferença não significativa entre os modelos m2 e m2.1: preferimos o m2.1 (mais simples), com efeitos aditivos de adubo e cor.

# o modelo mais simples é o m2.1

# o modelo minimo é o m2.1
rdiagnostic(m2.1)

Interpretação (diagnóstico): Resíduos sem padrões marcantes; pressupostos atendidos para o modelo final. Distribuição Gaussian usada adequadamente.

m2.1
## 
## Call:  glm(formula = peso ~ adubo + varied)
## 
## Coefficients:
##    (Intercept)           adubo  variedvermelha  
##         3.9080          0.8743         -3.8504  
## 
## Degrees of Freedom: 19 Total (i.e. Null);  17 Residual
## Null Deviance:       247.3 
## Residual Deviance: 47.02     AIC: 81.86
summary(m2.1)
## 
## Call:
## glm(formula = peso ~ adubo + varied)
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      3.9080     0.8853   4.414 0.000379 ***
## adubo            0.8743     0.1295   6.752 3.38e-06 ***
## variedvermelha  -3.8504     0.7438  -5.177 7.58e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 2.766119)
## 
##     Null deviance: 247.267  on 19  degrees of freedom
## Residual deviance:  47.024  on 17  degrees of freedom
## AIC: 81.856
## 
## Number of Fisher Scoring iterations: 2

Interpretação (coeficientes): Inclinação (b) = 0,874: cada unidade de adubo aumenta o peso em ~0,874 para ambas as cores/variedades (já que não há interação). O coeficiente variedvermelha é -3.8504, o que indica que a variedade vermelha tem peso médio menor. Para diagnóstico é possível fazer 47.024 (Residual deviance)/17 (degrees of freedom of residual deviance)= 2,766118 e compará-lo com o parametro de dispersão para família Gaussian -> 2.766119 (praticamente iguais), indica que o modelo está utilizando distribuição adequada.

plot(peso~adubo,pch=c(16),col=c("blue","red"),ylim=c(0,15), xlim=c(0,12),bty="l",las=1)
curve(3.9080+0.8743*x,add=T,col="blue",lty=2)
curve((3.9080-3.8504)+0.8743*x,add=T,col="red")
legend(0,15,legend=c("Azul","Vermelho"),lty=1,pch=16,col=c("blue","red"),bty="n")

Interpretação: As retas são paralelas (confirma ausência de interação): mesma inclinação para as duas cores. Ou seja, a medida que aumenta a quantidade de adubo, as plantulas vermelha e azul crescem com o mesmo efeito. Porém, a reta azul fica acima da vermelha em todo o eixo de adubo (intercepto maior): para o mesmo adubo, a cor azul resulta em maior peso médio.

# modelo nulo
m.n=glm(peso~1)
m2.1=glm(peso~adubo+varied)

anova(m.n,m2.1,test="F")
## Analysis of Deviance Table
## 
## Model 1: peso ~ 1
## Model 2: peso ~ adubo + varied
##   Resid. Df Resid. Dev Df Deviance      F    Pr(>F)    
## 1        19    247.267                                 
## 2        17     47.024  2   200.24 36.196 7.461e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Interpretação (ajuste global): Altamente significativo (p=7.461e-07). Rejeitamos a hipótese de que os preditores não explicam nada; o modelo final ajusta muito melhor que o nulo. O modelo final explica ~81% da deviance do nulo (200,24 (Deviance do m2.1)/247,27 (Deviance Residual do Modelo Nulo)). Decomposição (ANOVA): adubo ≈ 51% (126,115/247,267) e cor ≈ 30% (74,128/247,267). Ambos contribuem de forma relevante.

# Utilizando o resultado do modelo minimo adequado contra o modelo nulo, faga uma ANOVA e  divida o valor de deviance do moldeo test/Resid.Dev (Null)

200.24/247.267 
## [1] 0.8098129
# 0.8098129 = peseudo R2

# Utilizando o resultado do modelo minimo adequado do ANOVA,  dividir o valor de deviance/Resid.Dev  (Null)

anova(m2.1,test="F")
## Analysis of Deviance Table
## 
## Model: gaussian, link: identity
## 
## Response: peso
## 
## Terms added sequentially (first to last)
## 
## 
##        Df Deviance Resid. Df Resid. Dev      F    Pr(>F)    
## NULL                      19    247.267                     
## adubo   1  126.115        18    121.152 45.593 3.383e-06 ***
## varied  1   74.128        17     47.024 26.799 7.583e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# qual i a explicagco do adubo:
# adubo: 126.115/247.267

126.115/247.267
## [1] 0.5100357
# Adubo explica: 0.5100357

# varied: 74.128/247.267

74.128/247.267
## [1] 0.2997893
# Varied explica: 0.2997893


###################################################################

# Interações e crítica ao modelo

#1- Carregue o arquivo "a1.2txt" 

#2- Faça todos os procedimentos iniciais para verificar se os dados estão ok

dados= read.table("a1.2.txt",header=T)
attach(dados)
summary(dados)
##      Altura       Volume        Especie         
##  Min.   :63   Min.   :10.20   Length:31         
##  1st Qu.:72   1st Qu.:19.40   Class :character  
##  Median :76   Median :24.20   Mode  :character  
##  Mean   :76   Mean   :30.17                     
##  3rd Qu.:80   3rd Qu.:37.30                     
##  Max.   :87   Max.   :77.00
#3- Descubra se o volume da arvore é afetado pela altura e se a relação é a mesma para as duas espécies?

hist(Volume)

mn = glm(Volume~1) #modelo nulo
m.1= glm(Volume~Altura*Especie)
anova(mn,m.1,test="F")
## Analysis of Deviance Table
## 
## Model 1: Volume ~ 1
## Model 2: Volume ~ Altura * Especie
##   Resid. Df Resid. Dev Df Deviance      F    Pr(>F)    
## 1        30     8106.1                                 
## 2        27      672.1  3   7433.9 99.541 1.032e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Interpretação: O modelo com altura, espécie e interação é muito superior ao nulo, diferença altamente significativa (p= 1.032e-14); ou seja, há dependência clara do volume em relação à altura, espécie e interação entre altura e espécie.

# abandonamos o modelo nulo

m.1= glm(Volume~Altura*Especie)
anova(m.1,test="F")
## Analysis of Deviance Table
## 
## Model: gaussian, link: identity
## 
## Response: Volume
## 
## Terms added sequentially (first to last)
## 
## 
##                Df Deviance Resid. Df Resid. Dev       F    Pr(>F)    
## NULL                              30     8106.1                      
## Altura          1   2901.2        29     5204.9 116.541 2.681e-11 ***
## Especie         1   3798.3        28     1406.6 152.577 1.277e-12 ***
## Altura:Especie  1    734.5        27      672.1  29.504 9.568e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Interpretação: Altura (p=2.681e-11), espécie (p=1.277e-12) e a interação altura:espécie (p=9.568e-06) são altamente significativos. Todas as variáveis são importantes, modelo não deve mais ser simplificado. A interação significativa indica que o efeito da altura no volume depende da espécie.

library(RT4Bio)
rdiagnostic(m.1)

summary(m.1)
## 
## Call:
## glm(formula = Volume ~ Altura * Especie)
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     -48.0478    15.0013  -3.203 0.003474 ** 
## Altura            0.9079     0.1961   4.631 8.22e-05 ***
## Especieb        -95.8767    21.8622  -4.386 0.000159 ***
## Altura:Especieb   1.5586     0.2869   5.432 9.57e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 24.89414)
## 
##     Null deviance: 8106.08  on 30  degrees of freedom
## Residual deviance:  672.14  on 27  degrees of freedom
## AIC: 193.35
## 
## Number of Fisher Scoring iterations: 2

Interpretação (diagnóstico): Pressupostos razoavelmente atendidos; o modelo linear é adequado como primeira aproximação. A partir do summary do modelo, também é possível fazer o diagnóstico para saber se a distribuição escolhida foi adequada para o modelo, 24.894 (672.14/27) é igual à 24.894 (Dispersion parameter for gaussian family). Além disso, Deviance nula = 8106,08 e residual = 672,14, então R² = 1 − 672,14/8106,08 = 0,92 (92%, ótima explicação).

672.14/27 
## [1] 24.89407
# considero a distribuição ok!!!

#4- Construa o grafico

plot(Volume~Altura,ylab="Volume", xlab="Altura (m)",pch=c(16,17)[unclass(Especie)],col=c("green","blue")[unclass(Especie)],ylim=c(0,90),las=1,bty="l")
curve(-48.0478+0.9079*x,add=T,lty=2,col=("green"))
curve((-48.0478-95.8767)+(0.9079+1.5586)*x,add=T,lty=2,col=("blue"))
legend(63,90,c("Espicie A","Espicie B"),pch=c(16,17),col=c("green","blue"),bty="n")

Interpretação: Retas não paralelas (interação): na espécie A (verde) a inclinação = 0,908; na espécie B (azul) a inclinação é maior (= 2,468). Assim, para alturas maiores, a espécie B ganha volume mais rapidamente do que a A.