#----------------------------------------------------------------------------------
#                   EXERCÍCIO 2 - ANÁLISE PAS
#----------------------------------------------------------------------------------

library(car)
library(RcmdrMisc)

# Altere o diretório de trabalho:

setwd("")

library(readxl)
dados <- read_excel("C:/A DATA/PSP5103/2021/Dados/PAS.xlsx")
View(dados)

# Criando um fator a partir da variável smk:

class(dados$smk)

dados$fsmk <- factor(dados$smk, levels=c(0,1), labels=c("Não fumante", "Fumante"))

class(dados$fsmk)


#----------------------------------------------------------------------------------
#                  MODELO DE REGRESSÃO LINEAR MÚLTIPLA
#----------------------------------------------------------------------------------

scatterplotMatrix(~pas+age+quet, regLine=TRUE, smooth=FALSE, diagonal=list(method="none"),data=dados)

scatterplotMatrix(~pas+age+quet, regLine=FALSE, smooth=list(span=0.5, spread=FALSE), diagonal=list(method="none"),data=dados)

par(mfrow=c(1,2))
boxplot(pas~fsmk, ylab="Pressão arterial sistólica (mm/Hg)", xlab="", data=dados)
plotMeans(dados$pas, dados$fsmk, error.bars="conf.int", level=0.95, connect=FALSE, ylab="M?dia e IC para a Pressão arterial sistólica (mm/Hg)", xlab="", main="")

rcorr.adjust(dados[,c("pas", "quet", "age", "smk")], type="pearson", use="complete")


# ordem de entrada: age, quet, smk

mod.m1 <- lm(pas~age, data=dados) 
summary(mod.m1)  
confint(mod.m1)
anova(mod.m1)

mod.m2 <- lm(pas~age+quet, data=dados) # ou
mod.m2 <- update(mod.m1, .~. + quet) 
summary(mod.m2)
anova(mod.m2)
confint(mod.m2)
anova(mod.m1, mod.m2)    # Teste F-parcial

# A entrada de quet diminuiu o coeficiente de age em mais de 10%
# 1.6045/1.0387
# 1 - 1.0387/1.6045 

mod.m3 <- update(mod.m2, .~. + smk) 
summary(mod.m3)
anova(mod.m3)
confint(mod.m3)
anova(mod.m2, mod.m3)    


SQT <- 3861.6+268.0+772.8+1523.6
SQT

SQM <- 3861.6+268.0+772.8 
SQM

QMM <- SQM/3  
QMM

QMR <- 1523.6/28

F3.28 <- QMM/QMR
F3.28

pf(F3.28, 3, 28, lower.tail = FALSE)

summary(mod.m3)


# Testando interações
#----------------------------

mod.m4 <- update(mod.m3, .~. + age*quet) 
summary(mod.m4)
anova(mod.m3, mod.m4, test="F")

mod.m5 <- update(mod.m3, .~. + age*fsmk) 
summary(mod.m5)
anova(mod.m3, mod.m5, test="F")

mod.m6 <- update(mod.m3, .~. + quet*fsmk) 
summary(mod.m6)
anova(mod.m3, mod.m6, test="F")

# Modelo final: mod.m3


# Análise de resíduos do modelo final
#------------------------------------

# modm.3
#---------

par(mfrow=c(2,2))

plot(rstandard(mod.m3) ~ quet, ylab="Resíduo Padronizado", xlab="Índice de Quetelet (kg/m2)", data=dados)
abline(h=0, col="blue"); abline(h=-1.96, col="blue", lty=3); abline(h= 1.96, col="blue", lty=3)

plot(rstandard(mod.m3) ~ age, xlab="Idade (anos)", ylab="Resíduo Padronizado", data=dados)
abline(h=0, col="blue"); abline(h=-1.96, col="blue", lty=3); abline(h= 1.96, col="blue", lty=3)

plot(rstandard(mod.m3) ~ smk, xlab="smk", ylab="Resíduo Padronizado", data=dados) 
abline(h=0, col="blue"); abline(h=-1.96, col="blue", lty=3); abline(h= 1.96, col="blue", lty=3)

plot(rstandard(mod.m3) ~ fitted(mod.m3), xlab="Valores ajustados", ylab="Resíduo Padronizado"); 
abline(h=0, col="blue"); abline(h=-1.96, col="blue", lty=3); abline(h= 1.96, col="blue", lty=3)

boxplot(rstandard(mod.m3), ylab="Resíduo Padronizado")

hist(rstandard(mod.m3), scale="density",  xlab="Resíduo Padronizado", ylab="Densidade de frequência", main="")

qqPlot(rstandard(mod.m3), dist="norm", xlab="Percentis da N(0,1)", ylab="Resíduo Padronizado")


# Entendendo variáveis dummy
#-----------------------------

# Vamos fazer um experimento

mod.m7 <- lm(pas~age+smk, data=dados) 
summary(mod.m7)

# Note que:
# pas = 48.0496+1.7092age + 10.2944smk
# Não fumantes:  smk=0  pas = 48.0496+1.7092age
# Fumantes:      smk=1  pas = 48.0496+1.7092age + 10.2944 = 58.344+1.7092age

plot(fitted(mod.m7)~ dados$age)
abline(coef=c(48.0496, 1.7092), col="blue")   # Não fumantes
abline(coef=c(58.344, 1.7092), col="red")     # Fumantes