Control 3
Camila Diaz Orrego, William Awad
22-06-2021
Control 3 - Camila Diaz Orrego, William Awad
Predicción ataques cárdiacos Suponga que Ud. trabaja como analista en conjunto con médicos del área cardiovascular. Estos médicos le entregan una base de datos cuyas observaciones corresponde a personas que han tenido o no un ataque al corazón.
P0.
Cargar paquetes.
rm(list=ls())
library(data.table)
library(caret)## Loading required package: lattice## Loading required package: ggplot2library(jtools)
library(rpart.plot)## Loading required package: rpartlibrary(rpart)
library(factoextra)## Welcome! Want to learn more? See two factoextra-related books at https://goo.gl/ve3WBalibrary(dbscan)
library(ggplot2)
library(moderndive)## Registered S3 methods overwritten by 'broom':
## method from
## tidy.glht jtools
## tidy.summary.glht jtoolslibrary(knitr)
datads= fread("DATOS CONTROL 3.csv")
str(data)## function (..., list = character(), package = NULL, lib.loc = NULL, verbose = getOption("verbose"),
## envir = .GlobalEnv, overwrite = TRUE)datads[,sex:=as.factor(sex)]
datads[,cp:=as.factor(cp)]
datads[,fbs:=as.factor(fbs)]
datads[,restecg:=as.factor(restecg)]
datads[,exng:=as.factor(exng)]
datads[,output:=as.factor(output)]
str(data)## function (..., list = character(), package = NULL, lib.loc = NULL, verbose = getOption("verbose"),
## envir = .GlobalEnv, overwrite = TRUE)P1.
Realice dos modelos de regresión lineal multiple para predecir la Presión arterial en reposo ¿Cuál predice mejor dentro de muestra?. (8 puntos)
Modelo 1
reg1=lm(data=datads,formula =trtbps~sex + age +restecg)
summary(reg1)##
## Call:
## lm(formula = trtbps ~ sex + age + restecg, data = datads)
##
## Residuals:
## Min 1Q Median 3Q Max
## -39.322 -10.790 -0.734 9.438 65.019
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 107.0815 6.5219 16.419 < 2e-16 ***
## sex1 -1.1394 2.1019 -0.542 0.5882
## age 0.4982 0.1088 4.577 6.93e-06 ***
## restecg1 -3.6063 1.9727 -1.828 0.0685 .
## restecg2 3.3123 8.6000 0.385 0.7004
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 16.84 on 298 degrees of freedom
## Multiple R-squared: 0.09015, Adjusted R-squared: 0.07794
## F-statistic: 7.382 on 4 and 298 DF, p-value: 1.111e-05datads[,prediccion_2:=predict(reg1)]
datads[,residuo_2:=resid(reg1)]
pred1=predict(reg1)
datads[,pred1:=predict(reg1)]
predicciones1=data.table(RMSE=RMSE(pred1,datads$trtbps),
MAE=MAE(pred1,datads$trtbps))
predicciones1## RMSE MAE
## 1: 16.7013 12.88092Modelo 2
reg2=lm(data=datads,formula =trtbps~sex + age + restecg + output + chol + fbs + exng)
summary(reg2)##
## Call:
## lm(formula = trtbps ~ sex + age + restecg + output + chol + fbs +
## exng, data = datads)
##
## Residuals:
## Min 1Q Median 3Q Max
## -38.837 -11.028 -2.039 10.233 58.750
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 109.27671 8.45151 12.930 < 2e-16 ***
## sex1 -2.17577 2.24692 -0.968 0.333673
## age 0.40900 0.11317 3.614 0.000355 ***
## restecg1 -2.70668 1.99499 -1.357 0.175905
## restecg2 3.83047 8.55697 0.448 0.654740
## output1 -3.16227 2.30042 -1.375 0.170287
## chol 0.01436 0.01959 0.733 0.464214
## fbs1 7.25779 2.72741 2.661 0.008218 **
## exng1 0.04023 2.27386 0.018 0.985895
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 16.68 on 294 degrees of freedom
## Multiple R-squared: 0.1197, Adjusted R-squared: 0.09579
## F-statistic: 4.999 on 8 and 294 DF, p-value: 7.937e-06datads[,prediccion_2:=predict(reg2)]
datads[,residuo_2:=resid(reg2)]
pred2=predict(reg2)
datads[,pred2:=predict(reg2)]
predicciones2=data.table(RMSE=RMSE(pred2,datads$trtbps),
MAE=MAE(pred2,datads$trtbps))
predicciones2## RMSE MAE
## 1: 16.42744 13.0191Podemos ver que en el modelo 2, donde incorporamos más variables el error cuadrático medio es levemente menor que el del modelo 1. Esto quiere decir, que si bien no es muy preciso, es más preciso que el con menor cantidad de variables. Además, respecto al R cuadrado, vemos que en el modelo 2 es mayor en un poco menos de 0.01, pero aún así nos habla mejor del modelo 2.
P2.
Realice validación cruzada (CV) a los modelos de la pregunta anterior por el método K-folds con 5 folds. ¿Se mantienen las conclusiones anteriores?. (8 putos) Pista1: Recuerde setear la semilla set.seed(12345).
Pista2: Si existen variables con NA recuerde que puede excluirlas esas observaciones del análisis, pero no las elimine.
set.seed(12345)
setupKCV = trainControl(method = "cv" , number = 5)
predkfolds1=train(trtbps~sex + age +restecg,data=datads,method="lm",trControl= setupKCV, na.action = na.omit)
predkfolds2=train(trtbps~sex + age + restecg + output + chol + fbs + exng,data=datads,method="lm",trControl= setupKCV, na.action = na.omit)
print(predkfolds1)## Linear Regression
##
## 303 samples
## 3 predictor
##
## No pre-processing
## Resampling: Cross-Validated (5 fold)
## Summary of sample sizes: 242, 242, 242, 244, 242
## Resampling results:
##
## RMSE Rsquared MAE
## 17.12826 0.06884883 13.13225
##
## Tuning parameter 'intercept' was held constant at a value of TRUEprint(predkfolds2)## Linear Regression
##
## 303 samples
## 7 predictor
##
## No pre-processing
## Resampling: Cross-Validated (5 fold)
## Summary of sample sizes: 243, 241, 242, 243, 243
## Resampling results:
##
## RMSE Rsquared MAE
## 16.67958 0.09187585 13.3158
##
## Tuning parameter 'intercept' was held constant at a value of TRUEEn la predicción fuera de muestra, se ve una disminución en el R2 de ambos modelos pero aún así el modelo 2 sigue siendo mejor que el modelo 1. Por las mismas razones que en la predicción dentro de la muestra.
P3.
Proponga dos variables sobre las cuales segmentar la muestra, en orden de hacer análisis de clusters con el método de kmeans. Muestre sus resultados gráficamente. (5 puntos) Pista1: Si existen outliers recuerde eliminarlos, para esto cree un nuevo objeto donde se encuentren estas dos variables.
datads1=datads[!trtbps>180]
datads1=datads1[!age>73]
datads1=datads1[!age<30]
data_age=datads1[,.(age, trtbps)]
ggplot(data = datads1,aes(x=age,y=trtbps))+
geom_point()
k1=kmeans(x=data_age,centers=5 ,nstart=30)
fviz_cluster(k1,data=data_age,geom = "point")
P4.
Realice dos modelos de árboles de clasificación de la variable output. Pruebe cuál modelo clasifica mejor con validación cruzada. Entrene el modelo con un 80% de la muestra y testee con el 20% restante. Explicite qué modelo es mejor y porqué. (12 puntos) Importante: Recuerde setear la semilla set.seed(12345). Observación: No obtendrá puntaje si compara un modelo de clasificación de una variable.
Modelo 1
#Modelo 1: predecir output a través de sex, age, chol, thalachh, restecg y trtbps
arbol1 = rpart(output ~ sex + age + chol + thalachh + restecg + trtbps, data = datads , method = "class")
rpart.plot(arbol1, main = "Tipo de ataque cardiaco")
set.seed(12345)
div = createDataPartition(datads$output, times = 1, p = 0.8, list = F)
traindata = datads[div,]
testdata = datads[-div,]
arbol1.1 = rpart(output~sex + age + chol + thalachh + restecg + trtbps, data=traindata, method = "class")
rpart.plot(arbol1.1, main = "Tipo de ataque cardiaco")
### Modelo 2
#Predicción de output con variables sex, age, chol, thalachh , restecg , trtbps , cp , exng , fbs y oldpeak.
arbol2 = rpart(output ~sex + age + chol + thalachh + restecg + trtbps + cp + exng + fbs + oldpeak, data=datads, method = "class")
rpart.plot(arbol2, main = "Tipo de ataque cardiaco")
set.seed(12345)
div = createDataPartition(datads$output, times = 1, p = 0.8, list = F)
traindata = datads[div,]
testdata = datads[-div,]
arbol2.1 = rpart(output ~sex + age + chol + thalachh + restecg + trtbps + cp + exng + fbs + oldpeak , data=traindata, method = "class")
rpart.plot(arbol2.1, main = "Tipo de ataque cardiaco")
Conclusiones de ambos modelos
#Predicción y su precisión modelo 1
prediccion1 = predict(arbol1.1, newdata = testdata, type = "class")
matriz1 <- table(testdata$output, prediccion1)
matriz1## prediccion1
## 0 1
## 0 17 10
## 1 7 26precision1 <- sum(diag(matriz1))/sum(matriz1)
precision1## [1] 0.7166667#Predicción y su precisión modelo 2
prediccion2 <- predict(arbol2.1, newdata = testdata, type = "class")
matriz2 <- table(testdata$output, prediccion2)
matriz2## prediccion2
## 0 1
## 0 18 9
## 1 8 25precision2 <- sum(diag(matriz2))/sum(matriz2)
precision2## [1] 0.7166667