## corrplot 0.92 loaded
## Loading required package: MASS
## Loading required package: visdat
## Loading required package: car
## Loading required package: carData
## Loading required package: caret
## Loading required package: ggplot2
## Loading required package: lattice
## Loading required package: pROC
## Type 'citation("pROC")' for a citation.
## 
## Attaching package: 'pROC'
## The following objects are masked from 'package:stats':
## 
##     cov, smooth, var
## Loading required package: class
## All packages loaded successfully

Fallas cardiacas

(Proyecto Final)

Se trabajara con un dataset extraido de kaggle, del cual hace referencia a un paper en el que existe

##  [1] "age"                      "anaemia"                 
##  [3] "creatinine_phosphokinase" "diabetes"                
##  [5] "ejection_fraction"        "high_blood_pressure"     
##  [7] "platelets"                "serum_creatinine"        
##  [9] "serum_sodium"             "sex"                     
## [11] "smoking"                  "time"                    
## [13] "DEATH_EVENT"
## 'data.frame':    5000 obs. of  13 variables:
##  $ Edad     : num  55 65 45 60 95 70 63 70 50 53 ...
##  $ Anemia   : Factor w/ 2 levels "No","Yes": 1 1 1 2 2 1 2 2 1 2 ...
##  $ Creatina : int  748 56 582 754 582 232 122 171 482 446 ...
##  $ Diabetes : Factor w/ 2 levels "No","Yes": 1 1 2 2 1 2 2 1 2 1 ...
##  $ EjecFrec : int  45 25 38 40 30 30 60 50 30 45 ...
##  $ PresAlta : Factor w/ 2 levels "No","Yes": 1 1 1 2 1 1 1 2 1 2 ...
##  $ Plaquetas: num  263358 305000 319000 328000 461000 ...
##  $ SerumC   : num  1.3 5 0.9 1.2 2 1.2 1.2 0.9 0.9 1 ...
##  $ SerumS   : int  137 130 140 126 132 132 145 141 132 133 ...
##  $ Genero   : Factor w/ 2 levels "F","M": 2 2 1 2 2 2 1 1 2 2 ...
##  $ Fumador  : Factor w/ 2 levels "No","Yes": 2 1 1 1 1 1 1 1 1 1 ...
##  $ Tiempo   : int  88 207 244 90 50 210 147 196 109 215 ...
##  $ CausaM   : Factor w/ 2 levels "No","Yes": 1 1 1 1 2 1 1 1 1 1 ...
##       Edad       Anemia        Creatina      Diabetes      EjecFrec    
##  Min.   :40.00   No :2628   Min.   :  23.0   No :2803   Min.   :14.00  
##  1st Qu.:50.00   Yes:2372   1st Qu.: 121.0   Yes:2197   1st Qu.:30.00  
##  Median :60.00              Median : 248.0              Median :38.00  
##  Mean   :60.29              Mean   : 586.8              Mean   :37.73  
##  3rd Qu.:68.00              3rd Qu.: 582.0              3rd Qu.:45.00  
##  Max.   :95.00              Max.   :7861.0              Max.   :80.00  
##  PresAlta     Plaquetas          SerumC          SerumS      Genero  
##  No :3176   Min.   : 25100   Min.   :0.500   Min.   :113.0   F:1772  
##  Yes:1824   1st Qu.:215000   1st Qu.:0.900   1st Qu.:134.0   M:3228  
##             Median :263358   Median :1.100   Median :137.0           
##             Mean   :265075   Mean   :1.369   Mean   :136.8           
##             3rd Qu.:310000   3rd Qu.:1.400   3rd Qu.:140.0           
##             Max.   :850000   Max.   :9.400   Max.   :148.0           
##  Fumador        Tiempo      CausaM    
##  No :3441   Min.   :  4.0   No :3432  
##  Yes:1559   1st Qu.: 74.0   Yes:1568  
##             Median :113.0             
##             Mean   :130.7             
##             3rd Qu.:201.0             
##             Max.   :285.0
num_vars <- c('Edad', 'Creatina', 'EjecFrec', 
              'Plaquetas', 'SerumC', 'SerumS', 'Tiempo')

cat_vars <- c('Anemia', 'Diabetes', 'PresAlta', 'Genero', 'Fumador')
attach(data)
color=ifelse(CausaM=="No","black","red")
par(mfrow=c(3,7))

for (i in 1:(length(num_vars) - 1)) {
  for (j in (i + 1):length(num_vars)) {
    plot(data[[num_vars[i]]], data[[num_vars[j]]],
         col = color,
         pch = 20,
         xlab = num_vars[i],
         ylab = num_vars[j])
  }
}

detach(data)
par(mfrow=c(2,3))
lapply(names(data[,-c(1,3,5,7,8,9,12)]),function(y){
  barplot(table(data[,y]), ylab = y)
})

## [[1]]
##      [,1]
## [1,]  0.7
## [2,]  1.9
## 
## [[2]]
##      [,1]
## [1,]  0.7
## [2,]  1.9
## 
## [[3]]
##      [,1]
## [1,]  0.7
## [2,]  1.9
## 
## [[4]]
##      [,1]
## [1,]  0.7
## [2,]  1.9
## 
## [[5]]
##      [,1]
## [1,]  0.7
## [2,]  1.9
## 
## [[6]]
##      [,1]
## [1,]  0.7
## [2,]  1.9
attach(data)
par(mfrow = c(1, 2))
boxplot(EjecFrec ~ CausaM, data = data, main = "Ejection Fraction by DEATH_EVENT",
        xlab = "DEATH_EVENT", ylab = "Ejection Fraction", col = c("blue", "red"))

boxplot(SerumC ~ CausaM, data = data, main = "Serum Creatinine by DEATH_EVENT",
        xlab = "DEATH_EVENT", ylab = "Serum Creatinine", col = c("blue", "red"))

Modelo Con dos Variables a Evaluar

model <- glm(CausaM ~ EjecFrec + SerumC, data = data, family = binomial)

# Display model summary
summary(model)
## 
## Call:
## glm(formula = CausaM ~ EjecFrec + SerumC, family = binomial, 
##     data = data)
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  0.439468   0.140592   3.126  0.00177 ** 
## EjecFrec    -0.066428   0.003488 -19.046  < 2e-16 ***
## SerumC       0.843058   0.048138  17.513  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 6219.5  on 4999  degrees of freedom
## Residual deviance: 5247.8  on 4997  degrees of freedom
## AIC: 5253.8
## 
## Number of Fisher Scoring iterations: 5
Score_RL1 = predict(model)                    # Evalua los Score estimados para el modelo logistico
Prob_RL1  = predict(model,type = "response")


roc <- pROC::roc(data$CausaM,Prob_RL1, auc = TRUE, ci = TRUE)
## Setting levels: control = No, case = Yes
## Setting direction: controls < cases
print(roc)
## 
## Call:
## roc.default(response = data$CausaM, predictor = Prob_RL1, auc = TRUE,     ci = TRUE)
## 
## Data: Prob_RL1 in 3432 controls (data$CausaM No) < 1568 cases (data$CausaM Yes).
## Area under the curve: 0.7718
## 95% CI: 0.7565-0.7871 (DeLong)
pROC::plot.roc(roc, legacy.axes = TRUE, print.thres = "best", print.auc = TRUE,
          auc.polygon = FALSE, max.auc.polygon = FALSE, auc.polygon.col = "gainsboro",
          col = 2, grid = TRUE)

# Generar clasificaciones - Convertir probabilidades predichas en clasificaciones
pc=0.5                              # Seleccionamos un punto de corte pc=0.5
Class_RL1 = as.factor(ifelse(Prob_RL1>pc,"Yes","No")) 
#data.frame(data$CausaM,Class_RL1) Lo estoy comentando para que no saque la tabla de comparacion ya que es muy larga 5000 datos

# Evaluando la bondad de clasificación
table(data$CausaM,Class_RL1,dnn=c("observado","predicho"))          # Matriz de confusión
##          predicho
## observado   No  Yes
##       No  3239  193
##       Yes  977  591
caret::confusionMatrix(Class_RL1,data$CausaM,positive = "Yes")
## Confusion Matrix and Statistics
## 
##           Reference
## Prediction   No  Yes
##        No  3239  977
##        Yes  193  591
##                                          
##                Accuracy : 0.766          
##                  95% CI : (0.754, 0.7777)
##     No Information Rate : 0.6864         
##     P-Value [Acc > NIR] : < 2.2e-16      
##                                          
##                   Kappa : 0.3711         
##                                          
##  Mcnemar's Test P-Value : < 2.2e-16      
##                                          
##             Sensitivity : 0.3769         
##             Specificity : 0.9438         
##          Pos Pred Value : 0.7538         
##          Neg Pred Value : 0.7683         
##              Prevalence : 0.3136         
##          Detection Rate : 0.1182         
##    Detection Prevalence : 0.1568         
##       Balanced Accuracy : 0.6603         
##                                          
##        'Positive' Class : Yes            
## 
ICC_RL1=caret::confusionMatrix(Class_RL1,data$CausaM,positive = "Yes")
ICC_RL1$byClass  
##          Sensitivity          Specificity       Pos Pred Value 
##            0.3769133            0.9437646            0.7538265 
##       Neg Pred Value            Precision               Recall 
##            0.7682638            0.7538265            0.3769133 
##                   F1           Prevalence       Detection Rate 
##            0.5025510            0.3136000            0.1182000 
## Detection Prevalence    Balanced Accuracy 
##            0.1568000            0.6603389

Modelo Saturado y sin variables que no son significativas o que no aportan

modeloSat = glm(formula= CausaM~., data=data,family = "binomial")  # Ajusta un Modelo saturad --> todas las variables
summary(modeloSat)
## 
## Call:
## glm(formula = CausaM ~ ., family = "binomial", data = data)
## 
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  1.132e+01  1.382e+00   8.194 2.53e-16 ***
## Edad         4.651e-02  3.912e-03  11.888  < 2e-16 ***
## AnemiaYes    9.193e-02  8.834e-02   1.041  0.29804    
## Creatina     2.899e-04  4.654e-05   6.229 4.70e-10 ***
## DiabetesYes -4.953e-02  8.634e-02  -0.574  0.56619    
## EjecFrec    -7.734e-02  4.106e-03 -18.837  < 2e-16 ***
## PresAltaYes  2.462e-01  8.768e-02   2.808  0.00498 ** 
## Plaquetas   -5.879e-07  4.440e-07  -1.324  0.18553    
## SerumC       6.566e-01  4.714e-02  13.928  < 2e-16 ***
## SerumS      -8.018e-02  9.804e-03  -8.178 2.87e-16 ***
## GeneroM      3.973e-03  1.002e-01   0.040  0.96838    
## FumadorYes   5.482e-02  1.019e-01   0.538  0.59043    
## Tiempo      -2.042e-02  7.377e-04 -27.683  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 6219.5  on 4999  degrees of freedom
## Residual deviance: 3601.6  on 4987  degrees of freedom
## AIC: 3627.6
## 
## Number of Fisher Scoring iterations: 6
Modelo_RLM1=step(modeloSat,direction = "backward")
## Start:  AIC=3627.55
## CausaM ~ Edad + Anemia + Creatina + Diabetes + EjecFrec + PresAlta + 
##     Plaquetas + SerumC + SerumS + Genero + Fumador + Tiempo
## 
##             Df Deviance    AIC
## - Genero     1   3601.6 3625.6
## - Fumador    1   3601.8 3625.8
## - Diabetes   1   3601.9 3625.9
## - Anemia     1   3602.6 3626.6
## - Plaquetas  1   3603.3 3627.3
## <none>           3601.6 3627.6
## - PresAlta   1   3609.4 3633.4
## - Creatina   1   3647.5 3671.5
## - SerumS     1   3669.7 3693.7
## - Edad       1   3754.1 3778.1
## - SerumC     1   3806.9 3830.9
## - EjecFrec   1   4033.4 4057.4
## - Tiempo     1   4755.6 4779.6
## 
## Step:  AIC=3625.55
## CausaM ~ Edad + Anemia + Creatina + Diabetes + EjecFrec + PresAlta + 
##     Plaquetas + SerumC + SerumS + Fumador + Tiempo
## 
##             Df Deviance    AIC
## - Diabetes   1   3601.9 3623.9
## - Fumador    1   3601.9 3623.9
## - Anemia     1   3602.6 3624.6
## - Plaquetas  1   3603.3 3625.3
## <none>           3601.6 3625.6
## - PresAlta   1   3609.5 3631.5
## - Creatina   1   3648.5 3670.5
## - SerumS     1   3670.1 3692.1
## - Edad       1   3756.0 3778.0
## - SerumC     1   3806.9 3828.9
## - EjecFrec   1   4038.9 4060.9
## - Tiempo     1   4756.8 4778.8
## 
## Step:  AIC=3623.89
## CausaM ~ Edad + Anemia + Creatina + EjecFrec + PresAlta + Plaquetas + 
##     SerumC + SerumS + Fumador + Tiempo
## 
##             Df Deviance    AIC
## - Fumador    1   3602.3 3622.3
## - Anemia     1   3602.9 3622.9
## - Plaquetas  1   3603.9 3623.9
## <none>           3601.9 3623.9
## - PresAlta   1   3610.0 3630.0
## - Creatina   1   3649.1 3669.1
## - SerumS     1   3670.2 3690.2
## - Edad       1   3759.5 3779.5
## - SerumC     1   3807.3 3827.3
## - EjecFrec   1   4039.1 4059.1
## - Tiempo     1   4756.9 4776.9
## 
## Step:  AIC=3622.3
## CausaM ~ Edad + Anemia + Creatina + EjecFrec + PresAlta + Plaquetas + 
##     SerumC + SerumS + Tiempo
## 
##             Df Deviance    AIC
## - Anemia     1   3603.2 3621.2
## - Plaquetas  1   3604.1 3622.1
## <none>           3602.3 3622.3
## - PresAlta   1   3610.2 3628.2
## - Creatina   1   3649.3 3667.3
## - SerumS     1   3671.2 3689.2
## - Edad       1   3761.3 3779.3
## - SerumC     1   3809.1 3827.1
## - EjecFrec   1   4040.1 4058.1
## - Tiempo     1   4758.9 4776.9
## 
## Step:  AIC=3621.16
## CausaM ~ Edad + Creatina + EjecFrec + PresAlta + Plaquetas + 
##     SerumC + SerumS + Tiempo
## 
##             Df Deviance    AIC
## - Plaquetas  1   3605.1 3621.1
## <none>           3603.2 3621.2
## - PresAlta   1   3611.4 3627.4
## - Creatina   1   3649.5 3665.5
## - SerumS     1   3671.3 3687.3
## - Edad       1   3766.5 3782.5
## - SerumC     1   3809.7 3825.7
## - EjecFrec   1   4041.7 4057.7
## - Tiempo     1   4779.2 4795.2
## 
## Step:  AIC=3621.14
## CausaM ~ Edad + Creatina + EjecFrec + PresAlta + SerumC + SerumS + 
##     Tiempo
## 
##            Df Deviance    AIC
## <none>          3605.1 3621.1
## - PresAlta  1   3613.0 3627.0
## - Creatina  1   3650.9 3664.9
## - SerumS    1   3674.3 3688.3
## - Edad      1   3768.9 3782.9
## - SerumC    1   3812.4 3826.4
## - EjecFrec  1   4047.4 4061.4
## - Tiempo    1   4780.1 4794.1
summary(Modelo_RLM1) 
## 
## Call:
## glm(formula = CausaM ~ Edad + Creatina + EjecFrec + PresAlta + 
##     SerumC + SerumS + Tiempo, family = "binomial", data = data)
## 
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  1.112e+01  1.346e+00   8.261  < 2e-16 ***
## Edad         4.711e-02  3.834e-03  12.286  < 2e-16 ***
## Creatina     2.818e-04  4.554e-05   6.188  6.1e-10 ***
## EjecFrec    -7.761e-02  4.077e-03 -19.036  < 2e-16 ***
## PresAltaYes  2.417e-01  8.636e-02   2.799  0.00512 ** 
## SerumC       6.569e-01  4.692e-02  13.999  < 2e-16 ***
## SerumS      -7.967e-02  9.648e-03  -8.257  < 2e-16 ***
## Tiempo      -2.042e-02  7.299e-04 -27.975  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 6219.5  on 4999  degrees of freedom
## Residual deviance: 3605.1  on 4992  degrees of freedom
## AIC: 3621.1
## 
## Number of Fisher Scoring iterations: 6
Modelo_RLM1=update(modeloSat, . ~ . -Genero-Diabetes-Fumador-Anemia-Plaquetas)                 # Elimina los terminos no significativos
summary(Modelo_RLM1)                                              # Verifica la significancia individual en el nuevo modelo
## 
## Call:
## glm(formula = CausaM ~ Edad + Creatina + EjecFrec + PresAlta + 
##     SerumC + SerumS + Tiempo, family = "binomial", data = data)
## 
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  1.112e+01  1.346e+00   8.261  < 2e-16 ***
## Edad         4.711e-02  3.834e-03  12.286  < 2e-16 ***
## Creatina     2.818e-04  4.554e-05   6.188  6.1e-10 ***
## EjecFrec    -7.761e-02  4.077e-03 -19.036  < 2e-16 ***
## PresAltaYes  2.417e-01  8.636e-02   2.799  0.00512 ** 
## SerumC       6.569e-01  4.692e-02  13.999  < 2e-16 ***
## SerumS      -7.967e-02  9.648e-03  -8.257  < 2e-16 ***
## Tiempo      -2.042e-02  7.299e-04 -27.975  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 6219.5  on 4999  degrees of freedom
## Residual deviance: 3605.1  on 4992  degrees of freedom
## AIC: 3621.1
## 
## Number of Fisher Scoring iterations: 6
anova(modeloSat,Modelo_RLM1,test="Chisq")
## Analysis of Deviance Table
## 
## Model 1: CausaM ~ Edad + Anemia + Creatina + Diabetes + EjecFrec + PresAlta + 
##     Plaquetas + SerumC + SerumS + Genero + Fumador + Tiempo
## Model 2: CausaM ~ Edad + Creatina + EjecFrec + PresAlta + SerumC + SerumS + 
##     Tiempo
##   Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1      4987     3601.6                     
## 2      4992     3605.1 -5  -3.5848   0.6106
R2=function(Modelo){
   R2=(Modelo$null.deviance - Modelo$deviance)/Modelo$null.deviance
}
Bondad_Ajuste=data.frame(AIC=model$aic,deviance=model$deviance,R2=R2(model))
row.names(Bondad_Ajuste)="Modelo_2Var"
Bondad_Ajuste
##                  AIC deviance        R2
## Modelo_2Var 5253.775 5247.775 0.1562396
Bondad_Ajuste[2,]=data.frame(AIC=modeloSat$aic,deviance=modeloSat$deviance,R2=R2(modeloSat))
row.names(Bondad_Ajuste)[2]="Modelo_Sat"
Bondad_Ajuste
##                  AIC deviance        R2
## Modelo_2Var 5253.775 5247.775 0.1562396
## Modelo_Sat  3627.552 3601.552 0.4209266
Bondad_Ajuste[3,]=data.frame(AIC=Modelo_RLM1$aic,deviance=Modelo_RLM1$deviance,R2=R2(Modelo_RLM1))
row.names(Bondad_Ajuste)[3]="Modelo_Fin"
Bondad_Ajuste
##                  AIC deviance        R2
## Modelo_2Var 5253.775 5247.775 0.1562396
## Modelo_Sat  3627.552 3601.552 0.4209266
## Modelo_Fin  3621.137 3605.137 0.4203502

Evaluacion del modelo que no tiene las variables que no aportan

Prob_RLM1  = predict(Modelo_RLM1,type = "response")

# Buscando el punto de corte para el modelo múltiple
roc1 <- pROC::roc(data$CausaM,Prob_RLM1, auc = TRUE, ci = TRUE)
## Setting levels: control = No, case = Yes
## Setting direction: controls < cases
print(roc1)
## 
## Call:
## roc.default(response = data$CausaM, predictor = Prob_RLM1, auc = TRUE,     ci = TRUE)
## 
## Data: Prob_RLM1 in 3432 controls (data$CausaM No) < 1568 cases (data$CausaM Yes).
## Area under the curve: 0.901
## 95% CI: 0.8921-0.9099 (DeLong)
pROC::plot.roc(roc1, legacy.axes = TRUE, print.thres = "best", print.auc = TRUE,
          auc.polygon = FALSE, max.auc.polygon = FALSE, auc.polygon.col = "gainsboro",
          col = 2, grid = TRUE) 

# Calidad de la Clasificacion en el modelo múltiple
pc=0.397
Predict_RLM1 = as.factor(ifelse(Prob_RLM1>pc,"Yes","No")) 

# Indicadores de correcta clasificación.
caret::confusionMatrix(Predict_RLM1,data$CausaM,positive = "Yes")
## Confusion Matrix and Statistics
## 
##           Reference
## Prediction   No  Yes
##        No  3023  351
##        Yes  409 1217
##                                           
##                Accuracy : 0.848           
##                  95% CI : (0.8377, 0.8578)
##     No Information Rate : 0.6864          
##     P-Value [Acc > NIR] : < 2e-16         
##                                           
##                   Kappa : 0.6504          
##                                           
##  Mcnemar's Test P-Value : 0.03868         
##                                           
##             Sensitivity : 0.7761          
##             Specificity : 0.8808          
##          Pos Pred Value : 0.7485          
##          Neg Pred Value : 0.8960          
##              Prevalence : 0.3136          
##          Detection Rate : 0.2434          
##    Detection Prevalence : 0.3252          
##       Balanced Accuracy : 0.8285          
##                                           
##        'Positive' Class : Yes             
## 
ICC_RLM1=caret::confusionMatrix(Predict_RLM1,data$CausaM,positive = "Yes")
ICC_RLM1$byClass 
##          Sensitivity          Specificity       Pos Pred Value 
##            0.7761480            0.8808275            0.7484625 
##       Neg Pred Value            Precision               Recall 
##            0.8959692            0.7484625            0.7761480 
##                   F1           Prevalence       Detection Rate 
##            0.7620539            0.3136000            0.2434000 
## Detection Prevalence    Balanced Accuracy 
##            0.3252000            0.8284877