Actividad 1 - Módulo 3
Lisset Hernández A01284611
2024-02-29
Instrucciones
Seleccionar una de las dos opciones de bases de datos i) automobile_insurance_claims o ii) health_insurance. A partir de dicha selección realizar las instrucciones 1 – 5. En el desarrollo del archivo de R-Markdown, por favor incluir data storytelling de los resultados del análisis exploratorio de los datos (EDA) así como la interpretación de los resultados estimados.
datos <- read.csv("/Users/lishdz/Desktop/6to/R/Módulo 3/health_insurance.csv")library(DataExplorer)
library(corrplot)## corrplot 0.92 loadedlibrary(caret)## Loading required package: ggplot2## Loading required package: latticelibrary(car) #vif## Loading required package: carDatalibrary(lmtest) #bptest## Loading required package: zoo##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numericlibrary(dlookr) #plot_normality## Registered S3 methods overwritten by 'dlookr':
## method from
## plot.transform scales
## print.transform scales##
## Attaching package: 'dlookr'## The following object is masked from 'package:base':
##
## transformlibrary(spatialreg) #SAR## Loading required package: spData## To access larger datasets in this package, install the spDataLarge
## package with: `install.packages('spDataLarge',
## repos='https://nowosad.github.io/drat/', type='source')`## Loading required package: Matrix## Loading required package: sf## Linking to GEOS 3.11.0, GDAL 3.5.3, PROJ 9.1.0; sf_use_s2() is TRUElibrary(sp)
library(spdep) #moran test##
## Attaching package: 'spdep'## The following objects are masked from 'package:spatialreg':
##
## get.ClusterOption, get.coresOption, get.mcOption,
## get.VerboseOption, get.ZeroPolicyOption, set.ClusterOption,
## set.coresOption, set.mcOption, set.VerboseOption,
## set.ZeroPolicyOptionlibrary(xgboost)
library(Metrics) #rmse() function for xgboost##
## Attaching package: 'Metrics'## The following objects are masked from 'package:caret':
##
## precision, recalllibrary(rpart)
library(rpart.plot)
library(e1071)##
## Attaching package: 'e1071'## The following objects are masked from 'package:dlookr':
##
## kurtosis, skewnesslibrary(randomForest)## randomForest 4.7-1.1## Type rfNews() to see new features/changes/bug fixes.##
## Attaching package: 'randomForest'## The following object is masked from 'package:ggplot2':
##
## marginlibrary(neuralnet)
library(MASS) 1) Contexto
- Brevemente responder con tus propias palabras 2 de las siguientes 3 preguntas:
i) ¿Qué es Supervised Machine Learning y cuáles son algunas
de sus aplicaciones en Inteligencia de Negocios?
Es la rama de la inteligencia artificial que utiliza datas previamente
etiquetada para posteriormente poder hacer predicciones basada en el
aprendizaje que haya hecho. La meta es entrenar modelos con los datos
etiquetados y que después sean capazes de clasificar (o predecir)
correctamente datos nunca antes vistos. En el ámbito de negocios se
puede utilizar para predecir el comportamiento de clientes, por ejemplo
predecir si serán buenos candidatos para recibir un crédito o predecir
el precio de casas.
iii) ¿Qué es la R2 Ajustada? ¿Qué es la métrica RMSE? ¿Cuál
es la diferencia entre la R2 Ajustada y la métrica RMSE?
La R2 Ajustada sirve para explicar el % de variación en la variable
dependiente que es causada por las variables independientes usadas en el
modelo. La métrica RMSE — En general muestra la diferencia promedio
entre los datos reales y los valores predichos por un modelo. Por lo
tanto, entre más pequeño, mejores son las predicciones hechas por el
modelo. Entonces, la diferencia es que ambas métricas es su propósito.
El Adj R2 te dice qué tanto de la variable independiente te explica el
modelo, y el RMSE qué tan confiable es esa explicación.
2) Análisis Exploratorio de los Datos (EDA):
a. Identificación de NA’s
nas <- colSums(is.na(datos))
nas## age sex bmi children smoker region expenses
## 0 0 0 0 0 0 0b. Reemplazo de NA’s
No hay NAs – por lo tanto no se hace tratamiento.
c. Medidas descriptivas
datos$sex <- as.factor(datos$sex)
datos$smoker <- as.factor(datos$smoker)
datos$region <- as.factor(datos$region)
summary(datos)## age sex bmi children smoker
## Min. :18.00 female:662 Min. :16.00 Min. :0.000 no :1064
## 1st Qu.:27.00 male :676 1st Qu.:26.30 1st Qu.:0.000 yes: 274
## Median :39.00 Median :30.40 Median :1.000
## Mean :39.21 Mean :30.67 Mean :1.095
## 3rd Qu.:51.00 3rd Qu.:34.70 3rd Qu.:2.000
## Max. :64.00 Max. :53.10 Max. :5.000
## region expenses
## northeast:324 Min. : 1122
## northwest:325 1st Qu.: 4740
## southeast:364 Median : 9382
## southwest:325 Mean :13270
## 3rd Qu.:16640
## Max. :63770plot(datos$age, datos$expenses, main = "Expenses by age", xlab = "age", ylab = "expenses", col = "orange", pch = 16)
Parece haber tres grupos de expenses. La edad aparente tener poca influencia en expenses ya que, para los tres aparentes grupos de expenses hay de todas las edades en todos los rangos. Aunque conforme va aumentando la edad, el mínimo de expenses va aumantando.
plot(datos$region, datos$expenses, main = "Expenses por Region", xlab = "region", ylab = "expenses", col = "orange", pch = 16)
La región no paraece tener mucha influencia. La mediana de expenses para todas las regiones son cercanas, southeast muestra el mayor rango de costo.
b. Medidas de dispersión
plot_histogram(datos)
## Eliminar outliers
mean <- mean(datos$expenses)
sd <- sd(datos$expenses)
z_scores <- (datos$expenses - mean) / sd
threshold <- 3
outliers_index <- which(abs(z_scores) > threshold)
datos <- datos[-outliers_index, ]plot_histogram(datos$expenses)
plot_correlation(datos)
La mátriz de correlación muestra indicios de que la variable que más influye a expenses (variable dependiente) es si fuman o no.
plot_boxplot(datos, by = "smoker")
Al analizar por smoker, la unica diferencia significativa en cuanto a dispersión es en expenses. La mediana de expense de un no fumador es alrededor de 10,000 y sube hasta 30,000 cuando son fumadores.
plot_boxplot(datos, by = "expenses")
plot_normality(datos, age, sex, bmi, children, smoker, region, expenses)
BMI tiene distribución normal.
plot_histogram(datos$expenses)
3) Hipótesis
A partir de los resultados de EDA describir la especificación del modelo de regresión lineal a estimar. Brevemente, describir cómo es el posible impacto de cada una de las variables explicativas sobre la principal variable de estudio.
En al EDA, infiero que la variable más influyente será SMOKER. Porque tiene una fuerte correlación con expenses. Age en menor medida pero también.
4) Estimación de los modelos
#Crear la partición de datos
set.seed(123)
partition <- createDataPartition(y = datos$expenses, p=0.7, list=F)
train = datos[partition, ]
test = datos[-partition, ]1. OLS
#Crear el modelo
set.seed(123)
ols_model <- lm(expenses ~ age + sex + bmi + children + smoker + region, data = train)
#Calcular el RMSE del modelo
RMSE_ols_model_train <- sqrt(mean(ols_model$residuals^2))
RMSE_ols_model_train #5755.15## [1] 5755.16summary(ols_model)##
## Call:
## lm(formula = expenses ~ age + sex + bmi + children + smoker +
## region, data = train)
##
## Residuals:
## Min 1Q Median 3Q Max
## -10328.7 -2679.0 -962.7 1178.0 23722.4
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -11212.38 1108.30 -10.117 <2e-16 ***
## age 246.42 13.62 18.090 <2e-16 ***
## sexmale -68.04 380.67 -0.179 0.8582
## bmi 320.72 32.07 10.002 <2e-16 ***
## children 495.44 156.48 3.166 0.0016 **
## smokeryes 22869.14 470.05 48.653 <2e-16 ***
## regionnorthwest -69.93 554.54 -0.126 0.8997
## regionsoutheast -783.25 553.30 -1.416 0.1572
## regionsouthwest -1090.41 553.92 -1.969 0.0493 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5783 on 926 degrees of freedom
## Multiple R-squared: 0.7579, Adjusted R-squared: 0.7558
## F-statistic: 362.3 on 8 and 926 DF, p-value: < 2.2e-16Primero, entrené el modelo de OLS usando el subset de train. Según el output, el modelo tiene un p-value muy bajo, lo que indica que el modelo si es estadisticamente significativo. También tiene un Adj. R2 bastante alta. De acuerdo al modelo, las variables más influyentes son la edad, el BMI y si son fumadores. Este modelo muestra un error de 5755.16 en la base de train.
Test
#Crear predicciones
predicciones_ols_model <- predict(ols_model, newdata = test)
summary(predicciones_ols_model)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -1655 6296 9485 12560 14458 39434#Pasar las predicciones para poder utilizarla en calculos
predictions_df <- data.frame(fitted_values = predicciones_ols_model)
#Calcular el RMSE con los fitted values (predicciones) y los actual values (de la base de test)
RMSE_ols_model <-sqrt(mean((predictions_df$fitted_values-test$expenses)^2))
RMSE_ols_model## [1] 5976.548#RMSE = 5976.548Después, para comprobar la precisón del modelo en datos no vistos, hice una prueba con el subset de test. Utilizando las predicciones que se crearon con el modelo, el RMSE aumentó. Lo que quiere decir que el modelo es peor para predecir datos no antes visto.
2. LOG OLS
#Crear el modelo -- solo en LOG las variables relevantes
set.seed(123)
log_ols_model <- lm(log(expenses) ~ log(age) + sex + log(bmi) + children + smoker + region, data = train)
summary(log_ols_model)##
## Call:
## lm(formula = log(expenses) ~ log(age) + sex + log(bmi) + children +
## smoker + region, data = train)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.79387 -0.20917 -0.06750 0.07304 2.25338
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.91530 0.25987 11.218 < 2e-16 ***
## log(age) 1.25811 0.03683 34.158 < 2e-16 ***
## sexmale -0.08966 0.02850 -3.147 0.00170 **
## log(bmi) 0.40088 0.07189 5.576 3.22e-08 ***
## children 0.08237 0.01174 7.014 4.45e-12 ***
## smokeryes 1.49941 0.03520 42.602 < 2e-16 ***
## regionnorthwest -0.04865 0.04150 -1.172 0.24145
## regionsoutheast -0.12897 0.04137 -3.117 0.00188 **
## regionsouthwest -0.12714 0.04149 -3.064 0.00224 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4329 on 926 degrees of freedom
## Multiple R-squared: 0.7786, Adjusted R-squared: 0.7767
## F-statistic: 407 on 8 and 926 DF, p-value: < 2.2e-16De inicio, con los LOG practicamente todas las variables se vuelven influyentes al modelo. Pero edad y si son fumadores son las variables más signigicantes (tienen el p-value más chico). El Adj R2 aumenta 0.02 vs OLS original.
#Transformar los valores predichos a su magnitud original para poder compararlos con los valores reales
log_ols_model$fitted.values <- exp(log_ols_model$fitted.values)#Calcular el RMSE del modelo de LOG
RMSE_log_ols_model_train <-sqrt(mean((log_ols_model$fitted.values-train$expenses)^2))
RMSE_log_ols_model_train## [1] 7428.375#RMSE Train = 7428.375Test
#Crear predicciones
predicciones_log_ols_model <- predict(log_ols_model, newdata = test)
summary(predicciones_log_ols_model)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 7.610 8.552 9.048 9.045 9.370 11.055#Pasar las predicciones para poder utilizarla en calculos
predictions_df_log <- data.frame(fitted_values = exp(predicciones_log_ols_model))
#Calcular el RMSE con los fitted values (predicciones) y los actual values (de la base de test)
RMSE_log_ols_model <-sqrt(mean((predictions_df_log$fitted_values-test$expenses)^2))
RMSE_log_ols_model## [1] 7539.91Al probar el modelo de LOG OLS para predecir se puede ver que mejoro un poco su precisión con datos no vistos antes. Sin embargo, sigue teniendo un RMSE más alto que el modelo de OLS sin transformación logaritmica.
Pruebas de Diagnóstico
Prueba de multicolinealidad
vif(ols_model)## GVIF Df GVIF^(1/(2*Df))
## age 1.022151 1 1.011015
## sex 1.012787 1 1.006373
## bmi 1.107041 1 1.052160
## children 1.004820 1 1.002407
## smoker 1.011864 1 1.005914
## region 1.094123 3 1.015105No hay multicolinealidad – dado que ninguna variable tiene un valor de vif mayor a 5.
Prueba de heterocedasticidad
bptest(ols_model)##
## studentized Breusch-Pagan test
##
## data: ols_model
## BP = 89.623, df = 8, p-value = 5.547e-16Debido a que el valor p es considerablemente menor a 0.05, se rechaza
la hipótesis nula de homocedasticidad. Esto indica que existen
heterocedasticidad en los residuos del modelo.
Debido a la presencia de heterocedasticidad se modeló un modelo de
WLS.
Weighted Least Squared
Train
residuals <- residuals(ols_model)
rmse_ols_wls <- sqrt(mean(residuals^2))set.seed(123)
weights <- 1 / (abs(residuals) + rmse_ols_wls)
wls_model <- lm(expenses ~ age + sex + bmi + children + smoker + region, data = train, weights = weights)summary(wls_model)##
## Call:
## lm(formula = expenses ~ age + sex + bmi + children + smoker +
## region, data = train, weights = weights)
##
## Weighted Residuals:
## Min 1Q Median 3Q Max
## -90.302 -18.860 -5.726 14.527 142.262
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -9287.596 826.814 -11.233 < 2e-16 ***
## age 251.048 9.879 25.411 < 2e-16 ***
## sexmale -125.544 277.388 -0.453 0.6509
## bmi 233.156 24.973 9.336 < 2e-16 ***
## children 486.635 113.556 4.285 2.01e-05 ***
## smokeryes 23549.271 416.716 56.512 < 2e-16 ***
## regionnorthwest -208.297 401.706 -0.519 0.6042
## regionsoutheast -742.206 408.612 -1.816 0.0696 .
## regionsouthwest -937.453 399.706 -2.345 0.0192 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 45.57 on 926 degrees of freedom
## Multiple R-squared: 0.8163, Adjusted R-squared: 0.8147
## F-statistic: 514.4 on 8 and 926 DF, p-value: < 2.2e-16Test
set.seed(123)
wls_predict <- predict(wls_model, newdata = test)
RMSE_wls_model <-sqrt(mean((wls_predict-test$expenses)^2))
RMSE_wls_model## [1] 6031.981bptest(wls_model)##
## studentized Breusch-Pagan test
##
## data: wls_model
## BP = 0.014444, df = 8, p-value = 1Utilizando el método de WLS, dandole un peso proporcional al error que tiene cada observación se logra eliminar la tendencia de heterocedasticidad.
Autocorrelación serial y Espacial
Al no tener un componente de tiempo o espacio (concreto) en los datos, no hay riesgo de autocorrelación serial o espacial
Normalidad de los residuales
hist(wls_model$residuals, main = "Histogram of Residuals", xlab = "Residuals")
plot(ols_model$residuals, xlab= "Dependent Variable", ylab = "Residuals", main = 'OLS Regression Residuals')
abline(0,0)
3. SAR & 4. SEM
Los modelos de SAR (Spatial Autoregressive Model) consideran tanto el comportamiento histórico de una variable como el comportamiento de esa variable en ubicaciones vecinas para realizar predicciones. Sin embargo, con esta base de datos no se tiene suficiente información geografica para hacer un modelo de SAR.
5. XGBoost Regression
data_xgb <- datos %>% dplyr::select(expenses, age, sex, bmi, children, smoker, region)
data_xgb$sex <- ifelse(data_xgb$sex == "female", 1, 0)
# female = 1
# male = 0
data_xgb$smoker <- ifelse(data_xgb$smoker == "yes", 1, 0)
# smoker = 1
# no smoker = 0
data_xgb$region <- as.integer(factor(data_xgb$region, levels = c("northwest", "southeast", "southwest", "northeast"))) - 1
#northwest = 0
#southeast = 1
#soutwest = 2
#northeast = 3
data_xgb$age <- log(data_xgb$age)
data_xgb$bmi <- log(data_xgb$bmi)set.seed(123)
cv_data <- createDataPartition(y = data_xgb$expenses, p=0.7, list=F)
cv_train = data_xgb[cv_data, ]
cv_test = data_xgb[-cv_data, ]
# define explanatory variables (X's) and dependent variable (Y) in training set
train_x = data.matrix(cv_train[, -1]) #seleccionar todas las columnas menos la primera pq es la dependiente
train_y = cv_train[,1] #seleccionar solamente la primera pq es la dependiente
# define explanatory variables (X's) and dependent variable (Y) in testing set
test_x = data.matrix(cv_test[, -1])
test_y = cv_test[, 1]
# define final training and testing sets
xgb_train = xgb.DMatrix(data = train_x, label = train_y)
xgb_test = xgb.DMatrix(data = test_x, label = test_y)# Lets fit XGBoost regression model and display RMSE for both training and testing data at each round
watchlist = list(train=xgb_train, test=xgb_test)
model_xgb = xgb.train(data=xgb_train, max.depth=3, watchlist=watchlist, nrounds=70)## [1] train-rmse:12698.948008 test-rmse:12830.938891
## [2] train-rmse:9453.026194 test-rmse:9647.038046
## [3] train-rmse:7308.136508 test-rmse:7614.969459
## [4] train-rmse:5939.287404 test-rmse:6322.438367
## [5] train-rmse:5106.060187 test-rmse:5580.661249
## [6] train-rmse:4620.218506 test-rmse:5156.688296
## [7] train-rmse:4341.595139 test-rmse:4928.135124
## [8] train-rmse:4182.786050 test-rmse:4794.995494
## [9] train-rmse:4090.964227 test-rmse:4717.928616
## [10] train-rmse:4034.004734 test-rmse:4678.770793
## [11] train-rmse:3987.028983 test-rmse:4661.195458
## [12] train-rmse:3944.847082 test-rmse:4641.082752
## [13] train-rmse:3915.882903 test-rmse:4628.702088
## [14] train-rmse:3891.668431 test-rmse:4625.272617
## [15] train-rmse:3859.956647 test-rmse:4597.826432
## [16] train-rmse:3833.503030 test-rmse:4603.486348
## [17] train-rmse:3819.170199 test-rmse:4604.548695
## [18] train-rmse:3795.564395 test-rmse:4611.952679
## [19] train-rmse:3770.129916 test-rmse:4620.257176
## [20] train-rmse:3744.657857 test-rmse:4625.036175
## [21] train-rmse:3726.065892 test-rmse:4605.296091
## [22] train-rmse:3695.416217 test-rmse:4618.539794
## [23] train-rmse:3689.422434 test-rmse:4621.970822
## [24] train-rmse:3680.883096 test-rmse:4621.031111
## [25] train-rmse:3656.329281 test-rmse:4647.701420
## [26] train-rmse:3630.738857 test-rmse:4655.319650
## [27] train-rmse:3614.602771 test-rmse:4656.852225
## [28] train-rmse:3597.552535 test-rmse:4658.813095
## [29] train-rmse:3584.101633 test-rmse:4654.933999
## [30] train-rmse:3569.062711 test-rmse:4656.475085
## [31] train-rmse:3566.337150 test-rmse:4657.611154
## [32] train-rmse:3554.301959 test-rmse:4661.074501
## [33] train-rmse:3545.396841 test-rmse:4663.433794
## [34] train-rmse:3540.343443 test-rmse:4663.819863
## [35] train-rmse:3515.629475 test-rmse:4671.462321
## [36] train-rmse:3488.713019 test-rmse:4677.821777
## [37] train-rmse:3484.605703 test-rmse:4676.320020
## [38] train-rmse:3465.582213 test-rmse:4675.660260
## [39] train-rmse:3446.532350 test-rmse:4689.529080
## [40] train-rmse:3437.266563 test-rmse:4696.554977
## [41] train-rmse:3425.082849 test-rmse:4698.109859
## [42] train-rmse:3406.535735 test-rmse:4700.732617
## [43] train-rmse:3393.302257 test-rmse:4708.861979
## [44] train-rmse:3372.675808 test-rmse:4711.385226
## [45] train-rmse:3369.191495 test-rmse:4714.231994
## [46] train-rmse:3354.253598 test-rmse:4722.935440
## [47] train-rmse:3339.538087 test-rmse:4729.513663
## [48] train-rmse:3332.932077 test-rmse:4735.218592
## [49] train-rmse:3324.163932 test-rmse:4739.821060
## [50] train-rmse:3321.223821 test-rmse:4740.794787
## [51] train-rmse:3320.326252 test-rmse:4741.110612
## [52] train-rmse:3315.595921 test-rmse:4740.458508
## [53] train-rmse:3304.431218 test-rmse:4752.471478
## [54] train-rmse:3294.444166 test-rmse:4759.653225
## [55] train-rmse:3283.052539 test-rmse:4768.196457
## [56] train-rmse:3278.770235 test-rmse:4773.326782
## [57] train-rmse:3271.511622 test-rmse:4777.815396
## [58] train-rmse:3261.429279 test-rmse:4776.975397
## [59] train-rmse:3254.980231 test-rmse:4777.501234
## [60] train-rmse:3254.361950 test-rmse:4777.887334
## [61] train-rmse:3242.165437 test-rmse:4782.527513
## [62] train-rmse:3238.118449 test-rmse:4782.875334
## [63] train-rmse:3227.372039 test-rmse:4789.700039
## [64] train-rmse:3214.749719 test-rmse:4797.915797
## [65] train-rmse:3200.652279 test-rmse:4801.888020
## [66] train-rmse:3187.230930 test-rmse:4801.220100
## [67] train-rmse:3175.199375 test-rmse:4812.412603
## [68] train-rmse:3151.365418 test-rmse:4815.973371
## [69] train-rmse:3137.222644 test-rmse:4825.309412
## [70] train-rmse:3132.148084 test-rmse:4826.759749Basado en los resultado parece que el RMSE más bajo en ambos sets ocurre en la línea 15.
reg_xgb = xgboost(data = xgb_train, max.depth = 3, nrounds = 15, verbose = 0)
prediction_xgb_test<-predict(reg_xgb, xgb_test)
RMSE_XGB <- rmse(prediction_xgb_test, cv_test$expenses)
RMSE_XGB## [1] 4597.826# Plot first 3 trees of model
xgb.plot.tree(model=reg_xgb, trees=0:2)importance_matrix <- xgb.importance(model = reg_xgb)
xgb.plot.importance(importance_matrix, xlab = "Explanatory Variables X's Importance")
xgb_reg_residuals<-cv_test$expenses - prediction_xgb_test
plot(xgb_reg_residuals, xlab= "Dependent Variable", ylab = "Residuals", main = 'XGBoost Regression Residuals')
abline(0,0)
El modelo de XGBoost obtuvo un RMSE de 4597.826, hasta ahora el más bajo. Debido a que no es tan fácil su interpretacación se generó el árbol de decisión y el gráfico de barras para identificar que las variables que más influyen en este modelo son: smoker, age y bmi.
Decision Tree
set.seed(123)
#Crear el modelo con la base de train
decision_tree_model <- rpart(log(expenses) ~ log(age) + log(bmi) + region + sex + children + smoker, data = train)
# summary(decision_tree_regression)
plot(decision_tree_model, compress = TRUE)
text(decision_tree_model, use.n = TRUE)
rpart.plot(decision_tree_model)
### RMSE of DECISION TREE REGRESSION
decision_tree_prediction <- predict(decision_tree_model,test)
RMSE_decision_tree <- rmse(exp(decision_tree_prediction), test$expenses)
RMSE_decision_tree## [1] 5160.639El árbol de decisión es un modelo sumamente fácil de interpretar, en este caso también obtuvo un RMSE bastante bajo con los otros modelos. Nuevamente, smoker, edad y bmi son los factores principales.
Random Forest
rf_model <- randomForest(expenses ~ smoker + age + bmi + children + sex + region, data= train, proximity=TRUE)
print(rf_model) ### the train data set model accuracy is around 85%.##
## Call:
## randomForest(formula = expenses ~ smoker + age + bmi + children + sex + region, data = train, proximity = TRUE)
## Type of random forest: regression
## Number of trees: 500
## No. of variables tried at each split: 2
##
## Mean of squared residuals: 20116171
## % Var explained: 85.3 #Mean of squared residuals = MSE = benchmark
#% of Var explained = R-sq de los Random Forrest
# Prediction & Confusion Matrix – test data
rf_prediction <- predict(rf_model,test)
RMSE_rf <- rmse(rf_prediction, test$expenses)
RMSE_rf## [1] 4915.729# Evalute Variables' Importance
varImpPlot(rf_model, n.var = 5, main = "Top Variables")
importance(rf_model) ## IncNodePurity
## smoker 76278382456
## age 17455600981
## bmi 17645637950
## children 2753263528
## sex 1213982605
## region 2787283454El modelo de Random Forest combina muchos árboles de decisión para mejorar la predicción. En este caso el forest es capaz de explicar el 85% de la variabilid de la variable dependiente, lo cual es bastante alto. El RMSE también es bueno, es menor al del árbol de decisión individual.
Red Neuronal
library(dplyr)##
## Attaching package: 'dplyr'## The following object is masked from 'package:MASS':
##
## select## The following object is masked from 'package:neuralnet':
##
## compute## The following object is masked from 'package:randomForest':
##
## combine## The following object is masked from 'package:xgboost':
##
## slice## The following object is masked from 'package:car':
##
## recode## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionnn_datos <- dplyr::select(datos, expenses, age, sex, bmi, children, smoker, region)
nn_datos$sex <- ifelse(data_xgb$sex == "female", 1, 0)
nn_datos$smoker <- ifelse(nn_datos$smoker == "yes", 1, 0)
# smoker = 1
# no smoker = 0
nn_datos$region <- as.integer(factor(nn_datos$region, levels = c("northwest", "southeast", "southwest", "northeast"))) - 1
#northwest = 0
#southeast = 1
#soutwest = 2
#northeast = 3#Crear la partición de datos
set.seed(123)
partition <- createDataPartition(y = nn_datos$expenses, p=0.7, list=F)
nntrain = nn_datos[partition, ]
nntest = nn_datos[-partition, ]#Crear el modelo
nn_model <- neuralnet(expenses ~ age + bmi + smoker + region + sex + children, data = nntrain, hidden = c(5, 3), linear.output = TRUE)
# Plot the neural network
plot(nn_model)nn_test <- nntest %>% select(age, bmi, smoker, sex, region, children)
nnpred <- predict(nn_model, nn_test)
RMSE_nn <- rmse(nnpred, nntest$expenses)
RMSE_nn## [1] 11692.29El modelo de redes neuronales es bastante poderosos para crear modelos de regressión, sin embargo en este caso se obtuvo un RMSE bastante alto. Su interpretabilidad es baja y no es fácil identificar cuales son las variables que más influyen el modelo.
5) Evaluación y Selección de Modelo de Regresión
df_rmse <- data.frame( OLS = RMSE_ols_model, LOG = RMSE_log_ols_model, XGBoost = RMSE_XGB, DT = RMSE_decision_tree, RF = RMSE_rf, NN = RMSE_nn, WLS = RMSE_wls_model)
df_rmse## OLS LOG XGBoost DT RF NN WLS
## 1 5976.548 7539.91 4597.826 5160.639 4915.729 11692.29 6031.981df_rmse <- df_rmse[, order(colMeans(df_rmse))]
barplot(as.matrix(df_rmse), beside = TRUE, horiz = TRUE, xlim = c(0, 12000), xlab = "RMSE")
Con base unicamente en los resultados del RMSE, el modelo que
seleccionaría sería XGBoost ya que fue el más chico. Sin embargo, si
tomara en cuenta la interpretabilidad, seleccionaría Random Forest
porque se puede identificar las “decisiones” que toma el modelo para
llegar la las predicciones.
6) Hallazgos.
Desarrollar una breve descripción de los 6 – 10 principales hallazgos de: a. EDA:
- La edad aparente tener poca influencia en expenses ya que, para los
tres aparentes grupos de expenses hay de todas las edades en todos los
rangos. Aunque conforme va aumentando la edad, el mínimo de expenses va
aumantando.
- La región no paraece tener mucha influencia. La mediana de expenses
para todas las regiones son cercanas, southeast muestra el mayor rango
de costo.
- La mátriz de correlación muestra indicios de que la variable que más influye a expenses (variable dependiente) es si fuman o no. Al analizar por smoker, la unica diferencia significativa en cuanto a dispersión es en expenses. La mediana de expense de un no fumador es alrededor de 10,000 y sube hasta 30,000 cuando son fumadores.
- La cantidad de observaciones de los que si fuman a los que no fuman es muy despropocional. Esto influye en ocaciones a la exactitud del modelo.
- La variable de expenses tiene una distribución positiva, el fuerte
de los datos se concentran en la izquierda, menor a $15,000.
- Originalmente existian outliers en la variable de expenses que fueron eliminados bajo la métodología del z-score.
b. Modelo seleccionado: El modelo seleccionado fue XGBoost ya que tienen el RMSE más bajo de todos los modelos probados.
i. ¿Cuáles son las variables que contribuyen a explicar los cambios de la principal variable de estudio?
Las variables que más afectan son: Smoker, Age y BMI.
ii. ¿Cómo es el impacto de dichas variables explicativas sobre la variable dependiente?. Las variables de age, bmi y smoker, tienen un impacto positivo en la variable dependiente, es decir que si estos aumentan, expenses aumenta también.