Actividad 1 M3

Instalar Librerias

library(foreign)       
library(ggplot2)        
library(dplyr)          
library(mapview)        
library(naniar)         
library(tmap)           
library(RColorBrewer)   
library(dlookr)         
library(regclass)      
library(mctest)         
library(lmtest)        
library(spdep)          
library(sf)            
library(spData)         
library(spatialreg)     
library(caret)          
library(e1071)         
library(SparseM)       
library(Metrics)        
library(randomForest)   
library(jtools)        
library(xgboost)        
library(DiagrammeR)    
library(effects)        
library(shinyjs)
library(sp)
library(geoR)
library(gstat)
library(caret)
library(readtext)
library(kableExtra)
library(nortest) 
library(rpart.plot)
library(neuralnet)

1.Preguntas

#1)¿Qué es Supervised Machine Learning y cuáles son algunas de sus aplicaciones en Inteligencia de Negocios? Un modelo de aprendizaje supervisado es aquel a el que se le entregan datos etiquetados, para realizar predicciones o clasificaciones de estos datos de manera precisa, algunas de sus principales aplicaciones son las investigaciones de mercado, para saber cual es el mejor publico, segmento, cual es la mejor ubicacion para tu producto e incluso cual es el mejor momento para vender tu producto, tambien ppueden ayudar a fijar precios dinamicos de vuelos o viajes, a travez de datos que estan en constante cambio como la epoca del año, clima, demanda, entre otros

#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 y el RMSE son medidas o metodos de evaluacion, pero se centran en cosas distintas, la R2 evalua la calidad del modelo atraves de la varianza de las variables, nos ayuda a ver que variables no son significantes para el modelo, mientras que el RMSE evalua la precision del modelo ya que se centra en la diferencia entre los valores observados y los valores predichos por el modelo

#Obtener Base de datos

df <- read.csv("C:\\Users\\Luis Rodriguez\\Downloads\\health_insurance.csv")
#df

df[ , c("sex", "smoker", "region")] <- 
  lapply(df[ , 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.

2.Analisis Exploratorio de los datos (EDA):

Identificacion de NA’s:

# Identificación de NA's por columna
nas <- colSums(is.na(df))
nas

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

Remplazo de NA’s: En este caso no tenemos NA’s por lo cual no es necesario realizar un reemplazo

Medidas Descriptivas:

# Medidas descriptivas para todas las columnas numéricas
summary(df)

##       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

Medidas de Dispersion:

varianza <- apply(df,2,var)
varianza

##          age          sex          bmi     children       smoker       region 
## 1.974014e+02           NA 3.719027e+01 1.453213e+00           NA           NA 
##     expenses 
## 1.466524e+08

desviacion <- apply(df,2,sd)
desviacion

##          age          sex          bmi     children       smoker       region 
##    14.049960           NA     6.098382     1.205493           NA           NA 
##     expenses 
## 12110.011240

Identifiacion de patrones y/o tendencias

Ejemplos de gráficos

boxplot(expenses ~ sex, data = df, col = c("red","green","blue","purple", main = "Expenses por sexo"))

# Diagrama de caja para la variable 'bmi' por género
boxplot(bmi ~ region, data = df, col = c("blue","purple", "cyan", "lightpink"), main = "BMI por Region", xlab = "Region", ylab = "BMI")

boxplot(expenses ~ region, data = df, col = c("blue","purple", "cyan", "lightpink"), main = "Expenses por Region", xlab = "Region", ylab = "Expenses")

# Gráfico de dispersión para la relación entre 'age' y 'expenses'
plot(df$children, df$expenses, main = "Children vs. Gastos Médicos", xlab = "Children", ylab = "Gastos Médicos", col = "blue")

# Gráfico de dispersión para la relación entre 'smoke' y 'expenses'
plot(df$age, df$expenses, main = "Edad vs. Gastos Médicos", xlab = "Edad", ylab = "Gastos Médicos", col = "darkgreen")

# Gráfico de dispersión para la relación entre 'bmi' y 'expenses' tomando en cuenta 'Expenses'
ggplot(df, aes(x=bmi, y=expenses, shape=smoker, color=smoker)) +
  geom_point() +
  theme_minimal() +
  scale_color_brewer(palette = "Greens")

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

Basándonos en los resultados del Análisis Exploratorio de Datos (EDA), identificamos que las variables más relevantes para predecir los gastos médicos (expenses) son el índice de masa corporal (bmi), la edad (age) y si la persona es fumadora (smoker).

Al examinar los gráficos de dispersión, notamos que existe una tendencia clara: a medida que la edad aumenta, también lo hacen los gastos médicos. Además, observamos que las personas que fuman y tienen un índice de masa corporal elevado tienden a registrar gastos médicos más altos.

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

A.OLS Regresión

lets split data into training and test sets the training set is used to build the model and the test set to evaluate its predictive accuracy.

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

ols_model <- lm(expenses ~ age + bmi + smoker, data = df)
summary(ols_model)

## 
## Call:
## lm(formula = expenses ~ age + bmi + smoker, data = df)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -12413  -2970   -977   1476  28961 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -11679.05     937.53  -12.46   <2e-16 ***
## age            259.53      11.93   21.75   <2e-16 ***
## bmi            322.69      27.49   11.74   <2e-16 ***
## smokeryes    23822.61     412.86   57.70   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6092 on 1334 degrees of freedom
## Multiple R-squared:  0.7475, Adjusted R-squared:  0.7469 
## F-statistic:  1316 on 3 and 1334 DF,  p-value: < 2.2e-16

log_ols_model <- lm(log(expenses) ~ log(age) + log(bmi) + smoker, data = df)
summary(log_ols_model)

## 
## Call:
## lm(formula = log(expenses) ~ log(age) + log(bmi) + smoker, data = df)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.02436 -0.23021 -0.05241  0.09755  2.15745 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  2.98367    0.23094  12.920  < 2e-16 ***
## log(age)     1.27543    0.03232  39.468  < 2e-16 ***
## log(bmi)     0.35582    0.06233   5.708  1.4e-08 ***
## smokeryes    1.54369    0.03108  49.666  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4587 on 1334 degrees of freedom
## Multiple R-squared:  0.7517, Adjusted R-squared:  0.7512 
## F-statistic:  1346 on 3 and 1334 DF,  p-value: < 2.2e-16

AIC(ols_model)      # AIC = 27123.76

## [1] 27123.76

AIC(log_ols_model)  # AIC = 1717.374

## [1] 1717.374

RMSE_ols_model     <- sqrt(mean(ols_model$residuals^2)) #RMSE = 0.7475
RMSE_log_ols_model <- sqrt(mean(log_ols_model$residuals^2)) #RMSE = 7517

#Modelo OLS ABSOLUTO
residuos_ols <- resid(ols_model)
#residuos_ols
valores_ajustados_ols <- fitted(ols_model)

plot(valores_ajustados_ols, residuos_ols, main = "Residuos vs. Ajustes", xlab = "Valores Ajustados", ylab = "Residuos")

#Modelo LOG_OLS
residuos_log_ols <- resid(log_ols_model)
#residuos_log_ols
valores_ajustados_log_ols <- fitted(ols_model)

plot(valores_ajustados_log_ols, residuos_log_ols, main = "Residuos vs. Ajustes del LOG", xlab = "Valores Ajustados", ylab = "Residuos")

# Pruebas del modelo ols y ols log

# Prueba de multicolinealidad
vif_result <- car::vif(ols_model)# Factor de inflación de la varianza
vif_result

##      age      bmi   smoker 
## 1.012764 1.012146 1.000672

#No hay multicolinealidad entre la variable predictora y el resto de las variables predictoras.

# Prueba de heterocedasticidad (Test de Breusch-Pagan)
bp_test <- bptest(ols_model)
bp_test

## 
##  studentized Breusch-Pagan test
## 
## data:  ols_model
## BP = 112.59, df = 3, p-value < 2.2e-16

#El valor p es menor que el nivel de significancia (p < 0.05), hay evidencia de heterocedasticidad y se rechaza la hipótesis nula de homocedasticidad.

# Prueba de autocorrelación serial (Test de Durbin-Watson)
dw_test <- dwtest(ols_model)
dw_test

## 
##  Durbin-Watson test
## 
## data:  ols_model
## DW = 2.0765, p-value = 0.9191
## alternative hypothesis: true autocorrelation is greater than 0

#El valor del estadístico de prueba está cerca de 2, lo que sugiere que los residuos del modelo de regresión no muestran una autocorrelación serial significativa. Un valor de DW cercano a 2 indica que no hay autocorrelación serial.

# Prueba de normalidad de los residuos (Anderson-Darling)
ad_test <- ad.test(ols_model$residuals)
ad_test

## 
##  Anderson-Darling normality test
## 
## data:  ols_model$residuals
## A = 41.885, p-value < 2.2e-16

#Dado que el valor p es extremadamente bajo, se rechaza la hipótesis nula de que los residuos siguen una distribución normal. Los residuos del modelo de regresión no se ajustan a una distribución normal según el test de Anderson-Darling. 

# Imprimir resultados
print("Prueba de Multicolinealidad (VIF):")

## [1] "Prueba de Multicolinealidad (VIF):"

print(vif_result)

##      age      bmi   smoker 
## 1.012764 1.012146 1.000672

print("Prueba de Heterocedasticidad (Breusch-Pagan):")

## [1] "Prueba de Heterocedasticidad (Breusch-Pagan):"

print(bp_test)

## 
##  studentized Breusch-Pagan test
## 
## data:  ols_model
## BP = 112.59, df = 3, p-value < 2.2e-16

print("Prueba de Autocorrelación Serial (Durbin-Watson):")

## [1] "Prueba de Autocorrelación Serial (Durbin-Watson):"

print(dw_test)

## 
##  Durbin-Watson test
## 
## data:  ols_model
## DW = 2.0765, p-value = 0.9191
## alternative hypothesis: true autocorrelation is greater than 0

#print("Prueba de Autocorrelación Espacial (Moran's I):")
#print(moran_test)
print("Prueba de Normalidad de los Residuos (Anderson-Darling):")

## [1] "Prueba de Normalidad de los Residuos (Anderson-Darling):"

print(ad_test)

## 
##  Anderson-Darling normality test
## 
## data:  ols_model$residuals
## A = 41.885, p-value < 2.2e-16

# Prueba de multicolinealidad
vif_result <- car::vif(log_ols_model)# Factor de inflación de la varianza
vif_result

## log(age) log(bmi)   smoker 
## 1.012653 1.012136 1.000521

#Los valores de VIF cercanos a 1 para todas las variables indican que no hay evidencia sustancial de multicolinealidad entre las variables incluidas en el modelo. 

# Prueba de heterocedasticidad (Test de Breusch-Pagan)
bp_test <- bptest(log_ols_model)
bp_test

## 
##  studentized Breusch-Pagan test
## 
## data:  log_ols_model
## BP = 102.39, df = 3, p-value < 2.2e-16

#Dado que el valor p es muy pequeño, se rechaza la hipótesis nula de que no hay heterocedasticidad en los residuos. En otras palabras, hay evidencia de heterocedasticidad en los residuos del modelo de regresión. Podriamos utilizar métodos de ponderación para dar más peso a las observaciones con menor varianza y menos peso a las observaciones con mayor varianza. Esto puede ayudar a compensar la heterocedasticidad y mejorar la precisión de las estimaciones.

# Prueba de autocorrelación serial (Test de Durbin-Watson)
dw_test <- dwtest(log_ols_model)
dw_test

## 
##  Durbin-Watson test
## 
## data:  log_ols_model
## DW = 2.0265, p-value = 0.6861
## alternative hypothesis: true autocorrelation is greater than 0

#El valor del estadístico de prueba está cerca de 2, lo que sugiere que los residuos del modelo de regresión no muestran una autocorrelación serial significativa. Un valor de DW cercano a 2 indica que no hay autocorrelación serial.

# Prueba de normalidad de los residuos (Anderson-Darling)
ad_test <- ad.test(log_ols_model$residuals)
ad_test

## 
##  Anderson-Darling normality test
## 
## data:  log_ols_model$residuals
## A = 55.307, p-value < 2.2e-16

#Dado que el valor p es mayor que el nivel de significancia, no hay suficiente evidencia para rechazar la hipótesis nula de que no hay autocorrelación serial en los residuos del modelo. En este caso, la hipótesis nula sugiere que no hay autocorrelación positiva en los residuos del modelo.

# Imprimir resultados
print("Prueba de Multicolinealidad (VIF):")

## [1] "Prueba de Multicolinealidad (VIF):"

print(vif_result)

## log(age) log(bmi)   smoker 
## 1.012653 1.012136 1.000521

print("Prueba de Heterocedasticidad (Breusch-Pagan):")

## [1] "Prueba de Heterocedasticidad (Breusch-Pagan):"

print(bp_test)

## 
##  studentized Breusch-Pagan test
## 
## data:  log_ols_model
## BP = 102.39, df = 3, p-value < 2.2e-16

print("Prueba de Autocorrelación Serial (Durbin-Watson):")

## [1] "Prueba de Autocorrelación Serial (Durbin-Watson):"

print(dw_test)

## 
##  Durbin-Watson test
## 
## data:  log_ols_model
## DW = 2.0265, p-value = 0.6861
## alternative hypothesis: true autocorrelation is greater than 0

#print("Prueba de Autocorrelación Espacial (Moran's I):")
#print(moran_test)
print("Prueba de Normalidad de los Residuos (Anderson-Darling):")

## [1] "Prueba de Normalidad de los Residuos (Anderson-Darling):"

print(ad_test)

## 
##  Anderson-Darling normality test
## 
## data:  log_ols_model$residuals
## A = 55.307, p-value < 2.2e-16

B.SAR

No podemos realizar este modelo porque no tenemos todos los datos necesarios

C.SEM

No podemos realizar este modelo porque no tenemos todos los datos necesarios

D.XGBoost Regresión

df_alt <- df %>% select(age,bmi,smoker,expenses,children,sex, region)

df_alt$bmi <- as.numeric(df_alt$bmi)
df_alt$expenses <- as.numeric(df_alt$expenses)
df_alt$smoker_numeric <- as.numeric(df_alt$smoker == "yes")
df_alt$sex_numeric <- as.numeric(df_alt$sex == "yes")
df_alt$reg_numeric <- as.numeric(df_alt$region == "yes")

df_alt$age     <- log(df_alt$age)
df_alt$bmi <- log(df_alt$bmi)
df_alt$smoker_numeric <- log(df_alt$smoker_numeric) 
df_alt$expenses <- log(df_alt$expenses)
df_alt$children <- log(df_alt$children + 0.01)
df_alt$sex_numeric <- log(df_alt$sex_numeric) 
df_alt$reg_numeric     <- log(df_alt$reg_numeric)

summary(df_alt)

##       age             bmi        smoker        expenses         children       
##  Min.   :2.890   Min.   :2.773   no :1064   Min.   : 7.023   Min.   :-4.60517  
##  1st Qu.:3.296   1st Qu.:3.270   yes: 274   1st Qu.: 8.464   1st Qu.:-4.60517  
##  Median :3.664   Median :3.414              Median : 9.147   Median : 0.00995  
##  Mean   :3.597   Mean   :3.403              Mean   : 9.099   Mean   :-1.67105  
##  3rd Qu.:3.932   3rd Qu.:3.547              3rd Qu.: 9.720   3rd Qu.: 0.69813  
##  Max.   :4.159   Max.   :3.972              Max.   :11.063   Max.   : 1.61144  
##      sex            region    smoker_numeric  sex_numeric    reg_numeric  
##  female:662   northeast:324   Min.   :-Inf   Min.   :-Inf   Min.   :-Inf  
##  male  :676   northwest:325   1st Qu.:-Inf   1st Qu.:-Inf   1st Qu.:-Inf  
##               southeast:364   Median :-Inf   Median :-Inf   Median :-Inf  
##               southwest:325   Mean   :-Inf   Mean   :-Inf   Mean   :-Inf  
##                               3rd Qu.:-Inf   3rd Qu.:-Inf   3rd Qu.:-Inf  
##                               Max.   :   0   Max.   :-Inf   Max.   :-Inf

set.seed(123) # What is set.seed()? We want to make sure that we get the same results for randomization each time you run the script.   
cv_data   <- createDataPartition(y = df$expenses, p=0.7, list=F)
cv_train = df_alt[cv_data, ]
cv_test = df_alt[-cv_data, ]

# define explanatory variables (X's) and dependent variable (Y) in training set
train_x = data.matrix(cv_train[, -4])
train_y = cv_train[,4]

# define explanatory variables (X's) and dependent variable (Y) in testing set
test_x = data.matrix(cv_test[, -4])
test_y = cv_test[, 4]

any(is.infinite(test_y))

## [1] FALSE

train_x[is.infinite(train_x)] <- NA
test_x[is.infinite(test_x)] <- NA

# Calcular la media de cada columna
means <- colMeans(train_x, na.rm = TRUE)

# Imputar los valores NA con la media de cada columna
for (i in 1:ncol(train_x)) {
  train_x[is.na(train_x[, i]), i] <- means[i]
}

# Calcular la media de cada columna
means2 <- colMeans(test_x, na.rm = TRUE)

# Imputar los valores NA con la media de cada columna
for (i in 1:ncol(test_x)) {
  test_x[is.na(test_x[, i]), i] <- means2[i]
}

any(is.na(train_x))

## [1] TRUE

any(is.infinite(train_x))

## [1] FALSE

any(is.na(train_y))

## [1] FALSE

any(is.infinite(train_y))

## [1] FALSE

any(is.na(test_x))

## [1] TRUE

any(is.infinite(test_x))

## [1] FALSE

any(is.na(test_y))

## [1] FALSE

# 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) # 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

# Looks like the lowest RMSE for both training and test dataset is achieved at 59 round. 
# Lets estimate our final regression model
reg_xgb = xgboost(data = xgb_train, max.depth = 3, nrounds = 59, verbose = 0) # setting verbose = 0 avoids to display the training and testing error for each round. 
prediction_xgb_test<-predict(reg_xgb, xgb_test)
rmse(prediction_xgb_test, cv_test$expenses)

## [1] 0.4063379

# Lets do some diagnostic check of regression residuals 
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)

# 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")

# Calcular los residuos
residuals <- prediction_xgb_test - cv_test$expenses

# Calcular la suma de los cuadrados de los residuos
SS_residual <- sum(residuals^2)

# Calcular la suma total de los cuadrados
SS_total <- sum((cv_test$expenses - mean(cv_test$expenses))^2)

# Calcular R cuadrada
R_squared <- 1 - (SS_residual / SS_total)

# Mostrar el valor de R cuadrada
print(R_squared)

## [1] 0.7948389

E.Decision Trees

decision_tree_model <- rpart(log(expenses) ~ log(bmi) + log(age) + log(children + 0.01) + 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)

# Hacer predicciones con el modelo de árbol de decisión en los datos de entrenamiento
decision_tree_predictions <- predict(decision_tree_model, train)

# Calcular los residuos
residuals <- decision_tree_predictions - log(train$expenses)

# Calcular el RMSE
rmsedt <- sqrt(mean(residuals^2))

# Mostrar el valor del RMSE
print(rmsedt) #0.4010

## [1] 0.4010215

# Calcular la suma de los cuadrados de los residuos
SSE <- sum(residuals^2)

# Calcular la suma total de los cuadrados
SST <- sum((log(train$expenses) - mean(log(train$expenses)))^2)

# Calcular R cuadrada
R_squared <- 1 - (SSE / SST)

# Mostrar el valor de R cuadrada
print(R_squared)

## [1] 0.8102677

F.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: 21873431
##                     % Var explained: 85.31

# Prediction & Confusion Matrix – test data
rf_prediction <- predict(rf_model,cv_test)

# confusionMatrix(rf_prediction_train_data, train$MEDV) #Una matriz de confusión es esencialmente una tabla que clasifica las predicciones con respecto a los valores reales.
RMSE_rf <- rmse(rf_prediction, cv_test$expenses)
RMSE_rf #8521.414

## [1] 8303.481

#Evaluar la importancia de las variables
#¿Cómo interpretar varImpPlot()? Cuanto mayor sea el valor de la precisión de disminución media, mayor será la importancia de la variable en el modelo.
# En otras palabras, la disminución media de la precisión representa cuánto reduce la precisión del modelo la eliminación de cada variable.
varImpPlot(rf_model, n.var = 5, main = "Top 5 - Variable") # Muestra un gráfico de importancia de las variables del modelo de bosque aleatorio.

importance(rf_model)                                        # Vale la pena mencionar que IncNodePurity aumenta el error del modelo al eliminar cada una de las variables explicativas especificadas.

##          IncNodePurity
## age        19733986905
## bmi        21699766084
## children    3102584526
## smoker     98443051346
## sex         1307424273
## region      3137177082

                                                            # Brevemente, varImpPlot() indica la importancia de cada variable para explicar el desempeño de la variable dependiente (Y).

G.Neural Networks Regresión

# Lets estimate a Neural Network Regression
summary(train)

##       age            sex           bmi           children     smoker   
##  Min.   :18.00   female:530   Min.   :16.00   Min.   :0.000   no :850  
##  1st Qu.:27.00   male  :542   1st Qu.:26.40   1st Qu.:0.000   yes:222  
##  Median :39.00                Median :30.50   Median :1.000            
##  Mean   :39.19                Mean   :30.76   Mean   :1.081            
##  3rd Qu.:51.00                3rd Qu.:34.80   3rd Qu.:2.000            
##  Max.   :64.00                Max.   :53.10   Max.   :5.000            
##        region       expenses    
##  northeast:254   Min.   : 1122  
##  northwest:257   1st Qu.: 4745  
##  southeast:309   Median : 9382  
##  southwest:252   Mean   :13318  
##                  3rd Qu.:16604  
##                  Max.   :63770

train$bmi <- as.numeric(train$bmi)
train$expenses <- as.numeric(train$expenses)
train$smoker_numeric <- as.numeric(train$smoker == "yes")
train$sex_numeric <- as.numeric(train$sex == "yes")
train$reg_numeric <- as.numeric(train$region == "yes")

nn_model <- neuralnet(expenses ~ age + sex_numeric + bmi + children + reg_numeric + smoker_numeric, data = train, hidden = c(5, 3), linear.output = TRUE) 

# Plot the neural network 
plot(nn_model)

# Hacer predicciones con el modelo de red neuronal en los datos de entrenamiento
nn_predictions <- predict(nn_model, train)

# Calcular los residuos
residuals <- nn_predictions - train$expenses

# Calcular el RMSE
rmsenn <- sqrt(mean(residuals^2))

# Lets do some diagnostic check of regression residuals 
residuals<-train$expenses - nn_predictions
plot(residuals, xlab= "Dependent Variable", ylab = "Residuals", main = 'Neural Net Regression Residuals')
abline(0,0)

# Mostrar el valor del RMSE
print(rmsenn) #12201.58

## [1] 12201.58

Principales hallagazos encontrados durante el EDA:

A travez del EDA pudimos observar:

• Al examinar los gráficos de dispersión, notamos que existe una tendencia clara: a medida que la edad aumenta, también lo hacen los gastos médicos.

• Observamos que las personas que fuman y tienen un índice de masa corporal elevado tienden a registrar gastos médicos más altos.

• Observamos que la base de datos no presenta NA’s.

• Podemos observar que el gasto entre hombres y mujeres es casi el mismo porque sus medias se encuentran igual pero el grupo de los homres presenta una dispersion mayor

• Podemos observar que la region southeast es la que genera mayor expenses,a traves del analisis del boxplot

• Ademas que tambien podemos obversar que la region southeast es la que genera un BMI mas alto lo que puede ser un factor importante

RMSE de los diferentes modelos

# Datos 
model_names <- c("OLS", "OLS Log", "XGBoost", "DT", "RF", "NN")  # Nombres de los modelos
rmse_values <- c(0.7475, 7517,0.4063 , 0.4010, 8521.414, 12201.58)  # Valores de RMSE correspondientes a cada modelo

# Crear un dataframe con los datos
rmse_data <- data.frame(Modelo = model_names, RMSE = rmse_values)

# Crear el gráfico de barras
barplot(rmse_data$RMSE, names.arg = rmse_data$Modelo, col = "skyblue", 
        main = "RMSE de Modelos Estimados", xlab = "Modelo", ylab = "RMSE")

Seleccion del Modelo

El modelo que mejor se ajusta según el RMSE es el modelo de Decision Tree (DT) ya que tiene el valor de RMSE más bajo (0.4010).ya que un RMSE más bajo indica un mejor ajuste del modelo a los datos.otro modelo que pudiera ser buena opcion seria el XGBoost, ya que varia muy poco la diferencia entre sus RMSE, ademas de que DT muestra una R cuadrada mas alta lo que nos indica que tiene un ajuste mas adecuado que XGBoost.

• R cuadrada XGBoost: 0.7948.
• R cuadrada DT: 8102.