library(readxl)
df<- read_excel("C:/Users/wsand/Downloads/Procesamiento diciembre 2022 10311 registros.xls")
df$HIPOACUSIA <- ifelse(df$HIPOACUSIA == 2, 0, df$HIPOACUSIA)
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
df1=df %>% select(Hipoacusia_multi,SEXO,EDAD, TIPOZONAEXPOSICIÓNARUIDO, MOLESTIARUIDOENACTIVIDADESEXTRAMURALES,
              FUMAOFUMÓ,TRATAMIENTOCONALGUNMEDICAMENTOOTOTÓXICOPORMÁSDEDOSMESES, EXAMENFÍSICOMEMBRANATIMPANICAOD,
               EXAMENFÍSICOMEMBRANATIMPANICAOI, `@2HAPRESENTADONINGUNSÍNTOMA`, SÍNDROMEALTERACIÓNENSALUD)
df1$SEXO <- relevel(as.factor(df1$SEXO), ref = "Mujeres")
df1$TIPOZONAEXPOSICIÓNARUIDO<-relevel(as.factor(df1$TIPOZONAEXPOSICIÓNARUIDO), ref = "2")
df1$MOLESTIARUIDOENACTIVIDADESEXTRAMURALES <- relevel(as.factor(df1$MOLESTIARUIDOENACTIVIDADESEXTRAMURALES), ref="2")
df1$FUMAOFUMÓ <-relevel(as.factor(df1$FUMAOFUMÓ), ref="2")
df1$TRATAMIENTOCONALGUNMEDICAMENTOOTOTÓXICOPORMÁSDEDOSMESES<-relevel(as.factor(df1$TRATAMIENTOCONALGUNMEDICAMENTOOTOTÓXICOPORMÁSDEDOSMESES), ref="2" )
df1$EXAMENFÍSICOMEMBRANATIMPANICAOD<-relevel(as.factor(df1$EXAMENFÍSICOMEMBRANATIMPANICAOD), ref="2")
df1$EXAMENFÍSICOMEMBRANATIMPANICAOI <- relevel(as.factor(df1$EXAMENFÍSICOMEMBRANATIMPANICAOI), ref="2")
df1$`@2HAPRESENTADONINGUNSÍNTOMA`<-relevel(as.factor(df1$`@2HAPRESENTADONINGUNSÍNTOMA`), ref="2")
df1$SÍNDROMEALTERACIÓNENSALUD<-relevel(as.factor(df1$SÍNDROMEALTERACIÓNENSALUD), ref="2")
modelo1 <- glm(Hipoacusia_multi ~SEXO + EDAD+TIPOZONAEXPOSICIÓNARUIDO+ 
                 MOLESTIARUIDOENACTIVIDADESEXTRAMURALES+ 
                 FUMAOFUMÓ+TRATAMIENTOCONALGUNMEDICAMENTOOTOTÓXICOPORMÁSDEDOSMESES+
                 EXAMENFÍSICOMEMBRANATIMPANICAOD+
                EXAMENFÍSICOMEMBRANATIMPANICAOI+
                `@2HAPRESENTADONINGUNSÍNTOMA`+
                 SÍNDROMEALTERACIÓNENSALUD,
                 ,data = df1, family = binomial(logit))
summary(modelo1)
## 
## Call:
## glm(formula = Hipoacusia_multi ~ SEXO + EDAD + TIPOZONAEXPOSICIÓNARUIDO + 
##     MOLESTIARUIDOENACTIVIDADESEXTRAMURALES + FUMAOFUMÓ + TRATAMIENTOCONALGUNMEDICAMENTOOTOTÓXICOPORMÁSDEDOSMESES + 
##     EXAMENFÍSICOMEMBRANATIMPANICAOD + EXAMENFÍSICOMEMBRANATIMPANICAOI + 
##     `@2HAPRESENTADONINGUNSÍNTOMA` + SÍNDROMEALTERACIÓNENSALUD, 
##     family = binomial(logit), data = df1)
## 
## Coefficients:
##                                                           Estimate Std. Error
## (Intercept)                                              -4.556323   0.112410
## SEXOHombres                                               0.583303   0.049822
## EDAD                                                      0.062987   0.001918
## TIPOZONAEXPOSICIÓNARUIDO1                                 0.273610   0.047358
## MOLESTIARUIDOENACTIVIDADESEXTRAMURALES1                   0.327222   0.054593
## FUMAOFUMÓ1                                                0.181239   0.064022
## TRATAMIENTOCONALGUNMEDICAMENTOOTOTÓXICOPORMÁSDEDOSMESES1  0.354863   0.115301
## EXAMENFÍSICOMEMBRANATIMPANICAOD1                          0.405836   0.079481
## EXAMENFÍSICOMEMBRANATIMPANICAOI1                          0.361137   0.081185
## `@2HAPRESENTADONINGUNSÍNTOMA`1                            0.482225   0.046763
## SÍNDROMEALTERACIÓNENSALUD1                                0.264933   0.052331
##                                                          z value Pr(>|z|)    
## (Intercept)                                              -40.533  < 2e-16 ***
## SEXOHombres                                               11.708  < 2e-16 ***
## EDAD                                                      32.845  < 2e-16 ***
## TIPOZONAEXPOSICIÓNARUIDO1                                  5.778 7.58e-09 ***
## MOLESTIARUIDOENACTIVIDADESEXTRAMURALES1                    5.994 2.05e-09 ***
## FUMAOFUMÓ1                                                 2.831  0.00464 ** 
## TRATAMIENTOCONALGUNMEDICAMENTOOTOTÓXICOPORMÁSDEDOSMESES1   3.078  0.00209 ** 
## EXAMENFÍSICOMEMBRANATIMPANICAOD1                           5.106 3.29e-07 ***
## EXAMENFÍSICOMEMBRANATIMPANICAOI1                           4.448 8.66e-06 ***
## `@2HAPRESENTADONINGUNSÍNTOMA`1                            10.312  < 2e-16 ***
## SÍNDROMEALTERACIÓNENSALUD1                                 5.063 4.14e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 13407  on 10310  degrees of freedom
## Residual deviance: 11445  on 10300  degrees of freedom
## AIC: 11467
## 
## Number of Fisher Scoring iterations: 4

#Linealidad Uno de los supuestos más importantes en regresión logística es que la relación entre el logit o log-odds de la variable respuesta y cada variable predictora o variable independiente es lineal; este supuesto se verifica únicamente para las variables numéricas continuas que se tengan en el modelo.

Test de Box Tidwell

Para verificar este supuesto, es posible usar el Test Box-Tidwell, el cual se realiza por medio de la función boxTidwell El test se plantea tomando como hipótesis nula el cumplimineto del supuesto de linealidad.

library(car)
## Loading required package: carData
## 
## Attaching package: 'car'
## The following object is masked from 'package:dplyr':
## 
##     recode
logodds <- modelo1$linear.predictors
boxTidwell(logodds ~ df1$EDAD)
##  MLE of lambda Score Statistic (t) Pr(>|t|)
##         1.0425              1.0777   0.2812
## 
## iterations =  1

La salida anterior presenta el valor p asociado al estadístico de prueba del test de Box Tidwell, a partir del cual, no se rechaza el supuesto de linealidad entre el logit y la variable Edad.

library(ggplot2)
logEdad <- data.frame(logodds, Edad = df1$EDAD)
ggplot(data = logEdad, aes(x = Edad, y = logodds)) + 
  geom_point() +
  ggtitle("Edad vs Logaritmo del Odds") +
  theme_bw() +
  theme(plot.title = element_text(hjust = 0.5))

Independencia entre las observaciones

library(ggplot2)
Indice <- seq(1,10311,1)
Residuales <- modelo1$residuals
residuales <- data.frame(Indice, Residuales)
ggplot(data = residuales, aes(x = Indice, y = Residuales)) +
  geom_point() +
  ggtitle("Residuales del modelo ajustado") +
  theme_bw() +
  theme(plot.title = element_text(hjust = 0.5))

No hay un patrón obvio que indique una falta de independencia, pero hay algunas observaciones que podrían ser outliers o ejercer una gran influencia en el modelo.

Multicolinealidad

Aunque este fenómeno es particularmente relevante en la presencia de múltiples predictores continuos, en nuestro modelo solo se incluyó una variable independiente continua (Edas). Las demás variables en el estudio fueron binarias, reduciendo considerablemente la preocupación por la multicolinealidad. No obstante, se realizaron análisis para asegurar que no existiera colinealidad perfecta entre las variables categóricas binarias. Se inspeccionaron las tablas de contingencia para cada par de variables categóricas, lo que no reveló ninguna separación completa o alta correlación entre las categorías de dichas variables.

Bondad de ajuste basados en devianza

deviance(modelo1)
## [1] 11444.84

La devianza revela que el modelo completo logra una reducción significativa en la devianza de 13407 a 11445 en comparación con el modelo nulo, indicando un mejor ajuste a los datos observados. Esta mejora en el ajuste del modelo es apoyada por una disminución en los grados de libertad de 10310 a 10300, reflejando la inclusión de los predictores en el modelo.

Prueba de Hosmer-Lemeshow:

library(ResourceSelection)
## ResourceSelection 0.3-6   2023-06-27
hl <- hoslem.test(modelo1$y, fitted(modelo1), g=10)
hl
## 
##  Hosmer and Lemeshow goodness of fit (GOF) test
## 
## data:  modelo1$y, fitted(modelo1)
## X-squared = 11.877, df = 8, p-value = 0.1568

Esto da p=0,49, lo que indica que no hay evidencia de un ajuste deficiente