#----------------------------------------------------------------------------------
#                 PRÁTICA 4 -  ANÁLISE DE SATISFAÇÃO DO PACIENTE
#----------------------------------------------------------------------------------
library(RcmdrMisc)
library(epiDisplay)

setwd(" ")

# Importe o banco de dados "Nível de satisfaçao do paciente.xlsx" e dê o nome de "dados".

library(readxl)
dados <- read_excel("C:/A DATA/PSP5103/2021/Dados/Nível de satisfaçao do paciente.xlsx")
View(dados)

names(dados)


dados$fsetor <- factor(dados$setor, levels=c(0,1), labels=c("público","convênio")) 


#-------------------------------------------------------------------------------------------------------------------
#                                            Análise descritiva
#-------------------------------------------------------------------------------------------------------------------

# Descrevendo uma variável de cada vez

# Medidas descritivas numéricas

numSummary(dados$satisf)
numSummary(dados$idade)
numSummary(dados$gravidade)
numSummary(dados$ansiedade)

table(dados$satisf)  # o que este comando faz?


# Box plots

par(mfrow=c(2,2))

boxplot(dados$satisf, ylab="Nível de satisfação")
boxplot(dados$idade, ylab="Idade (anos)")
boxplot(dados$gravidade, ylab="Índice de gravidade da doença")
boxplot(dados$ansiedade, ylab="Escore de ansiedade")

# Histogramas

hist(dados$satisf)
hist(dados$idade)
hist(dados$gravidade)
hist(dados$ansiedade)

# Tabelas de frequencia

table(dados$fsetor)
prop.table(table(dados$fsetor))

# Faça outros gráficos ...

# Descrevendo associações entre variáveis
#-----------------------------------------------------------------------------------------------# 

# Descrevendo a associação entre Y e as variáveis explicativas quantitativas

# e a associação entre as variáveis explicativas quantitativas

rcorr.adjust(dados[,c("satisf", "ansiedade","gravidade","idade", "setor")], type="pearson", use="complete")

partial.cor(dados[,c("satisf", "ansiedade","gravidade","idade", "setor")], use="complete")


scatterplotMatrix(~satisf+ansiedade+gravidade+idade, regLine=TRUE, smooth=FALSE, 
                  diagonal=list(method="none"), data=dados)

scatterplotMatrix(~satisf+ansiedade+gravidade+idade, regLine=FALSE, smooth=TRUE, 
                  diagonal=list(method="none"), data=dados)

scatterplotMatrix(~satisf+ansiedade+gravidade+idade, regLine=FALSE, smooth=list(var=FALSE), 
                  diagonal=list(method="none"), data=dados)


# Descrevendo a associação entre Y e as variável explicativas qualitativas

par(mfrow=c(1,2))
boxplot(satisf~fsetor, ylab="Nível de satisfação", xlab="Setor de atendimento", data=dados, id=list(method="y"))
plotMeans(dados$satisf, dados$fsetor, error.bars="conf.int", level=0.95, connect=FALSE)

par(mfrow=c(1,1))
stripchart(satisf ~ fsetor, vertical=TRUE, method="stack", ylab="Nível de satisfação", xlab="Setor de atendimento", data=dados)

# Descrevendo a associação entre as variáveis explicativas quantitativas e a qualitativa (setor)

par(mfrow=c(1,3))
boxplot(idade~fsetor, ylab="Idade (anos)", xlab="Setor de atendimento", data=dados, id=list(method="y"))
boxplot(gravidade~fsetor, ylab="Índice de gravidade da doença", xlab="Setor de atendimento", data=dados, id=list(method="y"))
boxplot(ansiedade~fsetor, ylab="Escore de ansiedade", xlab="Setor de atendimento", data=dados, id=list(method="y"))

par(mfrow=c(1,3))
plotMeans(dados$idade, dados$fsetor, error.bars="conf.int", level=0.95, connect=FALSE, ylab="Idade (anos)", xlab="Setor de atendimento")
plotMeans(dados$gravidade, dados$fsetor, error.bars="conf.int", level=0.95, connect=FALSE, ylab="Índice de gravidade da doença", xlab="Setor de atendimento")
plotMeans(dados$ansiedade, dados$fsetor, error.bars="conf.int", level=0.95, connect=FALSE, ylab="Escore de ansiedade", xlab="Setor de atendimento")


# Descrevendo interações 
#----------------------------
# A associação entre Y e as variáveis explicativas quantitativas (idade, gravidade e ansiedade)
# dependem dos níveis da qualitativa (setor) ?

library(lattice)

scatterplot(satisf~idade | fsetor, regLine=FALSE, smooth=FALSE, boxplots=FALSE, by.groups=TRUE,
            ylab="Nível de satisfação", xlab="Idade (anos)", data=dados)

scatterplot(satisf~idade | fsetor, regLine=TRUE, smooth=FALSE, boxplots=FALSE, by.groups=TRUE,
            ylab="Nível de satisfação", xlab="Idade (anos)", data=dados)

scatterplot(satisf~gravidade | fsetor, regLine=TRUE, smooth=FALSE, boxplots=FALSE, by.groups=TRUE,
            ylab="Nível de satisfação", xlab="Índice de gravidade da doença", data=dados)

scatterplot(satisf~ansiedade | fsetor, regLine=TRUE, smooth=FALSE, boxplots=FALSE, by.groups=TRUE,
            ylab="Nível de satisfação", xlab="Escore de ansiedade", data=dados)


#-------------------------------------------------------------------------------------------------------------------
#                                     Análise univariada
#-------------------------------------------------------------------------------------------------------------------

# Verificando normalidade da variável resposta

par(mfrow=c(1,2))
hist(dados$satisf, scale="density", breaks="Sturges", col="darkgray")
qqPlot(dados$satisf, dist="norm")

names(dados)
shapiro.test(dados$satisf)
normalityTest(~satisf, test="lillie.test", data=dados)

warnings()

shapiro.test(dados$idade)
normalityTest(~idade, test="lillie.test", data=dados)

shapiro.test(dados$gravidade)
normalityTest(~gravidade, test="lillie.test", data=dados)

shapiro.test(dados$ansiedade)
normalityTest(~ansiedade, test="lillie.test", data=dados)


# Inicialmente, vamos avaliar a associação entre as variáveis explicativas e a resposta 
# através de modelos de regressão simples, uma variável de cada vez.

simples.1 <- lm(satisf~idade, data=dados)
summary(simples.1)

simples.2 <- lm(satisf~gravidade, data=dados)
summary(simples.2)

simples.3 <- lm(satisf~ansiedade, data=dados)
summary(simples.3)

simples.4 <- lm(satisf~setor, data=dados)         # Note que, se o setor fosse a única var explicativa, seria melhor um teste t-Student.
summary(simples.4)

simples.4 <- lm(satisf~fsetor, data=dados)         # Note que, se o setor fosse a única var explicativa, seria melhor um teste t-Student.
summary(simples.4)
#-------------------------------------------------------------------------------------------------------------------
#                                  Análise de Regressão Múltipla - Stepwise
#-------------------------------------------------------------------------------------------------------------------

# Agora vamos ajustar um modelo de regressão múltipla, com todas as variáveis conjuntamente.

rcorr.adjust(dados[,c("satisf","idade","gravidade","ansiedade","setor")], 
             type="pearson", use="complete")

# Ordem de entrada: idade, ansiedade, gravidade, setor)

mult.1 <- lm(satisf~idade, data=dados)
summary(mult.1)
anova(mult.1)

mult.2 <- update(mult.1, .~. + ansiedade)
summary(mult.2)
anova(mult.2)
anova(mult.1, mult.2)

mult.3 <- update(mult.2, .~. + gravidade)
summary(mult.3)
anova(mult.3)
anova(mult.2, mult.3)

mult.4 <- update(mult.2, .~. + setor)
summary(mult.4)
anova(mult.4)
anova(mult.2, mult.4, test="F")

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

mult.5 <- update(mult.4, .~. + idade*setor)
summary(mult.5)

mult.6 <- update(mult.4, .~. + ansiedade*setor)
summary(mult.6)

mult.7 <- update(mult.4, .~. + idade*ansiedade)
summary(mult.7)

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

par(mfrow=c(1,3))

boxplot(rstandard(mult.4), ylab="Residuo Padronizado")
hist(rstandard(mult.4), scale="density",  xlab="Residuo Padronizado", ylab="Densidade de frequ?ncia", main="")
qqPlot(rstandard(mult.4), dist="norm", xlab="Percentis da N(0,1)", ylab="Resíduo Padronizado")


par(mfrow=c(1,4))

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

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

plot(rstandard(mult.4) ~ dados$ansiedade, ylab="Resíduo Padronizado", xlab="Escore de ansiedade")
abline(h=0, col="blue")
abline(h=-1.96, col="blue", lty=2)
abline(h= 1.96, col="blue", lty=2)

plot(rstandard(mult.4) ~ dados$setor, ylab="Resíduo Padronizado", xlab="Setor de atendimento")
abline(h=0, col="blue")
abline(h=-1.96, col="blue", lty=2)
abline(h= 1.96, col="blue", lty=2)



#-------------------------------------------------------------------------------------------------------------------
#                           Entendendo a interpretação de variáveis categóricas
#-------------------------------------------------------------------------------------------------------------------



# Utilizando a gravidade da doença como uma variável categórica
#-------------------------------------------------------------------------------------------------------------------
par(mfrow=c(1,1))

plot(satisf~gravidade, data=dados)

table(dados$gravidade)
class(dados$gravidade)

dados$gravid <- cut(dados$gravidade, 3, right=FALSE, labels=FALSE)
table(dados$gravid)
class(dados$gravid)

dados$gravid1 <- cut(dados$gravidade, breaks=c(-Inf, 48, 52, Inf), right=TRUE)
table(dados$gravid1)
class(dados$gravid1)

dados$gravid2 <- Recode(dados$gravidade, "41:48=1; 49:52=2; 53:62=3")
table(dados$gravid2)
class(dados$gravid2)

dados$fgravid <- factor(dados$gravid2, levels=c(1,2,3), labels=c("baixa","média","alta"))
table(dados$fgravid)
class(dados$fgravid)

boxplot(satisf~fgravid, ylab="Nível de satisfação", xlab="Gravidade da doença", data=dados, id=list(method="y"))

# Criando variáveis dummy para gravidade

dados$g.baixa <- ifelse(dados$gravid2==1, 1, 0)
dados$g.media <- ifelse(dados$gravid2==2, 1, 0)
dados$g.alta <- ifelse(dados$gravid2==3, 1, 0)

# Inicialmente, vamos ver como fica o modelo só com gravidade:

tapply(dados$satisf, dados$fgravid, mean)

# baixa    média     alta 
# 73.35714 61.36842 49.15385 

73.35714 - 61.36842  
73.35714 - 49.15385 

mod.1 <- lm(satisf~g.media+g.alta, data=dados)
summary(mod.1)

73.35714 - 61.36842  
73.35714 - 49.15385 

mod.2 <- lm(satisf~g.baixa+g.media, data=dados)
summary(mod.2)

# Isto também funciona:

mod.1a <- lm(satisf~fgravid, data=dados)
summary(mod.1a)

mod.2a <- lm(satisf~relevel(fgravid, ref="alta"), data=dados)
summary(mod.2a)

# Agora um modelo com idade e gravidade

mod.3 <- lm(satisf~idade, data=dados)
summary(mod.3)
mod.4 <- lm(satisf~idade+g.media+g.alta, data=dados)
summary(mod.4)
anova(mod.3, mod.4)

# Valores ajustados

# grav baixa -> satisf = 118.8912 - 1.3739*idade
# grav media -> satisf = 118.8912 - 1.3739*idade - 6.3280 = 112.5632 - 1.3739*idade
# grav alta ->  satisf = 118.8912 - 1.3739*idade - 6.9619 = 111.9293 - 1.3739*dade

118.8912 - 6.3280
118.8912 - 6.9619

plot(fitted(mod.4)~idade, ylab="Nível de satisfação ajustado", xlab="Idade (anos)", data=dados)
abline(coef=c(118.8912, -1.3739), col="red")       # gravidade baixa
abline(coef=c(112.5632, -1.3739), col="green")     # gravidade média
abline(coef=c(111.9293, -1.3739), col="blue")      # gravidade alta


#-------------------------------------------------------------------------------------------------------------------
#                           Avaliando a presença de interação
#-------------------------------------------------------------------------------------------------------------------

# Testando interação entre idade (quantitativa) e gravidade (qualitativa)
#------------------------------------------------------------------------

scatterplot(satisf~idade | fgravid, regLine=TRUE, smooth=FALSE, boxplots=FALSE, by.groups=TRUE,
            ylab="Nível de satisfação", xlab="Idade (anos)", data=dados)

mod.4 <- lm(satisf~idade+g.media+g.alta, data=dados)
mod.5 <- lm(satisf~idade + g.media + g.alta + idade*g.media + idade*g.alta, data=dados)
summary(mod.5)
anova(mod.4, mod.5)
Anova(mod.5)

mod.5 <- lm(satisf~idade + fgravid + idade*fgravid, data=dados)
summary(mod.5)



# Supondo que as interações tivessem sido significativas
#---------------------------------------------------------

# Equação

# satisf = 137.4845 - 1.9349*idade - 37.1691*g.media - 11.7883*g.alta + 
#          0.8897*idade:g.media + 0.2597*idade:g.alta

# Valores ajustados

grav baixa -> satisf = 137.4845 - 1.9349*idade 
grav media -> satisf = 137.4845 - 1.9349*idade - 37.1691 + 0.8897*idade = 100.3154 - 1.0452*idade
grav alta  -> satisf = 137.4845 - 1.9349*idade - 11.7883 + 0.2597*idade = 125.6962 - 1.6752*idade

137.4845 - 37.1691
- 1.9349 + 0.8897

137.4845 - 11.7883
- 1.9349+ 0.2597

plot(dados$idade, fitted(mod.5))
abline(coef=c(137.4845, -1.9349), col="red")
abline(coef=c(100.3154, -1.0452), col="green")
abline(coef=c(125.6962, -1.6752), col="blue")


# Testando interação entre setor e gravidade (ambas qualitativas)
#------------------------------------------------------------------------

interaction.plot(dados$fsetor, dados$fgravid, dados$satisf)

interaction.plot(dados$fgravid, dados$fsetor, dados$satisf)

interaction.plot(dados$fgravid, dados$fsetor, dados$satisf, ylab="Satifação média", 
                 xlab="Gravidade", trace.label="Setor", lty=1, col=c("blue", "red", "green"), xpd=FALSE)

with(dados, plotMeans(satisf, fsetor, fgravid, error.bars="none", xlab="Gravidade", ylab="Satisfação média", main="", connect=TRUE, legend.pos="farright"))

with(dados, plotMeans(satisf, fgravid, fsetor, error.bars="none", xlab="Gravidade", ylab="Satisfação média", main="", connect=TRUE, legend.pos="farright"))


mod.6 <- lm(satisf~setor + g.media + g.alta, data=dados)
summary(mod.6)

mod.7 <- lm(satisf~setor + g.media + g.alta + setor*g.media + setor*g.alta, data=dados)
summary(mod.7)
anova(mod.6, mod.7)


# Supondo que as interações tivessem sido significativas
#---------------------------------------------------------

# Escreva a equação da reta ajustada

# Obtenha o valor ajustado para satisfação nos seguintes casos:


 Se setor = público  e grav = baixa -> satisf =  
 Se setor = público  e grav = media -> satisf = 
 Se setor = público  e grav = alta ->  satisf = 
  
 Se setor = convênio  e grav = baixa -> satisf =  
 Se setor = convênio  e grav = media -> satisf = 
 Se setor = convênio  e grav = alta ->  satisf =

# Interprete os valores que você encontrou