#
#  Exemplo de Sala de aula: Regressão Linear Simples
#      Pressão Intra-Ocular 
#
#==========================================================
#   Lendo o Banco de Dados de PIO e Idade
#
 reg<-read.table("dados_aula3.txt",h=T)
 reg
 dim(reg)
#
#    Explorando o Banco de Dados
#
 is.data.frame(reg)
 names(reg)
 is.numeric(reg$idade)
 is.numeric(reg$pio) # ambas variáveis são numéricas
 summary(reg)
 sd(reg$idade)
 sd(reg$pio)
#
#    Gráfico de Dispersão
# 
 require(ggplot2) # gráficos mais simpáticos
 x<-reg$idade
 y<-reg$pio
 plot(x,y,xlab="idade", ylab="PIO")
 qplot(x,y,xlab="idade", ylab="PIO") # forma alternativa
#
#   Correlação entre PIO e Idade
#----------------------------------------------------------------
 cor(x,y)
#
 cor.test(x,y) # teste e IC para \rho
#
# Regressão Linear Simples
#--------------------------------------------------------------------
 mod1<-lm(y~x) # ajuste do modelo regressão linear simples
 summary(mod1) # teste t para beta
#
#  Gráfico da reta ajustada
#
 plot(x,y,abline(mod1$coefficients[1], mod1$coefficients[2]),xlab="idade", ylab="PIO")
# lines(lowess(x,y),col=2)
 legend(60,17,"PIO = 6,64 + 0,23 idade",bty="n")
 legend(60.65,17.15,"^",bty="n")
#
#  Gráficos para verificar a adequação do ajuste
#
 par(mfrow=c(2,2))
 plot(mod1, which=c(1:4), add.smooth=T,pch=20)
# 
 par(mfrow=c(1,1))
 boxplot(mod1$resid)
# Teste de Normalidade
 shapiro.test(mod1$resid)
#
#  Um teste simples para verificar variância não-constante
#   Este teste não é completamente correto pois algum peso
#   deve ser utilizado e os graus de liberdade ajustados
#   Veja Faraway (2004)
#
 summary(lm(abs(residuals(mod1))~fitted(mod1)))
#  valor-p=0,064 para a hipótese de homocedasticidade
#
##################################################################################
#      TÓPICOS ESPECIAIS 
#
#  1- Excluir o  caso 25 (atípico)
#     Somente válido se houver uma justificativa para excluí-lo
#
 reg1<-reg[-25,]
 mod2<-lm(reg1$pio~reg1$idade)
 summary(mod2)
 par(mfrow=c(2,2))
 plot(mod2, which=c(1:4), add.smooth=T,pch=20)
 shapiro.test(mod2$resid)
 summary(lm(abs(residuals(mod2))~fitted(mod2)))
#
# Observe que a qualidade do ajuste melhorou em todos os aspectos.
# A conclusão continua mostrando uma associação entre PIO e idade.
# No entanto, observe que beta^ passou de 0,23 para 0,27.
#
#----------------------------------------------------------------------------
#   2- Modelar também a dispersão
#       Como existe uma leve evidência contra a suposição de homocedasticidade
#       vamos modelar também o desvio-padrão do erro em função de idade.
#       Ou seja, vamos testar se a dispersão aumenta com o aumento da idade.
#       Para fazer isso precisamos utilizar o pacote gamlss
#       Não vamos tirar o paciente 25.
#
 require(gamlss)
 mod3<-gamlss(reg$pio~reg$idade,sigma.formula= ~reg$idade,family=NO)
 summary(mod3)             
 wp(mod3)
 plot(mod3, which=c(1:4), add.smooth=T,pch=20)
 shapiro.test(mod3$resid)
 summary(lm(abs(residuals(mod3))~fitted(mod3)))
#
#  Não foi significativo o efeito de idade no modelo da dispersão. Os gráficos
#  confirmam não exister evidência contra a suposição de homocedasticidade.
#
####################################################################################