#--------------------------------------------------------------------------------------------------------
# EXERCÍCIO 3 - REGRESSÃO LOGÍSTICA
# BANCO LOW
#--------------------------------------------------------------------------------------------------------
# Altere o diretório de trabalho:
setwd("")
# Importe o banco de dados "LOW" e de o nome de "dados"
# Utilize o Menu
# File -> Import Dataset -> From Excel...
# Instale os pacotes:
install.packages("epiDisplay")
library(epiDisplay)
install.packages("MASS")
library(MASS)
install.packages("nnet")
library (nnet)
install.packages("ResourceSelection")
library (ResourceSelection)
install.packages("RcmdrMisc")
library(RcmdrMisc)
install.packages("epitools")
library(epitools)

# Conhecendo as variáveis
names(dados) 
names(dados) <- casefold(names(dados), upper=F) # convertendo os nomes das variaveis para letras minusculas
names(dados) # nomes das variaveis 
head(dados) # primeiras linhas do banco de dados

# Criando fatores 
#--------------------------------------------------------------------------------------------------------
class(dados$low)
dados$flow <- factor(dados$low, levels=c(0,1), labels=c("nao low", "sim low")) 
dados$flow
table(dados$flow)

class(dados$smoke)
dados$fsmoke <- factor(dados$smoke, levels=c(0,1), labels=c("nao smoke", "sim smoke"))
dados$fsmoke
table(dados$fsmoke)

class(dados$ht)
dados$fht <- factor(dados$ht, levels=c(0,1), labels=c("nao ht", "sim ht"))
dados$fht
table(dados$fht)

class(dados$ui)
dados$fui <- factor(dados$ui, levels=c(0,1), labels=c("nao ui", "sim ui"))
dados$fui
table(dados$fui)

class(dados$race)
dados$frace <- factor(dados$race, levels=c(1,2,3), labels=c("branca", "negra", "outras"))
dados$frace
table(dados$frace)

View(dados)

#--------------------------------------------------------------------------------------------------------
# LOW E SMOKE
#--------------------------------------------------------------------------------------------------------
# Criando uma tabela de contingencia
#--------------------------------------------------------------------------------------------------------
tabela <- xtabs(~ fsmoke + flow, data=dados)
tabela
(30/(30+44)) #RISCO EXPOSTOS
(29/(29+86)) #RISCO NÃO EXPOSTOS
(30/(30+44)) / (29/(29+86)) #RR

(30/(30+29)) #CHANCE EXPOSTOS
(44/(44+86)) #CHANCE NÃO EXPOSTOS 
(30*86) / (44*29) #OR

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

# Teste de associação pelo qui-quadrado
#--------------------------------------------------------------------------------------------------------
chisq.test(tabela, correct=FALSE) #QUI-QUADRADO SEM CORREÇÃO
chisq.test(tabela, correct=TRUE) #QUI-QUADRADO COM CORREÇÃO

# OR e RR
#--------------------------------------------------------------------------------------------------------
# Cálculo do Risco Relativo
result <- epitab(tabela, method = "riskratio", conf.level = 0.95)
print(result)

# Cálculo da Odds Ratio
result <- epitab(tabela, method = "oddsratio", conf.level = 0.95)
print(result)

# Modelo de regressão logística
#--------------------------------------------------------------------------------------------------------
# Só com o intercepto
fit.smoke0 <- glm(flow ~ 1, family=binomial(logit), data=dados)
summary(fit.smoke0)

# Modelo simples
fit.smoke1 <- glm(flow ~ fsmoke, family=binomial(logit), data=dados)
summary(fit.smoke1)
confint(fit.smoke1) # IC para Beta
exp(coef(fit.smoke1)) # OR = Exponencial do Beta1
exp(confint(fit.smoke1)) # IC para OR = Exponencial do IC para Beta1

# CÁLCULO DO LOGARITIMO DA FUNÇÃO DE VEROSSIMILHANÇA DO MODELO SÓ COM BETA 0:
59*log(59) + 130*log(130) - 189*log(189)
-117.336 * (-2)
-229.80 / 2
234.67-229.80 # OBTENHO O VALOR DA DEVIANCE

# AS DUAS DEVIANCES (234.67 e 229.81), A DIFERENÇA (4.8674) E O VALOR DO p DO QUI-QUADRADO DA RAZÃO DE VEROSSIMILHANÇA (0.02737)
anova(fit.smoke0, fit.smoke1, test="Chi") 

# A DIFERENÇA DAS DEVIANCES (4.8674) E AS DUAS DEVIANCES (234.67 e 229.81)
anova(fit.smoke1, test="Chi")

# VALOR DO p DO TESTE QUI-QUADRADO DA RAZÃO DE VEROSSIMILHANÇA - Um g.l. (duas categorias)
pchisq(4.8674, 1, lower.tail = FALSE) 

# Teste de Wald
# Encontrando o valor da z 
0.7041 / 0.3196
# Encontrando o p-valor do teste de Wald que aparece no modelo simples
2*pnorm(2.203066, lower.tail=FALSE)

# CALCULO DO MODELO
# Probabilidade doente e expostos
py1x1 <- 1 / (1+exp(1.0871-0.7041))
py1x1
0.4054035
# Probabilidade doente e não exposto 
py1x0 <- 1 / (1+exp(1.0871))
py1x0
0.2521648
# Risco relativo 
0.4054035 / 0.2521648
# OR 
exp(0.7041)

#--------------------------------------------------------------------------------------------------------
# LOW E HT
#--------------------------------------------------------------------------------------------------------
# Criando uma tabela de contingencia
#--------------------------------------------------------------------------------------------------------
tabela <- xtabs(~ fht + flow, data=dados)
tabela

(7/(7+5)) #RISCO EXPOSTOS 
(52/(52+125)) #RISCO NÃO EXPOSTOS 
(7/(7+5)) / (52/(52+125)) #RR

(7/(7+52)) #CHANCE EXPOSTOS
(5/(5+125)) #CHANCE NÃO EXPOSTOS
(7*125) / (5*52) #OR

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

# Teste de associação pelo qui-quadrado
#--------------------------------------------------------------------------------------------------------
chisq.test(tabela, correct=FALSE) #QUI-QUADRADO SEM CORREÇÃO
chisq.test(tabela, correct=TRUE) #QUI-QUADRADO COM CORREÇÃO

# OR e RR
#--------------------------------------------------------------------------------------------------------
# Cálculo do Risco Relativo
result <- epitab(tabela, method = "riskratio", conf.level = 0.95)
print(result)

# Cálculo da Odds Ratio
result <- epitab(tabela, method = "oddsratio", conf.level = 0.95)
print(result)

# Modelo de regressao logística
#--------------------------------------------------------------------------------------------------------
# Só com o intercepto
fit.ht0 <- glm(flow ~ 1, family=binomial(logit), data=dados)
summary(fit.ht0)

# Modelo simples
fit.ht1 <- glm(flow ~ fht, family=binomial(logit), data=dados)
summary(fit.ht1)
confint(fit.ht1) # IC para Beta

exp(coef(fit.ht1)) # OR = Exponencial do Beta1
exp(confint(fit.ht1)) # IC para OR = Exponencial do IC para Beta1

# CÁLCULO DO LOGARITIMO DA FUNÇÃO DE VEROSSIMILHANÇA DO MODELO SÓ COM BETA 0:
59*log(59) + 130*log(130) - 189*log(189)
-117.336 * (-2)
-230.65 / 2
234.672-230.65 #OBTENHO O VALOR DA DEVIANCE

# AS DUAS DEVIANCES (234.67 e 230.65), A DIFERENÇA (4.0221) E O VALOR DO p DO QUI-QUADRADO DA RAZÃO DE VEROSSIMILHANÇA (0.04491)
anova(fit.ht0, fit.ht1, test="Chi")

# A DIFERENÇA DAS DEVIANCES (4.0221) E AS DUAS DEVIANCES (234.67 e 230.65)
anova(fit.ht1, test="Chi")

#VALOR DO p DO TESTE QUI-QUADRADO DA RAZÃO DE VEROSSIMILHANÇA - Um g.l. (duas categorias)
pchisq(4.0221, 1, lower.tail = FALSE)

# Teste de Wald
# Encontrando o valor da z 
1.2135 / 0.6083
# Encontrando o p-valor do teste de Wald que aparece no modelo simples
2*pnorm(1.994904, lower.tail=FALSE) 

# CALCULO DO MODELO
# Probabilidade doente e expostos
py1x1 <- 1 / (1+exp(0.8771-1.2135))
py1x1
0.5833158
# Probabilidade doente e não exposto 
py1x0 <- 1 / (1+exp(0.8771))
py1x0
0.2937791
# Risco relativo 
0.5833158 / 0.2937791
# OR 
exp(1.2135)

#--------------------------------------------------------------------------------------------------------
# LOW E UI
#--------------------------------------------------------------------------------------------------------
# Criando uma tabela de contingencia
#--------------------------------------------------------------------------------------------------------
tabela <- xtabs(~ fui + flow, data=dados)
tabela

(14/(14+14)) #RISCO EXPOSTOS 
(45/(45+116)) #RISCO NÃO EXPOSTOS
(14/(14+14)) / (45/(45+116)) #RR

(14/(14+45)) #CHANCE EXPOSTOS
(14/(14+116)) #CHANCE NÃO EXPOSTOS
(14*116) / (14*45) #OR

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

# Teste de associação pelo qui-quadrado
#--------------------------------------------------------------------------------------------------------
chisq.test(tabela, correct=FALSE) #QUI-QUADRADO SEM CORREÇÃO
chisq.test(tabela, correct=TRUE) #QUI-QUADRADO COM CORREÇÃO

# OR e RR
#--------------------------------------------------------------------------------------------------------
# Cálculo do Risco Relativo
result <- epitab(tabela, method = "riskratio", conf.level = 0.95)
print(result)

# Cálculo da Odds Ratio
result <- epitab(tabela, method = "oddsratio", conf.level = 0.95)
print(result)

# Modelo de regressao logistica
#--------------------------------------------------------------------------------------------------------
# Só com o intercepto
fit.ui0 <- glm(flow ~ 1, family=binomial(logit), data=dados)
summary(fit.ui0)

# Modelo simples
fit.ui1 <- glm(flow ~ fui, family=binomial(logit), data=dados)
summary(fit.ui1)
confint(fit.ui1) # IC para Beta

exp(coef(fit.ui1)) # OR = Exponencial do Beta1
exp(confint(fit.ui1)) # IC para OR = Exponencial do IC para Beta1

# CÁLCULO DO LOGARÍTIMO DA FUNÇÃO DE VEROSSIMILHANÇA DO MODELO SÓ COM BETA 0:
59*log(59) + 130*log(130) - 189*log(189)
-117.336 * (-2)
-229.60 / 2
234.672-229.60 #OBTENHO O VALOR DA DEVIANCE

# AS DUAS DEVIANCES (234.67 e 229.60), A DIFERENÇA (5.0761) E O VALOR DO p DO QUI-QUADRADO DA RAZÃO DE VEROSSIMILHANÇA (0.02426)
anova(fit.ui0, fit.ui1, test="Chi")

# A DIFERENÇA DAS DEVIANCES (5.0761) E AS DUAS DEVIANCES (234.67 e 229.60)
anova(fit.ui1, test="Chi")

#VALOR DO p DO TESTE QUI-QUADRADO DA RAZÃO DE VEROSSIMILHANÇA - Um g.l. (duas categorias)
pchisq(5.0761, 1, lower.tail = FALSE)

# Teste de Wald
# Encontrando o valor da z 
0.9469 / 0.4168
# Encontrando o p-valor do teste de Wald que aparece no modelo simples
2*pnorm(2.271833, lower.tail=FALSE)

# CALCULO DO MODELO
# Probabilidade doente e expostos
py1x1 <- 1 / (1+exp(-0.9469 + 0.9469))
py1x1
0.5 
# Probabiliade doente e não exposto 
py1x0 <- 1 / (1+exp(0.9469))
py1x0
0.2795087
# Risco relativo 	
0.5 / 0.2795087
# OR 
exp(0.9469)

#--------------------------------------------------------------------------------------------------------
# LOW E RACE
#--------------------------------------------------------------------------------------------------------
# Criando uma tabela de contingencia
#--------------------------------------------------------------------------------------------------------
tabela <- xtabs(~ frace + flow, data=dados)
tabela

(23/(23+73)) # RISCO BRANCOS (REFERÊNCIA)
(11/(11+15)) # RISCO NEGROS
(25/(25+42)) # RISCO OUTROS
(11/(11+15)) / (23/(23+73)) # RR NEGROS
(25/(25+42)) / (23/(23+73)) # RR OUTROS

# CHANCE BRANCOS (REFERÊNCIA)
# CHANCE NEGROS
# CHANCE OUTROS
(11*73) / (15*23) # OR NEGROS
(25*73) / (42*23) # OR OUTROS

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

# Teste qui-quadrado
#--------------------------------------------------------------------------------------------------------
chisq.test(tabela, correct=FALSE) #QUI-QUADRADO SEM CORREÇÃO
chisq.test(tabela, correct=TRUE) #QUI-QUADRADO COM CORREÇÃO

# OR e RR
#--------------------------------------------------------------------------------------------------------
# Cálculo do Risco Relativo
result <- epitab(tabela, method = "riskratio", conf.level = 0.95)
print(result)

# Cálculo da Odds Ratio
result <- epitab(tabela, method = "oddsratio", conf.level = 0.95)
print(result)

# Modelo de regressao logistica
#--------------------------------------------------------------------------------------------------------
# Só com o intercepto
fit.race0 <- glm(flow ~ 1, family=binomial(logit), data=dados)
summary(fit.race0)

# Modelo simples
fit.race1 <- glm(flow ~ frace, family=binomial(logit), data=dados)
summary(fit.race1)
confint(fit.race1) # IC para Beta

exp(coef(fit.race1)) # OR = Exponencial do Beta1
exp(confint(fit.race1)) # IC para OR = Exponencial do IC para Beta1

# CÁLCULO DO LOGARÍTIMO DA FUNÇÃO DE VEROSSIMILHANÇA DO MODELO SÓ COM BETA 0:
59*log(59) + 130*log(130) - 189*log(189)
-117.336 * (-2)
-229.66 / 2
234.672-229.66 # OBTENHO O VALOR DA DEVIANCE

# AS DUAS DEVIANCES (234.67 e 229.66), A DIFERENÇA (5.0104) E O VALOR DO p DO QUI-QUADRADO DA RAZÃO DE VEROSSIMILHANÇA (0.08166)
anova(fit.race0, fit.race1, test="Chi")

# A DIFERENÇA DAS DEVIANCES (5.0104) E AS DUAS DEVIANCES (234.67 e 229.66)
anova(fit.race1, test="Chi")

#VALOR DO p DO TESTE QUI-QUADRADO DA RAZÃO DE VEROSSIMILHANÇA - Dois g.l. (três categorias)
pchisq(5.0104, 2, lower.tail = FALSE)

# Teste de Wald - RACE NEGRA
# Encontrando o valor da z 
0.8448 / 0.4634
# Encontrando o p-valor do teste de Wald que aparece no modelo simples
2*pnorm(1.823047, lower.tail=FALSE) 

# Teste de Wald - RACE OUTROS
# Encontrando o valor da z 
0.6362 / 0.3478
# Encontrando o p-valor do teste de Wald que aparece no modelo simples
2*pnorm (1.829212, lower.tail=FALSE)  

# CALCULO DO MODELO
# Probabilidade doente e expostos - NEGROS
py1x1 <- 1 / (1+exp(1.1550-0.8448))
py1x1
0.4230659
# Probabilidade doente e expostos - OUTROS
py1x1 <- 1 / (1+exp(1.1550-0.6362))
py1x1
0.3731329
# Probabilidade doente e não exposto 
py1x0 <- 1 / (1+exp(1.1550))
py1x0
0.239577
# Risco relativo - NEGROS
0.4230659 / 0.239577
# Risco relativo - OUTROS
0.3731329 / 0.239577
# OR - NEGROS
exp(0.8448)
# OR - OUTROS
exp(0.6362)

#--------------------------------------------------------------------------------------------------------
# LOW E AGE
#--------------------------------------------------------------------------------------------------------
table(dados$age)
xtabs(~age + flow, data=dados)

# CATEGORIZANDO <20 (baseline = 0) e 20-45
dados$age1 <- recode(dados$age, "14:19=0; 20:45=1") 
dados$fage1 <- factor(dados$age1, levels=c(0,1), labels=c("Menor que 20 anos", "20 a 45 anos"))
table(dados$fage1)
View(dados)

# Criando uma tabela de contingencia
#--------------------------------------------------------------------------------------------------------
# TABELA CATEGORIZADA
table(dados$fage1)
tabela <- xtabs(~fage1 + flow, data=dados)
tabela

(44/(44+94)) # RISCO MENOR QUE 20 ANOS (REFERÊNCIA)
(15/(15+36)) # RISCO DE 20 - 45 ANOS
(44/(44+94)) / (15/(15+36)) # RR 

(44/(44+15)) # CHANCE MENOR QUE 20 ANOS (REFERÊNCIA)
(94/(94+36)) # CHANCE DE 20 - 45 ANOS
(44*36) / (94*15) # OR 

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

# Teste de associação pelo qui-quadrado
#--------------------------------------------------------------------------------------------------------
chisq.test(tabela, correct=FALSE) #QUI-QUADRADO SEM CORREÇÃO
chisq.test(tabela, correct=TRUE) #QUI-QUADRADO COM CORREÇÃO

# OR e RR
#--------------------------------------------------------------------------------------------------------
# Cálculo do Risco Relativo
result <- epitab(tabela, method = "riskratio", conf.level = 0.95)
print(result)

# Cálculo da Odds Ratio
result <- epitab(tabela, method = "oddsratio", conf.level = 0.95)
print(result)
# Modelo de regressao logistica
#--------------------------------------------------------------------------------------------------------
# Só com o intercepto
fit.age0 <- glm(flow ~ 1, family=binomial(logit), data=dados)
summary(fit.age0)

# Modelo simples
fit.age1 <- glm(flow ~ fage1, family=binomial(logit), data=dados)
summary(fit.age1)
confint(fit.age1) # IC para Beta

exp(coef(fit.age1)) # OR = Exponencial do Beta1
exp(confint(fit.age1)) # IC para OR = Exponencial do IC para Beta1 

# CÁLCULO DO LOGARÍTIMO DA FUNÇÃO DE VEROSSIMILHANÇA DO MODELO SÓ COM BETA 0:
59*log(59) + 130*log(130) - 189*log(189)
-117.336 * (-2)
-234.57 / 2
234.672-234.57 # OBTENHO O VALOR DA DEVIANCE

# AS DUAS DEVIANCES (234.67 e 234.56), A DIFERENÇA (0.10674) E O VALOR DO p DO QUI-QUADRADO DA RAZÃO DE VEROSSIMILHANÇA (0.7439)
anova(fit.age0, fit.age1, test="Chi")

# A DIFERENÇA DAS DEVIANCES (0.10674) E AS DUAS DEVIANCES (234.67 e 234.56)
anova(fit.age1, test="Chi")

#VALOR DO p DO TESTE QUI-QUADRADO DA RAZÃO DE VEROSSIMILHANÇA - Um g.l. (duas categorias)
pchisq(0.10674, 1, lower.tail = FALSE)

# Teste de Wald
# Encontrando o valor da z 
0.1164 / 0.3575
# Encontrando o p-valor do teste de Wald que aparece no modelo simples
2*pnorm(0.3255944, lower.tail=FALSE)

# CALCULO DO MODELO
# Probabilidade doente e expostos
py1x1 <- 1 / (1+exp(0.8755-0.1164))
py1x1
0.3188417
# Probabilidade doente e não exposto 
py1x0 <- 1 / (1+exp(0.8755))
py1x0
0.2941112
# Risco relativo 
0.3188417 / 0.2941112
# OR 
exp(0.1164)

#--------------------------------------------------------------------------------------------------------
# LOW E LWT
#--------------------------------------------------------------------------------------------------------
table(dados$lwt)
xtabs(~lwt + flow, data=dados)

# CATEGORIZANDO <131, 131-150 e 151e+
dados$lwt1 <- recode(dados$lwt, "80:130=1; 131:150=2; 151:250=3") 
dados$flwt1 <- factor(dados$lwt1, levels=c(1,2,3), labels=c("Menor que 131", "131 a 150", "151 e mais"))
table(dados$flwt1)
View(dados)

# Criando uma tabela de contingencia
#--------------------------------------------------------------------------------------------------------
# TABELA DE CONTINGÊNCIA CATEGORIZADA
tabela <- xtabs(~ flwt1 + flow, data=dados)
tabela

(45/(45+78)) # RISCO <131 (REFERÊNCIA)
(7/(7+23)) # RISCO 131-150
(7/(7+29)) # RISCO 151 e + 
(7/(7+23)) / (45/(45+78)) # RR 131-150
(7/(7+29)) / (45/(45+78)) # RR 151 e +  

(7*78) / (23*45) # OR 131-150 
(7*78) / (29*45) # OR 151 e +

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

# Teste de associação pelo qui-quadrado
#--------------------------------------------------------------------------------------------------------
chisq.test(tabela, correct=FALSE) # QUI-QUADRADO SEM CORREÇÃO
chisq.test(tabela, correct=TRUE) # QUI-QUADRADO COM CORREÇÃO

# OR e RR
#--------------------------------------------------------------------------------------------------------
# Cálculo do Risco Relativo
result <- epitab(tabela, method = "riskratio", conf.level = 0.95)
print(result)

# Cálculo da Odds Ratio
result <- epitab(tabela, method = "oddsratio", conf.level = 0.95)
print(result)
# Modelo de regressao logística
#--------------------------------------------------------------------------------------------------------
# Só com o intercepto
fit.lwt0 <- glm(flow ~ 1, family=binomial(logit), data=dados)
summary(fit.lwt0)

# Modelo simples
fit.lwt1 <- glm(flow ~ flwt1, family=binomial(logit), data=dados)
summary(fit.lwt1)
confint(fit.lwt1) # IC para Beta

exp(coef(fit.lwt1)) # OR = Exponencial do Beta1
exp(confint(fit.lwt1)) # IC para OR = Exponencial do IC para Beta1

# CÁLCULO DO LOGARÍTIMO DA FUNÇÃO DE VEROSSIMILHANÇA DO MODELO SÓ COM BETA 0:
59*log(59) + 130*log(130) - 189*log(189)
-117.336 * (-2)
-229.61 / 2
234.672-229.61 # OBTENHO O VALOR DA DEVIANCE

# AS DUAS DEVIANCES (234.67 e 229.62), A DIFERENÇA (5.057) E O VALOR DO p DO QUI-QUADRADO DA RAZÃO DE VEROSSIMILHANÇA (0.07978)
anova(fit.lwt0, fit.lwt1, test="Chi")

# A DIFERENÇA DAS DEVIANCES (5.057) E AS DUAS DEVIANCES (234.67 e 229.62)
anova(fit.lwt1, test="Chi")

#VALOR DO p DO TESTE QUI-QUADRADO DA RAZÃO  DE VEROSSIMILHANÇA - Dois g.l. (três categorias)
pchisq(5.057, 2, lower.tail = FALSE)

# Teste de Wald - 131-150
# Encontrando o valor da z 
-0.6395 / 0.4705
# Encontrando o p-valor do teste de Wald que aparece no modelo simples
2*pnorm(1.359192, lower.tail=FALSE)  # MULTIPLICA POR DOIS PQ E A TABELA Z. SE FOSSE A TABELA T TAMBEM MULTIPLICA.

# Teste de Wald - 151 e +
# Encontrando o valor da z 
-0.8713 / 0.4608
# Encontrando o p-valor do teste de Wald que aparece no modelo simples
2*pnorm (1.890842, lower.tail=FALSE) # MULTIPLICA POR DOIS PQ E A TABELA Z. SE FOSSE A TABELA T TAMBEM MULTIPLICA.

# CALCULO DO MODELO
# Probabilidade doente e expostos - 131-150
py1x1 <- 1 / (1+exp(0.5500+0.6395))
py1x1
0.2333484
# Probabilidade doente e expostos - 151 e + 
py1x1 <- 1 / (1+exp(0.5500+0.8713))
py1x1
0.1944579
# Probabilidade doente e não exposto 
py1x0 <- 1 / (1+exp(0.5500))
py1x0
0.3658644
# Risco relativo - 131-150
0.2333484 / 0.3658644
# Risco relativo - 151 e +
0.1944579 / 0.3658644
# OR - 131-150
exp(-0.6395)
# OR - 151 e +
exp(-0.8713)

#--------------------------------------------------------------------------------------------------------
# LOW E PTL
#--------------------------------------------------------------------------------------------------------
table(dados$ptl)
xtabs(~ptl + flow, data=dados)

# CATEGORIZANDO 0; e 1 e +
dados$ptl1 <- recode(dados$ptl, "0=0; 1=1; 2=1; 3=1") 
dados$fptl1 <- factor(dados$ptl1, levels=c(0,1), labels=c("Nenhum", "1 ou mais"))
table(dados$fptl1)
View(dados)

# Criando uma tabela de contingencia
#--------------------------------------------------------------------------------------------------------
# TABELA CATEGORIZADA
table(dados$fptl1)
tabela <- xtabs(~ fptl1 + flow, data=dados)
tabela

(18/(18+12)) # RISCO NENHUM (REFERÊNCIA) 
(41/(41+118)) # RISCO 1 e + 
(18/(18+12)) / (41/(41+118)) # RR 

(12/(12+118)) # CHANCE NENHUM (REFERÊNCIA)
(18/(18+41)) # CHANCE 1 e + 
(18*118) / (12*41) # OR

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

# Teste de associação pelo qui-quadrado
#--------------------------------------------------------------------------------------------------------
chisq.test(tabela, correct=FALSE) #QUI-QUADRADO SEM CORREÇÃO
chisq.test(tabela, correct=TRUE) #QUI-QUADRADO COM CORREÇÃO

# OR e RR
#--------------------------------------------------------------------------------------------------------
# Cálculo do Risco Relativo
result <- epitab(tabela, method = "riskratio", conf.level = 0.95)
print(result)

# Cálculo da Odds Ratio
result <- epitab(tabela, method = "oddsratio", conf.level = 0.95)
print(result)
# Modelo de regressao logística
#--------------------------------------------------------------------------------------------------------
# Só com o intercepto
fit.ptl0 <- glm(flow ~ 1, family=binomial(logit), data=dados)
summary(fit.ptl0)

# Modelo simples
fit.ptl1 <- glm(flow ~ fptl1, family=binomial(logit), data=dados)
summary(fit.ptl1)
confint(fit.ptl1) # IC para Beta

exp(coef(fit.ptl1)) # OR = Exponencial do Beta1
exp(confint(fit.ptl1)) # IC para OR = Exponencial do IC para Beta1

# CÁLCULO DO LOGARÍTIMO DA FUNÇÃO DE VEROSSIMILHANÇA DO MODELO SÓ COM BETA 0:
59*log(59) + 130*log(130) - 189*log(189)
-117.336 * (-2)
-221.90/ 2
234.672-221.90 # OBTENHO O VALOR DA DEVIANCE

# AS DUAS DEVIANCES (234.67 e 221.90), A DIFERENÇA (12.772) E O VALOR DO p DO QUI-QUADRADO DA RAZÃO  DE VEROSSIMILHANÇA (0.0003514)
anova(fit.ptl0, fit.ptl1, test="Chi")

# A DIFERENÇA DAS DEVIANCES (12.772) E AS DUAS DEVIANCES (234.67 e 221.90)
anova(fit.ptl1, test="Chi")

# VALOR DO p DO TESTE QUI-QUADRADO DA RAZÃO DE VEROSSIMILHANÇA - Um g.l. (duas categorias)
pchisq(12.772, 1, lower.tail = FALSE)

# Teste de Wald
# Encontrando o valor da z 
1.4626 / 0.4144
# Encontrando o p-valor do teste de Wald que aparece no modelo simples
2*pnorm(3.52944, lower.tail=FALSE) 

# CALCULO DO MODELO
# Probabilidade doente e expostos
py1x1 <- 1 / (1+exp(1.0571-1.4626))
py1x1
0.6000084
# Probabilidade doente e não exposto 
py1x0 <- 1 / (1+exp(1.0571))
py1x0
0.257864
# Risco relativo 
0.6000084 / 0.257864
# OR 
exp(1.4626)

#--------------------------------------------------------------------------------------------------------
# LOW E FTV
#--------------------------------------------------------------------------------------------------------
table(dados$ftv)
xtabs(~ftv + flow, data=dados)

# CATEGORIZANDO 0; 1 a 3; e 4 e +
dados$ftv1 <- recode(dados$ftv, "0=0; 1=1; 2=1; 3=1; 4=2; 5=2; 6=2") 
dados$fftv1 <- factor(dados$ftv1, levels=c(0,1,2), labels=c("Nenhuma", "1 a 3", "4 e +"))
table(dados$fftv1)
View(dados)

# Criando uma tabela de contingencia
#--------------------------------------------------------------------------------------------------------
# TABELA DE CONTINGÊNCIA CATEGORIZADA
table(dados$fftv1)
tabela <- xtabs(~fftv1 + flow, data=dados)
tabela 

(36/(36+64)) # RISCO NENHUMA (REFERÊNCIA)
(22/(22+62)) # RISCO 1 a 3
(1/(1+4)) # RISCO 4 e MAIS
(22/(22+62)) / (36/(36+64)) # RR 
(1/(1+4)) / (36/(36+64)) # RR 

(22*64) / (62*36) # OR 
(1*64) / (4*36) # OR 

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

# Teste associação pelo qui-quadrado
#--------------------------------------------------------------------------------------------------------
chisq.test(tabela, correct=FALSE) #QUI-QUADRADO SEM CORREÇÃO
chisq.test(tabela, correct=TRUE) #QUI-QUADRADO COM CORREÇÃO

# OR e RR
#--------------------------------------------------------------------------------------------------------
# Cálculo do Risco Relativo
result <- epitab(tabela, method = "riskratio", conf.level = 0.95)
print(result)

# Cálculo da Odds Ratio
result <- epitab(tabela, method = "oddsratio", conf.level = 0.95)
print(result)

# Modelo de regressao logistica
#--------------------------------------------------------------------------------------------------------
# Só com o intercepto
fit.ftv0 <- glm(flow ~ 1, family=binomial(logit), data=dados)
summary(fit.ftv0)

# Modelo simples
fit.ftv1 <- glm(flow ~ fftv1, family=binomial(logit), data=dados)
summary(fit.ftv1)
confint(fit.ftv1) # IC para Beta

exp(coef(fit.ftv1)) # OR = Exponencial do Beta1
exp(confint(fit.ftv1)) # IC para OR = Exponencial do IC para Beta1

# CÁLCULO DO LOGARÍTIMO DA FUNÇÃO DE VEROSSIMILHANÇA DO MODELO SÓ COM BETA 0:
59*log(59) + 130*log(130) - 189*log(189)
-117.336 * (-2)
-232.29 / 2
234.672-232.29 # OBTENHO O VALOR DA DEVIANCE

# AS DUAS DEVIANCES (234.67 e 232.29), A DIFERENÇA (2.3776) E O VALOR DO p DO QUI-QUADRADO DA RAZÃO DE VEROSSIMILHANÇA (0.3046)
anova(fit.ftv0, fit.ftv1, test="Chi")

# A DIFERENÇA DAS DEVIANCES (2.3776) E AS DUAS DEVIANCES (234.67 e 232.29)
anova(fit.ftv1, test="Chi")

# VALOR DO p DO TESTE QUI-QUADRADO DA RAZÃO DE VEROSSIMILHANÇA - Dois g.l. (três categorias)
pchisq(2.3776, 2, lower.tail = FALSE)

# Teste de Wald – categoria 1 a 3
# Encontrando o valor da z 
-0.4607 / 0.3240
# Encontrando o p-valor do teste de Wald que aparece no modelo simples
2*pnorm(1.421914, lower.tail=FALSE) 

# Teste de Wald – categoria 4 e +
# Encontrando o valor da z 
-0.8109 / 1.1373
# Encontrando o p-valor do teste de Wald que aparece no modelo simples
2*pnorm (0.7130045, lower.tail=FALSE)

# CALCULO DO MODELO
# Probabilidade doente e expostos - 1 a 3
py1x1 <- 1 / (1+exp(0.5754+0.4607))
py1x1
0.2619032
# Probabilidade doente e expostos - 4 ou +
py1x1 <- 1 / (1+exp(0.5754+0.8109))
py1x1
0.1999991
# Probabilidade doente e não exposto 
py1x0 <- 1 / (1+exp(0.5754))
py1x0
0.3599917
# Risco relativo - 1 a 3 
0.2619032 / 0.3599917
# Risco relativo - 4 e +
0.1999991 / 0.3599917
# OR - 1 a 3
exp(-0.4607)
# OR - 4 e +
exp(-0.8109)
#--------------------------------------------------------------------------------------------------------
# MODELO MÚLTIPLO
#--------------------------------------------------------------------------------------------------------

mod.multip.1 <- glm(flow ~ fptl1, family=binomial(logit), data=dados)
summary(mod.multip.1)
exp(coef(mod.multip.1))

mod.multip.2 <- update(mod.multip.1, .~. + fui)
summary(mod.multip.2)
exp(coef(mod.multip.2))

anova (mod.multip.1, mod.multip.2, test="Chi")

mod.multip.3 <- update(mod.multip.2, .~. + fsmoke)
summary(mod.multip.3)
exp(coef(mod.multip.3))

anova (mod.multip.2, mod.multip.3, test="Chi")

mod.multip.4 <- update(mod.multip.3, .~. + fht)
summary(mod.multip.4)
exp(coef(mod.multip.4))

anova (mod.multip.3, mod.multip.4, test="Chi")

mod.multip.5 <- update(mod.multip.4, .~. + frace)
summary(mod.multip.5)
exp(coef(mod.multip.5))

anova (mod.multip.4, mod.multip.5, test="Chi") 

mod.multip.6 <- update(mod.multip.5, .~. + flwt1)
summary(mod.multip.6)
exp(coef(mod.multip.6))

anova (mod.multip.5, mod.multip.6, test="Chi")

mod.multip.7 <- update(mod.multip.6, .~. + fftv1)
summary(mod.multip.7)
exp(coef(mod.multip.7))

anova (mod.multip.6, mod.multip.7, test="Chi")

mod.multip.8 <- update(mod.multip.6, .~. + flwt1)
summary(mod.multip.8)
exp(coef(mod.multip.8))
exp(confint(mod.multip.8))

anova (mod.multip.7, mod.multip.8, test="Chi")

#--------------------------------------------------------------------------------------------------------
# Teste de aderência Hosmer-Lemeshow Goodness of fit (GOF)
#--------------------------------------------------------------------------------------------------------
library (ResourceSelection)
hoslem.test(mod.multip.5$y, fitted(mod.multip.5))

teste.hoslem5 <- hoslem.test(mod.multip.5$y, fitted(mod.multip.5))
teste.hoslem5$observed
teste.hoslem5$expected