library(epitools)
library(car)
library(RcmdrMisc)

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

#------------------------------------------------------------------------------------------------------

# Importe o banco de dados "LOW.xlsx" e dê o nome de "dados".

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

names(dados)

names(dados) <- casefold(names(dados), upper=F)      # convertendo os nomes das variáveis para letra minúscula

names(dados)                                         # nomes das variáveis 

head(dados)                                          # primeiras linhas do banco de dados


#------------------------------------------------------------------------------------------------------


# Criando fatores para as variáveis categóricas  
#-----------------------------------------------

dados$flow <- factor(dados$low, levels=c(0,1), labels=c("não", "sim")) 
dados$frace <- factor(dados$race, levels=c(1,2,3), labels=c("branca", "preta", "outras"))
dados$fsmoke <- factor(dados$smoke, levels=c(0,1), labels=c("não", "sim"))
dados$fht <- factor(dados$ht, levels=c(0,1), labels=c("não", "sim"))
dados$fui <- factor(dados$ui, levels=c(0,1), labels=c("não", "sim"))


# Categorizando as variáveis quantitativas
#----------------------------------------------------

# lwt
# ---

boxplot(dados$lwt)
hist(dados$lwt)
quantile(dados$lwt)

table(dados$lwt); table(dados$lwt, dados$flow); prop.table(table(dados$lwt, dados$flow))

# lwt1 com 3 categorias: <131; 131 a 150; 151+
# lwt2 com 2 categorias: <139; 140+
# lwt3 com 2 categorias: <121; 122+

dados$lwt1 <- recode(dados$lwt, "0:130=0; 131:150=1; 151:250=2") 
dados$lwt2 <- recode(dados$lwt, "0:150=0; 151:250=1")
dados$lwt3 <- recode(dados$lwt, "0:121=0; 122:250=1")

table(dados$lwt1); table(dados$lwt1, dados$flow); prop.table(table(dados$lwt1, dados$flow))
table(dados$lwt2); table(dados$lwt2, dados$flow); prop.table(table(dados$lwt2, dados$flow))
table(dados$lwt3); table(dados$lwt3, dados$flow); prop.table(table(dados$lwt3, dados$flow))

# ptl
# ---

table(dados$ptl); table(dados$ptl, dados$flow)
# ptl com 2 categorias: 0, 1+

dados$ptl1 <- recode(dados$ptl, "0=0; 1:3=1") 
table(dados$ptl1); table(dados$ptl1, dados$flow); prop.table(table(dados$ptl1, dados$flow))

# ftv
# ---

table(dados$ftv); table(dados$ftv, dados$flow)
# ftv1 com 3 categorias: 0, 1e2, 3+

dados$ftv1 <- recode(dados$ftv, "0=0; 1:2=1; 3:6=2") 
table(dados$ftv1); table(dados$ftv1, dados$flow); prop.table(table(dados$ftv1, dados$flow))

# age
# ---

boxplot(dados$age)
hist(dados$age)
quantile(dados$age)

table(dados$age); table(dados$age, dados$flow)
# age1 com 2 categorias: 14 a 19 e 20 a 45

dados$age1 <- recode(dados$age, "14:19=0; 20:45=1")
table(dados$age1); table(dados$age1, dados$flow); prop.table(table(dados$age1, dados$flow))

# Criando fatores para as novas variáveis qualitativas criadas acima

dados$flwt1 <- factor(dados$lwt1, levels=c(0,1,2), labels=c("0-130", "131-150", "151-250"))
dados$flwt2 <- factor(dados$lwt2, levels=c(0,1), labels=c("0-150", "151-250"))
dados$flwt3 <- factor(dados$lwt3, levels=c(0,1), labels=c("0-121", "122-250"))
dados$fptl1 <- factor(dados$ptl1, levels=c(0,1), labels=c("0", "1-3"))
dados$fftv1 <- factor(dados$ftv1, levels=c(0,1,2), labels=c("0", "1 a 2", "3 a 6"))
dados$fage1 <- factor(dados$age1, levels=c(0,1), labels=c("14-19", "20-45"))

#---------------------------------------------------------------------

#-----------------------------------
# variável explicativa: fsmoke
#-----------------------------------

# Tabela de contingência
#------------------------

tabela.fsmk <- xtabs(~ fsmoke + flow, data=dados)
tabela.fsmk

rowPercents(tabela.fsmk)      # percentuais na linha
colPercents(tabela.fsmk)      # percentuais na coluna

# Teste qui-quadrado
#------------------------

chisq.test(tabela.fsmk, correct=FALSE)
chisq.test(tabela.fsmk, correct=TRUE)

chisq.test(xtabs(~ fsmoke + flow, data=dados), correct=FALSE)

x2.fsmk <- chisq.test(tabela.fsmk, correct=FALSE)
x2.fsmk

x2.fsmk$observed    # valores observados
x2.fsmk$expected    # valores esperados

# Teste exato de fisher
#------------------------

fisher.test(tabela.fsmk, alternative = "two.sided")              
fisher.test(dados$fsmoke, dados$flow, alternative = "two.sided")  

# OR
#------------------------

tabela.fsmk

(30*86)/(44*29)     # OR

epitab(tabela.fsmk, method="oddsratio", verbose=TRUE)               # fornece a OR
epitab(dados$fsmoke, dados$flow, method="oddsratio", verbose=TRUE) # fornece a OR

  30/(29+30)      # probabilidade de EXPOSIÇÃO em doentes: p0, predictor = sim
# 0.5084746
  
  29/(29+30)      # probabilidade de NÃO EXPOSIÇÃO em doentes: p0, predictor = não
  1 - 0.5084746
# 0.4915254

  (30/(29+30))/(29/(29+30))  # chance de EXPOSIÇÃO em doentes - o epitab não mostra
   0.5084746/0.4915254
   30/29
# 1.034483

  44/(86+44)      # probabilidade de EXPOSIÇÃO em não doentes: p1, predictor = sim
# 0.3384615

  86/(86+44)      # probabilidade de NÃO EXPOSIÇÃO em não doentes: p1, predictor = não
# 0.3384615

 (44/(86+44))/(86/(86+44))  # chance de NÃO EXPOSIÇÃO em doentes - o epitab não mostra
  44/86
# 0.5116279

 (30/29)/(44/86)   # OR = chance de EXPOSIÇÃO em doentes / chance de NÃO EXPOSIÇÃO em doentes
  1.034483/0.5116279
# 2.021944

  fisher.test(dados$fsmoke, dados$flow, alternative = "two.sided")  # p.value
  

# RR
#------------------------

tabela.fsmk

(30/(30+44))/(29/(29+86))   # RR  

epitab(tabela.fsmk, method="riskratio", verbose=TRUE)                  # fornece o RR
epitab(x=dados$fsmoke, y=dados$flow, method="riskratio", verbose=TRUE) # fornece o RR


# RP
#------------------------

tabela.fsmk

(30/(30+44))/(29/(29+86))   # RP  

# Utilizar a mesma função do Risco Relativo: 

epitab(tabela.fsmk, method="riskratio", verbose=TRUE)                  # fornece o RR
epitab(x=dados$fsmoke, y=dados$flow, method="riskratio", verbose=TRUE) # fornece o RR

#-----------------------------------
# variável explicativa: fht
#-----------------------------------

table(dados$fht)

tabela.fht <- xtabs(~ fht + flow, data=dados)
tabela.fht

x2.fht <- chisq.test(tabela.fht, correct=FALSE)
x2.fht

x2.fht$observed    # valores observados
x2.fht$expected    # valores esperados 

# Note a Warning message - há um valor esperado menor do que 5
# Necessário usar o Teste Exato de Fisher

fisher.test(tabela.fht, alternative = "two.sided")  

epitab(tabela.fht, method="oddsratio", verbose=TRUE)              
epitab(tabela.fht, method="riskratio", verbose=TRUE)                   

#-----------------------------------
# variável explicativa: fui
#-----------------------------------

tabela.fui <- xtabs(~ fui + flow, data=dados)
tabela.fui

x2.fui <- chisq.test(tabela.fui, correct=FALSE)
x2.fui

x2.fui$observed    # valores observados
x2.fui$expected    # valores esperados

fisher.test(tabela.fui, alternative = "two.sided")  

epitab(tabela.fui, method="oddsratio", verbose=TRUE)                  
epitab(tabela.fui, method="riskratio", verbose=TRUE)    


#-----------------------------------
# variável explicativa: frace
#-----------------------------------

tabela.frace <- xtabs(~ frace + flow, data=dados)
tabela.frace

x2.frace <- chisq.test(tabela.frace, correct=FALSE)
x2.frace

x2.frace$observed    # valores observados
x2.frace$expected    # valores esperados

fisher.test(tabela.frace, alternative = "two.sided")  


# OR
#------------------------

(11*73)/(15*23)               # OR preta/branca
(25*73)/(42*23)               # OR outra/branca

epitab(tabela.frace, method="oddsratio", verbose=TRUE)    

# RR
#------------------------

(11/(11+15)) / (23/(23+73))   # RR preta/branca
(25/(25+42)) / (23/(23+73))   # RR outra/branca

epitab(tabela.frace, method="riskratio", verbose=TRUE)    


# Modelo de regressão logística
#---------------------------------------------------------------------


# modelo só com o intercepto
#---------------------------------------------------------------------

fit.0 <- glm(flow ~ 1, family=binomial(logit), data=dados)
summary(fit.0)

table(dados$flow)


59*log(59) + 130*log(130) - 189*log(189)    # =-117.336   # ln(L) = ln da verossimilhança do modelo só com o intercepto

-117.336*(-2)                               # = 234.672    # -2*ln(L) = Deviance do modelo só com o intercepto


# modelo simples
#---------------------------------------------------------------------

fit.1 <- glm(flow ~ fsmoke, family=binomial(logit), data=dados)
summary(fit.1)
anova(fit.0, fit.1, test="Chi")

#-----------------------------------------

coef(fit.1)
confint(fit.1)

exp(coef(fit.1))
exp(confint(fit.1))

# Compare isto com

epitab(tabela.fsmk, method="oddsratio")                      

epitab(tabela.fsmk, method="")                      

#-----------------------------------------

(30/(30+44))/(29/(29+86))   # RR  


# P(Y=1/X=1)
(exp(-1.0871 + 0.7041)) / (1+(exp(-1.0871 + 0.7041))) 


# P(Y=1/X=0)
(exp(-1.0871 )) / (1+(exp(-1.0871 ))) 


# Compare isto com

epitab(tabela.smk, method="riskratio")                     


# Teste da razão de verossimilhanças
#------------------------------------

234.67                  # Deviance do modelo só com bo    
229.81                  # Deviance do modelo com smoke
234.67 - 229.81         # Diferença de Deviances

234.67/(-2)        # = -117.335     # ln(L) = ln da verossimilhança para o modelo só com bo
229.80/(-2)        # = -114.9       # ln(L) = ln da verossimilhança para o modelo com smoke

anova(fit.0, fit.1, test="Chi")   # Teste da razão de verossimilhanças
anova(fit.1, test="Chi")          # Teste da razão de verossimilhanças

pchisq(4.8674, 1, lower.tail=FALSE)


# Teste de Wald    
#------------------------------------

# Ho: Beta1=0

0.7041 / 0.3196
2*pnorm(2.203066, lower.tail=FALSE)

# Intervalo de confiaça para Beta1    
#------------------------------------

0.7041-1.96*0.3196
0.7041+1.96*0.3196

confint(fit.1)