El hábito de lectura en estudiantes universitarios del Perú

Índice

1. Importar base de datos

2. Unión de base de datos

3. Recodificación de la variable dependiente

4. Eliminación de valores NA

5. Recodificación de variables independientes y creación de variables dummy

6. Renombre de variables

7. Descripción de variables mediante gráficos

8. Regresión logística binaria

1. Importamos las bases de datos

library(rio)
cap300400 = import("cap300-400.csv")
cap100600 = import("cap100-600.csv")

2. Unimos las datas según las variables en común mediante el método left_join

library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
# Unir datasets y seleccionar variables
data = left_join(cap300400, cap100600, 
          by = c("UBIGEO", "VIVIENDA", "HOGAR", "ESTRATOSOCIO")) %>%
  select(P401, P310_N, P601, P209, ESTRATOSOCIO, P435)
## Warning in left_join(cap300400, cap100600, by = c("UBIGEO", "VIVIENDA", : Detected an unexpected many-to-many relationship between `x` and `y`.
## ℹ Row 118 of `x` matches multiple rows in `y`.
## ℹ Row 1 of `y` matches multiple rows in `x`.
## ℹ If a many-to-many relationship is expected, set `relationship =
##   "many-to-many"` to silence this warning.

3. Recodificamos la variable dependiente

library(dplyr)

data <- data %>%
  mutate(lectorfrecuente = if_else(P401 %in% c(1, 2) & P310_N == 5, 1, 0))

4. Eliminamos los valores NA

data = na.omit(data)
library(dplyr)

data %>%
  summarise(across(everything(), ~sum(is.na(.)), .names = "NA_{.col}"))
##   NA_P401 NA_P310_N NA_P601 NA_P209 NA_ESTRATOSOCIO NA_P435 NA_lectorfrecuente
## 1       0         0       0       0               0       0                  0

5. Recodificamos las variables, cambiamos de nombre y creamos variables dummy

library(dplyr)

data <- data %>%
  # Eliminar columnas
  select(-P401, -P310_N) %>%
  
  # Transformaciones de variables
  mutate(
    P601 = as.numeric(P601),
    
    P209 = factor(P209, 
                  levels = c(1, 2), 
                  labels = c("Hombre", "Mujer")),
    
    ESTRATOSOCIO = factor(ESTRATOSOCIO,
                          levels = c(1, 2, 3, 4, 5, 6),
                          labels = c("Estrato Alto", "Estrato Medio Alto", "Estrato Medio",
                                     "Estrato Medio Bajo", "Estrato Bajo", "Rural"),
                          ordered = TRUE),
    
    P435 = factor(P435,
                  levels = c(1, 2),
                  labels = c("Sí", "No")),
    
    lectorfrecuente = factor(lectorfrecuente,
                             levels = c(1, 0),
                             labels = c("Sí", "No")),

    # Crear variables dummy (referencia: Estrato Alto)
    estrato_medio_alto = as.numeric(ESTRATOSOCIO == "Estrato Medio Alto"),
    estrato_medio      = as.numeric(ESTRATOSOCIO == "Estrato Medio"),
    estrato_medio_bajo = as.numeric(ESTRATOSOCIO == "Estrato Medio Bajo"),
    estrato_bajo       = as.numeric(ESTRATOSOCIO == "Estrato Bajo"),
    estrato_rural      = as.numeric(ESTRATOSOCIO == "Rural")
  )

6. Volvemos a renombrar algunas variables que quedaron sueltas

library(dplyr)

data <- data %>%
  rename(
    sexo = P209,
    librosxhogar = P601,
    asistenciafestival = P435
  )

7. Descripción de las variables mediante gráficos

7.1. Gráfico de libros en el hogar según si es lector frecuente

library(ggplot2)
library(ggthemes)

ggplot(data, aes(x = lectorfrecuente, y = librosxhogar)) +
  geom_boxplot(fill = "#6794a7") +
  labs(
    title = "Libros en el hogar según si es lector frecuente",
    x = "Lector frecuente",
    y = "Cantidad de libros en el hogar",
    caption = "Fuente: Encuesta Nacional de Lectura 2022"
  ) +
  theme_economist() +
  scale_y_continuous(expand = expansion(mult = c(0, 0.05)))

7.2. Gráfico de lectores frecuentes según sexo

library(ggplot2)
library(ggthemes)
library(scales)
library(dplyr)

# Paso 1: Calcular proporciones por grupo
data_plot <- data %>%
  group_by(sexo, lectorfrecuente) %>%
  summarise(n = n(), .groups = "drop") %>%
  group_by(sexo) %>%
  mutate(prop = n / sum(n),
         label = percent(prop, accuracy = 1))

# Paso 2: Graficar
ggplot(data_plot, aes(x = sexo, y = prop, fill = lectorfrecuente)) +
  geom_col() +
  geom_text(aes(label = label),
            position = position_stack(vjust = 0.5), color = "white", size = 4.2) +
  scale_fill_manual(values = c("#6794a7", "#014d64")) +
  scale_y_continuous(labels = percent_format()) +
  labs(
    title = "Proporción de lectores frecuentes según sexo",
    x = "Sexo",
    y = "Porcentaje",
    fill = "Lector frecuente",
    caption = "Fuente: Encuesta Nacional de Lectura 2022"
  ) +
  theme_economist() +
  theme(
    plot.caption = element_text(size = 9, face = "italic", hjust = 1)
  )

7.3. Recodificamos la variable estratosocio para que sea más amigable las etiquetas en el gráfico

levels(data$ESTRATOSOCIO) <- c(
  "Alto", "Medio Alto", "Medio", "Medio Bajo", "Bajo", "Rural"
)

7.4. Gráfico de lectores frecuentes según estrato socioeconómico

library(ggplot2)
library(ggthemes)
library(scales)
library(dplyr)

# Recalcular proporciones con etiquetas cortas
data_plot <- data %>%
  group_by(ESTRATOSOCIO, lectorfrecuente) %>%
  summarise(n = n(), .groups = "drop") %>%
  group_by(ESTRATOSOCIO) %>%
  mutate(prop = n / sum(n),
         label = percent(prop, accuracy = 1))

# Gráfico
ggplot(data_plot, aes(x = ESTRATOSOCIO, y = prop, fill = lectorfrecuente)) +
  geom_col() +
  geom_text(aes(label = label),
            position = position_stack(vjust = 0.5), color = "white", size = 4.2) +
  scale_fill_manual(values = c("#6794a7", "#014d64")) +
  scale_y_continuous(labels = percent_format()) +
  labs(
    title = "Proporción de lectores frecuentes según estrato socioeconómico",
    x = "Estrato socioeconómico",
    y = "Porcentaje",
    fill = "Lector frecuente",
    caption = "Fuente: Encuesta Nacional de Lectura 2022"
  ) +
  theme_economist() +
  theme(
    plot.caption = element_text(size = 9, face = "italic", hjust = 1),
    axis.text.x = element_text(angle = 30, hjust = 1)
  )

7.5. Gráfico de lectores frecuentes según la asistencia a festivales del libro

library(ggplot2)
library(ggthemes)
library(scales)
library(dplyr)

# Calcular proporciones por grupo
data_plot <- data %>%
  group_by(asistenciafestival, lectorfrecuente) %>%
  summarise(n = n(), .groups = "drop") %>%
  group_by(asistenciafestival) %>%
  mutate(prop = n / sum(n),
         label = percent(prop, accuracy = 1))

# Gráfico
ggplot(data_plot, aes(x = asistenciafestival, y = prop, fill = lectorfrecuente)) +
  geom_col() +
  geom_text(aes(label = label),
            position = position_stack(vjust = 0.5), color = "white", size = 4.2) +
  scale_fill_manual(values = c("#6794a7", "#014d64")) +
  scale_y_continuous(labels = percent_format()) +
  labs(
    title = "Proporción de lectores frecuentes según asistencia a festivales",
    x = "Asistió a festival literario",
    y = "Porcentaje",
    fill = "Lector frecuente",
    caption = "Fuente: Encuesta Nacional de Lectura 2022"
  ) +
  theme_economist() +
  theme(
    plot.caption = element_text(size = 9, face = "italic", hjust = 1),
    axis.text.x = element_text(angle = 0)
  )

8. Análisis con regresión logística binaria

8.1. Primera hipótesis

### semilla
set.seed(2019)

### first hypothesis
h1=formula(lectorfrecuente~sexo)

#regression
rlog1=glm(h1, data=data,family = binomial)
modelrl=list('Ser lector frecuente (I)'=rlog1)

#f <- function(x) format(x, digits = 4, scientific = FALSE)
library(modelsummary)
## `modelsummary` 2.0.0 now uses `tinytable` as its default table-drawing
##   backend. Learn more at: https://vincentarelbundock.github.io/tinytable/
## 
## Revert to `kableExtra` for one session:
## 
##   options(modelsummary_factory_default = 'kableExtra')
##   options(modelsummary_factory_latex = 'kableExtra')
##   options(modelsummary_factory_html = 'kableExtra')
## 
## Silence this message forever:
## 
##   config_modelsummary(startup_message = FALSE)
modelsummary(modelrl,
             title = "Regresión Logística",
             stars = TRUE,
             output = "kableExtra")
Regresión Logística
&nbsp;Ser lector frecuente (I)
(Intercept) 0.143***
(0.032)
sexoMujer 0.105*
(0.044)
Num.Obs. 8506
AIC 11706.5
BIC 11720.6
Log.Lik. -5851.259
F 5.815
RMSE 0.50
+ p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001
levels(data$sexo)
## [1] "Hombre" "Mujer"
data$sexo <- relevel(data$sexo, ref = "Mujer")

# Definir fórmula
h1 <- formula(lectorfrecuente ~ sexo)

# Reajustar modelo con nueva referencia
rlog1 <- glm(h1, data = data, family = binomial)

# Tabla con modelsummary
library(modelsummary)
modelrl <- list('Ser lector frecuente (I)' = rlog1)

modelsummary(modelrl,
             title = "Regresión Logística con nueva referencia en sexo",
             stars = TRUE,
             output = "kableExtra")
Regresión Logística con nueva referencia en sexo
&nbsp;Ser lector frecuente (I)
(Intercept) 0.248***
(0.030)
sexoHombre -0.105*
(0.044)
Num.Obs. 8506
AIC 11706.5
BIC 11720.6
Log.Lik. -5851.259
F 5.815
RMSE 0.50
+ p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001

8.2. Segunda hipótesis y comparación de modelos

### second hypothesis
h2=formula(lectorfrecuente~sexo + librosxhogar)
rlog2=glm(h2, data=data,family = binomial)

modelsrl=list('Ser lector frecuente (I)'=rlog1,
             'Ser lector frecuente (II)'=rlog2)

# formato creado para modelsummary
formatoNumero = function(x) format(x, digits = 4, scientific = FALSE)
modelsummary(modelsrl,
             fmt=formatoNumero, # usa función que creé antes
             exponentiate = T, # coeficientes sin logaritmo
             statistic = 'conf.int', # mostrar ICs
             title = "Regresión Logísticas (Coeficientes Exponenciados)",
             stars = TRUE,
             output = "kableExtra")
Regresión Logísticas (Coeficientes Exponenciados)
&nbsp;Ser lector frecuente (I) &nbsp;Ser lector frecuente (II)
(Intercept) 1.2817*** 1.3072***
[1.2086, 1.3595] [1.2297, 1.3899]
sexoHombre 0.9001* 0.9000*
[0.8262, 0.9805] [0.8262, 0.9805]
librosxhogar 0.9995*
[0.9990, 0.9999]
Num.Obs. 8506 8506
AIC 11706.5 11703.0
BIC 11720.6 11724.1
Log.Lik. -5851.259 -5848.498
F 5.815 5.448
RMSE 0.50 0.50
+ p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001

8.3. Tercera hipótesis y comparación de modelos

# Fórmula usando variables dummy
h3 <- formula(lectorfrecuente ~ sexo + librosxhogar +
              estrato_medio_alto + estrato_medio +
              estrato_medio_bajo + estrato_bajo + estrato_rural)

# Modelo logístico
rlog3 <- glm(h3, data = data, family = binomial)

# Resumen de resultados
summary(rlog3)
## 
## Call:
## glm(formula = h3, family = binomial, data = data)
## 
## Deviance Residuals: 
##    Min      1Q  Median      3Q     Max  
## -1.814  -1.263   1.046   1.093   1.661  
## 
## Coefficients:
##                      Estimate Std. Error z value Pr(>|z|)    
## (Intercept)         0.1869408  0.1194451   1.565   0.1176    
## sexoHombre         -0.1040511  0.0440731  -2.361   0.0182 *  
## librosxhogar       -0.0003648  0.0002201  -1.658   0.0974 .  
## estrato_medio_alto -0.1319821  0.1265570  -1.043   0.2970    
## estrato_medio       0.0613387  0.1231730   0.498   0.6185    
## estrato_medio_bajo  0.1371611  0.1246418   1.100   0.2711    
## estrato_bajo        0.1199025  0.1219569   0.983   0.3255    
## estrato_rural       1.2443133  0.2652748   4.691 2.72e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 11708  on 8505  degrees of freedom
## Residual deviance: 11650  on 8498  degrees of freedom
## AIC: 11666
## 
## Number of Fisher Scoring iterations: 4
modelsrl=list('Ser lector frecuente (I)'=rlog1,
             'Ser lector frecuente (II)'=rlog2,
             'Ser lector frecuente (III)'=rlog3)

f <- function(x) format(x, digits = 4, scientific = FALSE)
modelsummary(modelsrl,
             fmt=formatoNumero,
             exponentiate = T, 
             statistic = 'conf.int',
             title = "Comparando Regresión Logísticas (Coeficientes Exponenciados)",
             stars = TRUE,
             gof_map = c("nobs","aic","bic","rmse","logLik"), #comparar
             gof_omit = c("F"),
             output = "kableExtra")
Comparando Regresión Logísticas (Coeficientes Exponenciados)
&nbsp;Ser lector frecuente (I) &nbsp;Ser lector frecuente (II) &nbsp;Ser lector frecuente (III)
(Intercept) 1.2817*** 1.3072*** 1.2056
[1.2086, 1.3595] [1.2297, 1.3899] [0.9541, 1.5246]
sexoHombre 0.9001* 0.9000* 0.9012*
[0.8262, 0.9805] [0.8262, 0.9805] [0.8266, 0.9825]
librosxhogar 0.9995* 0.9996+
[0.9990, 0.9999] [0.9992, 1.0001]
estrato_medio_alto 0.8764
[0.6835, 1.1229]
estrato_medio 1.0633
[0.8347, 1.3534]
estrato_medio_bajo 1.1470
[0.8979, 1.4642]
estrato_bajo 1.1274
[0.8872, 1.4316]
estrato_rural 3.4706***
[2.0959, 5.9534]
Num.Obs. 8506 8506 8506
AIC 11706.5 11703.0 11665.5
BIC 11720.6 11724.1 11721.9
RMSE 0.50 0.50 0.50
Log.Lik. -5851.259 -5848.498 -5824.767
+ p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001

8.4. Cuarta hipótesis y comparación de modelos

# Estableces la categoría de referencia para asistencia al festival
data$asistenciafestival <- relevel(data$asistenciafestival, ref = "No")

# Fórmula usando las dummies creadas
h4 <- formula(lectorfrecuente ~ sexo + librosxhogar +
              estrato_medio_alto + estrato_medio +
              estrato_medio_bajo + estrato_bajo + estrato_rural +
              asistenciafestival)

# Modelo logístico
rlog4 <- glm(h4, data = data, family = binomial)

# Resultados
summary(rlog4)
## 
## Call:
## glm(formula = h4, family = binomial, data = data)
## 
## Deviance Residuals: 
##    Min      1Q  Median      3Q     Max  
## -1.849  -1.267   1.021   1.075   1.597  
## 
## Coefficients:
##                        Estimate Std. Error z value Pr(>|z|)    
## (Intercept)           0.2867678  0.1211303   2.367   0.0179 *  
## sexoHombre           -0.1075414  0.0441562  -2.435   0.0149 *  
## librosxhogar         -0.0003525  0.0002190  -1.610   0.1074    
## estrato_medio_alto   -0.1434360  0.1268587  -1.131   0.2582    
## estrato_medio         0.0472743  0.1234805   0.383   0.7018    
## estrato_medio_bajo    0.1010765  0.1250902   0.808   0.4191    
## estrato_bajo          0.0783527  0.1224660   0.640   0.5223    
## estrato_rural         1.2229273  0.2656991   4.603 4.17e-06 ***
## asistenciafestivalSí -0.2719966  0.0499243  -5.448 5.09e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 11708  on 8505  degrees of freedom
## Residual deviance: 11620  on 8497  degrees of freedom
## AIC: 11638
## 
## Number of Fisher Scoring iterations: 4
modelsrl=list('Ser lector frecuente (I)'=rlog1,
             'Ser lector frecuente (II)'=rlog2,
             'Ser lector frecuente (III)'=rlog3,
             'Ser lector frecuente (IV)'=rlog4)

f <- function(x) format(x, digits = 4, scientific = FALSE)
modelsummary(modelsrl,
             fmt=formatoNumero,
             exponentiate = T, 
             statistic = 'conf.int',
             title = "Comparando Regresión Logísticas (Coeficientes Exponenciados)",
             stars = TRUE,
             gof_map = c("nobs","aic","bic","rmse","logLik"), #comparar
             gof_omit = c("F"),
             output = "kableExtra")
Comparando Regresión Logísticas (Coeficientes Exponenciados)
&nbsp;Ser lector frecuente (I) &nbsp;Ser lector frecuente (II) &nbsp;Ser lector frecuente (III) &nbsp;Ser lector frecuente (IV)
(Intercept) 1.2817*** 1.3072*** 1.2056 1.3321*
[1.2086, 1.3595] [1.2297, 1.3899] [0.9541, 1.5246] [1.0509, 1.6903]
sexoHombre 0.9001* 0.9000* 0.9012* 0.8980*
[0.8262, 0.9805] [0.8262, 0.9805] [0.8266, 0.9825] [0.8236, 0.9792]
librosxhogar 0.9995* 0.9996+ 0.9996
[0.9990, 0.9999] [0.9992, 1.0001] [0.9992, 1.0001]
estrato_medio_alto 0.8764 0.8664
[0.6835, 1.1229] [0.6753, 1.1107]
estrato_medio 1.0633 1.0484
[0.8347, 1.3534] [0.8226, 1.3352]
estrato_medio_bajo 1.1470 1.1064
[0.8979, 1.4642] [0.8653, 1.4135]
estrato_bajo 1.1274 1.0815
[0.8872, 1.4316] [0.8502, 1.3746]
estrato_rural 3.4706*** 3.3971***
[2.0959, 5.9534] [2.0498, 5.8319]
asistenciafestivalSí 0.7619***
[0.6908, 0.8402]
Num.Obs. 8506 8506 8506 8506
AIC 11706.5 11703.0 11665.5 11637.9
BIC 11720.6 11724.1 11721.9 11701.3
RMSE 0.50 0.50 0.50 0.50
Log.Lik. -5851.259 -5848.498 -5824.767 -5809.931
+ p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001

8.5. Elección del mejor modelo

#Según el AIC y el valor Log.lik., el modelo 4 es el mejor modelo porque poseen los menores valores de AIC y Log.Lik. en comparación con los 3 modelos restantes.

8.6. Interpretación de los coeficientes:

#SEXOHombre = 0.8980 Este coeficiente indica que, manteniendo constante la edad, el estrato socioeconomico y la asistencia a festivales, ser hombre tienen en promedio una probabilidad un 10.2% menor de ser lectores frecuentes que las mujeres

#librosxhogar = esta variable no es significativa para la regresión, pues su p-valor es mayor a 0.05.

#estrato_medio_alto = esta variable no es significativa para la regresión, pues su p-valor es mayor a 0.05.

#estrato_medio = esta variable no es significativa para la regresión, pues su p-valor es mayor a 0.05.

#estrato_medio_bajo = esta variable no es significativa para la regresión, pues su p-valor es mayor a 0.05.

#estrato_bajo = esta variable no es significativa para la regresión, pues su p-valor es mayor a 0.05.

#estrato_rural = esta variable es altamente significativa. Su coeficiente exponenciado es 3.3971, lo que significa que las personas de estrato rural tienen en promedio 3.40 veces más probabilidades de ser lectoras frecuentes que los de estrato alto si se mantiene constantes las demás variables.

100 * ((3.3971) - 1)
## [1] 239.71

#Esto significa que los odds de ser lector frecuente en el estrato rural son 239.7% más altos que en el estrato alto.

#asistenciafestivalSí = 0.7619 / Las personas que asistieron a un festival literario tienen un 23.8% menos de odds de ser lectoras frecuentes que quienes no asistieron, mientras todo lo demás es constante.