Actividad 1 – Modelos de Regresión

1. Preguntas

i) ¿Qué es Supervised Machine Learning y cuáles son algunas de sus aplicaciones en Inteligencia de Negocios?

El aprendizaje automático supervisado (SML) es un proceso en el que se entrena un modelo predictivo utilizando datos etiquetados como entradas y buscando predecir o clasificar una salida, una vez comprendido esto podemos utilizarlo en el mundo de los negocios mediante predicciones como lo pueden ser de ventas, análisis de riestos, etc. o en el caso de segmentación con el uso de técnicas de clustering para identificar grupos de clientes con características similares.

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 (R cuadrada) es una medida que nos dice qué tan bien se ajusta un modelo a los datos, es decir que tanto porcentaje de ellos puede explicar, cuanto más alto sea el valor de R2, mejor se ajusta el modelo a los datos y más porcentaje de la variabilidad de la variable dependiente puede explicar La métrica RMSE es la medida que explica la magnitud promedio de los errores entre los valores predichos y los valores reales, es decir que entre menor sea el valor del RMSE, mejor será la capacidad predictiva del modelo, ya que indica que las predicciones del modelo están más cerca de los valores reales La principal diferencia entre estas 2 medidas es que la R2 evalúa la calidad general del modelo, mientras que el RMSE evalúa la precisión de las predicciones.

2. Análisis exploratorio (EDA)

2.1 Instalar paquetes y llamar librerías

library(dplyr)
library(ggplot2)
library(corrplot)
library(car)
library(lmtest)
library(lattice)
library(xgboost)
library(rpart)
library(randomForest)
library(neuralnet)

2.2 Importación base de datos

df<- read.csv("C:\\Users\\LuisD\\Documents\\Concentración\\MODULO 3\\health_insurance.csv")
dfrf<- read.csv("C:\\Users\\LuisD\\Documents\\Concentración\\MODULO 3\\health_insurance.csv")

#file.choose()

2.3 Primer acercamiento

summary(df)

##       age            sex                 bmi           children    
##  Min.   :18.00   Length:1338        Min.   :16.00   Min.   :0.000  
##  1st Qu.:27.00   Class :character   1st Qu.:26.30   1st Qu.:0.000  
##  Median :39.00   Mode  :character   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  
##     smoker             region             expenses    
##  Length:1338        Length:1338        Min.   : 1122  
##  Class :character   Class :character   1st Qu.: 4740  
##  Mode  :character   Mode  :character   Median : 9382  
##                                        Mean   :13270  
##                                        3rd Qu.:16640  
##                                        Max.   :63770

str(df)

## 'data.frame':    1338 obs. of  7 variables:
##  $ age     : int  19 18 28 33 32 31 46 37 37 60 ...
##  $ sex     : chr  "female" "male" "male" "male" ...
##  $ bmi     : num  27.9 33.8 33 22.7 28.9 25.7 33.4 27.7 29.8 25.8 ...
##  $ children: int  0 1 3 0 0 0 1 3 2 0 ...
##  $ smoker  : chr  "yes" "no" "no" "no" ...
##  $ region  : chr  "southwest" "southeast" "southeast" "northwest" ...
##  $ expenses: num  16885 1726 4449 21984 3867 ...

# df

2.4 Cambiar tipo de variables

Es necesario categorizar algunas variables para que pueda realizarse el primer método de SML (OLS)

df$sex <- ifelse(df$sex == "male", 1, 0)
df$smoker <- ifelse(df$smoker == "yes", 1, 0)
df$region <- as.numeric(factor(df$region, levels = c("southwest","southeast","northeast", "northwest"), labels = c(1,2,3,4)))
str(df)

## 'data.frame':    1338 obs. of  7 variables:
##  $ age     : int  19 18 28 33 32 31 46 37 37 60 ...
##  $ sex     : num  0 1 1 1 1 0 0 0 1 0 ...
##  $ bmi     : num  27.9 33.8 33 22.7 28.9 25.7 33.4 27.7 29.8 25.8 ...
##  $ children: int  0 1 3 0 0 0 1 3 2 0 ...
##  $ smoker  : num  1 0 0 0 0 0 0 0 0 0 ...
##  $ region  : num  1 2 2 4 4 2 2 4 3 4 ...
##  $ expenses: num  16885 1726 4449 21984 3867 ...

2.5 NA’s (Identificación y reemplazo)

na_count <- colSums(is.na(df))
na_count

##      age      sex      bmi children   smoker   region expenses 
##        0        0        0        0        0        0        0

#No contamos con NA's por lo que no es necesario remplazar ningun registro con mediana.

2.6 Medidas descriptivas

columnas_numericas <- sapply(df, is.numeric)

media <- colMeans(df[, columnas_numericas], na.rm = TRUE)
cat("Media:", media, "\n")

## Media: 39.20703 0.5052317 30.66547 1.094918 0.2047833 2.485052 13270.42

mediana <- sapply(df[, columnas_numericas], median, na.rm = TRUE) 
cat("Mediana:", mediana, "\n")

## Mediana: 39 1 30.4 1 0 2 9382.03

moda <- apply(df, 2, function(x) {
  unique_x <- unique(x[!is.na(x)])
  unique_x[which.max(tabulate(match(x, unique_x)))]})  
cat("Moda:", media, "\n")

## Moda: 39.20703 0.5052317 30.66547 1.094918 0.2047833 2.485052 13270.42

2.7 Medidas de dispersión

desviacion_estandar <- sapply(df, sd)
desviacion_estandar <- format(desviacion_estandar, scientific = FALSE)
desviacion_estandar

##             age             sex             bmi        children          smoker 
## "   14.0499604" "    0.5001596" "    6.0983822" "    1.2054927" "    0.4036940" 
##          region        expenses 
## "    1.1055720" "12110.0112397"

varianza <- sapply(df, var)
varianza <- format(varianza, scientific = FALSE)
varianza

##                 age                 sex                 bmi            children 
## "      197.4013867" "        0.2501596" "       37.1902653" "        1.4532127" 
##              smoker              region            expenses 
## "        0.1629689" "        1.2222895" "146652372.2258170"

range <- sapply(df, range)
range <- format(range, scientific = FALSE)
range

##      age        sex        bmi        children   smoker     region    
## [1,] "   18.00" "    0.00" "   16.00" "    0.00" "    0.00" "    1.00"
## [2,] "   64.00" "    1.00" "   53.10" "    5.00" "    1.00" "    4.00"
##      expenses  
## [1,] " 1121.87"
## [2,] "63770.43"

cuartiles <- apply(df, 2, quantile, probs = c(0, 0.25, 0.5, 0.75, 1))
rownames(cuartiles) <- c("0%", "25%", "50%", "75%", "100%")
colnames(cuartiles) <- colnames(df)
cuartiles

##      age sex  bmi children smoker region  expenses
## 0%    18   0 16.0        0      0      1  1121.870
## 25%   27   0 26.3        0      0      2  4740.288
## 50%   39   1 30.4        1      0      2  9382.030
## 75%   51   1 34.7        2      0      3 16639.915
## 100%  64   1 53.1        5      1      4 63770.430

2.8 Graficos para identificar patrones

# Distribución de géneros en gráfico de barras:
#barplot(table(df$sex), main="Distribución de género", xlab="Sexo", ylab="Frecuencia",ylim=c(0, 800))

# Distribución de fumadores en gráfico de barras:
barplot(table(df$smoker), main="Distribución de fumadores", xlab="Fumador o no fumador", ylab="Frecuencia",ylim=c(0, 1300))

# Histograma de la variable edad
hist(df$age, main="Histograma de Edad", xlab="Edad", ylab="Frecuencia")

# Histograma de la variable bmi
hist(df$bmi, main="Histograma de BMI", xlab="BMI", ylab="Frecuencia")

# Matriz de correlación para variables numéricas
corr_matrix <- cor(df[, c("age", "bmi", "children", "expenses")])
corrplot(corr_matrix, method="color")

# Box plot para 'expenses' vs. 'smoker'
boxplot(expenses ~ smoker, data=df, main="Box plot de Gastos por Fumador", xlab="Fumador", ylab="Gastos")

# Box plot para 'expenses' vs. 'region'
boxplot(expenses ~ region, data=df, main="Box plot de Gastos por Región", xlab="Región", ylab="Gastos")

#plot_normality(df, age, bmi, children)

3. Interpretación EDA

Tomando en cuenta los gráficos elaborados durante la anterior etapa podemos llegar a las siguientes conclusiones:
* Con respecto a las primeras 2 gráficas relacionadas a la distribución y cantidad de observaciones en las variables sex y smoker podemos ver que existe una caitdad similar de registros de hombres y mujeres, sin embargo dentro de la muestra contamos con un numero significante mayor de personas fumadoras a no fumadoras.
* De la misma forma podemos observar mediante el histograma de age que contamos con una gran variedad de registros de todas las edades, pues la diferencia entre las frecuencias del gráfico no es tan grande.
* En cuanto a la matriz de correlación entre cada variable, podemos observar una correlación mediana entre las viariables age y expenses y una correlación ligera entre las variables bmi y expenses, si tomamos en cuenta que expenses es nuestra variable dependiente puede ser que la correlación que presenta con estas dos variables independientes se vea reflejada mas adelante y las veamos como estadísticamente significantes para explicar el cambio en la variable Y.
* Finalmente mediante el uso de boxplots pudimos observar el comportamiento de la variable dependiente expenses cuando la relacionamos con la variable smoker, en este primer caso podemos observar una mediana bastante mayor dentro de las observaciones smoker-yes y también una mayor dispersión de observaciones, por el contrario del comportamiento presentado entre expenses y region pues no se observa una diferencia significativa entre las medias de las regiones, únicamente podemos observar un ligero aumento en cuanto a la variabilidad de las expenses relacionada a las regiones. * Tendremos que transformar la variable bmi para que no haya problemas de heterocedasticidad.

4. Modelos SML

4.0 Partición de datos y semilla

library(caret)
set.seed(123)
datos_partidos <- createDataPartition(y = df$expenses, p=0.8, list=F)
train = df[datos_partidos, ]
test <- df[-datos_partidos, ]

4.1 OLS Regresión

4.1.1 Creación del modelo

# Creación del modelo
ols_model <- lm(log(expenses)  ~ age + sex + log(bmi) + children + smoker + region, data = train)
summary(ols_model)

## 
## Call:
## lm(formula = log(expenses) ~ age + sex + log(bmi) + children + 
##     smoker + region, data = train)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.05636 -0.19007 -0.04657  0.07205  2.11194 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  5.9917571  0.2404797  24.916  < 2e-16 ***
## age          0.0339603  0.0009826  34.562  < 2e-16 ***
## sex         -0.0917478  0.0273668  -3.353 0.000829 ***
## log(bmi)     0.3891503  0.0690145   5.639 2.19e-08 ***
## children     0.0935763  0.0114271   8.189 7.48e-16 ***
## smoker       1.5683769  0.0337222  46.509  < 2e-16 ***
## region       0.0292215  0.0126078   2.318 0.020653 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.446 on 1065 degrees of freedom
## Multiple R-squared:  0.7669, Adjusted R-squared:  0.7656 
## F-statistic: 583.9 on 6 and 1065 DF,  p-value: < 2.2e-16

4.1.2 Prueba de Multicolinealidad

vif(ols_model)

##      age      sex log(bmi) children   smoker   region 
## 1.020795 1.009071 1.047174 1.003168 1.006477 1.027515

# Los valores de VIF son menores a 5 o 10 por lo que NO existe multicolinealidad

Tomando en cuenta los resultados de la prueba VIF podemos afirmar que NO EXISTE MULTICOLINEALIDAD

4.1.3 Prueba Heterocedasticidad

4.1.3.1 Visualización de residuales

# Como primer contacto hacemos un gráfico para ver los residuales del modelo
res<- residuals(ols_model)
plot(res, type="l")

El gráfico anterior parece ser prueba acerca de la existencia de Heterocedasticidad por lo que procederemos a realizar el BP-Test para confirmarlo

4.1.3.2 BP-Test

bptest(ols_model)

## 
##  studentized Breusch-Pagan test
## 
## data:  ols_model
## BP = 72.185, df = 6, p-value = 1.456e-13

Tomando en cuenta el resultado del Breusch-Pagan Test, observamos que este indica fuertemente la presencia de heterocedasticidad en el modelo OLS, por lo que es necesario realizar una estimación de mínimos cuadrados ponderados (WLS)…

4.1.3.3 WLS

wt <- 1 / lm(abs(ols_model$residuals) ~ ols_model$fitted.values)$fitted.values^2
wls_model<-lm(log(expenses) ~  age + sex + log(bmi) + children + smoker + region , data = train, weights=wt)
summary(wls_model)

## 
## Call:
## lm(formula = log(expenses) ~ age + sex + log(bmi) + children + 
##     smoker + region, data = train, weights = wt)
## 
## Weighted Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.0223 -0.7588 -0.1465  0.3821  6.3784 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  5.7753850  0.2345984  24.618  < 2e-16 ***
## age          0.0308326  0.0009478  32.530  < 2e-16 ***
## sex         -0.0721799  0.0264208  -2.732   0.0064 ** 
## log(bmi)     0.4990262  0.0671977   7.426 2.28e-13 ***
## children     0.0803374  0.0108946   7.374 3.32e-13 ***
## smoker       1.5032775  0.0291502  51.570  < 2e-16 ***
## region       0.0232378  0.0121946   1.906   0.0570 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.57 on 1065 degrees of freedom
## Multiple R-squared:  0.7853, Adjusted R-squared:  0.7841 
## F-statistic: 649.2 on 6 and 1065 DF,  p-value: < 2.2e-16

str(train)

## 'data.frame':    1072 obs. of  7 variables:
##  $ age     : int  19 18 28 33 31 46 37 60 25 62 ...
##  $ sex     : num  0 1 1 1 0 0 0 0 1 0 ...
##  $ bmi     : num  27.9 33.8 33 22.7 25.7 33.4 27.7 25.8 26.2 26.3 ...
##  $ children: int  0 1 3 0 0 1 3 0 0 0 ...
##  $ smoker  : num  1 0 0 0 0 0 0 0 0 1 ...
##  $ region  : num  1 2 2 4 2 2 4 4 3 2 ...
##  $ expenses: num  16885 1726 4449 21984 3757 ...

bptest(wls_model)

## 
##  studentized Breusch-Pagan test
## 
## data:  wls_model
## BP = 13971, df = 6, p-value < 2.2e-16

#  Podemos observar que se ha logrado eliminar la heterocedasticidad aumentando el p-value a 1

vif(wls_model)

##      age      sex log(bmi) children   smoker   region 
## 1.017365 1.012836 1.048003 1.001651 1.008322 1.027750

# Podemos observar también que a pesar de que los coeficientes de VIF aumentaron, NO tenemos multicolinealidad

4.1.4 Normalidad de residuales

# Residuales del WLS
res_wls<- residuals(wls_model)
plot(res_wls, type="l")

# Gráfico Q-Q de residuales
qqnorm(res_wls)
qqline(res_wls)

# Create residual vs. fitted plot

plot(fitted(wls_model), resid(wls_model), main="Linear Regression Residual vs. Fitted Values", xlab="Fitted Values", ylab="Residuals")
abline(0,0)

hist(wls_model$residuals, xlab="Estimated Regression Residuals", main='Distribution of OLS Estimated Regression Residuals', col='lightblue', border="white" )

# Shaphiro Test
shaphiro <- shapiro.test(res_wls)
shaphiro

## 
##  Shapiro-Wilk normality test
## 
## data:  res_wls
## W = 0.88152, p-value < 2.2e-16

A pesar de que aparentemente no contamos con normalidad de residuales, debemos recordar que estamos evaluando un modelo de WLS por lo que podemos continuar a la comparación de RMSE y R2 de los modelos OLS y WLS…

4.1.5 Comparación de desempeño RMSE

# Predicciones del modelo OLS
predicciones_ols <- predict(ols_model, newdata = test)
# Convertir las predicciones de logaritmo a su escala original
predicciones_ols <- exp(predicciones_ols)

# Predicciones del modelo WLS
predicciones_wls <- predict(wls_model, newdata = test)
# Convertir las predicciones de logaritmo a su escala original
predicciones_wls <- exp(predicciones_wls)

# Calcular el RMSE para cada modelo
rmse_ols <- sqrt(mean((predicciones_ols - test$expenses)^2))
rmse_wls <- sqrt(mean((predicciones_wls - test$expenses)^2))

# Imprimir los RMSE
cat("RMSE del modelo OLS:", rmse_ols, "\n")

## RMSE del modelo OLS: 8217.117

cat("RMSE del modelo WLS:", rmse_wls, "\n")

## RMSE del modelo WLS: 7302.724

4.2 XGBoost Regresion

reg_xgb <- xgboost(data = xgb_train, max_depth = 3, nrounds = 59, verbose = 0)
prediction_xgb_test <- predict(reg_xgb, xgb_test)
rmse_XGB <- sqrt(mean((exp(prediction_xgb_test) - exp(test_y))^2))
cat("RMSE del modelo XGBoost:", rmse_ols, "\n")

## RMSE del modelo XGBoost: 8217.117

4.3 Decision Tree Regresion

dt_regresion <- rpart(log(expenses) ~ age + sex + bmi + children + smoker + region, data = train)
plot(dt_regresion, compress = TRUE)
text(dt_regresion, use.n = TRUE)

library(rpart.plot)
rpart.plot(dt_regresion)

decision_tree_prediction <- predict(dt_regresion,test)
rmse_dt <- sqrt(mean((exp(decision_tree_prediction) - df$expenses)^2))
cat("RMSE de Decision Trees:", rmse_dt, "\n")

## RMSE de Decision Trees: 16132.53

4.4 Random Forest Regresion

# Crear el modelo de Random Forest
library(randomForest)
dfrf$sex<- factor(dfrf$sex, levels=c("male", "female"))
dfrf$smoker<- factor(dfrf$smoker, levels=c("yes", "no"))
dfrf$region <- factor(dfrf$region, levels = c("southwest", "southeast", "northeast", "northwest"))

library(caret)
set.seed(123)
datos_partidosrf <- createDataPartition(y = dfrf$expenses, p=0.8, list=F)
trainrf = dfrf[datos_partidos, ]
testrf <- dfrf[-datos_partidos, ]

random_forest <- randomForest(log(expenses) ~ age + sex + bmi + children + smoker + region, data = trainrf, na.action = na.exclude)
random_forest

## 
## Call:
##  randomForest(formula = log(expenses) ~ age + sex + bmi + children +      smoker + region, data = trainrf, na.action = na.exclude) 
##                Type of random forest: regression
##                      Number of trees: 500
## No. of variables tried at each split: 2
## 
##           Mean of squared residuals: 0.1504463
##                     % Var explained: 82.25

predictions_rf <- predict(random_forest, testrf)
rmse_rf <- sqrt(mean((exp(predictions_rf) - df$expenses)^2))
cat("RMSE Random Forest:", rmse_rf, "\n")

## RMSE Random Forest: 15421.87

4.5 Red Neuronal

library(neuralnet)
nn_model <- neuralnet(expenses ~ age + sex + bmi + children + smoker + region, data = train, hidden = c(5, 3), linear.output = TRUE)
plot(nn_model)
predictions_nnet <- predict(nn_model, test)
rmse_nn <- sqrt(mean((exp(predictions_nnet) - test$expenses)^2))
cat("RMSE de Redes Neuronales:", rmse_nn, "\n")

## RMSE de Redes Neuronales: Inf

6. Evaluación y selección de modelos

rmsesy<- c(rmse_ols,rmse_wls,rmse_XGB,rmse_dt,rmse_rf)

rmsesx<- c("OLS", "WLS", "XGBoost", "DT", "RF")

df_RMSES<- data.frame(rmsesx=rmsesx,rmsesy=rmsesy)


barplot(df_RMSES$rmsesy, names= df_RMSES$rmsesx, col = "skyblue", 
        main = "Valores de RMSE para diferentes modelos", 
        xlab = "RMSE", ylab = "Valor de RMSE")

Tomando en cuenta los resultados presentados anteriormente, en donde comparamos los RMSE de los diferentes modelos realizados (RMSE de Neural Networks no fue incluido debido a valor sumamente alto). Poemos observar que el modelo con un mejor desempeño comprobado por ser en RMSE con el valor más pequeño es el XGBoost model # 7. Conclusiones ## a. EDA Durante el Análisis Exploratorio de Datos tuve la oportunidad de comenzar el análisis realizando un acercamiento principal mediante str y summary para saber el tipo de variables y un poco acerca de su comportamiento, también incluimos medidas de dispersión y descriptivas, a la hora de buscar NA’s no logramos identificar la presencia de alguno, finalmente concluimos el EDA mediante la visualización de histogramas para poder visualizar las frecuencias de las variables numéricas, también usamos boxplots para poder observar un poco del comprotamiento de la variable dependiente y su relación con las variables categóricas, fnalmente podemos visualizar un heatmap en el que se muestran las relaciones entre variables.

b. Modelo seleccionado

Como conclusión y basándonos en la comparación de RMSE’s de los distintos modelos, podemos llegar a la conclusión de que el modelo XGBoost es aquel que muestra un mejor desempeño, en caso de tener alguna duda en cuanto a la comparación, se sugiere implementar algunas otras técnicas como AIC. ### i. ¿Cuáles son las variables que contribuyen a explicar los cambios de la principal variable de estudio?
Las principales variables que contribuyen a la explicación de los cambios son:
* age
* bmi
* smoker yes-no

ii. ¿Cómo es el impacto de dichas variables explicativas sobre la variable dependiente?

El impacto de todas las variables anteriores es positivo, es decir que si el valor de estas aumenta, también lo hará el valor de expenses, en el caso de smoker se justifica debido a que 0=“no” y 1=“yes”

iii. ¿Los resultados estimados del modelo seleccionado son similares a los otros modelos estimados? ¿Cuáles son las diferencias?

Si comparamos los rasultados obtenidos de XGBoost con los resultados de los demás modelos podemos observar que las vaiables significativas o con mayor impacto son las mismas en la mayoría de los modelos, así como su capacidad explicativa de los modelos, sin embargo logró diferenciarse en el RMSE que como su nombre lo dice es la diferencia entre los valores predichos por un modelo y los valores observados.