DATA SET: “COSTO SALUD”

El dataset fue extraido de Kaggle: https://www.kaggle.com/datasets/mirichoi0218/insurance

COLUMNAS

edad: Edad del beneficiario principal
sexo: Género femenino o masculino
bmi: Índice de masa corporal, que proporciona una comprensión del peso relativamente alto o bajo en relación con la altura, idealmente 18,5 a 24,9.
niños: Número de niños cubiertos por el seguro de salud / Número de dependientes
fumador: Fuma o no fuma
región: el área residencial del beneficiario en los EE. UU., noreste, sureste, suroeste, noroeste.
cargos: costos médicos individuales facturados por el seguro de salud

Cargar Librerías

library('GGally')

## Loading required package: ggplot2

## Registered S3 method overwritten by 'GGally':
##   method from   
##   +.gg   ggplot2

library('dplyr')

## 
## 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('lmtest')

## Loading required package: zoo

## 
## Attaching package: 'zoo'

## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric

Información de los Datos

df <- read.csv('insurance.csv')
head(df)

##   age    sex    bmi children smoker    region   charges
## 1  19 female 27.900        0    yes southwest 16884.924
## 2  18   male 33.770        1     no southeast  1725.552
## 3  28   male 33.000        3     no southeast  4449.462
## 4  33   male 22.705        0     no northwest 21984.471
## 5  32   male 28.880        0     no northwest  3866.855
## 6  31 female 25.740        0     no southeast  3756.622

A continuación vemos los tipos de datos y el tamaño del dataframe :

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 ...
##  $ children: int  0 1 3 0 0 0 1 3 2 0 ...
##  $ smoker  : chr  "yes" "no" "no" "no" ...
##  $ region  : chr  "southwest" "southeast" "southeast" "northwest" ...
##  $ charges : num  16885 1726 4449 21984 3867 ...

summary(df)

##       age            sex                 bmi           children    
##  Min.   :18.00   Length:1338        Min.   :15.96   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.66   Mean   :1.095  
##  3rd Qu.:51.00                      3rd Qu.:34.69   3rd Qu.:2.000  
##  Max.   :64.00                      Max.   :53.13   Max.   :5.000  
##     smoker             region             charges     
##  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

Gráfica de Dispersión

Utilizando la función nativa plot de R podemos visualizar un diagrama de dispersión para todas las variables del dataframe:

plot(df)

Del anterior gráfico podemos ver que hay una relación lineal subdividida en tres grupos para las variables edad y costo de salud, si pudieramos determinar porque estan ocurriendo estos grupos podríamos crear una relación lineal entre estas 2 variables para los subgrupos.

Visualizamos la relación lineal y los 3 subgrupos que aparecen entre la EDAD y el COSTO

plot(x = df$age, y=df$charges, xlab = 'AGE', ylab = 'CHARGES', main = 'charges vs age')

Visualizamos la distribución, la dispersión y los coeficientes de relación de los datos por colores según el sexo de la población

ggpairs(df, columns=1:7, ggplot2::aes(colour=sex), progress = FALSE)

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Podemos ver que tenemos unos coeficientes de relación de 0.2 y además vemos que la gráfica de edad vs costo no puede ser explicada por 1 solo modelo lineal.

Distribución, Dispersión y Coeficientes de relación de los datos por colores según ‘fumadores’ en la población

ggpairs(df, columns=1:7, ggplot2::aes(colour=smoker), progress = FALSE)

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Relación lineal aparente de los 2 subgrupos que aparecen entre BMI y el COSTO

plot(x = df$bmi, y=df$charges,  col = factor(df$smoker), main = 'Fumadores vs No fumadores', xlab = 'BMI', ylab = 'CHARGES')

Aquí podemos observar varias cosas útiles: - Existe una correlación muy alta para los ‘smoker=yes’ de 0.8 con respecto al ‘bmi’ - Existe una correlación de 0.6 para los ‘smokers=no’ con respecto a la edad

ggpairs(df, columns=1:7, ggplot2::aes(colour=as.factor(region)), progress = FALSE)

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

No encontramos mejores relaciones visualizando los datos según la ‘región’. No pudimos encontrar la razón de los 3 subgrupos del diagrama de dispersión edad vs costo así que decidímos realizar el análisis para los NO FUMADORES y aparte para los SI FUMADORES, de esta manera podremos explicar la variable ‘costo’ utilizando la variable más correlacionada en cada caso, que como se observo anteriormete es la edad para los no fumadores y el indice bmi para los fumadores.

Creamos un dataset para los fumadores y no fumadores

non<- select(filter(df, smoker == 'no'),c(age,sex,bmi, children, smoker, region, charges))
yes<- select(filter(df, smoker == 'yes'),c(age,sex,bmi, children, smoker, region, charges))

Observamos correlación de 0.6 entre ‘edad’ vs ‘costo’ para los NO FUMADORES

ggpairs(non, columns=1:7, ggplot2::aes(colour=as.factor(sex)), progress = FALSE)

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Observamos correlación de 0.8 entre ‘bmi’ vs ‘costo’ para los FUMADORES

ggpairs(yes, columns=1:7, ggplot2::aes(colour=as.factor(sex)), progress = FALSE)

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

MODELOS REGRESIÓN LINEAL

CREANDO LOS MODELOS DE REGRESIÓN LINEAL : 1 PARA FUMADORES Y OTRO PARA NO FUMADORES

Los siguientes modelos buscan explicar el ‘costo de salud’ para personas fumadoras y no fumadoras, ya que se encontro que esta es la variable más correlacionada con el ‘costo’. Además se hará un modelo de regresión lineal para los no fumadores basandose en la ‘edad’ , mientras que para los fumadores utilizaremos el índece de ‘bmi’ de las personas para predecir el costo de salud.

NO FUMADORES: COSTO vs EDAD

lm1 <- lm(non$charges~non$age)
summary(lm1)

## 
## Call:
## lm(formula = non$charges ~ non$age)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3182.9 -1948.6 -1363.8  -665.2 24470.7 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -2091.42     425.10   -4.92    1e-06 ***
## non$age       267.25      10.16   26.29   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4667 on 1062 degrees of freedom
## Multiple R-squared:  0.3943, Adjusted R-squared:  0.3937 
## F-statistic: 691.4 on 1 and 1062 DF,  p-value: < 2.2e-16

Verificamos los supuestos del modelo: Normalidad de residuiales y Homogeneidad de varianza

Test de Normalidad sobre los residuales:

Ho: Los errores son Normales Ha: Los errores No son Normales p-value <0.05 se rechaza la Ho

l1res<- residuals(lm1)
shapiro.test(l1res)

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

Los residuales NO pasa la prueba de NORMALIDAD para este modelo

Quitamos el intercepto del modelo

lm1.1 <- lm(non$charges~0 +non$age)
summary(lm1.1)

## 
## Call:
## lm(formula = non$charges ~ 0 + non$age)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3564.8 -2385.5 -1705.0  -571.1 23921.0 
## 
## Coefficients:
##         Estimate Std. Error t value Pr(>|t|)    
## non$age  220.162      3.458   63.67   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4718 on 1063 degrees of freedom
## Multiple R-squared:  0.7923, Adjusted R-squared:  0.7921 
## F-statistic:  4054 on 1 and 1063 DF,  p-value: < 2.2e-16

Verificamos supuestos

Homogeneidad de varianza

Breush-Pagan test: Ho: Las varianzas son homogeneas Ha: LAs varianzas no son homogeneas Necesitamos un p-value mayor a 0.05 para No rechazar la hipotesis nula

bptest(lm1)

## 
##  studentized Breusch-Pagan test
## 
## data:  lm1
## BP = 0.07064, df = 1, p-value = 0.7904

Normalidad Residuales

Test de Normalidad sobre los residuales: Ho: Los errores son Normales Ha: Los errores No son Normales p-value <0.05 se rechaza la Ho

l1res.1<- residuals(lm1.1)
shapiro.test(l1res.1)

## 
##  Shapiro-Wilk normality test
## 
## data:  l1res.1
## W = 0.53788, p-value < 2.2e-16

Por último vemos que los residuales son significativamente grandes, del orden de $1.600, esto se explica gráficamente.

plot(non$age, non$charges,main = 'Modelo SIN intercepto para NO FUMADORES', xlab = 'AGE', ylab = 'CHARGES')+
abline( 0, 267.25)

## integer(0)

plot(non$age, non$charges,main = 'Modelo CON intercepto para NO FUMADORES', xlab = 'AGE', ylab = 'CHARGES')+
abline( -2091.42, 267.25)

## integer(0)

A pesar de que en el modelo sin intercepto el R-cuadrado se eleva de 0.4 a 0.8, desafortunadamente los residuales incrementan $400 aún más por lo que no es tan conveniente en este caso, gráficamente se puede ver porque es preferible el modelo con intercepto.

La razón es que el modelo intenta ajustar TODOS LOS PUNTOS de la dispersión, y los valores muy alejados de la linea estan actuando como outlayers de los datos que tienen explicitamente una relación lineal.

RESUMEN DEL MODELO:

Y = 220.162xAGE - 2091.42

Este modelo sin intercepto el R-cuadrado ha subido de 0.4 a 0.8, sin embargo el error es mayor.

FUMADORES: COSTO vs BMI

lm2 <- lm(yes$charges~yes$bmi)
summary(lm2)

## 
## Call:
## lm(formula = yes$charges ~ yes$bmi)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -19768.0  -4487.9     34.4   3263.9  31055.9 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -13186.58    2052.88  -6.423 5.93e-10 ***
## yes$bmi       1473.11      65.48  22.496  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6837 on 272 degrees of freedom
## Multiple R-squared:  0.6504, Adjusted R-squared:  0.6491 
## F-statistic: 506.1 on 1 and 272 DF,  p-value: < 2.2e-16

Quitando el intercepto:

lm2.1 <- lm(yes$charges~ 0 + yes$bmi)
summary(lm2.1)

## 
## Call:
## lm(formula = yes$charges ~ 0 + yes$bmi)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -13492.8  -5983.4   -560.5   3936.2  30378.6 
## 
## Coefficients:
##         Estimate Std. Error t value Pr(>|t|)    
## yes$bmi  1061.08      14.11   75.19   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 7323 on 273 degrees of freedom
## Multiple R-squared:  0.9539, Adjusted R-squared:  0.9538 
## F-statistic:  5653 on 1 and 273 DF,  p-value: < 2.2e-16

Obtenemos un R-cuadrado del 95%¡

RESUMEN DEL MODELO:

Y = 1061.08 x BMI

plot(yes$bmi, yes$charges, main = "Modelo CON intercepto para FUMADORES", xlab = 'BMI', ylab = 'CHARGES')+
abline( -13186.58,  1473.11)

## integer(0)

plot(yes$bmi, yes$charges,main = "Modelo SIN intercepto para FUMADORES", xlab = 'BMI', ylab = 'CHARGES')+
abline( 0, 1061.08)

## integer(0)

HACEMOS VERIFICACIÓN DE LOS SUPUESTOS

l2res<- residuals(lm2.1)
shapiro.test(l2res)

## 
##  Shapiro-Wilk normality test
## 
## data:  l2res
## W = 0.97569, p-value = 0.0001279

bptest(lm2)

## 
##  studentized Breusch-Pagan test
## 
## data:  lm2
## BP = 2.3496, df = 1, p-value = 0.1253

Podemos decir para los SI FUMADORES que en el modelo con intercepto tiene uno residuales muy pequeños y adecuados, sin embargo el intercepto es negativo y NO EXPLICA la realidad de que existan cobros negativos, sin embargo, para el modelo sin intercepto los residuales son mayores.

ANÁLISIS DE VARIANZAS PARA LOS MODELOS EN GENERAL

summary.aov(lm1)

##               Df    Sum Sq   Mean Sq F value Pr(>F)    
## non$age        1 1.506e+10 1.506e+10   691.4 <2e-16 ***
## Residuals   1062 2.313e+10 2.178e+07                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

summary.aov(lm1.1)

##             Df    Sum Sq   Mean Sq F value Pr(>F)    
## non$age      1 9.022e+10 9.022e+10    4054 <2e-16 ***
## Residuals 1063 2.366e+10 2.226e+07                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

summary.aov(lm2)

##              Df    Sum Sq   Mean Sq F value Pr(>F)    
## yes$bmi       1 2.365e+10 2.365e+10   506.1 <2e-16 ***
## Residuals   272 1.271e+10 4.674e+07                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

summary.aov(lm2.1)

##            Df    Sum Sq   Mean Sq F value Pr(>F)    
## yes$bmi     1 3.032e+11 3.032e+11    5653 <2e-16 ***
## Residuals 273 1.464e+10 5.363e+07                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

CONCLUSIÓN

plot(yes$bmi,yes$charges,type="p",col="red", ylim = c(0, 60000), xlab="Edad y BMI ", ylab="Costo salud", main = 'Fumadores y No fumadores')+
lines(non$age, non$charges,type="p",col="green")+
  abline( -2091.42, 267.25)+
   abline( 0, 1061)

## integer(0)

Observando las relaciones de la variable de interés con las otras variables en los diagramas y analizando su comportamiento con respecto a las variables categóricas, pudímos darnos cuenta de relaciones lineales importantes que nos permiten análizar por separado el problema de regresión lineal de la variable ‘charges’.

El problema ahora sería verificar por qué los residuales de los modelos NO estan cumpliendo el supuesto de NORMALIDAD, tal vez haciendo una transformación adecuada de los datos para solucionarlo como una estandarización. Además de que tenemos el gran problema de los residuales, si observamos en cada modelo vemos que son diferentes de 0.

UNA SOLUCIÓN MÁS GENERAL sería sacar los ‘outlayers’ del modelo de AGE vs CHARGE con el fin de, obtener un mejor ajuste de los datos, que gráficamente se ven en la mayoría, suponiendo los otros datos como atípicos.

Regresión Lineal Simple

willandru

2022-09-04