1. Conjunto de datos

Para esta clase se utilizarán datos correspondiente a los pacientes egresados de todos los hospitales pertenecientes al Servicio de Salud Metropolitano Sur (SSMS) durante el año 2017. Estos datos fueron descargados desde el sitio web del Departamento de Estadísticas e Informaciones en Salud (DEIS). La base de datos en formato csv puede ser descargada de este vínculo.

Los datos contienen atributos como el código del establecimiento, su complejidad, el sexo del paciente egresado, la edad, el servicio de egreso, el diagnóstico, condición de egreso, días de estada y si éste fue intervenido o no.

library("knitr")

#Lectura de datos
datos=read.csv("BD_SSMS_2017.csv",sep=";",header=T,check.names = F)

#Organización de datos

  #Entradas
datos$ESTAB=as.factor(datos$ESTAB)                  # Código específico
datos$SEXO=as.factor(datos$SEXO)                    # Catogórica - 1: H, 2: M, 3: OTRO 9: DESCONOCIDO
datos$EDAD=as.numeric(datos$EDAD)                   # Numérica
datos$SERC_EGR=as.factor(datos$SERC_EGR)            # Catogórica - Código específico
datos$INTERV_Q=as.factor(datos$INTERV_Q)            # Catogórica - 1: NO INT 2: INT
datos$DIAG1=as.factor(datos$DIAG1)                  # Catogórica - Código específico
datos$DIAS_ESTAD=as.numeric(datos$DIAS_ESTAD)       # Numérica
datos$C_ESTAB=as.factor(datos$C_ESTAB)              # Complejidad - Alta, MEdia, Baja, Niños, Pediátrica

  #Clase - salida
datos$COND_EGR=as.factor(datos$COND_EGR)            # Catogórica - 1: Viv@ 2: Fallecid@

knitr::kable(head(datos))

2. Balace de datos

En este ejemplo estudiaremos usando una regresión logística, qué variables se relacionan con el fallecer en un hospital del SSMS. El balance de esta clase es la siguiente.

COND=table(datos$COND_EGR)
COND=round(rbind(COND,COND/sum(COND)),2)
colnames(COND)=c("Vivo/Viva","Fallecido/Fallecida")
rownames(COND)=c("Fa","Fr")
knitr::kable(COND,align = "r")

3. Revisión preliminar de variables

library("ggplot2")
library("plotly")
library("knitr")
library("gridExtra")


# Condiciónd de egreso según edad
grafico1=ggplot(datos, aes(x = EDAD, fill = COND_EGR)) + geom_histogram(aes(fill = COND_EGR), alpha = 0.4)
grafico1 = grafico1 + theme(legend.position="bottom") + scale_fill_discrete("Cond. de egreso",labels=c("Vivo/Viva", "Fallecido/Fallecida"))

grafico2 = ggplot(datos, aes(x = EDAD, fill = COND_EGR)) + geom_density(aes(fill = COND_EGR), alpha = 0.4)
grafico2 = grafico2 + theme(legend.position="bottom") + scale_fill_discrete("Cond. de egreso",labels=c("Vivo/Viva", "Fallecido/Fallecida"))

grid.arrange(grafico1, grafico2, ncol=2)

# Condiciónd de egreso según días de estada
grafico1=ggplot(datos, aes(x = log(DIAS_ESTAD), fill = COND_EGR)) + geom_histogram(aes(fill = COND_EGR), alpha = 0.4)
grafico1 = grafico1 + theme(legend.position="bottom") + scale_fill_discrete("Cond. de egreso",labels=c("Vivo/Viva", "Fallecido/Fallecida"))

grafico2 = ggplot(datos, aes(x = log(DIAS_ESTAD), fill = COND_EGR)) + geom_density(aes(fill = COND_EGR), alpha = 0.4)
grafico2 = grafico2 + theme(legend.position="bottom") + scale_fill_discrete("Cond. de egreso",labels=c("Vivo/Viva", "Fallecido/Fallecida"))

grid.arrange(grafico1, grafico2, ncol=2)

# Condición de egreso según IQ
grafico1=ggplot(datos, aes(x = INTERV_Q, fill = COND_EGR)) + geom_bar(aes(fill = COND_EGR), alpha = 0.4)
grafico1 = grafico1 + theme(legend.position="bottom") + scale_fill_discrete("Cond. de egreso",labels=c("Vivo/Viva", "Fallecido/Fallecida"))

# Condición de egreso según Complejidad del Hospital
grafico2 = ggplot(datos, aes(x = C_ESTAB, fill = COND_EGR)) + geom_bar(aes(fill = COND_EGR), alpha = 0.4)
grafico2 = grafico2 + theme(legend.position="bottom") + scale_fill_discrete("Cond. de egreso",labels=c("Vivo/Viva", "Fallecido/Fallecida"))
grafico2 = grafico2 + theme(axis.text.x = element_text(angle = -10, vjust = 0.5, hjust=0.4))


grid.arrange(grafico1, grafico2, ncol=2)

4. Regresión logística

Si bien no hemos visto aún regresión logística implementaremos un modelo de estas características para estudiar las variables relevantes en los pacientes fallecidos. Primero dividiremos el conjunto en un componente de entrenamiento (70%) y otra de prueba (30%), para luego configurar sobre el primero una valicación cruzada de n-fold.

library("caret")

set.seed(1985)
datos_idx = createDataPartition(datos$INTERV_Q, p = 0.7, list = FALSE)
datos_trn = datos[datos_idx,]
datos_tst = datos[-datos_idx,]

fitControl <- trainControl(## 10-fold CV
  method = "cv",
  number = 10,
  savePredictions = TRUE
)


modelo=train(COND_EGR~SEXO+EDAD+DIAS_ESTAD+ESTAB+INTERV_Q+C_ESTAB,data=datos_trn,method="glm",family=binomial(),trControl=fitControl)
print(modelo)

## Generalized Linear Model 
## 
## 45371 samples
##     6 predictor
##     2 classes: '1', '2' 
## 
## No pre-processing
## Resampling: Cross-Validated (10 fold) 
## Summary of sample sizes: 40833, 40835, 40833, 40834, 40835, 40834, ... 
## Resampling results:
## 
##   Accuracy   Kappa      
##   0.9634789  0.001030003

Los coeficientes de la regresión son los siguientes:

summary(modelo)

## 
## Call:
## NULL
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.3296  -0.2675  -0.1361  -0.0752   3.8575  
## 
## Coefficients: (1 not defined because of singularities)
##                                                Estimate Std. Error z value
## (Intercept)                                   -7.068714   0.128108 -55.178
## SEXO2                                         -0.324808   0.052664  -6.168
## SEXO3                                         -7.101184 882.743435  -0.008
## EDAD                                           0.051143   0.001481  34.532
## DIAS_ESTAD                                     0.008502   0.001281   6.635
## ESTAB113130                                   -0.506956   0.332305  -1.526
## ESTAB113150                                   -0.333961   0.097707  -3.418
## ESTAB113160                                   -0.458344   0.137746  -3.327
## ESTAB113170                                  -38.810118 151.327456  -0.256
## ESTAB113180                                   -0.036073   0.067022  -0.538
## ESTAB113190                                   -0.246125   0.101074  -2.435
## INTERV_Q2                                      1.463927   0.081143  18.041
## `C_ESTABEstablecimiento Mediana Complejidad`         NA         NA      NA
##                                              Pr(>|z|)    
## (Intercept)                                   < 2e-16 ***
## SEXO2                                        6.93e-10 ***
## SEXO3                                        0.993582    
## EDAD                                          < 2e-16 ***
## DIAS_ESTAD                                   3.24e-11 ***
## ESTAB113130                                  0.127116    
## ESTAB113150                                  0.000631 ***
## ESTAB113160                                  0.000876 ***
## ESTAB113170                                  0.797592    
## ESTAB113180                                  0.590422    
## ESTAB113190                                  0.014888 *  
## INTERV_Q2                                     < 2e-16 ***
## `C_ESTABEstablecimiento Mediana Complejidad`       NA    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 14209  on 45370  degrees of freedom
## Residual deviance: 11263  on 45359  degrees of freedom
## AIC: 11287
## 
## Number of Fisher Scoring iterations: 13

Mientras que la importancia de las variables se puede resumir como:

varImp(modelo)

## glm variable importance
## 
##              Overall
## EDAD        100.0000
## INTERV_Q2    52.2343
## DIAS_ESTAD   19.1965
## SEXO2        17.8413
## ESTAB113150   9.8770
## ESTAB113160   9.6148
## ESTAB113190   7.0300
## ESTAB113130   4.3956
## ESTAB113180   1.5357
## ESTAB113170   0.7196
## SEXO3         0.0000

mod_pred = predict(modelo,datos_tst)

## Warning in predict.lm(object, newdata, se.fit, scale = 1, type = if (type == :
## prediction from a rank-deficient fit may be misleading

confusionMatrix(mod_pred,datos_tst$COND_EGR)

## Confusion Matrix and Statistics
## 
##           Reference
## Prediction     1     2
##          1 18758   685
##          2     0     1
##                                           
##                Accuracy : 0.9648          
##                  95% CI : (0.9621, 0.9673)
##     No Information Rate : 0.9647          
##     P-Value [Acc > NIR] : 0.4946          
##                                           
##                   Kappa : 0.0028          
##                                           
##  Mcnemar's Test P-Value : <2e-16          
##                                           
##             Sensitivity : 1.000000        
##             Specificity : 0.001458        
##          Pos Pred Value : 0.964769        
##          Neg Pred Value : 1.000000        
##              Prevalence : 0.964719        
##          Detection Rate : 0.964719        
##    Detection Prevalence : 0.999949        
##       Balanced Accuracy : 0.500729        
##                                           
##        'Positive' Class : 1               
## 

Pregunta
¿Qué ocurrió con la predicción? ¿Qué se puede hacer en este caso?