Introducción

El consumo de cigarrillos ha sido objeto de amplio análisis en la economía de la salud, dado su impacto tanto en los individuos como en la sociedad. Comprender los factores que determinan el número de cigarrillos fumados por día es fundamental para el diseño de políticas públicas orientadas a reducir el tabaquismo y sus costos asociados.

El presente estudio busca responder la siguiente pregunta económica:

¿Qué factores socioeconómicos y demográficos influyen en la cantidad de cigarrillos que una persona fuma por día?

Para abordar esta cuestión, se utiliza la base de datos smoke del paquete wooldridge en R, que contiene información sobre el comportamiento de consumo de cigarrillos y características individuales. En particular, el modelo propuesto considera como variable dependiente el número de cigarrillos fumados por día (cigs), y como variables explicativas el nivel de ingreso (income), la edad del individuo (age), el precio de los cigarrillos (lcigpric) y el nivel educativo (educ).

\(cigs_i = \beta_0 + \beta_1 income_i + \beta_2 age_i + \beta_3 lcigpric_i + \beta_4 educ_i + u_i\).

A través de este modelo econométrico, se busca analizar cómo factores económicos, como el ingreso y el precio, junto con características personales como la edad y la educación, afectan el consumo de cigarrillos. Los resultados permitirán obtener evidencia empírica útil para comprender el comportamiento del fumador promedio y orientar políticas de salud pública más efectivas.

1. Preparación del entorno

# install.packages(c("wooldridge", "ggpubr", "performance", "car", "broom", "dplyr", "lmtest", "sandwich", "tidyverse"))

# Cargar librerías
library(wooldridge)
library(tidyverse)
library(car)
library(ggplot2)
library(broom)
library(dplyr)
library(ggpubr)
library(performance)
library(lmtest)
library(sandwich)



data("smoke", package = "wooldridge")
df <- as_tibble(smoke)
str(df)
## tibble [807 × 10] (S3: tbl_df/tbl/data.frame)
##  $ educ    : num [1:807] 16 16 12 13.5 10 6 12 15 12 12 ...
##  $ cigpric : num [1:807] 60.5 57.9 57.7 57.9 58.3 ...
##  $ white   : int [1:807] 1 1 1 1 1 1 1 1 1 1 ...
##  $ age     : int [1:807] 46 40 58 30 17 86 35 48 48 31 ...
##  $ income  : int [1:807] 20000 30000 30000 20000 20000 6500 20000 30000 20000 20000 ...
##  $ cigs    : int [1:807] 0 0 3 0 0 0 0 0 0 0 ...
##  $ restaurn: int [1:807] 0 0 0 0 0 0 0 0 0 0 ...
##  $ lincome : num [1:807] 9.9 10.3 10.3 9.9 9.9 ...
##  $ agesq   : int [1:807] 2116 1600 3364 900 289 7396 1225 2304 2304 961 ...
##  $ lcigpric: num [1:807] 4.1 4.06 4.05 4.06 4.07 ...
##  - attr(*, "time.stamp")= chr "25 Jun 2011 23:03"
summary(df)
##       educ          cigpric          white             age       
##  Min.   : 6.00   Min.   :44.00   Min.   :0.0000   Min.   :17.00  
##  1st Qu.:10.00   1st Qu.:58.14   1st Qu.:1.0000   1st Qu.:28.00  
##  Median :12.00   Median :61.05   Median :1.0000   Median :38.00  
##  Mean   :12.47   Mean   :60.30   Mean   :0.8786   Mean   :41.24  
##  3rd Qu.:13.50   3rd Qu.:63.18   3rd Qu.:1.0000   3rd Qu.:54.00  
##  Max.   :18.00   Max.   :70.13   Max.   :1.0000   Max.   :88.00  
##      income           cigs           restaurn         lincome      
##  Min.   :  500   Min.   : 0.000   Min.   :0.0000   Min.   : 6.215  
##  1st Qu.:12500   1st Qu.: 0.000   1st Qu.:0.0000   1st Qu.: 9.433  
##  Median :20000   Median : 0.000   Median :0.0000   Median : 9.903  
##  Mean   :19305   Mean   : 8.686   Mean   :0.2466   Mean   : 9.687  
##  3rd Qu.:30000   3rd Qu.:20.000   3rd Qu.:0.0000   3rd Qu.:10.309  
##  Max.   :30000   Max.   :80.000   Max.   :1.0000   Max.   :10.309  
##      agesq         lcigpric    
##  Min.   : 289   Min.   :3.784  
##  1st Qu.: 784   1st Qu.:4.063  
##  Median :1444   Median :4.112  
##  Mean   :1990   Mean   :4.096  
##  3rd Qu.:2916   3rd Qu.:4.146  
##  Max.   :7744   Max.   :4.250
#VERIFICAR NA Y CANTIDAD DE NA
any(is.na(df))
## [1] FALSE
sum(is.na(df))
## [1] 0
#VERIFICAS SI HAY VARIABLES CATEGORICAS
str(df)
## tibble [807 × 10] (S3: tbl_df/tbl/data.frame)
##  $ educ    : num [1:807] 16 16 12 13.5 10 6 12 15 12 12 ...
##  $ cigpric : num [1:807] 60.5 57.9 57.7 57.9 58.3 ...
##  $ white   : int [1:807] 1 1 1 1 1 1 1 1 1 1 ...
##  $ age     : int [1:807] 46 40 58 30 17 86 35 48 48 31 ...
##  $ income  : int [1:807] 20000 30000 30000 20000 20000 6500 20000 30000 20000 20000 ...
##  $ cigs    : int [1:807] 0 0 3 0 0 0 0 0 0 0 ...
##  $ restaurn: int [1:807] 0 0 0 0 0 0 0 0 0 0 ...
##  $ lincome : num [1:807] 9.9 10.3 10.3 9.9 9.9 ...
##  $ agesq   : int [1:807] 2116 1600 3364 900 289 7396 1225 2304 2304 961 ...
##  $ lcigpric: num [1:807] 4.1 4.06 4.05 4.06 4.07 ...
##  - attr(*, "time.stamp")= chr "25 Jun 2011 23:03"
# --- Definir variable dependiente (respuesta) ---
y_var <- "cigs"   # número de cigarrillos fumados por día

# --- Definir variables explicativas ---
x_vars <- c("lincome", "educ", "lcigpric", "age")

# --- Verificar ---
y_var
## [1] "cigs"
x_vars
## [1] "lincome"  "educ"     "lcigpric" "age"

2. Estimacion del modelo OLS

fml <- as.formula(paste(y_var, "~", paste(x_vars, collapse = " + ")))
m_ols <- lm(fml, data = df)

summary(m_ols)
## 
## Call:
## lm(formula = fml, data = df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -12.703  -9.027  -7.317   9.283  72.471 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)  
## (Intercept) 10.00443   24.26632   0.412   0.6802  
## lincome      1.77738    0.71323   2.492   0.0129 *
## educ        -0.38882    0.16851  -2.307   0.0213 *
## lcigpric    -2.91212    5.82577  -0.500   0.6173  
## age         -0.04265    0.02876  -1.483   0.1385  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 13.67 on 802 degrees of freedom
## Multiple R-squared:  0.01277,    Adjusted R-squared:  0.007843 
## F-statistic: 2.593 on 4 and 802 DF,  p-value: 0.03541

Resultados del modelo

Ingreso (lincome)
$ beta_1 = 1.777, ; p = 0.0129 $

El ingreso tiene un efecto positivo y estadísticamente significativo al 5%.
Esto quiere decir que, manteniendo todo lo demás igual, cuando el ingreso aumenta un 1%, el número de cigarrillos consumidos por día también aumenta ligeramente (alrededor de 0.0177).
En términos simples, las personas con mayores ingresos tienden a fumar un poco más, posiblemente porque tienen más recursos económicos para sostener ese hábito.

Educación (educ)
$ beta_2 = -0.389, ; p = 0.0213 $

La educación muestra un efecto negativo y significativo.
Cada año adicional de estudio se asocia con una reducción promedio de 0.39 cigarrillos al día.
En otras palabras, las personas con mayor nivel educativo tienden a fumar menos, lo cual coincide con la idea de que la educación aumenta la conciencia sobre los riesgos del tabaquismo y promueve hábitos más saludables.

Precio de los cigarrillos (lcigpric)
$ beta_3 = -2.912, ; p = 0.617 $

El coeficiente del precio es negativo, lo que sugiere que cuando los cigarrillos son más caros, el consumo tiende a disminuir.
Sin embargo, el resultado no es estadísticamente significativo, por lo que no se puede afirmar con certeza que el precio tenga un efecto claro sobre el consumo en esta muestra.

Edad (age)
$ beta_4 = -0.043, ; p = 0.138 $

La edad presenta un coeficiente negativo, lo que sugiere que, en promedio, las personas mayores tienden a fumar ligeramente menos que las más jóvenes.
Sin embargo, este efecto no es estadísticamente significativo, por lo que no se puede afirmar con certeza que la edad influya directamente en el consumo de cigarrillos en esta muestra.

De acuerdo con un estudio publicado en el International Journal of Environmental Research and Public Health, la relación entre edad y consumo de tabaco no es lineal: los fumadores de mediana edad tienden a presentar una mayor dependencia, mientras que los más jóvenes y los mayores muestran niveles más bajos de consumo (Lee et al., 2015).
Esto sugiere que el efecto de la edad sobre el consumo de cigarrillos es más complejo y puede variar según el grupo etario y otros factores sociales.

3. Linealidad (gráficos y pruebas)

cr_pl <- car::crPlots(m_ols)

# Residuos parciales con ggplot (para cada X)
hacer_grafico_pr <- function(modelo, x){
  aug <- augment(modelo)
  beta_x <- coef(modelo)[x]
  if(is.na(beta_x)) stop("La variable no está en el modelo.")
  aug |>
    mutate(partial = .resid + !!sym(x) * beta_x) |>
    ggplot(aes_string(x = x, y = "partial")) +
    geom_point(alpha = 0.6) +
    geom_smooth(method = "loess", se = FALSE, color = "magenta") +
    labs(
      title = paste("Gráfico de residuos parciales:", x),
      y = paste("e +", round(beta_x,4), "*", x),
      x = x
    ) +
    theme_minimal()
}
print(hacer_grafico_pr(m_ols, x_vars[1]))
## Warning: `aes_string()` was deprecated in ggplot2 3.0.0.
## ℹ Please use tidy evaluation idioms with `aes()`.
## ℹ See also `vignette("ggplot2-in-packages")` for more information.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
## `geom_smooth()` using formula = 'y ~ x'
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : pseudoinverse used at 10.329
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : neighborhood radius 0.89594
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : reciprocal condition number 3.8596e-15
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : There are other near singularities as well. 0.2209

print(hacer_grafico_pr(m_ols, x_vars[2]))
## `geom_smooth()` using formula = 'y ~ x'

print(hacer_grafico_pr(m_ols, x_vars[3]))
## `geom_smooth()` using formula = 'y ~ x'

print(hacer_grafico_pr(m_ols, x_vars[4]))
## `geom_smooth()` using formula = 'y ~ x'

lincome

La línea magenta tiene una tendencia creciente y suave, sin curvaturas visibles. La dispersión de los puntos aumenta, pero la relación parece monótona y estable.El supuesto de linealidad se cumple adecuadamente, Ya que lincome es el logaritmo del ingreso, la transformación logarítmica logra linearizar la relación.

edu

El gráfico muestra una línea magenta con ligera curvatura, también en forma de U invertida: el consumo tiende a disminuir con más años de educación, aunque no de forma estrictamente lineal. El efecto de la educación sobre el consumo presenta no linealidad leve, probablemente con rendimientos decrecientes del efecto educativo sobre la reducción del consumo.

lcigpric (log del precio del cigarrillo)

La línea magenta es casi plana, aunque con una leve curvatura en forma de U hacia los extremos. No se aprecia un patrón fuerte ni sistemático.La relación entre el log del precio y el consumo parece aproximadamente lineal. La ligera curvatura puede deberse a dispersión en los datos, no a un patrón real.

age (edad)

En el gráfico se nota que la relación entre la edad y el consumo de cigarrillos no es lineal. A medida que las personas se hacen adultas, el consumo tiende a subir, pero después de cierta edad empieza a bajar. Esto muestra una curva en forma de U invertida.

# RESET (funcionalidad específica)
reset_lm <- lmtest::resettest(m_ols, power = 2:3, type = "fitted")
reset_lm 
## 
##  RESET test
## 
## data:  m_ols
## RESET = 0.50538, df1 = 2, df2 = 800, p-value = 0.6035

Hipótesis nula (H₀)

el modelo está bien especificado (no faltan transformaciones ni variables relevantes).

Hipótesis alternativa (H₁)

el modelo está mal especificado (podrían faltar términos no lineales/interacciones o variables).

Regla

Si p < 0.05 → rechazar H₀ → hay indicios de mala especificación.

Si p ≥ 0.05 → no rechazar H₀ → sin evidencia de error funcional.

Resultado del test (modelo base, sin transformaciones)

RESET = 0.50538

p-valor = 0.6035

Como p = 0.6035 ≥ 0.05, no se rechaza H₀.con los datos disponibles, no hay evidencia de que el modelo lineal inicial esté mal especificado. En términos económicos, la relación entre cigs y los predictores (lincome, educ, lcigpric, age) puede considerarse aproximadamente lineal dentro del rango observado, sin embargo se analizarán los demas supuestos para validar.

Homocedasticidad

#  Breusch-Pagan
bp <- lmtest::bptest(m_ols) # H0: Var(u|X) = sigma^2 (homoced.)
bp
## 
##  studentized Breusch-Pagan test
## 
## data:  m_ols
## BP = 6.4463, df = 4, p-value = 0.1682

Hipótesis nula (H₀)

H₀: Var(u | X) = σ² → los residuos son homocedásticos (la varianza del error es constante).

Hipótesis alternativa (H₁)

H₁: Var(u | X) ≠ σ² → los residuos son heterocedásticos (la varianza del error depende de las variables explicativas).

Regla

Si p < 0.05 → rechazar H₀ → existe evidencia de heterocedasticidad.

Si p ≥ 0.05 → no rechazar H₀ → no hay evidencia de heterocedasticidad (el modelo cumple el supuesto de varianza constante).

Resultado del test (modelo base, sin transformaciones)

RESET = 6.4463

p-valor = 0.1682

Como el p-valor (0.1682) es mayor que 0.05, no se rechaza la hipótesis nula de homocedasticidad. Esto significa que, con base en el test de Breusch–Pagan, no hay evidencia estadística de que los errores del modelo presenten varianza no constante.

En términos económicos, esto sugiere que la dispersión del consumo de cigarrillos no cambia sistemáticamente con el ingreso, la educación, el precio o la edad de las personas. Por lo tanto, el supuesto de varianza constante de los errores se considera razonablemente cumplido en el modelo lineal inicial.

Test de White

#  Test de White  
X <- model.matrix(m_ols)
white_test <- lmtest::bptest(m_ols, ~ fitted(m_ols) + I(fitted(m_ols)^2))
white_test
## 
##  studentized Breusch-Pagan test
## 
## data:  m_ols
## BP = 2.9159, df = 2, p-value = 0.2327

Dado que el p-valor (0.2327) es mayor que 0.05, no se rechaza la hipótesis nula de homocedasticidad. Esto indica que no existe evidencia estadística de que la varianza de los errores dependa de los valores ajustados del modelo o de sus términos no lineales.

En otras palabras, el modelo lineal inicial cumple el supuesto de homocedasticidad, por lo que la dispersión del consumo de cigarrillos no varía sistemáticamente con los niveles de ingreso, educación, edad o precio.

Gráfico residuos vs ajustados

# 6.3 Gráfico residuos vs ajustados
augment(m_ols) |>
  ggplot(aes(.fitted, .resid)) +
  geom_hline(yintercept = 0, linetype = 2, color = "red") +
  geom_point(alpha = 0.6) +
  geom_smooth(method = "loess", se = FALSE, color = "blue") +
  labs(
    title = "Residuos vs Ajustados",
    x = "Valores ajustados (Ŷ)",
    y = "Residuos (ê)"
  )
## `geom_smooth()` using formula = 'y ~ x'

El gráfico no muestra evidencia visual de heterocedasticidad ni de error funcional. Esto sugiere que la varianza de los residuos se mantiene aproximadamente constante a lo largo de los valores ajustados del consumo de cigarrillos.

Normalidad de residuos Shapiro-Wilk

#  Shapiro-Wilk (n <= 5000 aprox.); H0: normalidad
if(nrow(df) <= 5000){
  sh <- shapiro.test(residuals(m_ols))
  print(sh)
}
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(m_ols)
## W = 0.75691, p-value < 2.2e-16

El resultado del test muestra un p-valor menor a 0.05, por lo tanto se rechaza la hipótesis nula de normalidad. Esto quiere decir que los residuos del modelo no siguen una distribución normal. En otras palabras, los errores no se reparten de manera completamente simétrica ni concentrada alrededor de cero; pueden tener colas más pesadas o cierta asimetría.

Desde el punto de vista económico, lo que sugiere que existen personas o grupos cuyo consumo de cigarrillos se desvía más de lo esperado según sus características (ingreso, educación, precio y edad).

Jarque-Bera

jb <- performance::check_normality(m_ols)
jb
## Warning: Non-normality of residuals detected (p < .001).

Como el p-valor es menor a 0.001, se rechaza la hipótesis nula.

En términos de diagnóstico, esto sugiere que el consumo de cigarrillos muestra comportamientos atípicos: hay individuos que fuman mucho más o mucho menos de lo que el modelo predice, lo cual introduce cierta irregularidad en los errores. Desde una perspectiva económica, la falta de normalidad refleja que el consumo de cigarrillos no es uniforme entre personas con características similares. Factores no observados como hábitos personales, condiciones de salud, entorno social o políticas locales de control del tabaco pueden hacer que algunos grupos se alejen del comportamiento promedio.

Grafico QQ

ggqqplot(residuals(m_ols)) + ggtitle("QQ-plot de residuos")

En este caso, los puntos se apartan de la línea de forma notoria, sobre todo en las colas (parte superior e inferior). Esto indica que los residuos tienen valores extremos o asimetrías.

Los errores no se reparten de manera perfectamente simétrica; algunos individuos fuman mucho más o mucho menos de lo que el modelo predice. Esto coincide con lo que mostró el test de Shapiro–Wilk y el Jarque–Bera: los residuos no son normales.

Multicolinealidad

vifs <- car::vif(m_ols) # VIF > 10 (regla práctica) sugiere problema
vifs
##  lincome     educ lcigpric      age 
## 1.114855 1.145114 1.006860 1.034822

Todos los valores de VIF son muy cercanos a 1, lo que indica que las variables explicativas no están correlacionadas entre sí.

Desde el punto de vista económico, esto es una buena señal: el modelo está identificando efectos independientes y claros de cada variable sobre el consumo de cigarrillos.

Condición (número de condición de X)

# 8.1 Condición (número de condición de X)
kappa_X <- kappa(model.matrix(m_ols))
kappa_X
## [1] 2758.374
# 8.2 Estandarizamos las variables explicativas
df_std <- df |> mutate(across(c(educ, lcigpric, age), scale))

# 8.3 Estimamos el modelo nuevamente
m_ols_std <- lm(lincome ~ educ + lcigpric + restaurn, data = df_std)

# 8.4 Calculamos el nuevo número de condición
kappa(model.matrix(m_ols_std))
## [1] 2.158011

Regla práctica

κ < 30 → sin problema.

30 ≤ κ ≤ 1000 → posible multicolinealidad moderada.

κ > 1000 → multicolinealidad severa o problema de escala

Número de condición (kappa_X = 2758.37) en la Primera estimacion

Un valor de κ muy alto indica que algunas variables están altamente correlacionadas o en escalas muy distintas, lo que puede generar inestabilidad numérica en la estimación.

Nuevo modelo y nuevo número de condición (κ = 2.16)

El número de condición del modelo original (κ = 2758.37) sugería un posible problema de escala entre las variables explicativas, más que una multicolinealidad real. Para corregirlo, se realizó una estandarización de las variables (educación, precio y edad), de manera que todas tuvieran media cero y desviación estándar uno.

Tras esta transformación, el número de condición se redujo a 2.16, lo que indica que la matriz de variables explicativas está bien condicionada y que no existen problemas de colinealidad ni inestabilidad numérica en el modelo.

Económicamente, esto significa que cada variable explica una parte distinta del comportamiento del consumo de cigarrillos y que las estimaciones obtenidas son confiables y precisas.

4. Soluciones Econometricas

# ---- 9) Soluciones econométricas si fallan supuestos: forma funcional ----
# Evidencia visual: curvatura en age y educ -> incluir términos cuadráticos

# 9.1 Especificación funcional (modelo final propuesto)
m_fun <- lm(cigs ~ lincome + educ + I(educ^2) + lcigpric + age + I(age^2), data = df)
summary(m_fun)
## 
## Call:
## lm(formula = cigs ~ lincome + educ + I(educ^2) + lcigpric + age + 
##     I(age^2), data = df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -14.142  -9.344  -5.837   7.557  73.774 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -14.377009  24.478747  -0.587 0.557150    
## lincome       0.534626   0.727244   0.735 0.462470    
## educ          2.825826   1.017585   2.777 0.005615 ** 
## I(educ^2)    -0.132965   0.039964  -3.327 0.000918 ***
## lcigpric     -2.664656   5.697764  -0.468 0.640150    
## age           0.826866   0.160222   5.161 3.11e-07 ***
## I(age^2)     -0.009356   0.001739  -5.379 9.85e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 13.37 on 800 degrees of freedom
## Multiple R-squared:  0.05812,    Adjusted R-squared:  0.05106 
## F-statistic: 8.228 on 6 and 800 DF,  p-value: 1.168e-08

Para corregir los problemas detectados de no linealidad en las variables edad y educación, se incluyeron los términos cuadráticos I(educ^2) y I(age^2) en el modelo.

R² ajustado: 0.051 → el modelo explica alrededor del 5% de la variación en el número de cigarrillos fumados. Aunque parece bajo, es común en datos individuales de comportamiento humano, donde influyen muchos factores no observados (gustos, salud, entorno, etc.).

F global (p < 0.001): el conjunto de variables es significativo en su totalidad, es decir, el modelo en general sí explica parte del consumo.

El modelo corregido mejora la especificación funcional: los signos y significancias de los términos cuadráticos confirman que había relaciones no lineales entre educación, edad y el consumo.

Las pendientes positivas de educ y age, combinadas con los coeficientes negativos de sus cuadrados, describen relaciones en forma de U invertida, es decir, las personas con baja educación o muy jóvenes fuman poco, el consumo crece en niveles intermedios,y luego vuelve a disminuir con mayor educación o mayor edad. Esto coincide con la lógica económica y social: a más educación o edad avanzada, suele haber mayor conciencia sobre los riesgos del tabaco.

Test de Reset modelo transformado

# 9.2 Test de especificación (RESET)
# H0: forma funcional correcta (no faltan potencias ni interacciones)
reset_fun <- lmtest::resettest(m_fun, power = 2:3, type = "fitted")
reset_fun
## 
##  RESET test
## 
## data:  m_fun
## RESET = 1.515, df1 = 2, df2 = 798, p-value = 0.2204

Como p = 0.2204 ≥ 0.05, no se rechaza H₀. Es decir, no hay evidencia de que falten transformaciones o interacciones después de incluir educ² y age².

En palabras simples: la forma funcional quedó adecuada; los términos cuadráticos corrigieron la no linealidad que veíamos antes.

Homocedasticidad (Modelo transformado)

# 9.3 Homocedasticidad
# Breusch–Pagan "clásico"
bp_fun <- lmtest::bptest(m_fun) 
bp_fun
## 
##  studentized Breusch-Pagan test
## 
## data:  m_fun
## BP = 25.111, df = 6, p-value = 0.0003257

El resultado indica que la dispersión de los errores cambia según los valores de las variables explicativas. En palabras simples, el modelo explica mejor el comportamiento de unos grupos de personas que de otros. Por ejemplo, es posible que en individuos con mayor ingreso o edad, las predicciones del consumo de cigarrillos sean más variables, lo que genera varianza desigual.

Inferencia Robusta

# 9.4 Inferencia robusta (si hay heterocedasticidad, HC1)
robust_fun <- lmtest::coeftest(m_fun, vcov = sandwich::vcovHC(m_fun, type = "HC1"))
robust_fun
## 
## t test of coefficients:
## 
##                Estimate  Std. Error t value  Pr(>|t|)    
## (Intercept) -14.3770092  25.1127489 -0.5725  0.567145    
## lincome       0.5346256   0.6033383  0.8861  0.375823    
## educ          2.8258260   0.9175552  3.0797  0.002143 ** 
## I(educ^2)    -0.1329647   0.0371086 -3.5831  0.000360 ***
## lcigpric     -2.6646557   5.8765200 -0.4534  0.650354    
## age           0.8268660   0.1410369  5.8628 6.655e-09 ***
## I(age^2)     -0.0093560   0.0014805 -6.3195 4.355e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Dado que el test de Breusch–Pagan detectó heterocedasticidad, se recalcularon los errores estándar de los coeficientes usando el estimador HC1 (White corregido), que ajusta las varianzas sin modificar los coeficientes del modelo. Esto permite que las pruebas t y p-valores sean válidas incluso cuando la varianza de los errores no es constante.

Tras aplicar los errores robustos, las variables educación y edad —junto con sus términos cuadráticos— se mantienen altamente significativas, lo que confirma que la forma funcional del modelo es adecuada. Las variables ingreso y precio siguen sin ser significativas, lo que indica que, dentro de esta muestra, el consumo de cigarrillos depende más de factores demográficos que económicos.

MIRADA ECONOMICA

Desde una mirada económica, la edad y la educación influyen en el consumo de manera no lineal: las personas tienden a fumar más en una etapa intermedia de su vida o nivel educativo, y menos cuando son más mayores o con mayor formación.

Ni el ingreso ni el precio parecen tener un efecto estadísticamente relevante en esta muestra, lo que podría explicarse porque el hábito de fumar está más asociado a factores sociales y de comportamiento que a las restricciones económicas.

Grafico de diagnóstico

# 9.5 Gráficos de diagnóstico
# 9.5.1 Componentes + residuos (parciales) para verificar curvaturas residuales
car::crPlots(m_fun)

Los gráficos muestran que las relaciones no lineales fueron corregidas correctamente.

Los ajustes son suaves y cercanas a la horizontal en la mayoría de las variables, sin patrones sistemáticos de curvatura residual. Esto significa que el modelo final cumple el supuesto de linealidad en los parámetros y representa adecuadamente la relación entre las variables explicativas y el consumo de cigarrillos.

Residuos VS los ajustados

# 9.5.2 Residuos vs ajustados
broom::augment(m_fun) |>
  ggplot2::ggplot(ggplot2::aes(.fitted, .resid)) +
  ggplot2::geom_hline(yintercept = 0, linetype = 2, color = "red") +
  ggplot2::geom_point(alpha = 0.6) +
  ggplot2::geom_smooth(method = "loess", se = FALSE) +
  ggplot2::labs(title = "Residuos vs Ajustados (m_fun)", x = "Ŷ", y = "ê")
## `geom_smooth()` using formula = 'y ~ x'

El hecho de que la línea azul esté relativamente plana y cercana al eje cero indica que el modelo ahora captura bien la tendencia media y no hay evidencia fuerte de error funcional o sesgo.

Normalidad de Residuos

# 9.6 Normalidad de residuos
if (nrow(df) <= 5000) {
  sh_fun <- shapiro.test(residuals(m_fun))
  print(sh_fun)
}
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(m_fun)
## W = 0.82357, p-value < 2.2e-16
jb_fun <- performance::check_normality(m_fun)
jb_fun
## Warning: Non-normality of residuals detected (p < .001).
ggpubr::ggqqplot(residuals(m_fun)) + ggplot2::ggtitle("QQ-plot residuos (m_fun)")

Desde un punto de vista económico, esta falta de normalidad es comprensible: el consumo de cigarrillos no es simétrico en la población —la mayoría fuma poco o nada, mientras que una minoría fuma mucho—. Por tanto, el modelo refleja un fenómeno real, más que un error técnico.

En resumen, aunque los residuos no sean perfectamente normales, el modelo sigue siendo útil y válido para inferencia al haber aplicado errores estándar robustos, lo que garantiza la confiabilidad de las conclusiones.

5 Informe ordenado (tabla de coeficientes)

# ---- 10 Informe ordenado (tabla de coeficientes) ----
tab_ols <- broom::tidy(m_ols, conf.int = TRUE)
tab_rob <- broom::tidy(coeftest(m_ols, vcov = sandwich::vcovHC(m_ols, type="HC1")))
list(
  resumen_ols = broom::glance(m_ols),
  coef_ci_ols = tab_ols,
  coef_se_rob = tab_rob,
  prueba_bp = bp,
  reset_test = reset_lm,
  vifs = vifs,
  kappa = kappa_X
)
## $resumen_ols
## # A tibble: 1 × 12
##   r.squared adj.r.squared sigma statistic p.value    df logLik   AIC   BIC
##       <dbl>         <dbl> <dbl>     <dbl>   <dbl> <dbl>  <dbl> <dbl> <dbl>
## 1    0.0128       0.00784  13.7      2.59  0.0354     4 -3253. 6518. 6546.
## # ℹ 3 more variables: deviance <dbl>, df.residual <int>, nobs <int>
## 
## $coef_ci_ols
## # A tibble: 5 × 7
##   term        estimate std.error statistic p.value conf.low conf.high
##   <chr>          <dbl>     <dbl>     <dbl>   <dbl>    <dbl>     <dbl>
## 1 (Intercept)  10.0      24.3        0.412  0.680  -37.6      57.6   
## 2 lincome       1.78      0.713      2.49   0.0129   0.377     3.18  
## 3 educ         -0.389     0.169     -2.31   0.0213  -0.720    -0.0580
## 4 lcigpric     -2.91      5.83      -0.500  0.617  -14.3       8.52  
## 5 age          -0.0427    0.0288    -1.48   0.138   -0.0991    0.0138
## 
## $coef_se_rob
## # A tibble: 5 × 5
##   term        estimate std.error statistic p.value
##   <chr>          <dbl>     <dbl>     <dbl>   <dbl>
## 1 (Intercept)  10.0      25.7        0.389 0.697  
## 2 lincome       1.78      0.592      3.00  0.00277
## 3 educ         -0.389     0.163     -2.38  0.0176 
## 4 lcigpric     -2.91      6.08      -0.479 0.632  
## 5 age          -0.0427    0.0238    -1.79  0.0731 
## 
## $prueba_bp
## 
##  studentized Breusch-Pagan test
## 
## data:  m_ols
## BP = 6.4463, df = 4, p-value = 0.1682
## 
## 
## $reset_test
## 
##  RESET test
## 
## data:  m_ols
## RESET = 0.50538, df1 = 2, df2 = 800, p-value = 0.6035
## 
## 
## $vifs
##  lincome     educ lcigpric      age 
## 1.114855 1.145114 1.006860 1.034822 
## 
## $kappa
## [1] 2758.374

Conclusión general del modelo

El modelo inicial OLS cumplió con la mayoría de los supuestos clásicos de MCO: no presentó problemas de multicolinealidad, la varianza de los errores fue aproximadamente constante y la especificación funcional general resultó adecuada según la prueba RESET. Sin embargo, los gráficos de residuos parciales revelaron cierta curvatura en las variables edad y educación, lo que sugería una relación no estrictamente lineal con el consumo de cigarrillos.

Por esta razón, se propuso un modelo ampliado con términos cuadráticos para capturar esas posibles no linealidades. Tras la transformación, el modelo mostró mejor ajuste visual en los gráficos de componentes + residuos y una representación más coherente del comportamiento económico: el consumo aumenta con la edad y educación hasta cierto punto, y luego disminuye, lo cual refleja un patrón realista del hábito de fumar en función del ciclo de vida y el nivel educativo.

Aunque el test de Breusch–Pagan en el modelo final detectó cierta heterocedasticidad, la corrección mediante errores estándar robustos (HC1) permitió conservar la validez estadística de las inferencias. Las variables educación y edad (junto con sus términos cuadráticos) se consolidaron como significativas, confirmando su efecto determinante sobre el consumo.

En conclusión, el modelo transformado representa una versión más fiel y econométricamente sólida del fenómeno analizado: mantiene los supuestos esenciales de MCO, mejora la forma funcional y ofrece interpretaciones económicas consistentes con la realidad del comportamiento de los consumidores de cigarrillos.