Introducción

En esta unidad se revisarán las funciones para la estimación de modelos lineales y algunos modelos generalizados. Para esto se trabajá con la librería stats y una primera base de datos que tiene datos simulados sobre: porcentaje de personas con enfermedad cardiaca (heart.disease), porcentaje de personas que usa la bicicleta como medio de transporte (biking) y porcentaje de personas que fuman (smoking) en una muestra de 498 ciudades independientes.

library(ggplot2)
library(dplyr)
library(stats)

Se exporta la base de datos que está en .csv

library(readr)
base_heart <- read_csv("~/Documents/Delia Ortega 2023/disco_d/DISCOE/CLASES/JAVERIANA_SALUD/REGISTRO_MEPI/CURSO_MÉTODOS_INFORMÁTICA/Material/Unidad2/heart_data.csv")

1. Modelo de Regresión Lineal

2. Aplicación del modelo de Regresión Lineal

Inicalmente se realizará un análisis exploratorio descriptivo de las variables en la base de datos.

summary(base_heart)
##        id            biking          smoking        heart.disease    
##  Min.   :  1.0   Min.   : 1.119   Min.   : 0.5259   Min.   : 0.5519  
##  1st Qu.:125.2   1st Qu.:20.205   1st Qu.: 8.2798   1st Qu.: 6.5137  
##  Median :249.5   Median :35.824   Median :15.8146   Median :10.3853  
##  Mean   :249.5   Mean   :37.788   Mean   :15.4350   Mean   :10.1745  
##  3rd Qu.:373.8   3rd Qu.:57.853   3rd Qu.:22.5689   3rd Qu.:13.7240  
##  Max.   :498.0   Max.   :74.907   Max.   :29.9467   Max.   :20.4535

Para observar la distribución de la variable respuesta: porcentaje de personas con enfermedad cardiaca, se realizará un histograma:

hist(base_heart$heart.disease, main="Histograma % personas con enfermedad cardíaca")

Para estimar la correlación entre las variables:

cor(base_heart$heart.disease, base_heart$biking, method = c("pearson"))
## [1] -0.9354555
# en método se puede seleccionar: pearson, kendall o spearman.
cor(base_heart, method = c("pearson"))
##                        id      biking    smoking heart.disease
## id             1.00000000  0.05708763 0.05267393   -0.05172509
## biking         0.05708763  1.00000000 0.01513618   -0.93545547
## smoking        0.05267393  0.01513618 1.00000000    0.30913098
## heart.disease -0.05172509 -0.93545547 0.30913098    1.00000000

También se pueden construir métodos gráficos para observar la relación entre variables de la base de datos.

#install.packages("PerformanceAnalytics")
library("PerformanceAnalytics")
chart.Correlation(base_heart, histogram=TRUE, pch=19)

Los asteriscos es el gráfico anterior muestran la significancia de la correlación a través de valores p que están asociados con los simbolos: p-values(0, 0.001, 0.01, 0.05, 0.1, 1) <=> simbolos(“”, “”, “”, “.”, ” “).

También se pueden realizar correlogramas usando la librería \(corrplot\).

#install.packages("corrplot")
library(corrplot)
## corrplot 0.92 loaded
cor=cor(base_heart, method = c("pearson"))
corrplot(cor, type = "upper", order = "hclust", 
         tl.col = "black", tl.srt = 45)

Aunque a través de la función chart.correlation se puede observar gráficamente la relación entre variables, también se puede hacer uso de la función \(plot\)

par(mfrow=c(1,2))  #permite tener una ventana con una fila y dos columnas para los gráficos.
plot(heart.disease ~ biking, data=base_heart)
plot(heart.disease ~ smoking, data=base_heart)

Para la estimación del modelo lineal se usa la función \(lm\). Ajustaremos el siguiente modelo múltiple: Heart_disease = \(\beta_0\) + \(\beta_1\) biking + \(\beta_2\) smoking

modelo<-lm(heart.disease ~ biking + smoking, data = base_heart)
summary(modelo)
## 
## Call:
## lm(formula = heart.disease ~ biking + smoking, data = base_heart)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.1789 -0.4463  0.0362  0.4422  1.9331 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 14.984658   0.080137  186.99   <2e-16 ***
## biking      -0.200133   0.001366 -146.53   <2e-16 ***
## smoking      0.178334   0.003539   50.39   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.654 on 495 degrees of freedom
## Multiple R-squared:  0.9796, Adjusted R-squared:  0.9795 
## F-statistic: 1.19e+04 on 2 and 495 DF,  p-value: < 2.2e-16

La salida anterior muestra un análisis descriptivo de los residuales del modelo, los coeficientes beta´s, el error estándar, la prueba t para cada coediciente y los valores p. Además la estimación de la prueba de bondad y ajuste del modelo \(R^2\) y la prueba F que plantea como hipótesis nula que todos los beta´s son iguales a cero.

Para el AIC y BIC del modelo:

AIC(modelo)
## [1] 995.3534
BIC(modelo)
## [1] 1012.196

3. Validación de supuestos

Para este modelo de regresión lineal se deben validar supuestos respecto a los residuales como la normalidad, homocedasticidad de la varianza e independencia.

par(mfrow=c(2,2))
plot(modelo)

El gráfico anterior permite evaluar la tendencia de residuales, la normalidad (QQ plot) y posibles valores atípicos.

Para determinar homocedasticidad de la varianza, se puede aplicar la prueba analítica Breush-Pagan (Hipótesis nula: residuales distribuidos con igual varianza).

library(lmtest)
lmtest::bptest(modelo)
## 
##  studentized Breusch-Pagan test
## 
## data:  modelo
## BP = 5.7775, df = 2, p-value = 0.05564

Para la normalidad de los residuales (Hipótesis nula:normalidad).

residuales=resid(modelo)
ks.test(residuales, 'pnorm')
## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  residuales
## D = 0.11967, p-value = 1.278e-06
## alternative hypothesis: two-sided
shapiro.test(residuales)
## 
##  Shapiro-Wilk normality test
## 
## data:  residuales
## W = 0.997, p-value = 0.4935