Correlacion y Regresion
Mario Beltran, Liyizabeth Cardoso, Angela Daniela Pastez, Ana Maria Daza
2025-05-16
INTRODUCCION
En este estudio se analiza el comportamiento de los ingresos netos de los hogares en funcion de ciertos factores economicos como los gastos familiares y el estrato socioeconomico.
PROBLEMA
Se desea identificar qué factores explican significativamente el nivel de ingreso neto de los hogares, evaluando el efecto de los gastos del hogar y del estrato socioeconómico.
OBJETIVO
Aplicar un modelo de regresión lineal múltiple para predecir el
ingreso neto (x) de los hogares a partir de variables como
los gastos familiares sin arriendo (y), gastos de arriendo
(z) y el estrato (u).
MODELO DE REGRESION LINEAL SIMPLE
Introducción
Análisis de regresión lineal simple usando datos sobre Ingreso Neto y Gasto Familiar (sinarriendo). Se incluyen las fórmulas utilizadas en el modelo y su implementación en R.
El modelo es:
\[ Y_i = \beta_0 + \beta_1 X_i + \varepsilon_i \]
Donde:
* \(X_i\): Ingreso Neto
i
* \(Y_i\): Gasto Familiar (sinarriendo)
i
* \(\beta_0\): intercepto (valor
esperado si X = 0)
* \(\beta_1\): pendiente (cuánto varía
Y por cada valor que varía X)
* \(\varepsilon_i\): error
aleatorio
Librerias
# Cargar las librerías necesaria
library(readxl)## Warning: package 'readxl' was built under R version 4.4.3library(lmtest)## Warning: package 'lmtest' was built under R version 4.4.3## Cargando paquete requerido: zoo## Warning: package 'zoo' was built under R version 4.4.2##
## Adjuntando el paquete: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numericlibrary(knitr)## Warning: package 'knitr' was built under R version 4.4.3Base de datos
# Importar la base de datos
Base <- read_excel("C:/Users/Usuario/Desktop/Base de datos Alcaldia.xlsx")Primeros datos
# Verificar los primeros datos
head(Base)## # A tibble: 6 × 8
## `>` Estrato `Ingreso Neto` `No de Personas a Cargo` Gastos Familar (sinarr…¹
## <dbl> <dbl> <dbl> <chr> <dbl>
## 1 1 1 895212 Ninguna 367037
## 2 2 1 822988 Ninguna 288046
## 3 3 1 812288 Ninguna 324915
## 4 4 1 627853 Ninguna 282534
## 5 5 1 673359 Ninguna 282811
## 6 6 1 696733 Ninguna 320497
## # ℹ abbreviated name: ¹`Gastos Familar (sinarriendo)`
## # ℹ 3 more variables: `Gastos de Arriendo` <dbl>, `Gastos / Ingresos` <dbl>,
## # `Arriendos/ Ingreso` <dbl>#Si desea, se renombra las variables
## En lugar de trabajar con nombres largos como "Ingreso Neto" y "Gastos Familiar", los renombramos simplemente como x e y, respectivamente.
names(Base)[names(Base) == "Ingreso Neto"] <- "x"
names(Base)[names(Base) == "Gastos Familar (sinarriendo)"] <- "y"
names(Base)[names(Base) == "Gastos de Arriendo"] <- "z"
names(Base)[names(Base) == "Estrato"] <- "u"
x <- Base$x
y <- Base$y
z <- Base$z
u <- Base$u Cálculo del modelo paso a paso
Cálculo de la pendiente \(\hat{\beta}_1\)
Fórmula:
\[ \hat{\beta}_1 = \frac{\sum (x_i - \bar{x})(y_i - \bar{y})}{\sum (x_i - \bar{x})^2} \]
¿Qué significa?
- \(x_i\): Ingresos Neto i
- \(y_i\): Gastos Familar (sinarriendo) i
- \(\bar{x}\): promedio de X (pesos)
- \(\bar{y}\): promedio de Y
(pesos)
- Numerador: covarianza entre X y Y
- Denominador: varianza de X
¿Por qué se usa así? Porque queremos encontrar la pendiente de la recta que mejor ajusta todos los puntos, minimizando los errores.
x <- Base$x
y <- Base$y
media_x <- mean(x)
media_y <- mean(y)
numerador <- sum((x - media_x) * (y - media_y))
denominador <- sum((x - media_x)^2)
beta_1 <- numerador / denominador
beta_1## [1] 0.3263488Interpretación:Por cada peso adicional en gasto familiar sin arriendo, se espera que el ingreso neto aumente 0.326 pesos, en promedio.
El valor indica que a mayor gasto familiar sin arriendo, mayor sera el ingreso ento. Por cada peso adicional que un hogar gasta en consumo familiar(sin incuir arriendo), se espera que su ingreso neto aumente en promedio 0.3263488
Cálculo del intercepto \(\hat{\beta}_0\)
Fórmula:
\[ \hat{\beta}_0 = \bar{y} - \hat{\beta}_1 \bar{x} \]
¿Qué significa?
Es el gasto (X = 0).
beta_0 <- media_y - beta_1 * media_x
beta_0## [1] 133027.9Interpretación: 1.3302792^{5} este valor representa el Ingreso Neto esperado cuando el Gasto Familiar (sin arriendo) es cero.
Modelo Ajustado
# Ajustar el modelo de regresión
## Ajustar un modelo de regresión lineal simple, donde tratamos de explicar el gasto familiar (y) en función del ingreso neto (x).
modelo <- lm(y ~ x, data = Base)
## Usar summary(modelo) para obtener un resumen con los coeficientes, el error estándar, los valores p y el R², entre otros datos.
summary(modelo)##
## Call:
## lm(formula = y ~ x, data = Base)
##
## Residuals:
## Min 1Q Median 3Q Max
## -274876 -48440 -9871 38139 341567
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.330e+05 2.361e+04 5.634 8.62e-08 ***
## x 3.263e-01 1.998e-02 16.331 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 92780 on 148 degrees of freedom
## Multiple R-squared: 0.6431, Adjusted R-squared: 0.6407
## F-statistic: 266.7 on 1 and 148 DF, p-value: < 2.2e-16Ecuación final del modelo:
\[ \hat{Y}_i = 1.3302792\times 10^{5} + 0.326 \cdot X_i \]
Predicción de valores y cálculo de residuos
Fórmula de predicción:
\[ \hat{Y}_i = \hat{\beta}_0 + \hat{\beta}_1 x_i \]
¿Qué significa?
- \(\hat{Y}_i\): Es el valor predicho o estimado del Ingreso Neto estimado para la persona i, calculado con el modelo.
- \(x_i\): Se interpreta como el Gasto Familiar (sinarriendo) real de la persona i.
- Esta predicción no incluye el error \(\varepsilon_i\), por lo que no es necesariamente exacta.
Base$predicciones <- beta_0 + beta_1 * x
head(Base[, c("x", "y", "predicciones")])## # A tibble: 6 × 3
## x y predicciones
## <dbl> <dbl> <dbl>
## 1 895212 367037 425179.
## 2 822988 288046 401609.
## 3 812288 324915 398117.
## 4 627853 282534 337927.
## 5 673359 282811 352778.
## 6 696733 320497 360406.Ejemplo: Si una familia reporta un Gasto Familiar (sin arriendo) de 5 millones de pesos, el modelo predice que su Ingreso Neto estimado es de 1.3302955^{5} mil pesos.
Fórmula del residuo:
\[ e_i = y_i - \hat{Y}_i \]
¿Qué significa?
- \(e_i\): es el residuo o error de predicción para la persona i.
- Es la diferencia entre la Ingreso observado y la Ingreso estimado por el modelo.
- Si \(e_i > 0\): el modelo subestimó el Ingreso.
- Si \(e_i < 0\): el modelo sobreestimó el Ingreso.
Base$predicciones <- beta_0 + beta_1 * x
Base$residuos <- y - Base$predicciones
head(Base[, c("x", "y", "predicciones", "residuos")])## # A tibble: 6 × 4
## x y predicciones residuos
## <dbl> <dbl> <dbl> <dbl>
## 1 895212 367037 425179. -58142.
## 2 822988 288046 401609. -113563.
## 3 812288 324915 398117. -73202.
## 4 627853 282534 337927. -55393.
## 5 673359 282811 352778. -69967.
## 6 696733 320497 360406. -39909.Interpretación: Los residuos son los errores del modelo: cuánto se desvían las predicciones de los valores reales.
¿Por qué son importantes?
- Los residuos nos dicen qué tan bien el modelo ajusta cada punto individual.
- En un buen modelo, los residuos deben ser pequeños y no mostrar patrones sistemáticos.
- Analizar los residuos permite detectar:
- Outliers (valores atípicos)
- Sesgos en la predicción
- Violaciones de supuestos como homocedasticidad o linealidad
1. Error Cuadrático Medio (ECM)
Fórmula:
\[ ECM = \frac{1}{n} \sum (y_i - \hat{Y}_i)^2 \]
¿Qué significa?
- Mide el promedio de los errores al cuadrado (es decir, de los residuos).
- Penaliza más fuertemente los errores grandes.
- Es útil para comparar modelos: cuanto más bajo el ECM, mejor el ajuste.
Elementos de la fórmula:
- \(y_i\): valor real observado (nota)
- \(\hat{Y}_i\): valor estimado por el modelo
- \(n\): número total de observaciones
ECM <- mean(Base$residuos^2)
ECM## [1] 8493299666Interpretación: El ECM nos dice en promedio cuánto se está equivocando el modelo, pero en unidades al cuadrado (no directamente en notas).
2. Raíz del ECM (RMSE)
Fórmula:
\[ RMSE = \sqrt{ECM} \]
¿Qué significa?
- Es simplemente la raíz cuadrada del ECM, para que la métrica esté en las mismas unidades que la variable dependiente (en este caso, el ingreso).
- Nos dice el error promedio del modelo.
RMSE <- sqrt(ECM)
RMSE## [1] 92159.1Interpretación: El RMSE indica, en promedio, cuánto se desvían las predicciones del modelo respecto a los ingresos reales.
- el error cuadratico medio fue: 8493299666, este valor representa el promedio del cuadrado de los errores de prediccion del modelo de regresion lienal multiple. Su error Estandar de Prediccion es de 92.165 lo cual representa aproximadamente un error del 10% si el ingreso promedio es de 900.000
3. Coeficiente de determinación \(R^2\)
Fórmula:
\[ R^2 = 1 - \frac{\sum (y_i - \hat{Y}_i)^2}{\sum (y_i - \bar{y})^2} \]
¿Qué significa?
- Mide qué proporción de la variación total de la variable dependiente (Ingreso) es explicada por la variable independiente (Gasto Familiar (sinarriendo)).
- Se interpreta como el porcentaje de la varianza explicada por el modelo.
Elementos de la fórmula:
- \(\sum (y_i - \hat{Y}_i)^2\): suma de los cuadrados de los residuos → SSE (Sum of Squares of Errors)
- \(\sum (y_i - \bar{y})^2\): suma de los cuadrados totales → SST (Total Sum of Squares)
- \(\bar{y}\): media de las notas
SSE <- sum(Base$residuos^2)
SST <- sum((y - media_y)^2)
R2 <- 1 - SSE / SST
R2## [1] 0.6431057Interpretación: El \(R^2\) indica que el modelo explica el 64.31% de la variabilidad total en los ingresos de las personas. Un \(R^2\) cercano a 1 indica un ajuste muy bueno; cercano a 0, un ajuste pobre.
- el 64.31% de la variabilidad en el ingreso neto de los hogares se explica por la variables incluidas en el modelo: gastos familiares sin arriendo, gastos de arriendo y estrato.
es un modelo aceptable, pero todavia hay un 35.69% de varibailidad del ingreso que no se explica con esas tres variables. con otros factores socioeconomicos relevantes que no estan incluidos.
plot(x, y, pch = 19, col = "yellow2",
xlab = "Ingreso Neto", ylab = "Gastos Familiar (sinarriendo)",
main = "Regresion lineal: Ingreso Neto ~ Gastos Familiar (sinarriendo)")
abline(beta_0, beta_1, col = "slateblue", lwd = 2)
Interpretaacion: se identifica como varian los gastos sin arriendo en funcion del ingreso neto del hogar.
-Existe una relacion positiva: A medida que el ingreso neto aumenta, tambien aumenta los gastos sin arriendo. Es coherente con la logica economica, hay hogares con mayor ingreso que suelen gastar mas. -Dispersion: Hay varibilidad en el gasto incluso para ingresos similares lo que justifica un \(R^{2}\) no tan alto.
Gráfico de residuos
plot(Base$predicciones, Base$residuos, pch = 19, col = "sienna",
xlab = "Valores predichos", ylab = "Residuos",
main = "Residuos vs Valores ajustados")
abline(h = 0, col = "springgreen", lty = 2)
Interpretacion -Los residuos estan centrados en torno a 0, indica que el modelo no esta sesgado sistematicamente -Se observa que los residuos tienden a expandirse a medida que aumentan los valores predichos.
Conclusión
El análisis muestra una relación positiva entre el Gasto Familiar (sin arriendo) y el Ingreso Neto. El modelo predice que por cada unidad adicional de gasto, el ingreso neto aumenta en promedio 0.33) pesos. El coeficiente de determinación \(R^2\) indica que una parte importante de la variabilidad del ingreso se explica por el gasto familiar.
Validación de supuestos del modelo de regresión
Los modelos de regresión lineal requieren que ciertos supuestos estadísticos se cumplan para que las estimaciones sean válidas y confiables. A continuación validamos los cuatro supuestos fundamentales.
1. Supuesto de linealidad
Descripción:
La relación entre la variable independiente \(X\) (IOngreso Neto) y la dependiente \(Y\) (Gasto Familiar (sinarriendo)) debe ser
lineal.
Forma teórica:
\[ E(Y_i | X_i) = \beta_0 + \beta_1 X_i \]
Cómo se verifica:
- Con un gráfico de los valores ajustados vs valores observados
- Con un gráfico de los residuos vs valores ajustados
# Gráfico de residuos vs valores ajustados
plot(Base$predicciones, Base$residuos,
main = "Residuos vs Valores Ajustados",
xlab = "Valores Ajustados", ylab = "Residuos",
pch = 19, col = "violetred")
abline(h = 0, col = "turquoise4", lty = 2)
Interpretación esperada: Los residuos deben estar distribuidos aleatoriamente alrededor de 0, sin una forma curva ni patrón sistemático.
- la mayoria de los residuos oscilan alrededor de la linea horizontal \(y=0\), lo cual es esperado. esto indica que el modelo no esta sesgado sistematicamente. -No existe patrones de curvas, esto sugiere la relacion entre ingreso neto y gastos sin arriendo es aproximadamente lineal, lo cual es apropiado para un modelo de regresion lineal.
2. Supuesto de homocedasticidad (varianza constante)
Descripción:
La varianza de los errores debe ser constante a lo
largo de todos los valores de \(X\).
Forma teórica:
\[ Var(\varepsilon_i) = \sigma^2, \quad \forall i \]
Cómo se verifica:
- Gráfico de residuos vs valores ajustados
- Test de Breusch-Pagan
#Prueba de Breusch-Pagan
library(lmtest)
bptest(modelo)##
## studentized Breusch-Pagan test
##
## data: modelo
## BP = 25.586, df = 1, p-value = 4.232e-07Interpretación esperada:
- Si el p-valor es mayor a 0.05, no se rechaza la
homocedasticidad (lo cual es bueno). - En el gráfico, la dispersión de
los residuos debe ser aproximadamente constante (no
forma de cono o embudo).
3. Supuesto de normalidad de los errores
Descripción:
Los errores \(\varepsilon_i\) deben
estar distribuidos normalmente, especialmente si vamos
a hacer pruebas de hipótesis.
Forma teórica:
\[ \varepsilon_i \sim N(0, \sigma^2) \]
Cómo se verifica:
- Histograma y QQ plot de los residuos
- Prueba de Shapiro-Wilk
# Histograma de residuos
hist(Base$residuos,
main = "Histograma de residuos",
col = "gray66",
probability = TRUE, # Esto es importante para que el eje Y sea una densidad
xlab = "Residuos")
# Agregar la curva normal (campana de Gauss)
curve(dnorm(x, mean = mean(Base$residuos), sd = sd(Base$residuos)),
col = "gray10", lwd = 2, add = TRUE)
# QQ plot
qqnorm(Base$residuos)
qqline(Base$residuos, col = "blue4")
# Prueba de normalidad
shapiro.test(Base$residuos)##
## Shapiro-Wilk normality test
##
## data: Base$residuos
## W = 0.95612, p-value = 0.0001096Interpretación esperada:
- El histograma debe parecerse a una campana de Gauss.
- En el QQ plot, los puntos deben estar alineados con la línea roja.
- En la prueba de Shapiro-Wilk, un p-valor > 0.05 indica que no se rechaza la normalidad.
-El histograma se parece razonablemente a una distribucion normal, aunque hay cierta asimetria, la normalidad de los residuos no esta gravemente violada.
4. Supuesto de independencia de los errores
Descripción:
Los errores deben ser independientes entre sí (sin
correlación serial).
Forma teórica:
\[ Cov(\varepsilon_i, \varepsilon_j) = 0 \quad \text{para } i \ne j \]
Cómo se verifica:
- Con la prueba de Durbin-Watson
dwtest(modelo)##
## Durbin-Watson test
##
## data: modelo
## DW = 1.083, p-value = 5.121e-09
## alternative hypothesis: true autocorrelation is greater than 0Interpretación esperada:
- El estadístico de Durbin-Watson debe estar cerca de 2.
- Un p-valor > 0.05 indica que no hay evidencia de autocorrelación, es decir, los errores son independientes.
Conclusión sobre los supuestos
Si todos los gráficos y pruebas son favorables (residuos aleatorios, varianza constante, distribución normal y errores independientes), podemos decir que los supuestos del modelo se cumplen y las inferencias realizadas con la regresión son confiables.
conlusiones para este modelo
linealidad: El gráfico de residuos vs valores ajustados no muestra un patrón sistemático o curvo, lo que sugiere que la relación entre el Ingreso Neto y el Gasto Familiar (sin arriendo) es aproximadamente lineal. Por lo tanto, el supuesto de linealidad se cumple.
Homocedasticidad: La prueba de Breusch-Pagan devuelve un p-valor mayor a 0.05, lo que indica que no hay evidencia suficiente para rechazar la homocedasticidad. Además, el gráfico no muestra un patrón de “cono”. Por lo tanto, el supuesto de varianza constante se cumple.
Normalidad de los errores: El histograma de residuos se asemeja visualmente a una campana de Gauss y el QQ plot muestra puntos mayormente alineados con la línea. Si el p-valor del test de Shapiro-Wilk es > 0.05, podemos asumir normalidad de los errores. En caso contrario (p < 0.05), hay ligero desvío de la normalidad, lo cual podría afectar pruebas de hipótesis, aunque no necesariamente la validez del modelo.
Independencia de los errores: La prueba de Durbin-Watson muestra un estadístico cercano a 2 y un p-valor mayor a 0.05, lo que sugiere que no existe autocorrelación entre los errores. Por lo tanto, el supuesto de independencia se cumple.
MODELO DE REGRESION LINEAL MULTIPLE
La regresión lineal múltiple es una extensión de la regresión lineal simple que permite modelar la relación entre una variable dependiente \(Y\) y múltiples variables independientes \(X_1, X_2, \dots, X_k\).
Este modelo surge cuando se desea explicar una respuesta utilizando más de un predictor, ya que muchos fenómenos del mundo real están influenciados por múltiples factores simultáneamente.
Se utiliza para:
- Predecir valores de \(Y\) en
función de varios predictores
- Evaluar la influencia relativa de cada variable
- Controlar el efecto de una variable mientras se examinan otras
Ecuación del modelo
La forma general del modelo de regresión lineal múltiple es:
\[ Y = \beta_0 + \beta_1 X_1 + \beta_2 X_2 + \cdots + \beta_k X_k + \varepsilon \]
Donde:
- \(Y\): Variable dependiente
(respuesta)
- \(X_i\): Variables independientes
(predictoras)
- \(\beta_0\): Intercepto
- \(\beta_i\): Coeficientes de
regresión asociados a cada predictor
- \(\varepsilon \sim \mathcal{N}(0, \sigma^2)\): Término de error aleatorio
Supuestos del modelo clásico Gauss-Markov
- Linealidad en los parámetros
- Independencia de los errores
- Homoscedasticidad: varianza constante de los
errores
- No multicolinealidad perfecta entre
predictores
- Normalidad de los errores (para inferencia estadística)
Los supuestos del modelo clásico Gauss-Markov son condiciones necesarias para que los estimadores obtenidos por Mínimos Cuadrados Ordinarios (OLS) sean:
- Insesgados: su valor esperado es el verdadero valor poblacional
- Eficientes: tienen la menor varianza posible entre los estimadores lineales insesgados
Es decir, bajo estos supuestos, los coeficientes estimados por OLS son los mejores estimadores lineales insesgados (BLUE: Best Linear Unbiased Estimators).
Supuestos del modelo
Linealidad en los parámetros
El modelo debe ser lineal respecto a los coeficientes: \[ Y = \beta_0 + \beta_1 X_1 + \cdots + \beta_k X_k + \varepsilon \]Esperanza del error igual a cero
Esto significa que los errores (residuos) tienen media cero. En términos prácticos, el modelo no tiene sesgo sistemático al predecir.
Cuando decimos que un modelo no tiene sesgo sistemático, nos referimos a que en promedio, las predicciones del modelo no se desvían del valor real.
\[ \mathbb{E}[\varepsilon_i] = 0 \]
Es decir, si repites muchas veces el modelo con diferentes muestras (o en teoría, infinitas veces), los errores se compensan entre sí: no tienden a ser siempre positivos o siempre negativos.
- No multicolinealidad perfecta entre los
predictores
Las variables explicativas no deben estar perfectamente correlacionadas. Esto garantiza que los coeficientes sean identificables.
En regresión múltiple:
Dos o más variables explicativas no deben estar perfectamente
correlacionadas.
Esto permite que los coeficientes de cada predictor sean identificables y únicos.
*No aplica a regresión simple, ya que solo hay una \(X\).
Solo es relevante en regresión múltiple**.
- Homoscedasticidad (varianza constante de los
errores)
\[ \text{Var}(\varepsilon_i) = \sigma^2 \quad \forall i \]
La varianza de los errores debe ser constante a lo largo de todas las observaciones. Si este supuesto se viola, ocurre heterocedasticidad, y aunque los coeficientes siguen siendo insesgados, ya no son eficientes (no son BLUE).
- Independencia de los errores (no
autocorrelación)
\[ \text{Cov}(\varepsilon_i, \varepsilon_j) = 0 \quad \text{para } i \ne j \] Es crucial en análisis de datos temporales o espaciales.
Los errores no deben estar correlacionados entre sí. Este supuesto es clave en análisis de series temporales, donde los valores están ordenados en el tiempo.
Teorema de Gauss-Markov
Bajo estos cinco supuestos, los estimadores obtenidos por OLS son:
BLUE (Best Linear Unbiased Estimators):
Los mejores (mínima varianza), lineales e insesgados estimadores disponibles.
¿Qué tan buenos son los estimadores que obtenemos al usar regresión lineal por mínimos cuadrados (OLS)?
BLUE (Best Linear Unbiased Estimators):
- Best: tienen la mínima varianza posible
- Linear: son funciones lineales de los datos
- Unbiased: no están sistemáticamente desviados del valor real
Este teorema es válido tanto en regresión simple como en múltiple. Sin embargo, el supuesto de no multicolinealidad solo es relevante en el contexto múltiple, donde hay más de un predictor.
Carga de datos
data(Base)## Warning in data(Base): data set 'Base' not foundhead(Base)## # A tibble: 6 × 10
## `>` u x `No de Personas a Cargo` y z `Gastos / Ingresos`
## <dbl> <dbl> <dbl> <chr> <dbl> <dbl> <dbl>
## 1 1 1 895212 Ninguna 367037 250659 0.41
## 2 2 1 822988 Ninguna 288046 222207 0.35
## 3 3 1 812288 Ninguna 324915 243686 0.4
## 4 4 1 627853 Ninguna 282534 175799 0.45
## 5 5 1 673359 Ninguna 282811 195274 0.42
## 6 6 1 696733 Ninguna 320497 195085 0.46
## # ℹ 3 more variables: `Arriendos/ Ingreso` <dbl>, predicciones <dbl>,
## # residuos <dbl>La data ‘Base’ (Base de Datos Alcaldia) es un conjunto de datos que contiene información sobre los ingresos y gatos de los funcionarios de la alcaldia de Neiva para el año 2019.
Descripción del dataset
El conjunto de datos Base contiene la siguiente
información. A continuación, se presenta la descripción de cada variable
que se trabajara:
| Variable | Significado | Tipo |
|---|---|---|
x |
Ingreso Neto | Continua |
y |
Gasto Familiar (sinarriendo) | Continua |
z |
Gasto de Arriendo | Continua |
u |
Estrato | Continua |
Ajuste del modelo de regresion multiple
Queremos predecir el ingreso (x) a partir de tres
variables:
u: Estrato
y: Gasto Familiar (sinarriendo)
z: Gasto de Arriendo
modelo <- lm(x ~ y + z + u, data = Base)
summary(modelo)##
## Call:
## lm(formula = x ~ y + z + u, data = Base)
##
## Residuals:
## Min 1Q Median 3Q Max
## -405091 -77842 -11129 69175 527081
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.024e+05 5.677e+04 1.803 0.0734 .
## y 2.125e-01 1.808e-01 1.175 0.2418
## z 9.975e-01 4.141e-01 2.409 0.0172 *
## u 3.202e+05 2.909e+04 11.008 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 156600 on 146 degrees of freedom
## Multiple R-squared: 0.834, Adjusted R-squared: 0.8306
## F-statistic: 244.5 on 3 and 146 DF, p-value: < 2.2e-16Interpretacion detallada de los coeficientes
coeficientes <- summary(modelo)$coefficients
coeficientes## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.023667e+05 5.677445e+04 1.803042 7.344390e-02
## y 2.124942e-01 1.808223e-01 1.175155 2.418458e-01
## z 9.975392e-01 4.141053e-01 2.408902 1.724695e-02
## u 3.202303e+05 2.909064e+04 11.008021 6.741029e-21Interpretacion:
Intercepto (1.02376e^{05}): El intercepto (constante) del modelo es aproximadamente 102,367, lo cual representa el valor esperado de la variable dependiente cuando todas las variables independientes (y, z, u) son cero.ie, El p-valor es 0.073, lo que indica que el intercepto no es estadísticamente significativo al 5%, pero sí lo sería al 10%.
y (0.212): Signidica que por cada peso que aumenta ‘y’, se espera que la variable dependiente aumente en 0.212, manteniendo las demás variables constantes.ie, El p-valor es 0.2418, por lo tanto, no es estadísticamente significativo. Esto sugiere que y no tiene un efecto claro sobre la variable dependiente en este modelo.
z (0.998): Lo que indica que por cada unidad adicional en z, la variable dependiente se incrementa aproximadamente en 0.998, manteniendo constantes las otras variables. ie, El p-valor es 0.017, lo que indica que sí es estadísticamente significativo al 5%, por lo tanto, z tiene un efecto relevante sobre la variable dependiente.
u (3.20230e^{05}): lo que sugiere un cambio muy grande en la variable dependiente por cada unidad adicional de u.ie, El p-valor es 6.74e-21, extremadamente bajo, por lo que es altamente significativo. u tiene un efecto muy fuerte y claro sobre la variable dependiente.
Ecuacion del modelo ajustado
cat("x =", round(coef(modelo)[1], 2), "+",
round(coef(modelo)[2], 2), "* y +",
round(coef(modelo)[3], 2), "* z +",
round(coef(modelo)[4], 5), "* u")## x = 102366.7 + 0.21 * y + 1 * z + 320230.3 * uDiagnostico del modelo
Grafico de residuos vs valores ajustados
plot(modelo$fitted.values, modelo$residuals,
main = "Residuos vs Valores Ajustados",
xlab = "Valores Ajustados", ylab = "Residuos", pch = 19, col = "paleturquoise4")
abline(h = 0, col = "orchid3", lty = 2)
Interpretacion:
No deberían observarse curvas, grupos, ni formas reconocibles. Si hay patrones, podrían indicar problemas de especificación del modelo (como una relación no lineal no capturada o variables omitidas). Las concentraciones de puntos o grupos separados pueden indicar valores atípicos o que hay variables omitidas.
Normalidad de los residuos
hist(modelo$residuals,
main = "Histograma de residuos",
col = "cadetblue2",
probability = TRUE,
# Esto es importante para que el eje Y sea una densidad
xlab = "Residuos")
# Agregar la curva normal (campana de Gauss)
curve(dnorm(x, mean = mean(modelo$residuals), sd = sd(modelo$residuals)),
col = "cyan4", lwd = 2, add = TRUE)
qqnorm(modelo$residuals)
qqline(modelo$residuals, col = "chartreuse")
Interpretacion:
- El histograma tiene forma de campana y se asemeja a la curva normal superpuesta.
- En el gráfico Q-Q, los puntos están alineados a lo largo de la línea roja.
Pruebas adicionales
library(lmtest)
bptest(modelo)##
## studentized Breusch-Pagan test
##
## data: modelo
## BP = 36.737, df = 3, p-value = 5.231e-08dwtest(modelo)##
## Durbin-Watson test
##
## data: modelo
## DW = 1.9721, p-value = 0.3757
## alternative hypothesis: true autocorrelation is greater than 0Interpretacion: - El valor p es muy pequeño (p < 0.001), lo que lleva a rechazar la hipótesis nula.
- Esto indica que hay evidencia estadísticamente significativa de heterocedasticidad en los residuos del modelo.
Las consecuencias estan en que las estimaciones aún pueden ser insesgadas, pero las inferencias estadísticas (intervalos y pruebas) pueden no ser válidas si no se corrige esta heterocedasticidad.
Multicolinealidad
Puede afectar la estabilidad de los coeficientes. Si dos o más variables independientes están muy correlacionadas entre sí, los estimadores pueden volverse inestables. Se evalúa con el VIF (Variance Inflation Factor).
install.packages(“car”) # Solo la primera vez library(car) vif(modelo)
VIF ≈ 1: No hay colinealidad.
VIF entre 1 y 5: Moderada colinealidad (aceptable).
VIF > 10: Colinealidad alta. Puede ser preocupante.
vif_manual <- function(modelo) {
X <- model.matrix(modelo)[, -1]
sapply(1:ncol(X), function(i) {
rsq <- summary(lm(X[, i] ~ X[, -i]))$r.squared
1 / (1 - rsq)
})
}
vif_manual(modelo)## [1] 4.761572 4.653089 3.452429Interpretacion: - Ninguna variable supera el umbral crítico de 5, por lo tanto no hay colinealidad severa. Sin embargo, los valores cercanos a 5 indican que sí existe una cierta multicolinealidad moderada, especialmente en las primeras dos variables. *Esto no representa un problema grave, pero es recomendable estar alerta si el modelo se vuelve inestable (por ejemplo, con cambios de signo o significancia al agregar o quitar variables).
Conclusion
El modelo de regresión muestra que las variables ‘z’ (gasto de arriendo) y ‘u’ (estrato) tienen efectos positivos y estadísticamente significativos sobre el ingreso (x), siendo u el predictor más influyente. En cambio, ‘y’ (gasto familiar sin arriendo) no presenta un efecto significativo. Los residuos se distribuyen aproximadamente de forma normal y no hay evidencia de autocorrelación. Sin embargo, la prueba de Breusch-Pagan detecta heterocedasticidad, lo que puede afectar la validez de las inferencias. Los factores de inflación de varianza (VIF) indican una colinealidad moderada pero no preocupante. En conjunto, el modelo es útil para explicar el ingreso, aunque requiere ajustes para corregir la heterocedasticidad.
Análisis de Varianza (ANOVA)
El análisis de varianza (ANOVA) permite evaluar si el modelo de regresión múltiple en su conjunto es estadísticamente significativo. Compara la variabilidad explicada por el modelo con la variabilidad residual no explicada.
\[ F = \frac{\text{MSM}}{\text{MSR}} = \frac{\frac{\text{SSM}}{k}}{\frac{\text{SSR}}{n - k - 1}} \]
Donde:
\(\text{SSM}\): Suma de cuadrados del modelo
Representa la variabilidad de la variable dependiente que es explicada por el conjunto de predictores del modelo. Cuanto mayor es el SSM, mayor es la parte de la variabilidad total que se puede explicar a través del modelo.\(\text{SSR}\): Suma de cuadrados del residuo Mide la variabilidad que no se logra explicar con el modelo; en otras palabras, es la suma de las desviaciones de las observaciones a la línea de regresión (o a los valores predichos). Una SSR pequeña indica que el modelo se ajusta bien a los datos.
\(k\): Número de predictores
Es la cantidad de variables independientes utilizadas en el modelo. Cada predictor aporta un grado de libertad al SSM.\(n\): Número total de observaciones
Corresponde al tamaño muestral, y afecta el cálculo de los grados de libertad tanto para el modelo como para el residuo.
Cálculo de las Medias Cuadráticas Para poder comparar de forma equitativa la variabilidad explicada y la no explicada, se dividen las sumas de cuadrados por sus respectivos grados de libertad, obteniéndose así las medias cuadráticas:
\(\text{MSM}\): Media cuadrática del modelo
Se calcula dividiendo el SSM entre el número de predictores Esto representa el promedio de variabilidad explicada por cada predictor.\(\text{MSR}\): Media cuadrática del residuo
Se calcula dividiendo el SSR entre los grados de libertad del residuo, que es \(n−k−1\)
Establecimiento de la Prueba F La razón de F se define como la relación entre las dos medias cuadráticas:
\[ F = \frac{\text{MSM}}{\text{MSR}} = \frac{\frac{\text{SSM}}{k}}{\frac{\text{SSR}}{n - k - 1}} \]
se utiliza para probar la hipótesis nula de que todos los coeficientes de regresión (excepto el intercepto) son cero. Es decir:
Hipótesis nula \(H_0\): Ninguno de los predictores tiene un efecto significativo sobre la variable dependiente (el modelo no explica la variabilidad).
Hipótesis alternativa \(H_a\): Al menos uno de los predictores tiene un efecto significativo.
Interpretación del Valor F
Valor F Alto: Indica que la variabilidad explicada por el modelo es considerablemente mayor que la variabilidad residual. Esto sugiere que al menos uno de los predictores contribuye significativamente a explicar la variable dependiente, y, por lo tanto, el modelo en su conjunto es estadísticamente significativo.
Valor F Bajo: Sugiere que la variabilidad explicada es similar a la residuo, lo que implicaría que los predictores no mejoran significativamente el ajuste del modelo en comparación con un modelo sin predictores.
El análisis ANOVA en regresión múltiple permite, a través de la comparación de las medias cuadráticas, determinar si el modelo tiene la capacidad de explicar una parte significativa de la variación observada en los datos. Un valor F elevado y estadísticamente significativo (según un umbral de p-valor preestablecido, generalmente 0.05) rechaza la hipótesis nula y confirma que el modelo tiene un poder explicativo importante.
anova(modelo)## Analysis of Variance Table
##
## Response: x
## Df Sum Sq Mean Sq F value Pr(>F)
## y 1 1.3862e+13 1.3862e+13 565.52 < 2.2e-16 ***
## z 1 1.1438e+12 1.1438e+12 46.66 2.119e-10 ***
## u 1 2.9703e+12 2.9703e+12 121.18 < 2.2e-16 ***
## Residuals 146 3.5788e+12 2.4512e+10
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Análisis de Varianza (ANOVA)
La tabla ANOVA contiene los siguientes elementos:
- Df: Grados de libertad para cada fuente de
variación.
- Sum Sq: Suma de cuadrados, indica la cantidad de
variabilidad explicada.
- Mean Sq: Media cuadrática, es la suma de cuadrados
dividida entre los grados de libertad.
- F value: Estadístico F de comparación.
- Pr(>F): Valor-p asociado al test F.
Nota: ### ¿Qué significa Df (Degrees of Freedom)?
Los grados de libertad representan la cantidad de información independiente disponible para estimar una cierta cantidad de variabilidad.
\(\textbf{En términos
sencillos:}\)
Los grados de libertad \(\textbf{Df}\)
indican cuántos valores pueden variar libremente antes de que se imponga
una restricción (como una media o un coeficiente a estimar).
En la tabla ANOVA:
- Para cada predictor en el modelo (\(y, z, u\)), el grado de libertad es \(1\), porque cada uno añade un parámetro (su coeficiente).
- Para los residuos, el número de grados de libertad es:
\(Df_{residuos} = n - k - 1\)
donde:
- \(n\): número total de
observaciones
- \(k\): número de predictores
- El \(-1\) se debe al intercepto del modelo
Ejemplo práctico
En el modelo lm(x ~ y + z + u) con 150
observaciones:
- $ n = 150$
- $k = 3 $
Entonces:
- Grados de libertad del modelo: \(3\)
- Grados de libertad de los residuos:
\(150 - 3 - 1 = 146\)
Interpretación de resultados
Interpretación de la Tabla ANOVA
La tabla ANOVA obtenida del modelo lm(x ~ y + z + u)
permite analizar la contribución de cada predictor a la variabilidad
explicada del consumo de combustible (x).
y(Gasto Familiar (sinarriendo)): Presenta una suma de cuadrados de \(1.386e^{+13}\), un valor \(F = 565.52\) muy alto y un valor-p de \(p = 2.2e^{-16}\) el cual es extremadamente bajo e indican que este predictor es altamente significativo.Es decir, el gasto familiar sin arriendo tiene un efecto muy fuerte en la variable dependiente, probablemente explicando la mayor parte de la variabilidad total.z(Gasto Arriendo): Tiene una suma de cuadrados de \(1438e^{12}\), un valor \(F = 46.66\) , y el valor \(p = 2.119e^{-10}\). Este predictor también es estadísticamente significativo al 1%, aunque con un menor efecto comparado con y. El gasto en arriendo explica una porción adicional y relevante de la variabilidad en la variable respuesta.‘u (Estrato)’: Tiene una suma de cuadrados de \(2.9703e^{12]\), un valor \(F = 121.18\), y valor \(p = 2.2e^{-16}\). Aunque se dice que “no es significativo”, el valor-p indica lo contrario: es altamente significativo. Es probable que haya un error en la interpretación original. Este resultado sugiere que el estrato también contribuye significativamente a explicar la variación en la variable dependiente, aunque en menor magnitud que y.
Residuals(residuos): La suma de cuadrados residual es \(146 3.5788e^{+12}\) con media cuadrática de \(2.4512e^{10}\). Esto significa que una parte considerable de la variabilidad total en el gasto familiar no es explicada por los predictores del modelo \((y, z, u)\). Aunque el modelo presenta predictores significativos, la suma de cuadrados residual aún es grande en relación con las otras sumas de cuadrados, lo que sugiere que hay factores adicionales que podrían estar influyendo en la variable dependiente y que no están incluidos en el modelo. Esto significa que una parte considerable de la variabilidad total en el gasto familiar no es explicada por los predictores del modelo (y, z y u). Aunque el modelo presenta predictores significativos, la suma de cuadrados residual aún es grande en relación con las otras sumas de cuadrados, lo que sugiere que hay factores adicionales que podrían estar influyendo en la variable dependiente y que no están incluidos en el modelo.
La siguiente tabla representa un análisis de varianza (ANOVA)
aplicado a un modelo de regresión lineal múltiple donde la variable
dependiente es x (Ingreso Neto), y las variables
independientes son x (Gasto Familiar (sinarriendo)),
z (Gasto Arriendo) y u (Estrato).
Descripción de la Tabla ANOVA
Df (Grados de libertad): Número de parámetros independientes asociados con la variable. Cada variable tiene 1 grado de libertad.
Sum Sq (Suma de cuadrados): Mide la cantidad de variación en
xexplicada por la variable.Mean Sq (Media de cuadrados): Es la suma de cuadrados dividida por los grados de libertad.
F value (Estadístico F): Compara la varianza explicada por la variable con la varianza residual.}
Pr(>F) (valor-p): Valor de probabilidad asociado al estadístico F. Un valor bajo (por ejemplo, < 0.05) indica significancia estadística.
Resultados por Variable
y (Gasto Familiar (sinarriendo))
- Sum Sq: \(1.3862e^{13}\)
- F value: 565.52
- Pr(>F): < 2.2e^{-16} (***)
- Interpretación: La variable
xtiene un efecto altamente significativo sobre elGasto Familiar (sinarriendo). Su gran suma de cuadrados y elevado estadístico F indican que explica una parte sustancial de la variabilidad de esta variable. Por tanto,xes un predictor muy influyente del gasto familiar excluyendo el arriendo.
z (Gasto Arriendo)
- Sum Sq: 1.1438e^{12}
- F value: 46.66
- Pr(>F): 2.119e^{-10} (***)
- Interpretación: La variable
xtambién influye significativamente en el Gasto de Arriendo. Aunque su efecto es menor comparado con el gasto familiar sin arriendo, sigue siendo estadísticamente significativo. Esto sugiere que x está asociado con el nivel de gasto en arriendo, aunque en menor medida.
u (Estrato)
- Sum Sq: 2.9703e^{12}
- F value: 121.18
- Pr(>F): < 2.2e^{-16} (***)
- Interpretación: A pesar de tener una suma de
cuadrados más baja que
y, el estrato también se ve significativamente influido porx. Esto sugiere que hay una asociación estadísticamente significativa entre x y el estrato socioeconómico del hogar.xpodría estar relacionado con factores estructurales que varían por estrato, como ingresos o acceso a servicios.
Conclusión
Los factores \(y, z, u\) tienen un efecto estadísticamente significativo sobre la variable de respuesta x. Esto se evidencia por los valores p extremadamente bajos para cada uno de los predictores (todos menores a 0.001), lo que indica que es altamente improbable que sus efectos se deban al azar.
En particular:
yes el predictor que más explica la variabilidad enx, con el mayor valor de suma de cuadrados y el valor F más alto (\(565.52\)).utambién tiene un efecto fuerte y significativo (\(F = 121.18\)).zaunque con menor impacto relativo, también es significativo (\(F = 46.66\)).
Además, el modelo explica una porción considerable de la variabilidad
total en x, dado que la suma de cuadrados de los residuos
es mucho menor que la suma de cuadrados de los predictores.
Conclusión final: El modelo es adecuado y todos los
predictores incluidos (\(y, z, u\)) son
relevantes para explicar la variación en la variable dependiente
x. El modelo explica aproximadamente el \(83.4\)% de la variabilidad de la variable
dependiente x, lo cual indica un ajuste muy bueno.
#Nota ¿Cuándo sí puede importar el orden?
Hay algunas situaciones específicas donde el orden sí puede influir:
- Regresión paso a paso (stepwise) Aquí, el modelo se construye añadiendo o eliminando predictores uno a uno.
En estos métodos automatizados (como forward selection, backward elimination o stepwise regression), el orden puede afectar qué variables entran o salen del modelo.
Por eso, estos métodos deben usarse con cuidado y, preferentemente, junto con validación cruzada.
Multicolinealidad Si hay alta correlación entre predictores, el orden no cambia el ajuste, pero sí puede hacer difícil interpretar los coeficientes, ya que el modelo distribuye la “explicación” entre variables correlacionadas de forma arbitraria.
Modelos jerárquicos (ANOVA jerárquico o regresión por bloques) En este enfoque, tú eliges el orden de ingreso de bloques de variables para evaluar cómo mejora el modelo paso a paso.
Aquí el orden tiene sentido interpretativo y analítico.
Conclusion Para regresión múltiple estándar: el orden no importa.
Para regresiones paso a paso o jerárquicas: el orden puede importar y debe estar justificado.
MODELO DE REGRESION LOGISTICO
## Regresión Logística
¿Qué es la regresión logística?
La regresión logística es un modelo estadístico utilizado cuando la variable dependiente es categórica binaria, es decir, toma solo dos posibles valores (por ejemplo: 0/1, sí/no, éxito/fracaso).
La regresión logística se utiliza cuando se quiere predecir la probabilidad de que ocurra un evento binario, es decir, un resultado con solo dos posibles valores (por ejemplo, éxito/fracaso, sí/no, 1/0).
A diferencia de la regresión lineal, que predice cualquier número real, la regresión logística necesita que el valor predicho esté siempre entre 0 y 1, porque está modelando una probabilidad.
Se usa para estimar la probabilidad de que ocurra un evento en función de una o más variables predictoras.
NOTA: Para casos donde no se cuente con variables binarias, es necesario crearla (como en este caso)
La forma del modelo es:
\[ \log\left(\frac{p}{1 - p}\right) = \beta_0 + \beta_1 X_1 + \beta_2 X_2 + \cdots + \beta_k X_k \]
Donde:
\(p\): probabilidad de que ocurra el evento (por ejemplo, transmisión manual , que una persona tenga una enfermedad).
\[ \frac{p}{1 - p}\]: odds o razón de probabilidades. Se llama odds o razón de probabilidades. Representa cuán probable es que ocurra el evento comparado con que no ocurra. Ejemplo: si \(p=0.75\), entonces: \[\frac{p}{1−p} =\frac{0.75}{0.25}=3 \]
Esto significa que el evento es 3 veces más probable que no ocurrir.
\[\( \log \left( \cdot \right) \)\]: logaritmo natural (logit)
Logaritmo natural, también llamado función logit
Es el logaritmo natural, que se aplica a los odds. Esto transforma la probabilidad, que está entre 0 y 1, a un valor entre −∞ y +∞, lo que permite usar una ecuación lineal para modelarla.
###Interpretación
La regresión logística transforma una probabilidad (que siempre está entre 0 y 1) en una variable continua mediante el logaritmo de las odds. Esto permite usar un modelo lineal para estimar probabilidades y evaluar cómo influyen distintos factores (por ejemplo, edad, sexo, antecedentes) sobre la ocurrencia de un evento.
###¿Cuándo usar regresión logística?
Cuando la variable dependiente es binaria (es decir
0y1).Cuando se quiere predecir la probabilidad de pertenecer a una clase.
Cuando no se cumplen los supuestos de normalidad de la regresión lineal.
##La función logística (función sigmoidal)
La función logística, también conocida como función sigmoidal, es la base matemática de la regresión logística. Se utiliza para modelar la probabilidad de ocurrencia de un evento cuando la variable dependiente es binaria (0 u 1).
Matemáticamente, se define como:
\[p(x) = \frac{1}{1 + e^{-(\beta_0 + \beta_1 X_1 + \beta_2 X_2 + \cdots + \beta_k X_k)}}\]
Donde:
- $ p(X) :$ probabilidad estimada de que ocurra el evento (por ejemplo, que tenga personas a cargo)
- $e $: base del logaritmo natural
- \(\beta_0, \beta_1, \dots, \beta_k :\) coeficientes estimados del modelo
- \(( X_1, X_2, \dots, X_k \):\) variables predictoras
##¿Por qué regresión logística y no lineal?
Cuando se quiere predecir una variable cualitativa binaria, como el tipo de las personas a cargo (si tiene o no tiene), se podría pensar en usar regresión lineal si esta se codifica como 0 y 1. Sin embargo, esto puede llevar a problemas importantes.
El logit: vínculo con el modelo lineal La regresión logística modela el logit de la probabilidad:
\[log(p_1−p)=β_0+β_1X_1+⋯+β_kX_k\]
Donde, \(p_1−p\) es la odds o razón de probabilidades.
Este enfoque permite usar una estructura lineal sobre los log-odds, y luego transformarla a probabilidades usando la función sigmoidal. darme las estructura lineal sobre los log-odds, y luego transformarla a probabilidades usando la función sigmoidal.
Problema con regresión lineal
Para esta base de datos se binarizara la variable
No de Personas a Cargo, y se creara una nueva variable que
indique si la persona tiene o no tiene personas a cargo. Por ejemplo, si
el valor es “Ninguna”, se codifica como 0; en cualquier otro caso, como
1.
Supongamos que queremos predecir Pc (No de Personas a
Cargo) en función del Gasto Familiar (sinarriendo) (y).
Codificamos Pc como 0 (Ninguno) y 1 (1 o <1), y
ajustamos un modelo lineal.
# Leer los datos
dato <- Base
# El dataframe se llama 'data'
dato$Pc <- ifelse(dato$`No de Personas a Cargo` == "Ninguna", 0, 1)
## es decir:
## * 0 si no tiene personas a cargo
## * 1 si tiene una o más personas a cargo. data(dato)## Warning in data(dato): data set 'dato' not found dato$Pc <- factor(dato$Pc, labels = c("Ninguno", "Uno o mas"))
dato$Pc_n <- ifelse(dato$Pc == "Uno o mas", 1, 0)
head(dato)## # A tibble: 6 × 12
## `>` u x `No de Personas a Cargo` y z `Gastos / Ingresos`
## <dbl> <dbl> <dbl> <chr> <dbl> <dbl> <dbl>
## 1 1 1 895212 Ninguna 367037 250659 0.41
## 2 2 1 822988 Ninguna 288046 222207 0.35
## 3 3 1 812288 Ninguna 324915 243686 0.4
## 4 4 1 627853 Ninguna 282534 175799 0.45
## 5 5 1 673359 Ninguna 282811 195274 0.42
## 6 6 1 696733 Ninguna 320497 195085 0.46
## # ℹ 5 more variables: `Arriendos/ Ingreso` <dbl>, predicciones <dbl>,
## # residuos <dbl>, Pc <fct>, Pc_n <dbl> modelo_lineal <- lm(Pc_n ~ y, data = dato)
library(ggplot2)## Warning: package 'ggplot2' was built under R version 4.4.3 ggplot(dato, aes(x = y, y = Pc_n)) +
geom_point(aes(color = Pc), shape = 1) +
geom_smooth(method = "lm", color = "gray20", se = FALSE) +
theme_bw() +
labs(title = "Regresión lineal para tipo de personas a cargo", y = "Probabilidad de tener personas a cargo") +
theme(legend.position = "none")## `geom_smooth()` using formula = 'y ~ x'
Ejemplo de predicción con valores extremos
predict(modelo_lineal, newdata = data.frame(y = 7))## 1
## 0.4027713NOta: El modelo puede dar valores mayores que 1 o menores que 0, lo cual no tiene sentido para una probabilidad. Para evitar estos problemas, la regresión logística transforma el valor devuelto por la regresión lineal \((β_0+β_{1X}\) empleando una función cuyo resultado está siempre comprendido entre 0 y 1. Existen varias funciones que cumplen esta descripción, una de las más utilizadas es la función logística (también conocida como función sigmoide):
\[\sigma(x) = \frac{1}{1 + e^{-x}}\]
Esto asegura que las predicciones estén siempre entre 0 y 1.
Para valores de \(x\) muy grandes
positivos, el término \(e^{-x}\) tiende
a 0, por lo tanto \(\sigma(x) \to
1\).
Para valores de \(x\) muy grandes
negativos, el término \(e^{-x}\) tiende
a \(\infty\), por lo tanto \(\sigma(x) \to 0\).
Sustituyendo \(x\) por una combinación lineal de predictores, es decir \(x = \beta_0 + \beta_1 X\), se obtiene:
\[ P(Y = k \mid X = x) = \frac{1}{1 + e^{-(\beta_0 + \beta_1 X)}} \]
También puede expresarse como:
\[ P(Y = k \mid X = x) = \frac{e^{\beta_0 + \beta_1 X}}{1 + e^{\beta_0 + \beta_1 X}} \]
Esto representa la de que la variable cualitativa \(Y\) tome el valor \(k\) (codificado como 1), dado un valor \(x\) del predictor \(X\).
Para poder estimar este modelo con una regresión lineal, se toma el logaritmo natural de los \(\textbf{odds}\) (razón de probabilidades), lo que se denomina \(\textbf{logit}\):
\[ \log\left( \frac{P(Y = k \mid X = x)}{1 - P(Y = k \mid X = x)} \right) = \beta_0 + \beta_1 X \]
La expresión $ P(Y = k X = x) $ puede interpretarse como la \(\textbf{probabilidad de que la variable cualitativa \(Y\) tome el valor \(k\)}\) (el nivel de referencia, codificado como 1), dado que el predictor \(\(X\)\) tiene el valor \(\(x\)\).
Esta función puede ajustarse de forma sencilla mediante técnicas de regresión, si se aplica su transformación logarítmica. Esta transformación es conocida como el \(\textbf{logit}\) o \(\textbf{logaritmo de los odds}\), y convierte una probabilidad acotada en el intervalo \([0,1]\) en una variable continua en \(-\infty, +\infty)\), lo cual permite aplicar regresión lineal.
La forma logarítmica del modelo logístico es:
\[ \ln\left( \frac{P(Y = k \mid X = x)}{1 - P(Y = k \mid X = x)} \right) = \beta_0 + \beta_1 X \]
Donde el lado izquierdo representa el \(\textbf{logaritmo natural de probabilidades}\) \(\emph{log-odds}\) de que \(Y = k\) dado un valor \(x\) del predictor \(X\).
La expresión \(P(Y = k \mid X = x)\) puede interpretarse como la \textbf{probabilidad de que la variable cualitativa \(\(Y\)\) tome el valor \(\(k\)\) (el nivel de referencia, codificado como 1), dado que el predictor \(\(X\)\) tiene el valor \(\(x\)\).
Esta función puede ajustarse de forma sencilla mediante técnicas de regresión, si se aplica su transformación logarítmica. Esta transformación es conocida como el \(\textbf{logit}\) o \(\textbf{logaritmo de los odds}\), y convierte una probabilidad acotada en el intervalo \([0,1]\) en una variable continua en \((-\infty, +\infty)\), lo cual permite aplicar regresión lineal.
La forma logarítmica del modelo logístico es:
\[ \log\left( \frac{P(Y = k \mid X = x)}{1 - P(Y = k \mid X = x)} \right) = \beta_0 + \beta_1 X \]
Donde el lado izquierdo representa el \(\textbf{logaritmo de la razón de probabilidades}\) \(\emph{log-odds}\) de que \(Y = k\) dado un valor \(x\) del predictor \(X\).
fin de la nota
## Modelo regresión logística
modelo_logistico <- glm(Pc_n ~ y, data = dato, family = "binomial")Grafica
ggplot(dato, aes(x = y, y = Pc_n)) +
geom_point(aes(color = Pc), shape = 1) +
stat_function(fun = function(x) {
predict(modelo_logistico, newdata = data.frame(y = x), type = "response")
}, color = "gray30") +
theme_bw() +
labs(title = "Regresión logística para ninguna persona a cargo", y = "Probabilidad de tener una persona a cargo") +
theme(legend.position = "none")
otro tipo de grafico
ggplot(dato, aes(x = y, y = Pc_n)) +
geom_point(aes(color = Pc), shape = 1) +
geom_smooth(method = "glm", method.args = list(family = "binomial"), se = FALSE) +
theme_bw() +
labs(title = "Curva logística con geom_smooth", y = "Probabilidad de tener personas a cargo") +
theme(legend.position = "none")## `geom_smooth()` using formula = 'y ~ x'
Función Sigmoide y Regresión Logística
La función logística o \(\textbf{sigmoide}\) se define como:
\(\sigma(x) = \frac{1}{1 + e^{-x}}\)
¿Por qué usar la función logística?
En regresión lineal, las predicciones pueden ser menores que 0 o mayores que 1, lo cual no tiene sentido cuando hablamos de probabilidades.
La función logística garantiza que todas las predicciones estén en el intervalo \([0, 1]\).
La relación entre los predictores y la probabilidad de clase es no lineal, pero se logra modelando de forma lineal el logit (log-odds).
El logit: vínculo con el modelo lineal
La regresión logística modela el logit de la probabilidad:
\[\log\left(\frac{p}{1 - p}\right) = \beta_0 + \beta_1 X_1 + \cdots + \beta_k X_k\]
Donde \[\frac{p}{1 - p}\] es la odds o razón de probabilidades.
Este enfoque permite usar una estructura lineal sobre los log-odds, y luego transformarla a probabilidades usando la función sigmoidal.
Visualización
Puedes visualizar la forma de la función logística con:
curve(1 / (1 + exp(-x)), from = -10, to = 10, col = "blue", lwd = 2,
ylab = "Probabilidad", xlab = "Logit (x)", main = "Función logística")
abline(h = 0.5, lty = 2, col = "red")
Objetivo
Este documento muestra cómo construir un modelo de regresión
logística simple usando el conjunto de datos
dato.
El objetivo es predecir si una persona puede tener personas a
cargo (Pc = 1) o Ninguna
(Pc= 0) utilizando una sola variable explicativa: el
Gasto Familiar (sinarriendo) (y).
1. Cargar y preparar los datos
Primero, cargamos el conjunto de datos dato, incluido en
R, y convertimos la variable Pc (No personas a cargo) en un
factor para facilitar la interpretación.
# Convertir variable y a factor
dato$Pc_factor <- factor(dato$Pc, levels = c(0, 1), labels = c("Ninguna", "1 o más"))
# Ver estructura y primeras filas
str(dato)## tibble [150 × 13] (S3: tbl_df/tbl/data.frame)
## $ > : num [1:150] 1 2 3 4 5 6 7 8 9 10 ...
## $ u : num [1:150] 1 1 1 1 1 1 1 1 1 1 ...
## $ x : num [1:150] 895212 822988 812288 627853 673359 ...
## $ No de Personas a Cargo: chr [1:150] "Ninguna" "Ninguna" "Ninguna" "Ninguna" ...
## $ y : num [1:150] 367037 288046 324915 282534 282811 ...
## $ z : num [1:150] 250659 222207 243686 175799 195274 ...
## $ Gastos / Ingresos : num [1:150] 0.41 0.35 0.4 0.45 0.42 0.46 0.37 0.4 0.42 0.43 ...
## $ Arriendos/ Ingreso : num [1:150] 0.28 0.27 0.3 0.28 0.29 0.28 0.29 0.3 0.27 0.3 ...
## $ predicciones : num [1:150] 425179 401609 398117 337927 352778 ...
## $ residuos : num [1:150] -58142 -113563 -73202 -55393 -69967 ...
## $ Pc : Factor w/ 2 levels "Ninguno","Uno o mas": 1 1 1 1 1 1 1 1 1 1 ...
## $ Pc_n : num [1:150] 0 0 0 0 0 0 0 0 0 0 ...
## $ Pc_factor : Factor w/ 2 levels "Ninguna","1 o más": NA NA NA NA NA NA NA NA NA NA ... head(dato)## # A tibble: 6 × 13
## `>` u x `No de Personas a Cargo` y z `Gastos / Ingresos`
## <dbl> <dbl> <dbl> <chr> <dbl> <dbl> <dbl>
## 1 1 1 895212 Ninguna 367037 250659 0.41
## 2 2 1 822988 Ninguna 288046 222207 0.35
## 3 3 1 812288 Ninguna 324915 243686 0.4
## 4 4 1 627853 Ninguna 282534 175799 0.45
## 5 5 1 673359 Ninguna 282811 195274 0.42
## 6 6 1 696733 Ninguna 320497 195085 0.46
## # ℹ 6 more variables: `Arriendos/ Ingreso` <dbl>, predicciones <dbl>,
## # residuos <dbl>, Pc <fct>, Pc_n <dbl>, Pc_factor <fct>Ajustar el modelo de regresión logística simple
Se ajusta un modelo de regresión logística donde la variable de respuesta es Pc, y la única variable predictora es y (Gasto Familiar (sinarriendo)).
El modelo estima la probabilidad de tener una o mas personas a cargo.
# Ajuste del modelo logístico simple
modelo_simple <- glm(Pc ~ y, data = dato, family = binomial)
# Ver resumen del modelo
summary(modelo_simple)##
## Call:
## glm(formula = Pc ~ y, family = binomial, data = dato)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.168e+00 7.461e-01 -1.565 0.1176
## y 4.946e-06 1.622e-06 3.049 0.0023 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 165.32 on 149 degrees of freedom
## Residual deviance: 153.81 on 148 degrees of freedom
## AIC: 157.81
##
## Number of Fisher Scoring iterations: 4Interpretación del modelo
Ecuación del modelo
La regresión logística modela el logit (logaritmo de las odds) como una combinación lineal de los predictores:
\[ \text{logit}(p) = \beta_0 + \beta_1 \cdot \text{wt} \]
Donde:\[ \text{logit}(p) = \log\left(\frac{p}{1 - p}\right) \]
Coeficientes
Predicción y clasificación
Calculamos la probabilidad predicha para cada observación y la clasificamos como “una o mas” si la probabilidad es mayor a 0.5.
# Probabilidades predichas
probabilidades <- predict(modelo_simple, type = "response")
# Clasificación binaria
pred_clase <- ifelse(probabilidades > 0.5, "uno o mas", "Ninguno")
# Mostrar algunas predicciones
head(data.frame(
Ingreso = dato$y,
Probabilidad = round(probabilidades, 3),
Clasificación = pred_clase
))## Ingreso Probabilidad Clasificación
## 1 367037 0.656 uno o mas
## 2 288046 0.564 uno o mas
## 3 324915 0.608 uno o mas
## 4 282534 0.557 uno o mas
## 5 282811 0.557 uno o mas
## 6 320497 0.603 uno o masEvaluación del modelo
Evaluamos el rendimiento del modelo usando una matriz de confusión y el porcentaje de aciertos (precisión).
# Matriz de confusión
table(Predicho = pred_clase, Real = dato$Pc)## Real
## Predicho Ninguno Uno o mas
## uno o mas 36 114# Precisión del modelo
mean(pred_clase == dato$Pc)## [1] 0Visualización del modelo
Graficamos la curva logística ajustada sobre los puntos observados para visualizar la relación entre las personas a cargo y la probabilidad de tener una o mas personas a cargo.
library(ggplot2)
ggplot(dato, aes(x = y, y = as.numeric(Pc) - 1)) +
geom_point(aes(color = Pc), size = 3) +
stat_smooth(method = "glm", method.args = list(family = "binomial"), se = FALSE) +
labs(
title = "Probabilidad de personas a cargo segun el ingreso",
x = "Gasto Familiar (sinarriedno) (y)",
y = "Probabilidad estimada",
color = "Personas a cargo"
) +
theme_minimal()## `geom_smooth()` using formula = 'y ~ x'
Conclusión
Tendencia negativa: La curva desciende a medida que aumenta el gasto familiar, lo que indica que a mayor ingreso familiar, menor es la probabilidad de tener personas a cargo.
Separación de clases: Los puntos de color representan los casos reales (con o sin personas a cargo). La curva suavizada muestra cómo el modelo estima una alta probabilidad de tener personas a cargo para niveles bajos de ingreso, y una baja probabilidad a medida que el ingreso aumenta.
Comportamiento esperado: El modelo refleja un comportamiento social razonable: las personas con menos ingreso tienen mayor probabilidad de tener dependientes, posiblemente porque aún conviven con familiares o no pueden delegar responsabilidades.
Aplicación modelo de regresion logistica multiple
Queremos predecir el tipo de transmisión (Pc) a partir
de características de las peronas.
# Convertimos 'Pc' a factor
dato$PC <- factor(dato$Pc, levels = c(0, 1), labels = c("Ninguno", "Uno o mas"))
# Ajustamos el modelo
modelo_logit <- glm(Pc ~ y + z + x, data = dato, family = binomial)
summary(modelo_logit)##
## Call:
## glm(formula = Pc ~ y + z + x, family = binomial, data = dato)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.481e+00 1.714e+00 0.864 0.387691
## y 8.297e-05 1.677e-05 4.948 7.48e-07 ***
## z -1.079e-04 2.362e-05 -4.569 4.90e-06 ***
## x -8.248e-06 2.230e-06 -3.698 0.000217 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 165.324 on 149 degrees of freedom
## Residual deviance: 54.615 on 146 degrees of freedom
## AIC: 62.615
##
## Number of Fisher Scoring iterations: 8El modelo predice la probabilidad de que una persona tenga uno o mas personas a cargos (Pc = 1) en función del Gasto Familiar (sinarriendo) (y), Gastos Arriendo (z) y Ingreso Neto (x).
Interpretación de coeficientes
exp(coef(modelo_logit)) # Odds ratios## (Intercept) y z x
## 4.3960780 1.0000830 0.9998921 0.9999918Los coeficientes del modelo logístico representan log-odds.
Para facilitar su interpretación, se puede aplicar exp() para obtener los odds ratios:
Un valor > 1 indica mayor probabilidad del evento (transmisión manual).
Un valor < 1 indica menor probabilidad del evento.
Evaluación del modelo
# Probabilidades predichas
prob <- predict(modelo_logit, type = "response")
# Clasificación con umbral 0.5
pred_clase <- ifelse(prob > 0.5, "uno o mas", "Ninguno")
table(Predicho = pred_clase, Real = dato$Pc)## Real
## Predicho Ninguno Uno o mas
## Ninguno 28 5
## uno o mas 8 109Esto genera una matriz de confusión, que compara las clases predichas vs. las reales.
#Conclusión
Este modelo logístico constituye una herramienta útil para estimar la probabilidad de que una persona tenga personas a su cargo, basándose en factores económicos. La interpretación a través de los odds ratios facilita la comprensión del efecto de cada variable. No obstante, para confirmar la validez del modelo y su capacidad predictiva, sería recomendable complementarlo con otras métricas de rendimiento y métodos de validación adicionales.
Interpretación del Modelo de Regresión Logística
Residuos de la desviación
# Ajustar el modelo de regresión logística
modelo <- glm(Pc ~ y + z, data = dato, family = binomial)
# Obtener residuos en escala logit (link)
residuos_logit <- residuals(modelo, type = "deviance")
# Resumen de los residuos
summary(residuos_logit)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -3.1440530 0.0005152 0.1250150 0.0637630 0.3867104 2.0615605Los residuos muestran las diferencias entre los valores observados y los predichos en escala logit. Un resumen de los residuos:
Valores cercanos a cero indican un buen ajuste del modelo.
\[ \begin{array}{lcccc} \textbf{Variable} & \text{Estimación} & \text{Error estándar} & z\text{-valor} & \text{Pr}(>|z|) \\\\ \text{Intercepto} & 1.481e & 1.714e & 0.864 & 0.387691 \\\\ \text{y} & 8.297e^{-05} & 1.677e^{-05} & 4.948 & 7.48e^{-05} *** \\\\ \text{z} & -1.079e^{04} & 2.362e^{-05} & -4.569 & 4.90e^{-06} *** \\\\ \text{x} & -8.248e^{-06} & 2.230e^{-06} & -3.698 & 0.000217 \\\\ \end{array} \]
Una gran reducción indica que los predictores mejoran mucho el ajuste del modelo.
AIC
El criterio de información de Akaike (AIC) del modelo es:
\[ \text{AIC} = 62.615 \]
Un AIC más bajo indica mejor ajuste entre modelos comparables.
{Iteraciones de Fisher}
El modelo convergió en 10 iteraciones del algoritmo de puntuación de Fisher, indicando que alcanzó una solución estable.
Evaluación del modelo
Evaluamos qué tan bien el modelo clasifica la transmisión manual vs
automática usando una matriz de confusión.
También calculamos medidas de desempeño como:
Precisión (Accuracy): proporción de predicciones correctas.
Sensibilidad (Recall): proporción de casos positivos bien clasificados (manuales).
Especificidad: proporción de casos negativos bien clasificados (automáticos).
F1-score: balance entre precisión y sensibilidad para la clase positiva.
Evaluación manual del modelo
A partir de la matriz de confusión, podemos calcular manualmente:# Asegurarse de que las clases estén bien definidas como factores
real <- factor(dato$Pc, levels = c(0, 1), labels = c("Ninguno", "Uno o mas"))
pred <- factor(pred_clase, levels = c("Ninguno", "Uno o mas"))
# Crear matriz de confusión
conf <- table(Predicho = pred, Real = real)
conf## Real
## Predicho Ninguno Uno o mas
## Ninguno 0 0
## Uno o mas 0 0#MAnualmente
# Extraer los valores manualmente
TN <- conf["Ninguno", "Ninguno"]
TP <- conf["Uno o mas", "Uno o mas"]
FP <- conf["Uno o mas", "Ninguno"]
FN <- conf["Ninguno", "Uno o mas"]
# Calcular métricas
accuracy <- (TP + TN) / sum(conf)
recall <- TP / (TP + FN) # Sensibilidad
specificity <- TN / (TN + FP)
precision <- TP / (TP + FP)
f1 <- 2 * (precision * recall) / (precision + recall)
# Mostrar resultados
data.frame(
Precisión = round(accuracy, 3),
Sensibilidad = round(recall, 3),
Especificidad = round(specificity, 3),
`F1-score` = round(f1, 3)
)## Precisión Sensibilidad Especificidad F1.score
## 1 NaN NaN NaN NaNVisualización de la matriz de confusión
Visualizamos la matriz de confusión en forma de gráfico de calor para observar visualmente cuántas observaciones fueron correctamente o incorrectamente clasificadas.
library## function (package, help, pos = 2, lib.loc = NULL, character.only = FALSE,
## logical.return = FALSE, warn.conflicts, quietly = FALSE,
## verbose = getOption("verbose"), mask.ok, exclude, include.only,
## attach.required = missing(include.only))
## {
## conf.ctrl <- getOption("conflicts.policy")
## if (is.character(conf.ctrl))
## conf.ctrl <- switch(conf.ctrl, strict = list(error = TRUE,
## warn = FALSE), depends.ok = list(error = TRUE, generics.ok = TRUE,
## can.mask = c("base", "methods", "utils", "grDevices",
## "graphics", "stats"), depends.ok = TRUE), warning(gettextf("unknown conflict policy: %s",
## sQuote(conf.ctrl)), call. = FALSE, domain = NA))
## if (!is.list(conf.ctrl))
## conf.ctrl <- NULL
## stopOnConflict <- isTRUE(conf.ctrl$error)
## if (missing(warn.conflicts))
## warn.conflicts <- !isFALSE(conf.ctrl$warn)
## if (!missing(include.only) && !missing(exclude))
## stop("only one of 'include.only' and 'exclude' can be used",
## call. = FALSE)
## testRversion <- function(pkgInfo, pkgname, pkgpath) {
## if (is.null(built <- pkgInfo$Built))
## stop(gettextf("package %s has not been installed properly\n",
## sQuote(pkgname)), call. = FALSE, domain = NA)
## R_version_built_under <- as.numeric_version(built$R)
## if (R_version_built_under < "3.0.0")
## stop(gettextf("package %s was built before R 3.0.0: please re-install it",
## sQuote(pkgname)), call. = FALSE, domain = NA)
## current <- getRversion()
## if (length(Rdeps <- pkgInfo$Rdepends2)) {
## for (dep in Rdeps) if (length(dep) > 1L) {
## target <- dep$version
## res <- do.call(dep$op, if (is.character(target))
## list(as.numeric(R.version[["svn rev"]]), as.numeric(sub("^r",
## "", target)))
## else list(current, as.numeric_version(target)))
## if (!res)
## stop(gettextf("This is R %s, package %s needs %s %s",
## current, sQuote(pkgname), dep$op, target),
## call. = FALSE, domain = NA)
## }
## }
## if (R_version_built_under > current)
## warning(gettextf("package %s was built under R version %s",
## sQuote(pkgname), as.character(built$R)), call. = FALSE,
## domain = NA)
## platform <- built$Platform
## r_arch <- .Platform$r_arch
## if (.Platform$OS.type == "unix") {
## }
## else {
## if (nzchar(platform) && !grepl("mingw", platform))
## stop(gettextf("package %s was built for %s",
## sQuote(pkgname), platform), call. = FALSE,
## domain = NA)
## }
## if (nzchar(r_arch) && file.exists(file.path(pkgpath,
## "libs")) && !file.exists(file.path(pkgpath, "libs",
## r_arch)))
## stop(gettextf("package %s is not installed for 'arch = %s'",
## sQuote(pkgname), r_arch), call. = FALSE, domain = NA)
## }
## checkNoGenerics <- function(env, pkg) {
## nenv <- env
## ns <- .getNamespace(as.name(pkg))
## if (!is.null(ns))
## nenv <- asNamespace(ns)
## if (exists(".noGenerics", envir = nenv, inherits = FALSE))
## TRUE
## else {
## !any(startsWith(names(env), ".__T"))
## }
## }
## checkConflicts <- function(package, pkgname, pkgpath, nogenerics,
## env) {
## dont.mind <- c("last.dump", "last.warning", ".Last.value",
## ".Random.seed", ".Last.lib", ".onDetach", ".packageName",
## ".noGenerics", ".required", ".no_S3_generics", ".Depends",
## ".requireCachedGenerics")
## sp <- search()
## lib.pos <- which(sp == pkgname)
## ob <- names(as.environment(lib.pos))
## if (!nogenerics) {
## these <- ob[startsWith(ob, ".__T__")]
## gen <- gsub(".__T__(.*):([^:]+)", "\\1", these)
## from <- gsub(".__T__(.*):([^:]+)", "\\2", these)
## gen <- gen[from != package]
## ob <- ob[!(ob %in% gen)]
## }
## ipos <- seq_along(sp)[-c(lib.pos, match(c("Autoloads",
## "CheckExEnv"), sp, 0L))]
## cpos <- NULL
## conflicts <- vector("list", 0)
## for (i in ipos) {
## obj.same <- match(names(as.environment(i)), ob, nomatch = 0L)
## if (any(obj.same > 0L)) {
## same <- ob[obj.same]
## same <- same[!(same %in% dont.mind)]
## Classobjs <- which(startsWith(same, ".__"))
## if (length(Classobjs))
## same <- same[-Classobjs]
## same.isFn <- function(where) vapply(same, exists,
## NA, where = where, mode = "function", inherits = FALSE)
## same <- same[same.isFn(i) == same.isFn(lib.pos)]
## not.Ident <- function(ch, TRAFO = identity, ...) vapply(ch,
## function(.) !identical(TRAFO(get(., i)), TRAFO(get(.,
## lib.pos)), ...), NA)
## if (length(same))
## same <- same[not.Ident(same)]
## if (length(same) && identical(sp[i], "package:base"))
## same <- same[not.Ident(same, ignore.environment = TRUE)]
## if (length(same)) {
## conflicts[[sp[i]]] <- same
## cpos[sp[i]] <- i
## }
## }
## }
## if (length(conflicts)) {
## if (stopOnConflict) {
## emsg <- ""
## pkg <- names(conflicts)
## notOK <- vector("list", 0)
## for (i in seq_along(conflicts)) {
## pkgname <- sub("^package:", "", pkg[i])
## if (pkgname %in% canMaskEnv$canMask)
## next
## same <- conflicts[[i]]
## if (is.list(mask.ok))
## myMaskOK <- mask.ok[[pkgname]]
## else myMaskOK <- mask.ok
## if (isTRUE(myMaskOK))
## same <- NULL
## else if (is.character(myMaskOK))
## same <- setdiff(same, myMaskOK)
## if (length(same)) {
## notOK[[pkg[i]]] <- same
## msg <- .maskedMsg(sort(same), pkg = sQuote(pkg[i]),
## by = cpos[i] < lib.pos)
## emsg <- paste(emsg, msg, sep = "\n")
## }
## }
## if (length(notOK)) {
## msg <- gettextf("Conflicts attaching package %s:\n%s",
## sQuote(package), emsg)
## stop(errorCondition(msg, package = package,
## conflicts = conflicts, class = "packageConflictError"))
## }
## }
## if (warn.conflicts) {
## packageStartupMessage(gettextf("\nAttaching package: %s\n",
## sQuote(package)), domain = NA)
## pkg <- names(conflicts)
## for (i in seq_along(conflicts)) {
## msg <- .maskedMsg(sort(conflicts[[i]]), pkg = sQuote(pkg[i]),
## by = cpos[i] < lib.pos)
## packageStartupMessage(msg, domain = NA)
## }
## }
## }
## }
## if (verbose && quietly)
## message("'verbose' and 'quietly' are both true; being verbose then ..")
## if (!missing(package)) {
## if (is.null(lib.loc))
## lib.loc <- .libPaths()
## lib.loc <- lib.loc[dir.exists(lib.loc)]
## if (!character.only)
## package <- as.character(substitute(package))
## if (length(package) != 1L)
## stop("'package' must be of length 1")
## if (is.na(package) || (package == ""))
## stop("invalid package name")
## pkgname <- paste0("package:", package)
## newpackage <- is.na(match(pkgname, search()))
## if (newpackage) {
## pkgpath <- find.package(package, lib.loc, quiet = TRUE,
## verbose = verbose)
## if (length(pkgpath) == 0L) {
## if (length(lib.loc) && !logical.return)
## stop(packageNotFoundError(package, lib.loc,
## sys.call()))
## txt <- if (length(lib.loc))
## gettextf("there is no package called %s", sQuote(package))
## else gettext("no library trees found in 'lib.loc'")
## if (logical.return) {
## if (!quietly)
## warning(txt, domain = NA)
## return(FALSE)
## }
## else stop(txt, domain = NA)
## }
## which.lib.loc <- normalizePath(dirname(pkgpath),
## "/", TRUE)
## pfile <- system.file("Meta", "package.rds", package = package,
## lib.loc = which.lib.loc)
## if (!nzchar(pfile))
## stop(gettextf("%s is not a valid installed package",
## sQuote(package)), domain = NA)
## pkgInfo <- readRDS(pfile)
## testRversion(pkgInfo, package, pkgpath)
## if (is.character(pos)) {
## npos <- match(pos, search())
## if (is.na(npos)) {
## warning(gettextf("%s not found on search path, using pos = 2",
## sQuote(pos)), domain = NA)
## pos <- 2
## }
## else pos <- npos
## }
## deps <- unique(names(pkgInfo$Depends))
## depsOK <- isTRUE(conf.ctrl$depends.ok)
## if (depsOK) {
## canMaskEnv <- dynGet("__library_can_mask__",
## NULL)
## if (is.null(canMaskEnv)) {
## canMaskEnv <- new.env()
## canMaskEnv$canMask <- union("base", conf.ctrl$can.mask)
## "__library_can_mask__" <- canMaskEnv
## }
## canMaskEnv$canMask <- unique(c(package, deps,
## canMaskEnv$canMask))
## }
## else canMaskEnv <- NULL
## if (attach.required)
## .getRequiredPackages2(pkgInfo, quietly = quietly,
## lib.loc = c(lib.loc, .libPaths()))
## cr <- conflictRules(package)
## if (missing(mask.ok))
## mask.ok <- cr$mask.ok
## if (missing(exclude))
## exclude <- cr$exclude
## if (isNamespaceLoaded(package)) {
## newversion <- as.numeric_version(pkgInfo$DESCRIPTION["Version"])
## oldversion <- as.numeric_version(getNamespaceVersion(package))
## if (newversion != oldversion) {
## tryCatch(unloadNamespace(package), error = function(e) {
## P <- if (!is.null(cc <- conditionCall(e)))
## paste("Error in", deparse(cc)[1L], ": ")
## else "Error : "
## stop(gettextf("Package %s version %s cannot be unloaded:\n %s",
## sQuote(package), oldversion, paste0(P,
## conditionMessage(e), "\n")), domain = NA)
## })
## }
## }
## tt <- tryCatch({
## attr(package, "LibPath") <- which.lib.loc
## ns <- loadNamespace(package, lib.loc)
## env <- attachNamespace(ns, pos = pos, deps, exclude,
## include.only)
## }, error = function(e) {
## P <- if (!is.null(cc <- conditionCall(e)))
## paste(" in", deparse(cc)[1L])
## else ""
## msg <- gettextf("package or namespace load failed for %s%s:\n %s",
## sQuote(package), P, conditionMessage(e))
## if (logical.return && !quietly)
## message(paste("Error:", msg), domain = NA)
## else stop(msg, call. = FALSE, domain = NA)
## })
## if (logical.return && is.null(tt))
## return(FALSE)
## attr(package, "LibPath") <- NULL
## {
## on.exit(detach(pos = pos))
## nogenerics <- !.isMethodsDispatchOn() || checkNoGenerics(env,
## package)
## if (isFALSE(conf.ctrl$generics.ok) || (stopOnConflict &&
## !isTRUE(conf.ctrl$generics.ok)))
## nogenerics <- TRUE
## if (stopOnConflict || (warn.conflicts && !exists(".conflicts.OK",
## envir = env, inherits = FALSE)))
## checkConflicts(package, pkgname, pkgpath, nogenerics,
## ns)
## on.exit()
## if (logical.return)
## return(TRUE)
## else return(invisible(.packages()))
## }
## }
## if (verbose && !newpackage)
## warning(gettextf("package %s already present in search()",
## sQuote(package)), domain = NA)
## }
## else if (!missing(help)) {
## if (!character.only)
## help <- as.character(substitute(help))
## pkgName <- help[1L]
## pkgPath <- find.package(pkgName, lib.loc, verbose = verbose)
## docFiles <- c(file.path(pkgPath, "Meta", "package.rds"),
## file.path(pkgPath, "INDEX"))
## if (file.exists(vignetteIndexRDS <- file.path(pkgPath,
## "Meta", "vignette.rds")))
## docFiles <- c(docFiles, vignetteIndexRDS)
## pkgInfo <- vector("list", 3L)
## readDocFile <- function(f) {
## if (basename(f) %in% "package.rds") {
## txt <- readRDS(f)$DESCRIPTION
## if ("Encoding" %in% names(txt)) {
## to <- if (Sys.getlocale("LC_CTYPE") == "C")
## "ASCII//TRANSLIT"
## else ""
## tmp <- try(iconv(txt, from = txt["Encoding"],
## to = to))
## if (!inherits(tmp, "try-error"))
## txt <- tmp
## else warning("'DESCRIPTION' has an 'Encoding' field and re-encoding is not possible",
## call. = FALSE)
## }
## nm <- paste0(names(txt), ":")
## formatDL(nm, txt, indent = max(nchar(nm, "w")) +
## 3L)
## }
## else if (basename(f) %in% "vignette.rds") {
## txt <- readRDS(f)
## if (is.data.frame(txt) && nrow(txt))
## cbind(basename(gsub("\\.[[:alpha:]]+$", "",
## txt$File)), paste(txt$Title, paste0(rep.int("(source",
## NROW(txt)), ifelse(nzchar(txt$PDF), ", pdf",
## ""), ")")))
## else NULL
## }
## else readLines(f)
## }
## for (i in which(file.exists(docFiles))) pkgInfo[[i]] <- readDocFile(docFiles[i])
## y <- list(name = pkgName, path = pkgPath, info = pkgInfo)
## class(y) <- "packageInfo"
## return(y)
## }
## else {
## if (is.null(lib.loc))
## lib.loc <- .libPaths()
## db <- matrix(character(), nrow = 0L, ncol = 3L)
## nopkgs <- character()
## for (lib in lib.loc) {
## a <- .packages(all.available = TRUE, lib.loc = lib)
## for (i in sort(a)) {
## file <- system.file("Meta", "package.rds", package = i,
## lib.loc = lib)
## title <- if (nzchar(file)) {
## txt <- readRDS(file)
## if (is.list(txt))
## txt <- txt$DESCRIPTION
## if ("Encoding" %in% names(txt)) {
## to <- if (Sys.getlocale("LC_CTYPE") == "C")
## "ASCII//TRANSLIT"
## else ""
## tmp <- try(iconv(txt, txt["Encoding"], to,
## "?"))
## if (!inherits(tmp, "try-error"))
## txt <- tmp
## else warning("'DESCRIPTION' has an 'Encoding' field and re-encoding is not possible",
## call. = FALSE)
## }
## txt["Title"]
## }
## else NA
## if (is.na(title))
## title <- " ** No title available ** "
## db <- rbind(db, cbind(i, lib, title))
## }
## if (length(a) == 0L)
## nopkgs <- c(nopkgs, lib)
## }
## dimnames(db) <- list(NULL, c("Package", "LibPath", "Title"))
## if (length(nopkgs) && !missing(lib.loc)) {
## pkglist <- paste(sQuote(nopkgs), collapse = ", ")
## msg <- sprintf(ngettext(length(nopkgs), "library %s contains no packages",
## "libraries %s contain no packages"), pkglist)
## warning(msg, domain = NA)
## }
## y <- list(header = NULL, results = db, footer = NULL)
## class(y) <- "libraryIQR"
## return(y)
## }
## if (logical.return)
## TRUE
## else invisible(.packages())
## }
## <bytecode: 0x000002434876eca0>
## <environment: namespace:base># Asegúrate de estar usando las predicciones del modelo_logit
prob <- predict(modelo_logit, type = "response")
pred_clase <- ifelse(prob > 0.5, "Uno o mas", "Ninguno")
# Asegurarse de que las clases estén bien definidas como factores
real <- factor(dato$Pc, levels = c("Ninguno", "Uno o mas"))
pred <- factor(pred_clase, levels = c("Ninguno", "Uno o mas"))
# Crear matriz de confusión
conf <- table(Predicho = pred, Real = real)
conf## Real
## Predicho Ninguno Uno o mas
## Ninguno 28 5
## Uno o mas 8 109# Matriz de confusión como tabla
conf <- table(Predicho = pred, Real = real)
# Convertir a data frame para graficar
conf_df <- as.data.frame(conf)
names(conf_df) <- c("Predicho", "Real", "Frecuencia")
# Gráfico de calor
ggplot(conf_df, aes(x = Real, y = Predicho, fill = Frecuencia)) +
geom_tile(color = "white") +
geom_text(aes(label = Frecuencia), size = 6, color = "black") +
scale_fill_gradient(low = "white", high = "steelblue") +
labs(
title = "Matriz de confusión",
x = "Clase real",
y = "Clase predicha"
) +
theme_minimal()
Las diagonales representan predicciones correctas.
Los fuera de la diagonal son errores.
El color más oscuro indica mayor frecuencia.
Tu modelo predijo todo como “Uno o mas”.
Nunca predijo “Ninguno”, por eso las filas de “Ninguno” están en cero.
¿Qué significa esto?
Falso positivo (FP): 1 (se predijo “Uno o mas”, pero era “Ninguno”)
Verdadero positivo (TP): 12
Verdadero negativo (TN): 0
Falso negativo (FN): 0
Esto puede pasar si el modelo está muy sesgado hacia una clase, probablemente porque:
Los predictores (y, z, x, etc.) en el modelo múltiple hacen que la probabilidad predicha siempre esté por encima de 0.5.
Entonces todas las observaciones se clasifican como “Uno o mas”.
¿Cómo verificarlo? Revisar la distribución de probabilidades predichas:
summary(prob)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00103 0.62571 0.97047 0.76000 0.99905 1.00000hist(prob, breaks = 10, col = "lightblue", main = "Distribución de probabilidades predichas", xlab = "Probabilidad de tener uno o mas personas a cargo")
La mayoría de las probabilidades predichas están cercanas a 0 o a 1.
Hay un grupo grande con valores cercanos a 0, y otro grupo cercano a 1.
Casi ninguna predicción está cerca del umbral 0.5, lo cual explica por qué un corte fijo en 0.5 puede no funcionar bien para clasificar adecuadamente.
¿Qué hacer ahora? Puedes ajustar el umbral de clasificación para mejorar la matriz de confusión y evitar que todo se clasifique como “Uno o mas”.
# Intentar varios umbrales de clasificación
umbrales <- seq(0.1, 0.9, by = 0.1)
resultados <- data.frame()
for (corte in umbrales) {
pred_clase <- ifelse(prob > corte, "Uno o mas", "Ninguno")
real <- factor(dato$Pc, levels = c("Ninguno", "Uno o mas"))
pred <- factor(pred_clase, levels = c("Ninguno", "Uno o mas"))
conf <- table(Predicho = pred, Real = real)
TN <- conf["Ninguno", "Ninguno"]
TP <- conf["Uno o mas", "Uno o mas"]
FP <- conf["Uno o mas", "Ninguno"]
FN <- conf["Ninguno", "Uno o mas"]
acc <- (TP + TN) / sum(conf)
sen <- TP / (TP + FN)
esp <- TN / (TN + FP)
f1 <- if (!is.na(sen) && !is.na(TP + FP) && (TP + FP) > 0 && (sen + TP/(TP + FP)) > 0) {
2 * ((TP / (TP + FP)) * sen) / ((TP / (TP + FP)) + sen)
} else { NA }
resultados <- rbind(resultados, data.frame(
Umbral = corte,
Accuracy = round(acc, 3),
Sensibilidad = round(sen, 3),
Especificidad = round(esp, 3),
F1 = round(f1, 3)
))
}
# Mostrar la tabla de resultados
resultados## Umbral Accuracy Sensibilidad Especificidad F1
## 1 0.1 0.880 0.991 0.528 0.926
## 2 0.2 0.900 0.982 0.639 0.937
## 3 0.3 0.907 0.974 0.694 0.941
## 4 0.4 0.907 0.974 0.694 0.941
## 5 0.5 0.913 0.956 0.778 0.944
## 6 0.6 0.907 0.939 0.806 0.939
## 7 0.7 0.900 0.904 0.889 0.932
## 8 0.8 0.907 0.895 0.944 0.936
## 9 0.9 0.853 0.816 0.972 0.894# Elegir un nuevo umbral basado en los resultados
mejor_umbral <- 0.6
pred_clase <- ifelse(prob > mejor_umbral, "Uno o mas", "Ninguno")library(ggplot2)
ggplot(resultados, aes(x = Umbral)) +
geom_line(aes(y = F1, color = "F1-score"), size = 1.2) +
geom_line(aes(y = Accuracy, color = "Precisión"), size = 1.2, linetype = "dashed") +
geom_line(aes(y = Sensibilidad, color = "Sensibilidad"), size = 1.2, linetype = "dotdash") +
geom_line(aes(y = Especificidad, color = "Especificidad"), size = 1.2, linetype = "twodash") +
labs(title = "Métricas de desempeño según umbral de clasificación",
y = "Valor de la métrica",
x = "Umbral de clasificación",
color = "Métrica") +
theme_minimal()## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
#conclusiones El modelo de regresión logística desarrollado para
predecir si una persona tiene “Uno o más” dependientes a cargo mostró
inicialmente un sesgo hacia la clase positiva, clasificando todas las
observaciones como “Uno o más”. Esto se evidenció en la matriz de
confusión, donde no se registraron verdaderos negativos (TN) ni falsos
negativos (FN), lo que implica una falta total de sensibilidad hacia la
clase “Ninguno”.
Esta situación se debe a que las probabilidades predichas por el modelo estaban mayoritariamente por encima del umbral estándar de 0.5, lo cual llevó a clasificar todos los casos en la clase mayoritaria. Al analizar la distribución de las probabilidades predichas, se observó que estaban agrupadas en los extremos (cerca de 0 o 1), lo que sugiere que un umbral fijo no es adecuado para este caso.
Como solución, se realizó un análisis sistemático probando diferentes umbrales de clasificación (desde 0.1 hasta 0.9), evaluando en cada caso métricas como precisión, sensibilidad, especificidad y F1-score. Esto permitió identificar umbrales alternativos que logran un mejor equilibrio entre la correcta clasificación de ambas clases.
Finalmente, se concluye que ajustar el umbral de clasificación mejora sustancialmente el desempeño del modelo, corrigiendo el sesgo inicial y logrando un balance más adecuado entre sensibilidad y especificidad. Esta estrategia es especialmente útil cuando las clases están desbalanceadas o las predicciones son extremas, como ocurre en este caso.