#==========================================================
#
#  Regressão Linear Múltipla
#    Resposta: PIO
#    Covariáveis: idade e 
#                 gênero (0- mulher e 1-homem)
#===========================================================
#  
#   Lendo Banco de Dados
  ex1<-read.table("dados_aula4.txt",h=T)
  ex1
#  Acertando a classe das variáveis
  Y<-ex1$pio
  X1<-ex1$idade
  X2<-ex1$sexo
  class(Y)
  Y<-as.numeric(Y)
  class(X1)
  X1<-as.numeric(X1)
  class(X2)
  X2<-factor(X2)
#
#  Regressão Linear Simples para Idade
#
  mod1 <- lm(Y~X1)
  summary(mod1)
#
#  Regressão Linear Simples para Gênero
#
  mod2 <- lm(Y~X2)
  summary(mod2)
#
#  Regressão Linear Múltipla: Idade + Gênero
#
  mod3 <- lm(Y~X1+X2)
  summary(mod3)
#
#  Verificando as suposções do Modelo
#  Gráfico 1: variância constante e linearidade
#  Gráfico 2: normalidade
#  Gráfico 3: somente variância constante
#  Gráfico 4: observações atípicas
#
  par(mfrow=c(2,2))
  plot(mod3, which=c(1:4),pch=16, add.smooth=T)
#
#  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(mod3))~fitted(mod3)))
#
#  Teste para verificar normalidade 
   shapiro.test(mod3$resid)
#
#  Gráfico para verificar independência (na sequência de 
#    coleta dos dados), se fizer sentido.
#
   par(mfrow=c(1,1))
   plot(mod2$residuals,pch=16,ylab="Residuals")
#
#  Como aparentemente as suposições foram validadas,
#   Porque gênero perdeu a significância na regressão múltipla?
# 
  levels(X2) <- c("mulher","homem")
  boxplot(X1~X2, xlab="Gênero", ylab="idade")
# 
# Homens são mais velhos que as mulheres. Isto significa que 
# idade é um fator de confusão para a associação entre Pio e gênero.
#
#
##################################################################################
#      TÓPICOS ESPECIAIS
#
#  Excluindo o paciente 25
#
  ex2<-ex1[-25,]
  mod7<-gamlss(ex2$pio~ex2$idade+ex2$sexo,sigma.formula= ~ex2$idade,family=NO)
  plot(mod7)
  wp(mod7)
  summary(mod7)
  mod8<-lm(ex2$pio~ex2$idade+ex2$sexo)
  summary(mod8)
  par(mfrow=c(2,2))
  plot(mod8)
  shapiro.test(mod8$resid)
#
#=========================================================================================