#-------------------------------------------------------------------------------------------------------- # EXERCÍCIO 4 - REGRESSÃO LOGÍSTICA # BANCO STA #-------------------------------------------------------------------------------------------------------- # Altere o diretório de trabalho: setwd("") # Importe o banco de dados "STA" 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 #-------------------------------------------------------------------------------------------------------- dados$fsta <- factor(dados$sta, levels=c(0,1), labels=c("vivo", "morto")) table(dados$fsta) dados$fsex <- factor(dados$sex, levels=c(0,1), labels=c("homem", "mulher")) table(dados$fsex) dados$fser <- factor(dados$ser, levels=c(0,1), labels=c("médico", "cirúrgico")) table(dados$fser) dados$fcan <- factor(dados$can, levels=c(0,1), labels=c("não", "sim")) table(dados$fcan) dados$fcrn <- factor(dados$crn, levels=c(0,1), labels=c("não", "sim")) table(dados$fcrn) dados$finf <- factor(dados$inf, levels=c(0,1), labels=c("não", "sim")) table(dados$finf) dados$fcpr <- factor(dados$cpr, levels=c(0,1), labels=c("não", "sim")) table(dados$fcpr) dados$ftyp <- factor(dados$typ, levels=c(0,1), labels=c("normal", "emergência")) table(dados$ftyp) dados$ffra <- factor(dados$fra, levels=c(0,1), labels=c("não", "sim")) table(dados$ffra) dados$fpo2 <- factor(dados$po2, levels=c(0,1), labels=c(">60","<=60")) table(dados$fpo2) dados$fph <- factor(dados$ph, levels=c(0,1), labels=c(">=7,25","<7,25")) table(dados$fph) dados$fpco <- factor(dados$pco, levels=c(0,1), labels=c("<=45",">45")) table(dados$fpco) dados$fbic <- factor(dados$bic, levels=c(0,1), labels=c(">=18","<18")) table(dados$fbic) dados$fcre <- factor(dados$cre, levels=c(0,1), labels=c("<=2",">2")) table(dados$fcre) dados$fpre <- factor(dados$pre, levels=c(0,1), labels=c("não", "sim")) table(dados$fpre) dados$fage1 <- as.numeric(as.character(dados$age)) dados$fage1 <- recode(dados$fage1, "16:59=0; 60:92=1") dados$fage1 <- factor(dados$fage1, levels=c(0,1), labels=c("16-59","60-92")) table(dados$fage1) dados$floc <- as.numeric(as.character(dados$loc)) dados$floc <- factor(dados$floc, levels=c(0,1,2), labels=c("ñ coma ou inconsicência","coma profundo","coma")) table(dados$floc) dados$floc1 <- as.numeric(as.character(dados$loc)) dados$floc1 <- recode(dados$floc1, "0=0; 1:2=1") dados$floc1 <- factor(dados$floc1, levels=c(0,1), labels=c("ñ coma ou inconsicência","coma")) table(dados$floc1) dados$fsys <- as.numeric(as.character(dados$sys)) dados$fsys <- recode(dados$fsys, "0:120=0; 121:900=1") dados$fsys <- factor(dados$fsys, levels=c(0,1), labels=c("<=120 mmhg",">120 mmhg")) table(dados$fsys) dados$frace <- factor(dados$race, levels=c(1,2,3), labels=c("branco", "preto", "outras")) table(dados$frace) #-------------------------------------------------------------------------------------------------------- # AGE #-------------------------------------------------------------------------------------------------------- tabelaage1 <- xtabs(~ fage1 + fsta, data=dados) tabelaage1 rowPercents(tabelaage1) # percentuais na linha colPercents(tabelaage1) # percentuais na coluna # Teste de associação pelo qui-quadrado #-------------------------------------------------------------------------------------------------------- chisq.test(tabelaage1, correct=FALSE) #melhor esse, pois pega mais var # OR e RR #-------------------------------------------------------------------------------------------------------- # Cálculo do Risco Relativo result <- epitab(tabelaage1, method = "riskratio", conf.level = 0.95) print(result) # Cálculo da Odds Ratio result <- epitab(tabelaage1, method = "oddsratio", conf.level = 0.95) print(result) #-------------------------------------------------------------------------------------------------------- # SEX #-------------------------------------------------------------------------------------------------------- tabelasex <- xtabs(~ fsex + fsta, data=dados) tabelasex rowPercents(tabelasex) # percentuais na linha colPercents(tabelasex) # percentuais na coluna # Teste de associação pelo qui-quadrado #-------------------------------------------------------------------------------------------------------- chisq.test(tabelasex, correct=FALSE) #sem correção chisq.test(tabelasex, correct=TRUE) #com correção # OR e RR #-------------------------------------------------------------------------------------------------------- # Cálculo do Risco Relativo result <- epitab(tabelasex, method = "riskratio", conf.level = 0.95) print(result) # Cálculo da Odds Ratio result <- epitab(tabelasex, method = "oddsratio", conf.level = 0.95) print(result) #-------------------------------------------------------------------------------------------------------- # RACE #-------------------------------------------------------------------------------------------------------- tabelarace <- xtabs(~ frace + fsta, data=dados) tabelarace rowPercents(tabelarace) # percentuais na linha colPercents(tabelarace) # percentuais na coluna # Teste de associação pelo qui-quadrado #-------------------------------------------------------------------------------------------------------- chisq.test(tabelarace, correct=FALSE) #sem correção chisq.test(tabelarace, correct=TRUE) #com correção # OR e RR #-------------------------------------------------------------------------------------------------------- # Cálculo do Risco Relativo result <- epitab(tabelarace, method = "riskratio", conf.level = 0.95) print(result) # Cálculo da Odds Ratio result <- epitab(tabelarace, method = "oddsratio", conf.level = 0.95) print(result) #-------------------------------------------------------------------------------------------------------- # SER #-------------------------------------------------------------------------------------------------------- tabelaser <- xtabs(~ fser + fsta, data=dados) tabelaser rowPercents(tabelaser) # percentuais na linha colPercents(tabelaser) # percentuais na coluna # Teste de associação pelo qui-quadrado #-------------------------------------------------------------------------------------------------------- chisq.test(tabelaser, correct=FALSE) #sem correção chisq.test(tabelaser, correct=TRUE) #com correção # OR e RR com IC 95% #-------------------------------------------------------------------------------------------------------- # Cálculo do Risco Relativo result <- epitab(tabelaser, method = "riskratio", conf.level = 0.95) print(result) # Cálculo da Odds Ratio result <- epitab(tabelaser, method = "oddsratio", conf.level = 0.95) print(result) #-------------------------------------------------------------------------------------------------------- # CAN #-------------------------------------------------------------------------------------------------------- tabelacan <- xtabs(~ fcan + fsta, data=dados) tabelacan rowPercents(tabelacan) # percentuais na linha colPercents(tabelacan) # percentuais na coluna # Teste de associação pelo qui-quadrado #-------------------------------------------------------------------------------------------------------- chisq.test(tabelacan, correct=FALSE) #sem correção chisq.test(tabelacan, correct=TRUE) #com correção # OR e RR #-------------------------------------------------------------------------------------------------------- # Cálculo do Risco Relativo result <- epitab(tabelacan, method = "riskratio", conf.level = 0.95) print(result) # Cálculo da Odds Ratio result <- epitab(tabelacan, method = "oddsratio", conf.level = 0.95) print(result) #-------------------------------------------------------------------------------------------------------- # CRN #-------------------------------------------------------------------------------------------------------- tabelacrn <- xtabs(~ fcrn + fsta, data=dados) #Histórico de falha renal cronica tabelacrn rowPercents(tabelacrn) # percentuais na linha colPercents(tabelacrn) # percentuais na coluna # Teste de associação pelo qui-quadrado #-------------------------------------------------------------------------------------------------------- chisq.test(tabelacrn, correct=FALSE) #sem correção chisq.test(tabelacrn, correct=TRUE) #com correção # OR e RR #-------------------------------------------------------------------------------------------------------- # Cálculo do Risco Relativo result <- epitab(tabelacrn, method = "riskratio", conf.level = 0.95) print(result) # Cálculo da Odds Ratio result <- epitab(tabelacrn, method = "oddsratio", conf.level = 0.95) print(result) #-------------------------------------------------------------------------------------------------------- # INF #-------------------------------------------------------------------------------------------------------- tabelainf <- xtabs(~ finf + fsta, data=dados) tabelainf rowPercents(tabelainf) # percentuais na linha colPercents(tabelainf) # percentuais na coluna # Teste de associação pelo qui-quadrado #-------------------------------------------------------------------------------------------------------- chisq.test(tabelainf, correct=FALSE) chisq.test(tabelainf, correct=TRUE) # OR e RR #-------------------------------------------------------------------------------------------------------- # Cálculo do Risco Relativo result <- epitab(tabelainf, method = "riskratio", conf.level = 0.95) print(result) # Cálculo da Odds Ratio result <- epitab(tabelainf, method = "oddsratio", conf.level = 0.95) print(result) #-------------------------------------------------------------------------------------------------------- # CPR #-------------------------------------------------------------------------------------------------------- tabelacpr <- xtabs(~ fcpr + fsta, data=dados) tabelacpr rowPercents(tabelacpr) # percentuais na linha colPercents(tabelacpr) # percentuais na coluna # Teste de associação pelo qui-quadrado #-------------------------------------------------------------------------------------------------------- chisq.test(tabelacpr, correct=FALSE) #sem correção chisq.test(tabelacpr, correct=TRUE) #com correção # OR e RR #-------------------------------------------------------------------------------------------------------- # Cálculo do Risco Relativo result <- epitab(tabelacpr, method = "riskratio", conf.level = 0.95) print(result) # Cálculo da Odds Ratio result <- epitab(tabelacpr, method = "oddsratio", conf.level = 0.95) print(result) #-------------------------------------------------------------------------------------------------------- # SYS - PAS #-------------------------------------------------------------------------------------------------------- tabelapas <- xtabs(~ fsys + fsta, data=dados) tabelapas rowPercents(tabelapas) # percentuais na linha colPercents(tabelapas) # percentuais na coluna # Teste de associação pelo qui-quadrado #-------------------------------------------------------------------------------------------------------- chisq.test(tabelapas, correct=FALSE) #sem correção chisq.test(tabelapas, correct=TRUE) #com correção # OR e RR #-------------------------------------------------------------------------------------------------------- # Cálculo do Risco Relativo result <- epitab(tabelapas, method = "riskratio", conf.level = 0.95) print(result) # Cálculo da Odds Ratio result <- epitab(tabelapas, method = "oddsratio", conf.level = 0.95) print(result) #-------------------------------------------------------------------------------------------------------- # PRE #-------------------------------------------------------------------------------------------------------- tabelapre <- xtabs(~ fpre + fsta, data=dados) tabelapre rowPercents(tabelapre) # percentuais na linha colPercents(tabelapre) # percentuais na coluna # Teste de associação pelo qui-quadrado #-------------------------------------------------------------------------------------------------------- chisq.test(tabelapre, correct=FALSE) #sem correção chisq.test(tabelapre, correct=TRUE) #com correção # OR e RR #-------------------------------------------------------------------------------------------------------- # Cálculo do Risco Relativo result <- epitab(tabelapre, method = "riskratio", conf.level = 0.95) print(result) # Cálculo da Odds Ratio result <- epitab(tabelapre, method = "oddsratio", conf.level = 0.95) print(result) #-------------------------------------------------------------------------------------------------------- # TYP #-------------------------------------------------------------------------------------------------------- tabelatyp <- xtabs(~ ftyp + fsta, data=dados) tabelatyp rowPercents(tabelatyp) # percentuais na linha colPercents(tabelatyp) # percentuais na coluna # Teste de associação pelo qui-quadrado #-------------------------------------------------------------------------------------------------------- chisq.test(tabelatyp, correct=FALSE) #sem correção chisq.test(tabelatyp, correct=TRUE) #com correção # OR e RR #-------------------------------------------------------------------------------------------------------- # Cálculo do Risco Relativo result <- epitab(tabelatyp, method = "riskratio", conf.level = 0.95) print(result) # Cálculo da Odds Ratio result <- epitab(tabelatyp, method = "oddsratio", conf.level = 0.95) print(result) #-------------------------------------------------------------------------------------------------------- # FRA #-------------------------------------------------------------------------------------------------------- tabelafra <- xtabs(~ ffra + fsta, data=dados) tabelafra rowPercents(tabelafra) # percentuais na linha colPercents(tabelafra) # percentuais na coluna # Teste de associação pelo qui-quadrado #-------------------------------------------------------------------------------------------------------- chisq.test(tabelafra, correct=FALSE) #sem correção chisq.test(tabelafra, correct=TRUE) #com correção # OR e RR #-------------------------------------------------------------------------------------------------------- # Cálculo do Risco Relativo result <- epitab(tabelafra, method = "riskratio", conf.level = 0.95) print(result) # Cálculo da Odds Ratio result <- epitab(tabelafra, method = "oddsratio", conf.level = 0.95) print(result) #-------------------------------------------------------------------------------------------------------- # PO2 #-------------------------------------------------------------------------------------------------------- tabelapo2 <- xtabs(~ fpo2 + fsta, data=dados) tabelapo2 rowPercents(tabelapo2) # percentuais na linha colPercents(tabelapo2) # percentuais na coluna # Teste de associação pelo qui-quadrado #-------------------------------------------------------------------------------------------------------- chisq.test(tabelapo2, correct=FALSE) #sem correção chisq.test(tabelapo2, correct=TRUE) #com correção # OR e RR #-------------------------------------------------------------------------------------------------------- # Cálculo do Risco Relativo result <- epitab(tabelapo2, method = "riskratio", conf.level = 0.95) print(result) # Cálculo da Odds Ratio result <- epitab(tabelapo2, method = "oddsratio", conf.level = 0.95) print(result) #-------------------------------------------------------------------------------------------------------- # PH #-------------------------------------------------------------------------------------------------------- tabelaph <- xtabs(~ fph + fsta, data=dados) tabelaph rowPercents(tabelaph) # percentuais na linha colPercents(tabelaph) # percentuais na coluna # Teste de associação pelo qui-quadrado #-------------------------------------------------------------------------------------------------------- chisq.test(tabelaph, correct=FALSE) #sem correção chisq.test(tabelaph, correct=TRUE) #com correção # OR e RR #-------------------------------------------------------------------------------------------------------- # Cálculo do Risco Relativo result <- epitab(tabelaph, method = "riskratio", conf.level = 0.95) print(result) # Cálculo da Odds Ratio result <- epitab(tabelaph, method = "oddsratio", conf.level = 0.95) print(result) #-------------------------------------------------------------------------------------------------------- # PCO #-------------------------------------------------------------------------------------------------------- tabelapco <- xtabs(~ fpco + fsta, data=dados) tabelapco rowPercents(tabelapco) # percentuais na linha colPercents(tabelapco) # percentuais na coluna # Teste qui-quadrado #-------------------------------------------------------------------------------------------------------- chisq.test(tabelapco, correct=FALSE) #sem correção chisq.test(tabelapco, correct=TRUE) #com correção # OR e RR #-------------------------------------------------------------------------------------------------------- # Cálculo do Risco Relativo result <- epitab(tabelapco, method = "riskratio", conf.level = 0.95) print(result) # Cálculo da Odds Ratio result <- epitab(tabelapco, method = "oddsratio", conf.level = 0.95) print(result) #-------------------------------------------------------------------------------------------------------- # BIC #-------------------------------------------------------------------------------------------------------- tabelabic <- xtabs(~ fbic + fsta, data=dados) tabelabic rowPercents(tabelabic) # percentuais na linha colPercents(tabelabic) # percentuais na coluna # Teste de associação pelo qui-quadrado #-------------------------------------------------------------------------------------------------------- chisq.test(tabelabic, correct=FALSE) #sem correção chisq.test(tabelabic, correct=TRUE) #com correção # OR e RR no R #-------------------------------------------------------------------------------------------------------- # Cálculo do Risco Relativo result <- epitab(tabelabic, method = "riskratio", conf.level = 0.95) print(result) # Cálculo da Odds Ratio result <- epitab(tabelabic, method = "oddsratio", conf.level = 0.95) print(result) #-------------------------------------------------------------------------------------------------------- # CRE #-------------------------------------------------------------------------------------------------------- tabelacre <- xtabs(~ fcre + fsta, data=dados) tabelacre rowPercents(tabelacre) # percentuais na linha colPercents(tabelacre) # percentuais na coluna # Teste de associação pelo qui-quadrado #-------------------------------------------------------------------------------------------------------- chisq.test(tabelacre, correct=FALSE) #sem correção chisq.test(tabelacre, correct=TRUE) #com correção # OR e RR #-------------------------------------------------------------------------------------------------------- # Cálculo do Risco Relativo result <- epitab(tabelacre, method = "riskratio", conf.level = 0.95) print(result) # Cálculo da Odds Ratio result <- epitab(tabelacre, method = "oddsratio", conf.level = 0.95) print(result) #-------------------------------------------------------------------------------------------------------- # LOC #-------------------------------------------------------------------------------------------------------- tabelaloc <- xtabs(~ floc + fsta, data=dados) tabelaloc tabelaloc1 <- xtabs(~ floc1 + fsta, data=dados) tabelaloc1 rowPercents(tabelaloc) # percentuais na linha colPercents(tabelaloc) # percentuais na coluna rowPercents(tabelaloc1) # percentuais na linha colPercents(tabelaloc1) # percentuais na coluna # Teste de associação pelo qui-quadrado #-------------------------------------------------------------------------------------------------------- chisq.test(tabelaloc, correct=FALSE) #sem correção chisq.test(tabelaloc, correct=TRUE) #com correção chisq.test(tabelaloc1, correct=FALSE) #sem correção chisq.test(tabelaloc1, correct=TRUE) #com correção # OR e RR no R #-------------------------------------------------------------------------------------------------------- # Cálculo do Risco Relativo result <- epitab(tabelaloc, method = "riskratio", conf.level = 0.95) print(result) result <- epitab(tabelaloc1, method = "riskratio", conf.level = 0.95) print(result) # Cálculo da Odds Ratio result <- epitab(tabelaloc, method = "oddsratio", conf.level = 0.95) print(result) result <- epitab(tabelaloc1, method = "riskratio", conf.level = 0.95) print(result) #-------------------------------------------------------------------------------------------------------- # MODELO MÚLTIPLO #-------------------------------------------------------------------------------------------------------- mod.1 <- glm(fsta ~ floc1, family=binomial(logit), data=dados) summary(mod.1) exp(confint(mod.1)) exp(coef(mod.1)) mod.2 <- glm(fsta ~ floc1 + ftyp, family=binomial(logit), data=dados) summary(mod.2) anova(mod.1, mod.2, test="Chi") exp(coef(mod.2)) mod.3 <- glm(fsta ~ floc1 + ftyp + fcpr , family=binomial(logit), data=dados) summary(mod.3) anova(mod.2, mod.3, test="Chi") exp(coef(mod.3)) mod.4 <- glm(fsta ~ floc1 + ftyp + fcpr + fser , family=binomial(logit), data=dados) summary(mod.4) anova(mod.3, mod.4, test="Chi") exp(coef(mod.4)) mod.5 <- glm(fsta ~ floc1 + ftyp + fcpr + finf, family=binomial(logit), data=dados) summary(mod.5) anova(mod.3, mod.5, test="Chi") #teste de verossimilhança exp(coef(mod.5)) mod.6 <- glm(fsta ~ floc1 + ftyp + fcpr + fcrn, family=binomial(logit), data=dados) summary(mod.6) anova(mod.3, mod.6, test="Chi") exp(coef(mod.6)) mod.7 <- glm(fsta ~ floc1 + ftyp + fcpr + fcre, family=binomial(logit), data=dados) #family fala qual a função (chance, OR, prob) summary(mod.7) anova(mod.3, mod.7, test="Chi") exp(coef(mod.7)) mod.8 <- glm(fsta ~ floc1 + ftyp + fcpr + fage1, family=binomial(logit), data=dados) summary(mod.8) anova(mod.3, mod.8, test="Chi") exp(coef(mod.8)) mod.9 <- glm(fsta ~ floc1 + ftyp + fcpr + fage1 + fbic , family=binomial(logit), data=dados) summary(mod.9) anova(mod.8, mod.9, test="Chi") exp(coef(mod.8)) #-------------------------------------------------------------------------------------------------------- # Teste de aderência Hosmer-Lemeshow Goodness of fit (GOF) #-------------------------------------------------------------------------------------------------------- library (ResourceSelection) hoslem.test(mod.8$y, fitted(mod.8))