Librerías cargadas

library(foreign)        # Read Data Stored by 'Minitab', 'S', 'SAS', 'SPSS', 'Stata', 'Systat', 'Weka', 'dBase'
library(ggplot2)        # It is a system for creating graphics
library(lattice)
library(DataExplorer)
library(dplyr)          # A fast, consistent tool for working with data frame like objects
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ lubridate 1.9.3     ✔ tibble    3.2.1
## ✔ purrr     1.0.2     ✔ tidyr     1.3.1
## ✔ readr     2.1.4
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(tidyr)
library(mapview)        # Quickly and conveniently create interactive visualizations of spatial data with or without background maps
library(naniar)         # Provides data structures and functions that facilitate the plotting of missing values and examination of imputations.
library(tmaptools)
#library(maptools)       # A collection of functions to create spatial weights matrix objects from polygon 'contiguities', for summarizing these objects, and for permitting their use in spatial data analysis
library(tmap)           # For drawing thematic maps
## Breaking News: tmap 3.x is retiring. Please test v4, e.g. with
## remotes::install_github('r-tmap/tmap')
library(RColorBrewer)   # It offers several color palettes 
library(dlookr)         # A collection of tools that support data diagnosis, exploration, and transformation
## 
## Attaching package: 'dlookr'
## 
## The following object is masked from 'package:tidyr':
## 
##     extract
## 
## The following object is masked from 'package:base':
## 
##     transform
# predictive modeling
library(regclass)       # Contains basic tools for visualizing, interpreting, and building regression models
## Loading required package: bestglm
## Loading required package: leaps
## Loading required package: VGAM
## Loading required package: stats4
## Loading required package: splines
## Loading required package: rpart
## Loading required package: 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:dplyr':
## 
##     combine
## 
## The following object is masked from 'package:ggplot2':
## 
##     margin
## 
## Important regclass change from 1.3:
## All functions that had a . in the name now have an _
## all.correlations -> all_correlations, cor.demo -> cor_demo, etc.
## 
## 
## Attaching package: 'regclass'
## 
## The following object is masked from 'package:lattice':
## 
##     qq
library(mctest)         # Multicollinearity diagnostics
library(lmtest)         # Testing linear regression models
## Loading required package: zoo
## 
## Attaching package: 'zoo'
## 
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
## 
## 
## Attaching package: 'lmtest'
## 
## The following object is masked from 'package:VGAM':
## 
##     lrtest
library(spdep)          # A collection of functions to create spatial weights matrix objects from polygon 'contiguities', for summarizing these objects, and for permitting their use in spatial data analysis
## 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: sf
## Linking to GEOS 3.11.0, GDAL 3.5.3, PROJ 9.1.0; sf_use_s2() is TRUE
library(sf)             # A standardized way to encode spatial vector data
library(spData)         # Diverse spatial datasets for demonstrating, benchmarking and teaching spatial data analysis
library(spatialreg)     # A collection of all the estimation functions for spatial cross-sectional models
## Loading required package: Matrix
## 
## Attaching package: 'Matrix'
## 
## The following objects are masked from 'package:tidyr':
## 
##     expand, pack, unpack
## 
## 
## Attaching package: 'spatialreg'
## 
## The following objects are masked from 'package:spdep':
## 
##     get.ClusterOption, get.coresOption, get.mcOption,
##     get.VerboseOption, get.ZeroPolicyOption, set.ClusterOption,
##     set.coresOption, set.mcOption, set.VerboseOption,
##     set.ZeroPolicyOption
library(caret)          # The caret package (short for Classification And Rgression Training) contains functions to streamline the model training process for complex regression and classification problems.
## 
## Attaching package: 'caret'
## 
## The following object is masked from 'package:VGAM':
## 
##     predictors
## 
## The following object is masked from 'package:purrr':
## 
##     lift
library(e1071)          # Functions for latent class analysis, short time Fourier transform, fuzzy clustering, support vector machines, shortest path computation, bagged clustering, naive Bayes classifier, generalized k-nearest neighbor.
## 
## Attaching package: 'e1071'
## 
## The following objects are masked from 'package:dlookr':
## 
##     kurtosis, skewness
library(SparseM)        # Provides some basic R functionality for linear algebra with sparse matrices
## 
## Attaching package: 'SparseM'
## 
## The following object is masked from 'package:base':
## 
##     backsolve
library(Metrics)        # An implementation of evaluation metrics in R that are commonly used in supervised machine learning
## 
## Attaching package: 'Metrics'
## 
## The following objects are masked from 'package:caret':
## 
##     precision, recall
library(randomForest)   # Classification and regression based on a forest of trees using random inputs 
library(jtools)         # This is a collection of tools for more efficiently understanding and sharing the results of (primarily) regression analyses 
library(xgboost)        # The package includes efficient linear model solver and tree learning algorithms
## 
## Attaching package: 'xgboost'
## 
## The following object is masked from 'package:dplyr':
## 
##     slice
library(DiagrammeR)     # Build graph/network structures using functions for stepwise addition and deletion of nodes and edges 
library(effects) 
## Loading required package: carData
## Registered S3 method overwritten by 'survey':
##   method      from  
##   summary.pps dlookr
## Use the command
##     lattice::trellis.par.set(effectsTheme())
##   to customize lattice options for effects plots.
## See ?effectsTheme for details.
library(neuralnet)      
## 
## Attaching package: 'neuralnet'
## 
## The following object is masked from 'package:dplyr':
## 
##     compute

Carga de la Base de datos Health Insurance

health0<- read.csv("/Users/jazmincortez/Desktop/Concentración IA/MOD3/Act 1/health_insurance.csv")

1) Brevemente responder con tus propias palabras 2 de las siguientes 3 preguntas:

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

Este modelo, proporciona al algoritmo ejemplos de entrada junto con las respuestas deseadas (etiquetas), y el modelo aprende a mapear las entradas a las salidas, a diferencia del Unsupervised Machine Learning.

Algunos de las aplicaciones que pueden ser en Inteligencia de Negocios son: * Análisis de Sentimientos * Previsión de Ventas * Segmentación de Clientes * Detección de Fraudes

1.2) ¿Cuáles son los principales algoritmos de Supervised Machine Learning? Brevemente describir con tus propias palarbas 5 – 7 de los principales algoritmos de Supervised Machine Learning. Los que principalmente vimos en el Módulo son : * Regresión Lineal, que es de tipo Regresión, el cual ayuda a predecir valores continuos y establece una relación lineal entre la variable de entrada y la de salida. Regresión Logística, es de tipo clasificación binaria, permite predecir la probabilidad de que una instancia pertenezca a una clase particular y se utiliza comúnmente en problemas de clasificación binaria. Árboles de Decisión, que es de tipo clasificación y regresión, el que permite tomar decisiones basadas en reglas de “si-entonces” “what if”, dividiendo el conjunto de datos en ramas según las características más informativas. *Random Forest, que es de timpo ensamblaje de árboles de decisión, permite combinar múltiples árboles de decisión para mejorar la precisión y generalización del modelo.

1.3) ¿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 es una métrica que ajusta la R2 (coeficiente de determinación) para tener en cuenta el número de predictores (variables independientes) en el modelo. La R2 mide la proporción de la varianza en la variable dependiente que es explicada por el modelo. * El RMSE mide la raíz cuadrada del promedio de los cuadrados de los errores entre las predicciones del modelo y los valores reales. Es una medida de la precisión de las predicciones, donde valores más bajos indican un mejor ajuste del modelo a los datos observados.

2) Desarrollar Análisis Exploratorio de los Datos (EDA) que incluye los siguientes elementos:

a. Identificación de NA’s

health0[ , c("sex", "smoker", "region")] <- 
  lapply(health0[ , c("sex", "smoker", "region")], as.factor) #Es necesario cambiar variables de tipo factor a una escala adecuada para poder manejar los datos de dichas variable.

str(health0) #Esta salida proporciona información sobre el tipo de objeto de cada variable (columna)
## 'data.frame':    1338 obs. of  7 variables:
##  $ age     : int  19 18 28 33 32 31 46 37 37 60 ...
##  $ sex     : Factor w/ 2 levels "female","male": 1 2 2 2 2 1 1 1 2 1 ...
##  $ 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  : Factor w/ 2 levels "no","yes": 2 1 1 1 1 1 1 1 1 1 ...
##  $ region  : Factor w/ 4 levels "northeast","northwest",..: 4 3 3 2 2 3 3 2 1 2 ...
##  $ expenses: num  16885 1726 4449 21984 3867 ...
nas_columnas<- colSums(is.na(health0)) #Con esta función podemos obtener cuantos NA´s hay en total por columna
nas_columnas
##      age      sex      bmi children   smoker   region expenses 
##        0        0        0        0        0        0        0

b. Reemplazo de NA’s

#En este caso se observó que no hubo presencia de NA´s , en caso de haber existido, lo óptimo sería reemplazar esos NA´s con la mediana, ya que así conservamos la escala ordinal de los datos. Cabe recalcar que esto es importante cuando los datos son ordinales o categóricos, ya que la media podría no tener sentido en estos casos, y puede afectar la varianza de los datos y, por lo tanto, la dispersión general de la distribución.
#health <- health0 %>% mutate_all(~ifelse(is.na(.), median(., na.rm = TRUE), .))

c. Medidas descriptivas

summary(health0)
##       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.   :63770
desc <- describe(health0)
print(desc)
## # A tibble: 4 × 26
##   described_variables     n    na   mean     sd se_mean    IQR skewness kurtosis
##   <chr>               <int> <int>  <dbl>  <dbl>   <dbl>  <dbl>    <dbl>    <dbl>
## 1 age                  1338     0 3.92e1 1.40e1 3.84e-1 2.4 e1   0.0557  -1.25  
## 2 bmi                  1338     0 3.07e1 6.10e0 1.67e-1 8.4 e0   0.285   -0.0534
## 3 children             1338     0 1.09e0 1.21e0 3.30e-2 2   e0   0.938    0.202 
## 4 expenses             1338     0 1.33e4 1.21e4 3.31e+2 1.19e4   1.52     1.61  
## # ℹ 17 more variables: p00 <dbl>, p01 <dbl>, p05 <dbl>, p10 <dbl>, p20 <dbl>,
## #   p25 <dbl>, p30 <dbl>, p40 <dbl>, p50 <dbl>, p60 <dbl>, p70 <dbl>,
## #   p75 <dbl>, p80 <dbl>, p90 <dbl>, p95 <dbl>, p99 <dbl>, p100 <dbl>

b. Medidas de dispersión

health1 <- as.data.frame(sapply(health0, as.numeric))

#Desviación Estándar

desviacion <- sd(health1$bmi)
desviacion
## [1] 6.098382
#Varianza

#varianza <- var(health0)
#varianza
var <- sapply(health1,var) 
var
##          age          sex          bmi     children       smoker       region 
## 1.974014e+02 2.501596e-01 3.719027e+01 1.453213e+00 1.629689e-01 1.220771e+00 
##     expenses 
## 1.466524e+08
#Rango
rango <- range(health1)
#rango <- sapply(health0, function(x) max(x, na.rm = TRUE) - min(x, na.rm = TRUE))
rango
## [1]     0.00 63770.43
# Calcular el coeficiente de variación
media <- mean(health1)
## Warning in mean.default(health1): argument is not numeric or logical: returning
## NA
coef_variacion <- (desviacion / media) * 100
print(paste("Coeficiente de variación (%):", coef_variacion))
## [1] "Coeficiente de variación (%): NA"

c. Identificación de patrones y/o tendencias en los datos mediante el uso de gráficos incluyendo bar plots, line plots, pie plots, histogramas, matriz de correlación, box plot, scatter plot, qq- plot, etc Mostrar al menos 4 – 6 gráficos.

# Calcular la matriz de correlación
correlation_matrix <- cor(health1)
# Visualizar la matriz de correlación usando corrplot
library(corrplot)
## corrplot 0.92 loaded
corrplot(correlation_matrix, type = 'upper', order = 'hclust', addCoef.col = 'black')

plot_boxplot(health0, by = "expenses")

# Distribución de las variables con mi dataset numérico
plot_histogram(health1)

plot_missing(health0)

# Crear un scatter plot con BMI en el eje x y Expenses en el eje y
ggplot(health0, aes(x = bmi, y = expenses)) +
  geom_point() +
  labs(title = "Scatter Plot", x = "BMI", y = "Expenses") +
    scale_color_brewer(palette = "Paired")

ggplot(health0, aes(x = bmi, y = expenses, shape = smoker, color = smoker)) +
  geom_point() +
  theme_minimal() +
   scale_color_brewer(palette = "Paired")

# Crear un scatter plot con Age en el eje x y Expenses en el eje y
ggplot(health0, aes(x = age, y = expenses)) +
  geom_point() +
  labs(title = "Scatter Plot", x = "Age", y = "Expenses")

3) 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.

4) Estimación de cada uno de los siguientes modelos de Supervised Machine Learning (SML):

Cross Validation

set.seed(123)
partition<- createDataPartition(y=health0$expenses, p=0.7, list = F)
train= health0[partition, ]
test= health0[-partition,]

a. OLS Regresión

#Primer modelo OLS ordinal
ols_model_1 <- lm(expenses ~ age + bmi + region + smoker + sex + children, data = train)
summary(ols_model_1)
## 
## Call:
## lm(formula = expenses ~ age + bmi + region + smoker + sex + children, 
##     data = train)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -11592.2  -2737.3   -874.7   1476.4  29975.4 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     -11356.39    1152.06  -9.857  < 2e-16 ***
## age                248.93      14.12  17.626  < 2e-16 ***
## bmi                332.80      33.70   9.877  < 2e-16 ***
## regionnorthwest   -322.52     566.75  -0.569  0.56945    
## regionsoutheast  -1109.94     562.02  -1.975  0.04857 *  
## regionsouthwest   -878.02     565.53  -1.553  0.12087    
## smokeryes        24014.62     486.11  49.401  < 2e-16 ***
## sexmale           -249.64     395.68  -0.631  0.52825    
## children           441.91     162.28   2.723  0.00659 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6006 on 929 degrees of freedom
## Multiple R-squared:  0.7603, Adjusted R-squared:  0.7582 
## F-statistic: 368.2 on 8 and 929 DF,  p-value: < 2.2e-16
#Modelo Base OLS al haber seleccionado las variables más relevantes y significativas con impacto en expenses.
ols_model <- lm(expenses ~ age + bmi  +  children + smoker + sex + region, data = train)
summary(ols_model)
## 
## Call:
## lm(formula = expenses ~ age + bmi + children + smoker + sex + 
##     region, data = train)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -11592.2  -2737.3   -874.7   1476.4  29975.4 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     -11356.39    1152.06  -9.857  < 2e-16 ***
## age                248.93      14.12  17.626  < 2e-16 ***
## bmi                332.80      33.70   9.877  < 2e-16 ***
## children           441.91     162.28   2.723  0.00659 ** 
## smokeryes        24014.62     486.11  49.401  < 2e-16 ***
## sexmale           -249.64     395.68  -0.631  0.52825    
## regionnorthwest   -322.52     566.75  -0.569  0.56945    
## regionsoutheast  -1109.94     562.02  -1.975  0.04857 *  
## regionsouthwest   -878.02     565.53  -1.553  0.12087    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6006 on 929 degrees of freedom
## Multiple R-squared:  0.7603, Adjusted R-squared:  0.7582 
## F-statistic: 368.2 on 8 and 929 DF,  p-value: < 2.2e-16
# Predicción | Test 
pred_ols_model  <-predict(ols_model, test)
# RMSE OLS base
RMSE_ols_model  <- rmse(pred_ols_model, test$expenses)
#Ajuste del modelo con las variables más significativas seleccionadas,
#Modelo método LOG OLS
log_ols_model <- lm(log(expenses) ~ log(age)+ log(bmi)  + log( children + 0.01) + smoker + sex + region , data = train)
summary(log_ols_model)  
## 
## Call:
## lm(formula = log(expenses) ~ log(age) + log(bmi) + log(children + 
##     0.01) + smoker + sex + region, data = train)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.7896 -0.2404 -0.0796  0.1017  2.2109 
## 
## Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)           3.070300   0.267498  11.478  < 2e-16 ***
## log(age)              1.232184   0.037551  32.814  < 2e-16 ***
## log(bmi)              0.432876   0.073903   5.857 6.52e-09 ***
## log(children + 0.01)  0.034900   0.005671   6.154 1.12e-09 ***
## smokeryes             1.549217   0.035830  43.238  < 2e-16 ***
## sexmale              -0.093984   0.029154  -3.224  0.00131 ** 
## regionnorthwest      -0.071042   0.041794  -1.700  0.08950 .  
## regionsoutheast      -0.168645   0.041371  -4.076 4.96e-05 ***
## regionsouthwest      -0.122549   0.041712  -2.938  0.00339 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4426 on 929 degrees of freedom
## Multiple R-squared:  0.775,  Adjusted R-squared:  0.773 
## F-statistic: 399.9 on 8 and 929 DF,  p-value: < 2.2e-16
AIC(ols_model)     
## [1] 18995.15
AIC(log_ols_model)  
## [1] 1143.579

En ambos modelos tanto OLS como el ajuste lOG OLS , como una fuerte evidencia tenemos un p-value apto para rechazar la hipótesis nula , es decir al menos una de las variables predictoras que seleccioné tiene un efecto significativo en cuanto a afectar la variable de expenses.

# Predicción | Test | Modelo ajustado método LOG
pred_log_ols <-predict(log_ols_model, test)
# RMSE del ajuste con LOG
RMSE_log_ols    <- rmse(exp(pred_log_ols), test$expenses)
# Modelo Weighted Least Squares otro método para poder arreglar heterocedasticidad
wt <- 1/lm(abs(ols_model$residuals) ~ ols_model$fitted.values)$fitted.values^2

wls_model <- lm(log(expenses) ~ log(age)+ log(bmi)  + log( children + 0.01) + smoker + sex + region , data = train, weights =wt)

summary(wls_model)
## 
## Call:
## lm(formula = log(expenses) ~ log(age) + log(bmi) + log(children + 
##     0.01) + smoker + sex + region, data = train, weights = wt)
## 
## Weighted Residuals:
##        Min         1Q     Median         3Q        Max 
## -2.615e-04 -7.583e-05 -3.772e-05  2.198e-05  1.336e-03 
## 
## Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)           3.262157   0.257653  12.661  < 2e-16 ***
## log(age)              1.413024   0.041974  33.664  < 2e-16 ***
## log(bmi)              0.228507   0.075299   3.035  0.00248 ** 
## log(children + 0.01)  0.066577   0.006209  10.723  < 2e-16 ***
## smokeryes             1.615290   0.092300  17.501  < 2e-16 ***
## sexmale              -0.180134   0.031186  -5.776 1.04e-08 ***
## regionnorthwest      -0.119607   0.044744  -2.673  0.00765 ** 
## regionsoutheast      -0.220736   0.045493  -4.852 1.43e-06 ***
## regionsouthwest      -0.260738   0.044253  -5.892 5.33e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.0001618 on 929 degrees of freedom
## Multiple R-squared:  0.7046, Adjusted R-squared:  0.702 
## F-statistic:   277 on 8 and 929 DF,  p-value: < 2.2e-16
# Predicción - cv test
pred_wls_model <-predict(wls_model, test)
# RMSE
RMSE_wls       <- rmse(exp(pred_wls_model), test$expenses)

b. SAR

c. SEM

d. XGBoost Regresión

# Extracción de variables determinadas para el modelo
df_xboost <- health0 %>% 
  dplyr::select(  expenses,age, bmi, children, smoker, sex, region )


# Transformaciones previas a la construcción del modelo
df_xboost$expenses                       <- log(df_xboost$expenses)
df_xboost$bmi                               <- log(df_xboost$bmi)
df_xboost$children                           <- log(df_xboost$children + 0.01)
df_xboost$age                                <- log(df_xboost$age)


summary(df_xboost)
##     expenses           age             bmi           children        smoker    
##  Min.   : 7.023   Min.   :2.890   Min.   :2.773   Min.   :-4.60517   no :1064  
##  1st Qu.: 8.464   1st Qu.:3.296   1st Qu.:3.270   1st Qu.:-4.60517   yes: 274  
##  Median : 9.147   Median :3.664   Median :3.414   Median : 0.00995             
##  Mean   : 9.099   Mean   :3.597   Mean   :3.403   Mean   :-1.67105             
##  3rd Qu.: 9.720   3rd Qu.:3.932   3rd Qu.:3.547   3rd Qu.: 0.69813             
##  Max.   :11.063   Max.   :4.159   Max.   :3.972   Max.   : 1.61144             
##      sex            region   
##  female:662   northeast:324  
##  male  :676   northwest:325  
##               southeast:364  
##               southwest:325  
##                              
## 
#df_xboost
#Cross Validation para XGBOOST
set.seed(123)
partitionxboost<- createDataPartition(y=df_xboost$expenses, p=0.7, list = F)
train2= df_xboost[partitionxboost, ]
test2= df_xboost[-partitionxboost,]


# Espeficar  la ubicación de la variable dependiente en el set de train .
train_x = data.matrix(train2[, -1])
train_y = train2[,1]

# Espeficar  la ubicación de la variable dependiente en el set de test .
test_x = data.matrix(test2[, -1])
test_y = test2[, 1]

# Creación de la matriz para los tests de train y test del modelo XGBoost
xgb_train = xgb.DMatrix(data = train_x, label = train_y)
xgb_test  = xgb.DMatrix(data = test_x, label = test_y)

# Evaluación del modelo
watchlist = list(train=xgb_train, test=xgb_test)
model_xgb = xgb.train(data=xgb_train, max.depth=3, watchlist=watchlist, nrounds=70) # the more the number of rounds selected, the longer the time to display the results. 
## [1]  train-rmse:6.072340 test-rmse:6.080231 
## [2]  train-rmse:4.269586 test-rmse:4.273616 
## [3]  train-rmse:3.009975 test-rmse:3.016854 
## [4]  train-rmse:2.130089 test-rmse:2.138267 
## [5]  train-rmse:1.519448 test-rmse:1.532367 
## [6]  train-rmse:1.099072 test-rmse:1.116282 
## [7]  train-rmse:0.813936 test-rmse:0.837907 
## [8]  train-rmse:0.626154 test-rmse:0.651870 
## [9]  train-rmse:0.506887 test-rmse:0.540346 
## [10] train-rmse:0.432948 test-rmse:0.473401 
## [11] train-rmse:0.389128 test-rmse:0.435379 
## [12] train-rmse:0.363638 test-rmse:0.414419 
## [13] train-rmse:0.350199 test-rmse:0.401414 
## [14] train-rmse:0.342434 test-rmse:0.396992 
## [15] train-rmse:0.337626 test-rmse:0.392986 
## [16] train-rmse:0.333211 test-rmse:0.391730 
## [17] train-rmse:0.330166 test-rmse:0.392686 
## [18] train-rmse:0.328498 test-rmse:0.393087 
## [19] train-rmse:0.327184 test-rmse:0.391769 
## [20] train-rmse:0.325128 test-rmse:0.392075 
## [21] train-rmse:0.324032 test-rmse:0.391605 
## [22] train-rmse:0.322020 test-rmse:0.391637 
## [23] train-rmse:0.319798 test-rmse:0.393234 
## [24] train-rmse:0.318123 test-rmse:0.392708 
## [25] train-rmse:0.317147 test-rmse:0.392356 
## [26] train-rmse:0.316132 test-rmse:0.392560 
## [27] train-rmse:0.314582 test-rmse:0.392418 
## [28] train-rmse:0.311536 test-rmse:0.394391 
## [29] train-rmse:0.310501 test-rmse:0.395032 
## [30] train-rmse:0.308835 test-rmse:0.395847 
## [31] train-rmse:0.307523 test-rmse:0.396483 
## [32] train-rmse:0.305734 test-rmse:0.397351 
## [33] train-rmse:0.305119 test-rmse:0.396509 
## [34] train-rmse:0.304878 test-rmse:0.396479 
## [35] train-rmse:0.302952 test-rmse:0.397100 
## [36] train-rmse:0.302601 test-rmse:0.397160 
## [37] train-rmse:0.300799 test-rmse:0.399331 
## [38] train-rmse:0.300183 test-rmse:0.399595 
## [39] train-rmse:0.299513 test-rmse:0.399475 
## [40] train-rmse:0.299039 test-rmse:0.399634 
## [41] train-rmse:0.298510 test-rmse:0.399970 
## [42] train-rmse:0.297043 test-rmse:0.401154 
## [43] train-rmse:0.296887 test-rmse:0.401448 
## [44] train-rmse:0.295576 test-rmse:0.401904 
## [45] train-rmse:0.294417 test-rmse:0.401840 
## [46] train-rmse:0.293584 test-rmse:0.401712 
## [47] train-rmse:0.292662 test-rmse:0.401116 
## [48] train-rmse:0.291352 test-rmse:0.401927 
## [49] train-rmse:0.291132 test-rmse:0.402094 
## [50] train-rmse:0.290560 test-rmse:0.402451 
## [51] train-rmse:0.289503 test-rmse:0.401954 
## [52] train-rmse:0.288683 test-rmse:0.402227 
## [53] train-rmse:0.286285 test-rmse:0.403134 
## [54] train-rmse:0.284967 test-rmse:0.403995 
## [55] train-rmse:0.283773 test-rmse:0.404111 
## [56] train-rmse:0.282890 test-rmse:0.403945 
## [57] train-rmse:0.280574 test-rmse:0.404495 
## [58] train-rmse:0.279666 test-rmse:0.405907 
## [59] train-rmse:0.279033 test-rmse:0.406338 
## [60] train-rmse:0.278440 test-rmse:0.406274 
## [61] train-rmse:0.277746 test-rmse:0.406726 
## [62] train-rmse:0.277380 test-rmse:0.406652 
## [63] train-rmse:0.276981 test-rmse:0.406254 
## [64] train-rmse:0.276834 test-rmse:0.406287 
## [65] train-rmse:0.276475 test-rmse:0.406133 
## [66] train-rmse:0.275757 test-rmse:0.406308 
## [67] train-rmse:0.274788 test-rmse:0.406433 
## [68] train-rmse:0.271534 test-rmse:0.407788 
## [69] train-rmse:0.270561 test-rmse:0.407387 
## [70] train-rmse:0.269874 test-rmse:0.407535
#Creación del Modelo ajustado

reg_xgb = xgboost(data = xgb_train, max.depth = 3, nrounds = 100, verbose = 0) # setting verbose = 0 avoids to display the training and testing error for each round. 


# Predicción del modelo - Test (con el test2 )
prediction_xgb_test<-predict(reg_xgb, xgb_test)

#RMSE XGBOOST 
RMSE_XGBOOST<- rmse(prediction_xgb_test, test2$expenses)

# Diagnóstico de los residuales del modelo
xgb_reg_residuals<-test2$expenses - prediction_xgb_test
plot(xgb_reg_residuals, xlab= "Dependent Variable", ylab = "Residuals", main = 'XGBoost Regression Residuals')
abline(0,0)

* Los puntos se distribuyen aproximadamente en una línea recta, lo que indica que los residuales son normales. * Hay algunos puntos que se alejan de la línea recta, lo que podría indicar que hay algunos valores atípicos en los datos.

# Gráfico de los 3 primeros árboles del modelo
xgb.plot.tree(model=reg_xgb, trees=0:2)
#Importancia de las variables explicativas de acuerdo al modelo.
importance_matrix <- xgb.importance(model = reg_xgb)
xgb.plot.importance(importance_matrix, xlab = "Explanatory Variables X's Importance")

e. Decision Trees

library(rpart)
library(rpart.plot)
decision_tree <- rpart(log(expenses) ~ log(age) + log(bmi)  + log( children + 0.01) + smoker + sex + region, data = train)

rpart.plot(decision_tree)

# Predicción - test | Para Árboles de Decisión
decision_tree_prediction <- predict(decision_tree,test)

# RMSE 
RMSE_decision_tree <- rmse(exp(decision_tree_prediction), test$expenses)

f. Random Forest

# Creación del modelo Random Forest
rf_model <- randomForest(expenses ~ age + bmi  +  children + smoker + sex + region, data = train, proximity=TRUE)

print(rf_model)
## 
## Call:
##  randomForest(formula = expenses ~ age + bmi + children + smoker +      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: 21415171
##                     % Var explained: 85.63
# Predicción – cv test data
rf_prediction <- predict(rf_model,test)

# RMSE
RMSE_rf <- rmse(rf_prediction, test$expenses)


# Importancia de las variables
varImpPlot(rf_model, n.var = 4, main = "Top 4 - Variable")

importance(rf_model)                                       
##          IncNodePurity
## age        17854853081
## bmi        19079473372
## children    2931171353
## smoker     85901684084
## sex         1111177849
## region      3043017719

g. Neural Networks Regresión

# Semilla para reproductibilidad
set.seed(123)

# Transformación a variables numéricas
train_numeric <- as.data.frame(lapply(train2, as.numeric))
test_numeric <- as.data.frame(lapply(test2, as.numeric))


# Creación del modelo con un mayor número de pasos
#nn_model <- neuralnet(expenses ~ ., data = train_numeric, hidden = c(3,3), linear.output = TRUE, stepmax = 1000000) 
# Gráfico de la Red Neuronal 
#summary(nn_model)
# Predicción -  test numeric
#nnet_prediction <- predict(nn_model, test_numeric)

# RMSE
#RMSE_nnet <- rmse(exp(nnet_prediction), exp(test_numeric$expenses))

#plot(nn_model)

5) Pruebas de Diagnóstico de los Resultados Obtenidos de la Estmación de Modelos de Regresión

a. Multicolinealidad

#Un VIF alto (generalmente mayor que 5 o 10) puede indicar multicolinealidad.
#VIF Modelo OLS
VIF(ols_model)
##              GVIF Df GVIF^(1/(2*Df))
## age      1.026148  1        1.012989
## bmi      1.096808  1        1.047286
## children 1.006354  1        1.003172
## smoker   1.019387  1        1.009647
## sex      1.017577  1        1.008750
## region   1.095389  3        1.015301
#VIF Modelo ajustado con LOG OLS 
VIF(log_ols_model)
##                          GVIF Df GVIF^(1/(2*Df))
## log(age)             1.035111  1        1.017404
## log(bmi)             1.092813  1        1.045377
## log(children + 0.01) 1.014228  1        1.007089
## smoker               1.020115  1        1.010007
## sex                  1.017556  1        1.008740
## region               1.091687  3        1.014728
#VIF Modelo WLS
VIF(wls_model)
##                          GVIF Df GVIF^(1/(2*Df))
## log(age)             1.145330  1        1.070201
## log(bmi)             1.157654  1        1.075943
## log(children + 0.01) 1.066923  1        1.032920
## smoker               1.009111  1        1.004545
## sex                  1.018886  1        1.009399
## region               1.107699  3        1.017194
# No se observa Multicolinealidad en el modelo absoluto, logaritmico ni con la alternativa de wls, ya que , un VIF mayor que 5 o 10 a menudo se considera una señal de multicolinealidad, y no se observaron

### Nota : La Multicolinealidad se debe detectar desde el modelo de Regresión Lineal Base, cuando se tiene los valores absolutos

b. Heterocedasticidad

En este caso para poder identificar si hay o no presencia de heterocedasticisad, debemos de considerar los siguientes criterios: * Si el p-value es mayor que el nivel de significancia elegido (por ejemplo, 0.05), no hay evidencia suficiente para rechazar la hipótesis nula. En este caso, se interpretaría que no hay indicios de heterocedasticidad. H0: No hay presencia de heterocedasticidad p-value> 0.05

  • Si el p-value es menor que el nivel de significancia, se considera evidencia suficiente para rechazar la hipótesis nula. En este caso, se interpretaría que hay evidencia de heterocedasticidad. Ha: Hay presencia de Heterocedasticidad p-value< 0.05
#BP test Modelo Base OLS
bptest(ols_model)
## 
##  studentized Breusch-Pagan test
## 
## data:  ols_model
## BP = 97.012, df = 8, p-value < 2.2e-16
#En el modelo base de OLS se observa heterocedasticidad, ya que obtuvimos un p-value menor a 0.05

#BP test Modelo ajustado LOG OLS
bptest(log_ols_model) 
## 
##  studentized Breusch-Pagan test
## 
## data:  log_ols_model
## BP = 61.174, df = 8, p-value = 2.742e-10
#En el modelo LOG OLS que es en dónde se trató de arreglar la heterocedasticidad, si mejoró el p-value que era lo que se buscaba, sin embargo sigue habiendo presencia de heterocedasticidad.

#Alternativa para poder seguir mejorando la Heterocedasticidad
#BP test Modelo WLS
bptest(wls_model)
## 
##  studentized Breusch-Pagan test
## 
## data:  wls_model
## BP = 1.9488e-05, df = 8, p-value = 1
#Con esta alternativa, se observa que ya se pudo arreglar la presencia de heterocedasticidad, es decir que con esta prueba , presentamos ya una homocedasticidad de los datos.

c. Autocorrelación Serial

No es posible ya que no tenemos la información(datos) necesarios para realizar este modelo. ## d. Autocorrelación Espacial No es posible ya que no tenemos la información(datos) necesarios para realizar este modelo.

e. Normalidad de los Residuales

Evaluación de la muestra de datos para identificar si sigue o no una distribución normal, en dónde, si el p-value del test de Shapiro-Wilk es menor que un nivel de significancia de 0.05, se rechaza la hipótesis nula, es decir que estaremos observando que los datos no siguen una distribución normal. H0: Hay una distribución normal p-value> 0.05 Ha: No hay una distribución normal p-value< 0.05

#Modelo OLS
norm_ols <- resid(ols_model)
shapiro_ols <- shapiro.test(norm_ols)
print(shapiro_ols) 
## 
##  Shapiro-Wilk normality test
## 
## data:  norm_ols
## W = 0.90396, p-value < 2.2e-16
#Modelo LOG_OLS
norm_log_ols <- resid(log_ols_model)
shapiro_log_ols <- shapiro.test(norm_log_ols)
print(shapiro_log_ols) 
## 
##  Shapiro-Wilk normality test
## 
## data:  norm_log_ols
## W = 0.84456, p-value < 2.2e-16
#Alternativa
# Modelo wls
norm_wls <- resid(wls_model)
shapiro_wls <- shapiro.test(norm_wls)
print(shapiro_wls)
## 
##  Shapiro-Wilk normality test
## 
## data:  norm_wls
## W = 0.82565, p-value < 2.2e-16

Mediante la prueba de shapiro , la normalidad de los residuos generados por tres modelos distintos: el modelo de Regresión Lineal Ordinaria (OLS), el modelo Log-OLS y el modelo de Regresión Ponderada por Mínimos Cuadrados (WLS), se puede observar que los 3 modelos no están siguiendo una distribución normal de los datos, ya que todos presentan un p-value menor a 0.05. Siendo así que los resultados de las pruebas en los tres modelos están indicando que los residuos no siguen una distribución normal, ya que los valores de p asociados son significativamente bajos.

Nota: En caso de que las pruebas de diagnóstico identifiquen cualquiera de los anteriores a) – e) plantear una solución para mejorar la estimación de la especificiación del modelo.

6) Evaluación y Selección de Modelo de Regresión

a. Mediante el cálculo de la métrica RMSE para cada uno de los modelos estimados en 4) seleccionar el modelo que muestra los mejores resultados estimados.

RMSE_ols_model
## [1] 6196.778
RMSE_log_ols
## [1] 7535.944
RMSE_wls
## [1] 8890.723
RMSE_XGBOOST
## [1] 0.4147628
RMSE_decision_tree
## [1] 5232.392
RMSE_rf
## [1] 4989.633
#RMSE_nnet

En el proceso de construcción y evaluación de diversos modelos predictivos para nuestro conjunto de datos de health insurance, se hizo la exploración una variedad de enfoques, desde métodos lineales tradicionales hasta técnicas más avanzadas de aprendizaje automático. Nuestro objetivo principal fue identificar el modelo que proporciona las predicciones más precisas posibles para nuestro problema específico que en este caso es conocer como se afectan los gastos médicos.

Después de realizar una comparación de los modelos , medimos el rendimiento utilizando el Root Mean Squared Error (RMSE), una métrica que cuantifica la precisión de las predicciones en términos de la raíz cuadrada de las diferencias entre las predicciones y los valores reales, viendo el valor de errores lo menormente posible.

Por lo que con estos resultados observamos que el modelo de XGBOOST proporcionó el menor RMSE de 0.4147628, indicando una mayor precisión en las predicciones en comparación con los demás modelos. Esto sugiere que el modelo Random Forest es el más adecuado para nuestro conjunto de datos y podría ser la elección más adecuada para futuras predicciones.

b. Presentar los valores de la métrica RMSE de cada uno de los modelos estimados en 4) en un gráfico de barras.

rmse_data <- data.frame(
  Model = c("OLS", "LOG OLS","WLS" ,"XGBoost", "RPART", "RF"),
  RMSE = c(RMSE_ols_model, RMSE_log_ols,RMSE_wls, RMSE_XGBOOST,RMSE_decision_tree, RMSE_rf))
rmse_data <- rmse_data[order(-rmse_data$RMSE), ]

rmse_data$RMSE <- round(rmse_data$RMSE, 4)

ggplot(rmse_data, aes(x = reorder(Model, -RMSE), y = RMSE)) +
  geom_bar(stat = "identity", fill = "purple") +
  labs(title = "RMSE por Modelo", x = "Modelo", y = "RMSE") +
  coord_flip() +
  geom_text(aes(label = RMSE), hjust = 1.2)

7) Desarrollar una breve descripción de los 6 – 10 principales hallazgos de:

a. EDA

  • La variable Edad tiene una distribución normal, con la mayoría de las personas en el rango de 20 a 60 años.
  • La variable BMI tiene una distribución ligeramente sesgada a la derecha, con un mayor número de personas con un IMC alto.
  • La variable Niños tiene una distribución bimodal, con un pico en “Sí” y otro en “No”.
  • La variable Gastos tiene una distribución muy sesgada a la derecha, con un mayor número de personas con gastos bajos y un pequeño número de personas con gastos muy altos.
  • La edad es un factor importante que influye en los gastos del seguro de salud, por lo que las personas mayores tienen un mayor riesgo de tener gastos de seguro más altos.
  • Se observa una tendencia general al alza en los gastos del seguro a medida que aumenta el IMC.
  • Los personas con hábitos sin tabaco (no fumadores) tienden a estar por debajo de la línea de tendencia general de gastos médicos mayores, mientras que los que sí consumen tabaco (fumadores) tienden a estar por encima.
  • Hay una gran dispersión en los datos, lo que significa que hay mucha variabilidad en los gastos del seguro incluso dentro de cada categoría de BMI y hábito de fumar.
  • El BMI y el hábito de fumar son dos factores importantes que influyen en los gastos del seguro de salud.

b. Modelo seleccionado:

i. ¿Cuáles son las variables que contribuyen a explicar los cambios de la principal variable de estudio?

Las variables que principalmente ayudana explicar mejor el modelos, son age, bmi, children , smoker , sex, region

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

  • La variable smoker tiene una importancia significativamente alta en comparación con otras variables. Esto sugiere que la condición de si una persona fuma o no es crucial para predecir la variable objetivo en el modelo.
  • Las variables age y bmi también tienen una importancia considerable, aunque no tan alta como smoker. Esto indica que la edad y el índice de masa corporal (IMC) también son factores relevantes en la predicción del modelo, pero en menor medida que el hábito de fumar.
  • La variable children tiene una importancia más baja en comparación con las anteriores, pero aún sigue contribuyendo al modelo. Puede sugerir que el número de niños asegurados (hijos) tiene un impacto, aunque menos impactante, en las predicciones del modelo.
  • Las variables sex (género) y region tienen importancias relativamente más bajas en comparación con las otras variables. Esto podría indicar que estas características tienen menos influencia en las predicciones del modelo en el conjunto de datos particular.

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

  • Se podría inferir decir que el Random Forest también está presentando un RMSE relativamente bajo, por lo que podria ser otro modelo similar a considerar,