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

library(car)
library(RcmdrMisc)

# Altere o diretório de trabalho:

setwd("")
setwd("C:/A DATA/PSP5103/2021/Praticas")

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

names(dados)
dim(dados)
head(dados)
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)

#----------------------------------------------------------------------------------
#                            ANÁLISE DESCRITIVA
#----------------------------------------------------------------------------------

# DESCREVENDO UMA VARIÁVEL DE CADA VEZ
#----------------------------------------------------------------------------------

# Variáveis quantitativas:
#-------------------------

# Pas
#-----

par(mfrow=c(1,3))
boxplot(dados$pas, ylab="Pressão arterial sistólica (mm/Hg)", col="orange")
hist(dados$pas, scale="density", col="lightblue", xlab="Pressão arterial sistólica (mm/Hg)", ylab="Densidade de frequência", main="")
qqPlot(dados$pas, dist="norm", ylab="Pressão arterial sistólica (mm/Hg)")                                   

par(mfrow=c(1,1))

shapiro.test(dados$pas)              
ks.test(dados$pas,"pnorm", mean(dados$pas), sd(dados$pas))    

# Quet
#-----

par(mfrow=c(1,3))
boxplot(dados$quet, ylab="Indice de Quetelet (kg/m2)", col="orange")
hist(dados$quet, scale="density", col="lightblue", xlab="Índice de Quetelet (kg/m2)", ylab="Densidade de frequência", main="")
qqPlot(dados$quet, dist="norm", ylab="Índice de Quetelet (kg/m2)")                                   

# Age
#-----

par(mfrow=c(1,3))
boxplot(dados$age, ylab="Idade (anos)", col="orange")
hist(dados$age, scale="density", col="lightblue", xlab="Idade (anos)", ylab="Densidade de frequência", main="")
qqPlot(dados$age, dist="norm", ylab="Idade (anos)")                                   


# Variáveis qualitativas:
#------------------------

#Smk
#---

table(dados$smk)
prop.table(table(dados$smk))


table(dados$fsmk)
prop.table(table(dados$fsmk))


# DESCREVENDO ASSOCIAÇÕES
#----------------------------------------------------------------------------------

# Variáveis explicativas quantitativas:
#---------------------------------------

# Diagramas de dispersão 
#------------------------

par(mfrow=c(1,2))
plot(pas~quet, ylab="Pressão arterial sistólica (mm/Hg)", xlab="Índice de Quetelet (kg/m2)", col="blue", data=dados)
plot(pas~age, ylab="Pressão arterial sistólica (mm/Hg)", xlab="Idade (anos)", col="blue", data=dados)

# Matriz de dispersão
#---------------------

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


# Matriz de correlação
#---------------------

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


# Variáveis explicativas quantitativas:
#---------------------------------------

# estranho!
par(mfrow=c(1,1))
plot(pas~smk, ylab="Pressão arterial sistólica (mm/Hg)", col="blue", data=dados)  


# Box plot e IC 
#---------------------

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="")

# Medidas descritivas
#--------------------------

tapply(dados$pas, dados$fsmk,  ci.numeric, na.rm=TRUE)

# Teste t para duas médias
#--------------------------

# Primeiro o teste para variâncias
tapply(dados$pas, dados$fsmk,  var, na.rm=TRUE)
var.test(pas ~ fsmk, alternative='two.sided', conf.level=.95, data=dados)
bartlett.test(pas ~ fsmk, data=dados)

# Depois o teste para médias
t.test(pas~fsmk, alternative="two.sided", conf.level=0.95, var.equal=TRUE, data=dados)


#----------------------------------------------------------------------------------
#                                  ANÁLISE UNIVARIADA
#----------------------------------------------------------------------------------


mod.1 <- lm(pas~quet, data=dados); anova(mod.1); summary(mod.1); confint(mod.1)
mod.2 <- lm(pas~age, data=dados);  anova(mod.2); summary(mod.2); confint(mod.2)
mod.3 <- lm(pas~smk, data=dados);  anova(mod.3); summary(mod.3); confint(mod.3)

mod.3a <- lm(pas~fsmk, data=dados); anova(mod.3a); summary(mod.3a); confint(mod.3a)


# Análise de resíduos
#--------------------

# mod.1
#-------

par(mfrow=c(2,3))

plot(pas~quet, ylab="Pressão arterial sistólica (mm/Hg)", xlab="Índice de Quetelet (kg/m2)", col="blue", data=dados)
abline(mod.1, col="red")

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

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

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

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

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


# mod.2
#-------

par(mfrow=c(2,3))

plot(pas~age, ylab="Pressão arterial sistólica (mm/Hg)", xlab="Idade (anos)", col="blue", data=dados)
abline(mod.2, col="red")

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

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

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

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

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

# mod.3
#-------

par(mfrow=c(2,3))

plot(pas~smk, ylab="Pressão arterial sistólica (mm/Hg)", xlab="smk", col="blue", data=dados)
abline(mod.3, col="red")

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

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

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

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

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


# Intervalo de confiança para E[Yi]=ui quando quet=30
#-----------------------------------------------------

anova(mod.1)

qmr <- anova(mod.1)[2,3] 
qmr

yh <- coef(mod.1)[1] + coef(mod.1)[2]*30
yh

dp.yh <- sqrt(qmr *  ( 1/32 + (30 - mean(dados$quet))^2  /  (31*var(dados$quet)) )  )
dp.yh

t <- qt(c(0.025), df=30, lower.tail=FALSE)
t

yh - t*dp.yh
yh + t*dp.yh

# ou utilizando a função predict do R

nd1 <- data.frame(quet=30); nd1

mod1.conf <- predict.lm( mod.1, newdata=nd1, interval="conf", se.fit=TRUE ); mod1.conf
mod1.conf

# Teste de hipóteses para ui quando quet=30
#-----------------------------------------------------

# Ho: ui=150
# Ha: ui<>150

to <- (mean(dados$pas) + coef(mod.1)[2]*(30 - mean(dados$quet))-150)/dp.yh
to

# ou

to <- (yh - 150)/dp.yh
to

pt(to, 30, lower.tail=TRUE)


# Intervalo de predição para Yi quando quet=30
#-----------------------------------------------------

nd1 <- data.frame(quet=30); nd1

mod1.pred <- predict.lm( mod.1, newdata=nd1, interval="pred", se.fit=TRUE ); mod1.pr
mod1.pred


# Intervalos de estimação e predição para Yi quando age=50
#---------------------------------------------------------

nd2 <- data.frame(age=50); nd2
mod2.conf <- predict.lm( mod.2, newdata=nd2, interval="conf", se.fit=TRUE ); mod2.conf
mod2.pr <- predict.lm( mod.2, newdata=nd2, interval="pred", se.fit=TRUE ); mod2.pr

# Teste de hipóteses para ui quando age=50
#---------------------------------------------------------

# Ho: ui=150
# Ha: ui<>150

yh <- 139.3166
dp.yh <- 1.809169

to <- (yh - 150)/dp.yh
to

pt(to, 30, lower.tail=TRUE)


# Intervalos de estimação e predição para Yi em fumantes e não fumantes
#----------------------------------------------------------------------

nd3 <- data.frame(fsmk=c("Não fumante", "Fumante")); nd3
mod3.conf <- predict.lm( mod.3, newdata=nd3, interval="conf", se.fit=TRUE ); mod3.conf
mod3.pr <- predict.lm( mod.3, newdata=nd3, interval="pred", se.fit=TRUE ); mod3.pr


# Teste de hipóteses para ui em fumantes e não fumantes
#----------------------------------------------------------------------

# Não Fumantes
#-----------------------------

# Ho: ui=150
# Ha: ui<>150

yh <- 140.8000
dp.yh <- 3.661472

to <- (yh - 150)/dp.yh
to

pt(to, 30, lower.tail=TRUE)

# Fumante
#-----------------------------

# Ho: ui=150
# Ha: ui<>150

yh <- 147.8235
dp.yh <- 3.439354

to <- (yh - 150)/dp.yh
to

pt(to, 30, lower.tail=TRUE)