07 DE JUNIO 2025

Quarto

install.packages(c(“MASS”, “msm”, “car”, “margins”, “modelsummary”, “Rchoice”)) install.packages(“modelsummary”) install.packages(“VGAM”) install.packages(“pROC”) install.packages(“nnet”) install.packages(“pandoc”) install.packages(“officer”) install.packages(“flextable”) install.packages(“broom”)

library(MASS)
library(VGAM)

Cargando paquete requerido: stats4

Cargando paquete requerido: splines

library(pROC)

Type 'citation("pROC")' for a citation.


Adjuntando el paquete: 'pROC'

The following objects are masked from 'package:stats':

    cov, smooth, var

library(ggplot2)
library(nnet)
library(officer)
library(flextable)
library(broom)
library(foreign)
library(msm)
library(car)

Cargando paquete requerido: carData


Adjuntando el paquete: 'car'

The following object is masked from 'package:VGAM':

    logit

library(margins)
library(modelsummary)


Adjuntando el paquete: 'modelsummary'

The following object is masked from 'package:VGAM':

    Max

library(haven)
library(Rchoice)

Cargando paquete requerido: Formula

Cargando paquete requerido: maxLik

Cargando paquete requerido: miscTools


Please cite the 'maxLik' package as:
Henningsen, Arne and Toomet, Ott (2011). maxLik: A package for maximum likelihood estimation in R. Computational Statistics 26(3), 443-458. DOI 10.1007/s00180-010-0217-1.

If you have questions, suggestions, or comments regarding the 'maxLik' package, please use a forum or 'tracker' at maxLik's R-Forge site:
https://r-forge.r-project.org/projects/maxlik/

setwd("C:/Users/alexd/Desktop/TAREA07_06_25/07 DE JUNIO 2025")

##Leer los datos (ajustar la ruta del archivo)

library(haven)
data <- read_dta("Data1_R.dta")
View(data)

##Ver las primeras filas

head(data)

# A tibble: 6 × 50
  area       empleo          region   edad t_hijos nac_vivo_murieron mortinato_2
  <dbl+lbl>  <dbl+lbl>       <dbl+l> <dbl>   <dbl> <dbl+lbl>         <dbl+lbl>  
1 1 [Urbano] 1 [Trabajó al … 1 [Sie…    19       1 0 [No]            0 [No]     
2 1 [Urbano] 0 [No trabajó]  1 [Sie…    23       1 0 [No]            0 [No]     
3 1 [Urbano] 1 [Trabajó al … 1 [Sie…    38       5 0 [No]            0 [No]     
4 1 [Urbano] 0 [No trabajó]  1 [Sie…    18       1 0 [No]            0 [No]     
5 1 [Urbano] 0 [No trabajó]  1 [Sie…    21       1 0 [No]            0 [No]     
6 1 [Urbano] 1 [Trabajó al … 1 [Sie…    22       1 0 [No]            0 [No]     
# ℹ 43 more variables: depresion_pp <dbl+lbl>, intensidad_dpp <dbl+lbl>,
#   etnia <dbl+lbl>, f2_s2_216_1 <dbl+lbl>, f2_s2_216_2 <dbl>,
#   f2_s2_218_1_a <dbl+lbl>, tiempo_dpp <dbl+lbl>, f2_s5_504a_1 <dbl+lbl>,
#   f2_s5_504b_1 <dbl+lbl>, f2_s5_504c_1 <dbl+lbl>, f2_s5_504d_1 <dbl+lbl>,
#   f2_s5_504e_1 <dbl+lbl>, f2_s5_504f_1 <dbl+lbl>, f2_s5_504g_1 <dbl+lbl>,
#   f2_s5_504h_1 <dbl+lbl>, f2_s5_504i_1 <dbl+lbl>, f2_s5_504j_1 <dbl+lbl>,
#   f2_s5_504k_1 <dbl+lbl>, est_civil <dbl+lbl>, q_usted <dbl+lbl>, …

Revisar estructura de los datos

str(data)

tibble [16,451 × 50] (S3: tbl_df/tbl/data.frame)
 $ area             : dbl+lbl [1:16451] 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1...
   ..@ label       : chr "Área"
   ..@ format.stata: chr "%8.0g"
   ..@ labels      : Named num [1:2] 0 1
   .. ..- attr(*, "names")= chr [1:2] "Rural" "Urbano"
 $ empleo           : dbl+lbl [1:16451] 1, 0, 1, 0, 0, 1, 1, 1, 0, 0, 0, 1, 0, 0, 0, 1, 1, 1...
   ..@ label       : chr "1.  ¿Qué hizo (…) la semana pasada:"
   ..@ format.stata: chr "%26.0g"
   ..@ labels      : Named num [1:2] 0 1
   .. ..- attr(*, "names")= chr [1:2] "No trabajó" "Trabajó al menos una hora"
 $ region           : dbl+lbl [1:16451] 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1...
   ..@ label       : chr "Región"
   ..@ format.stata: chr "%9.0g"
   ..@ labels      : Named num [1:4] 1 2 3 4
   .. ..- attr(*, "names")= chr [1:4] "Sierra" "Costa" "Amazonía" "Insular"
 $ edad             : num [1:16451] 19 23 38 18 21 22 40 31 19 34 ...
  ..- attr(*, "label")= chr "101.  Entonces, ¿Cuántos años cumplidos tiene?"
  ..- attr(*, "format.stata")= chr "%8.0g"
 $ t_hijos          : num [1:16451] 1 1 5 1 1 1 3 2 1 2 ...
  ..- attr(*, "label")= chr "208.c  Total hijos/as en casa"
  ..- attr(*, "format.stata")= chr "%8.0g"
 $ nac_vivo_murieron: dbl+lbl [1:16451] 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0...
   ..@ label       : chr " 210.  ¿Tuvo usted hijos o hijas que nacieron vivos/as y que murieron"
   ..@ format.stata: chr "%8.0g"
   ..@ labels      : Named num [1:2] 0 1
   .. ..- attr(*, "names")= chr [1:2] "No" "Si"
 $ mortinato_2      : dbl+lbl [1:16451] 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0...
   ..@ label       : chr "212.1  De todos los embarazos que usted ha tenido en su vida"
   ..@ format.stata: chr "%8.0g"
   ..@ labels      : Named num [1:2] 0 1
   .. ..- attr(*, "names")= chr [1:2] "No" "Si"
 $ depresion_pp     : dbl+lbl [1:16451] 1, 1, 0, 0, 0, 0, 1, 0, 0, 1, 0, 1, 1, 0, 0, 0, 0, 0...
   ..@ label       : chr "219.  ¿Después de alguno de sus partos, sintió tristeza?"
   ..@ format.stata: chr "%8.0g"
   ..@ labels      : Named num [1:2] 0 1
   .. ..- attr(*, "names")= chr [1:2] "No" "Si"
 $ intensidad_dpp   : dbl+lbl [1:16451] 2, 3, 1, 1, 1, 1, 3, 1, 1, 3, 1, 3, 2, 1, 1, 1, 1, 1...
   ..@ label       : chr "Intensidad DPP"
   ..@ format.stata: chr "%12.0g"
   ..@ labels      : Named num [1:3] 0 1 2
   .. ..- attr(*, "names")= chr [1:3] "No tiene dpp" "Poca dpp" "Mucha dpp"
 $ etnia            : dbl+lbl [1:16451] 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0...
   ..@ label       : chr "9.1  ¿Cómo se Identifica (…) según su cultura y costumbres:"
   ..@ format.stata: chr "%8.0g"
   ..@ labels      : Named num [1:2] 0 1
   .. ..- attr(*, "names")= chr [1:2] "Otra" "Indigena"
 $ f2_s2_216_1      : dbl+lbl [1:16451] 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 1, 2, 2, 2, 1, 2, 2...
   ..@ label       : chr "216.1  Hay mujeres que pierden sus embarazos antes de cumplir el quinto mes"
   ..@ format.stata: chr "%8.0g"
   ..@ labels      : Named num [1:2] 1 2
   .. ..- attr(*, "names")= chr [1:2] "si" "no"
 $ f2_s2_216_2      : num [1:16451] NA NA NA NA NA NA NA NA NA NA ...
  ..- attr(*, "label")= chr "216.2  ¿Cuántos?"
  ..- attr(*, "format.stata")= chr "%8.0g"
 $ f2_s2_218_1_a    : dbl+lbl [1:16451] 2, 2, 2, 2, 1, 2, 1, 2, 2, 1, 2, 1, 1, 2, 2, 1, 1, 2...
   ..@ label       : chr "218.a  ¿Hombre o mujer?"
   ..@ format.stata: chr "%8.0g"
   ..@ labels      : Named num [1:2] 1 2
   .. ..- attr(*, "names")= chr [1:2] "hombre" "mujer"
 $ tiempo_dpp       : dbl+lbl [1:16451] 1, 1, 0, 0, 0, 0, 3, 0, 0, 1, 0, 1, 1, 0, 0, 0, 0, 0...
   ..@ label       : chr "221.  ¿Durante cuánto tiempo se sintió triste?"
   ..@ format.stata: chr "%8.0g"
   ..@ labels      : Named num [1:5] 1 2 3 4 88
   .. ..- attr(*, "names")= chr [1:5] "menos de 3 meses" "3 meses" "6 meses" "1año o más" ...
 $ f2_s5_504a_1     : dbl+lbl [1:16451] 2, 3, 3, 3, 3, 3, 3, 3, 2, 3, 3, 3, 3, 3, 3, 3, 3, 3...
   ..@ label       : chr "504.a  El trato que recibió del personal de salud que le atendió?PN"
   ..@ format.stata: chr "%8.0g"
   ..@ labels      : Named num [1:3] 1 2 3
   .. ..- attr(*, "names")= chr [1:3] "Malo" "Regular" "Bueno"
 $ f2_s5_504b_1     : dbl+lbl [1:16451] 3, 3, 3, 3, 3, 3, 3, 3, 2, 3, 3, 3, 3, 3, 3, 3, 3, 3...
   ..@ label       : chr "504.b  La experiencia/conocimiento de la persona que le atendió?PN"
   ..@ format.stata: chr "%8.0g"
   ..@ labels      : Named num [1:3] 1 2 3
   .. ..- attr(*, "names")= chr [1:3] "Malo" "Regular" "Bueno"
 $ f2_s5_504c_1     : dbl+lbl [1:16451] 3, 3, 2, 3, 3, 3, 3, 3, 2, 3, 3, 3, 3, 3, 3, 3, 3, 3...
   ..@ label       : chr "504.c  La disponibilidad del médico, enfermera u otro personal"
   ..@ format.stata: chr "%8.0g"
   ..@ labels      : Named num [1:3] 1 2 3
   .. ..- attr(*, "names")= chr [1:3] "Malo" "Regular" "Bueno"
 $ f2_s5_504d_1     : dbl+lbl [1:16451] 3, 3, 1, 3, 3, 3, 3, 3, 2, 3, 3, 3, 3, 3, 3, 3, 3, 3...
   ..@ label       : chr "504.d  La privacidad durante su atención? - PNoC"
   ..@ format.stata: chr "%8.0g"
   ..@ labels      : Named num [1:3] 1 2 3
   .. ..- attr(*, "names")= chr [1:3] "Malo" "Regular" "Bueno"
 $ f2_s5_504e_1     : dbl+lbl [1:16451] 2, 3, 2, 3, 3, 3, 3, 3, 2, 3, 3, 2, 3, 3, 3, 3, 3, 3...
   ..@ label       : chr "504.e  La comodidad del establecimiento? - Parto normal o cesárea"
   ..@ format.stata: chr "%8.0g"
   ..@ labels      : Named num [1:3] 1 2 3
   .. ..- attr(*, "names")= chr [1:3] "Malo" "Regular" "Bueno"
 $ f2_s5_504f_1     : dbl+lbl [1:16451] 3, 3, 3, 3, 3, 3, 3, 3, 2, 3, 3, 2, 3, 3, 3, 3, 3, 3...
   ..@ label       : chr "504.f  El tiempo que tardaron en atenderla? - Parto normal o cesárea"
   ..@ format.stata: chr "%8.0g"
   ..@ labels      : Named num [1:3] 1 2 3
   .. ..- attr(*, "names")= chr [1:3] "Malo" "Regular" "Bueno"
 $ f2_s5_504g_1     : dbl+lbl [1:16451] 1, 3, 3, 3, 3, 3, 3, 3, 2, 3, 3, 3, 3, 3, 3, 3, 3, 3...
   ..@ label       : chr "504.g  El horario de atención? - Parto normal o cesárea"
   ..@ format.stata: chr "%8.0g"
   ..@ labels      : Named num [1:3] 1 2 3
   .. ..- attr(*, "names")= chr [1:3] "Malo" "Regular" "Bueno"
 $ f2_s5_504h_1     : dbl+lbl [1:16451] 2, 3, 3, 3, 3, 3, 3, 3, 2, 3, 3, 3, 3, 3, 3, 3, 3, 3...
   ..@ label       : chr "504.h  La limpieza/aseo del establecimiento? - Parto normal o cesárea"
   ..@ format.stata: chr "%8.0g"
   ..@ labels      : Named num [1:3] 1 2 3
   .. ..- attr(*, "names")= chr [1:3] "Malo" "Regular" "Bueno"
 $ f2_s5_504i_1     : dbl+lbl [1:16451] 2, 3, 1, 3, 3, 3, 3, 3, 2, 3, 3, 3, 3, 3, 3, 3, 3, 3...
   ..@ label       : chr "504.i  La ropa/vestimenta utilizada en el parto? - Parto normal o cesárea"
   ..@ format.stata: chr "%8.0g"
   ..@ labels      : Named num [1:3] 1 2 3
   .. ..- attr(*, "names")= chr [1:3] "Malo" "Regular" "Bueno"
 $ f2_s5_504j_1     : dbl+lbl [1:16451] 2, 3, 3, 3, 3, 3, 3, 3, 2, 3, 3, 3, 3, 3, 3, 3, 3, 3...
   ..@ label       : chr "504.j  La información que le dio el personal de salud sobre su salud"
   ..@ format.stata: chr "%8.0g"
   ..@ labels      : Named num [1:3] 1 2 3
   .. ..- attr(*, "names")= chr [1:3] "Malo" "Regular" "Bueno"
 $ f2_s5_504k_1     : dbl+lbl [1:16451] 3, 3, 2, 3, 3, 3, 3, 3, 2, 3, 3, 3, 3, 3, 3, 3, 3, 3...
   ..@ label       : chr "504.k  El respeto a su cultura y costumbres relacionadas con el parto? - Parto n"
   ..@ format.stata: chr "%8.0g"
   ..@ labels      : Named num [1:3] 1 2 3
   .. ..- attr(*, "names")= chr [1:3] "Malo" "Regular" "Bueno"
 $ est_civil        : dbl+lbl [1:16451] 0, 0, 0, 0, 0, 1, 0, 1, 0, 1, 0, 0, 1, 1, 0, 0, 1, 1...
   ..@ label       : chr "900.  ¿Actualmente usted esta:"
   ..@ format.stata: chr "%12.0g"
   ..@ labels      : Named num [1:2] 0 1
   .. ..- attr(*, "names")= chr [1:2] "En otro caso" "Casada"
 $ q_usted          : dbl+lbl [1:16451] 0, 0, 0, 1, 0, 1, 1, 1, 0, 1, 0, 1, 1, 1, 1, 1, 0, 1...
   ..@ label       : chr "403.  En la época en la que quedó embarazada de (…), quería Usted:"
   ..@ format.stata: chr "%17.0g"
   ..@ labels      : Named num [1:2] 0 1
   .. ..- attr(*, "names")= chr [1:2] "No tener ese hijo" "Tener ese hijo"
 $ q_pareja         : dbl+lbl [1:16451] 0, 1, 0, 1, 1, 1, 1, 1, 0, 1, 0, 0, 1, 1, 1, 1, 0, 1...
   ..@ label       : chr "405.  ¿Quería su pareja:"
   ..@ format.stata: chr "%17.0g"
   ..@ labels      : Named num [1:2] 0 1
   .. ..- attr(*, "names")= chr [1:2] "No tener ese hijo" "Tener ese hijo"
 $ f2_s4b_406_      : dbl+lbl [1:16451] 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 1, 1, 1, 1, 1, 1, 1, 1...
   ..@ label       : chr "406.  ¿Tuvo algún control prenatal cuando estaba embarazada de (...)?"
   ..@ format.stata: chr "%8.0g"
   ..@ labels      : Named num [1:2] 1 2
   .. ..- attr(*, "names")= chr [1:2] "si" "no"
 $ lugar            : dbl+lbl [1:16451] 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 0, 1, 0, 1, 0, 0, 0, 0...
   ..@ label       : chr "423.a  ¿En que lugar tuvo el parto de (..).?"
   ..@ format.stata: chr "%32.0g"
   ..@ labels      : Named num [1:2] 0 1
   .. ..- attr(*, "names")= chr [1:2] "Establecimiento privado" "Establecimiento de salud del msp"
 $ form_parto       : dbl+lbl [1:16451] 0, 1, 1, 1, 1, 1, 1, 0, 0, 0, 1, 1, 0, 1, 0, 0, 0, 0...
   ..@ label       : chr "425.  El parto de (...) fue:"
   ..@ format.stata: chr "%15.0g"
   ..@ labels      : Named num [1:2] 0 1
   .. ..- attr(*, "names")= chr [1:2] "Cesárea" "Normal"
 $ c_posparto       : dbl+lbl [1:16451] 1, 0, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1...
   ..@ label       : chr "440.  ¿Tuvo usted algún control después del parto de (...)?"
   ..@ format.stata: chr "%8.0g"
   ..@ labels      : Named num [1:2] 0 1
   .. ..- attr(*, "names")= chr [1:2] "No" "Si"
 $ anios_esc        : num [1:16451] 14 10 6 8 16 9 10 13 14 6 ...
  ..- attr(*, "label")= chr "Años escolaridad"
  ..- attr(*, "format.stata")= chr "%9.0g"
 $ ingrl            : num [1:16451] 0 0 240 0 0 0 0 60 0 0 ...
  ..- attr(*, "label")= chr "Ingreso laboral"
  ..- attr(*, "format.stata")= chr "%9.0g"
 $ planf            : dbl+lbl [1:16451] 0, 0, 0, 1, 0, 1, 1, 1, 0, 1, 0, 0, 1, 1, 1, 1, 0, 1...
   ..@ label       : chr "Planificación familiar"
   ..@ format.stata: chr "%13.0g"
   ..@ labels      : Named num [1:2] 0 1
   .. ..- attr(*, "names")= chr [1:2] "No planifico" "Si planificó"
 $ lingrl           : num [1:16451] 0 0 5.48 0 0 ...
  ..- attr(*, "label")= chr "Logaritmo ingreso laboral"
  ..- attr(*, "format.stata")= chr "%9.0g"
 $ lugarx1          : num [1:16451] 0 0 0 0 0 0 0 1 0 0 ...
  ..- attr(*, "label")= chr "lugar==Establecimiento privado"
  ..- attr(*, "format.stata")= chr "%8.0g"
 $ lugarx2          : num [1:16451] 1 1 1 1 1 1 1 0 1 1 ...
  ..- attr(*, "label")= chr "lugar==Establecimiento de salud del msp"
  ..- attr(*, "format.stata")= chr "%8.0g"
 $ regionx1         : num [1:16451] 1 1 1 1 1 1 1 1 1 1 ...
  ..- attr(*, "label")= chr "region==Sierra"
  ..- attr(*, "format.stata")= chr "%8.0g"
 $ regionx2         : num [1:16451] 0 0 0 0 0 0 0 0 0 0 ...
  ..- attr(*, "label")= chr "region==Costa"
  ..- attr(*, "format.stata")= chr "%8.0g"
 $ regionx3         : num [1:16451] 0 0 0 0 0 0 0 0 0 0 ...
  ..- attr(*, "label")= chr "region==Amazonía"
  ..- attr(*, "format.stata")= chr "%8.0g"
 $ regionx4         : num [1:16451] 0 0 0 0 0 0 0 0 0 0 ...
  ..- attr(*, "label")= chr "region==Insular"
  ..- attr(*, "format.stata")= chr "%8.0g"
 $ intensidad_dppx1 : num [1:16451] 0 0 1 1 1 1 0 1 1 0 ...
  ..- attr(*, "label")= chr "intensidad_dpp==No tiene dpp"
  ..- attr(*, "format.stata")= chr "%8.0g"
 $ intensidad_dppx2 : num [1:16451] 1 0 0 0 0 0 0 0 0 0 ...
  ..- attr(*, "label")= chr "intensidad_dpp==Poca dpp"
  ..- attr(*, "format.stata")= chr "%8.0g"
 $ intensidad_dppx3 : num [1:16451] 0 1 0 0 0 0 1 0 0 1 ...
  ..- attr(*, "label")= chr "intensidad_dpp==Mucha dpp"
  ..- attr(*, "format.stata")= chr "%8.0g"
 $ tiempo_dppx1     : num [1:16451] 0 0 1 1 1 1 0 1 1 0 ...
  ..- attr(*, "label")= chr "tiempo_dpp==     0.0000"
  ..- attr(*, "format.stata")= chr "%8.0g"
 $ tiempo_dppx2     : num [1:16451] 1 1 0 0 0 0 0 0 0 1 ...
  ..- attr(*, "label")= chr "tiempo_dpp==menos de 3 meses"
  ..- attr(*, "format.stata")= chr "%8.0g"
 $ tiempo_dppx3     : num [1:16451] 0 0 0 0 0 0 0 0 0 0 ...
  ..- attr(*, "label")= chr "tiempo_dpp==3 meses"
  ..- attr(*, "format.stata")= chr "%8.0g"
 $ tiempo_dppx4     : num [1:16451] 0 0 0 0 0 0 1 0 0 0 ...
  ..- attr(*, "label")= chr "tiempo_dpp==6 meses"
  ..- attr(*, "format.stata")= chr "%8.0g"
 $ tiempo_dppx5     : num [1:16451] 0 0 0 0 0 0 0 0 0 0 ...
  ..- attr(*, "label")= chr "tiempo_dpp==1año o más"
  ..- attr(*, "format.stata")= chr "%8.0g"

str(data)

Ver distribución de la variable dependiente

table(data$intensidad_dpp) #(ordenada)


    1     2     3 
12828  1790  1833

table(data$depresion_pp)  #(dicótomica)


    0     1 
12828  3623

##Ejemplo 1: Modelo con variable dependiente odenada####

library("Rchoice")
logitord <- Rchoice(intensidad_dpp ~ lingrl + anios_esc + edad + t_hijos + etnia + area, 
                  data = data, 
                  na.action = na.omit, 
                  family = ordinal('logit'))

summary(logitord)


Model: ordinal
Model estimated on: sáb. jun. 07 23:47:39 2025 

Call:
Rchoice(formula = intensidad_dpp ~ lingrl + anios_esc + edad + 
    t_hijos + etnia + area, data = data, na.action = na.omit, 
    family = ordinal("logit"), method = "bfgs")


Frequencies of categories:
y
     1      2      3 
0.7798 0.1088 0.1114 
The estimation took: 0h:0m:2s 

Coefficients:
           Estimate Std. Error z-value Pr(>|z|)    
kappa.1    0.822194   0.018816  43.697  < 2e-16 ***
constant  -2.352638   0.100572 -23.393  < 2e-16 ***
lingrl     0.001703   0.007114   0.239  0.81077    
anios_esc -0.009924   0.004865  -2.040  0.04136 *  
edad       0.034483   0.003190  10.810  < 2e-16 ***
t_hijos    0.039105   0.018717   2.089  0.03669 *  
etnia      0.344926   0.059842   5.764 8.22e-09 ***
area       0.118337   0.042259   2.800  0.00511 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Optimization of log-likelihood by BFGS maximization
Log Likelihood: -11050
Number of observations: 16451
Number of iterations: 129
Exit of MLE: successful convergence

Análisis

El ingreso es una variable que no presenta significancia estadística, por lo tanto, no permite conocer su efecto sobre la intensidad de depresión. A medida que aumenta un año de escolaridad reduce la intensidad de depresión post parto, dado que, es estadísticamente significativo. Cuando las mujeres aumentan un año de edad, de igual forma aumenta la depresion post parto A medida que tienen un hijo más, es más probable que sufran depresión post parto Las mujeres de la etnia indígena tienen mayor probabilidad de padecer depresión post parto. Las mujeres que habitan en zonas rurales tienen mayor probabilidad de sufrir de depresión.

Exportar información a word

library(modelsummary)
modelsummary(logitord, output = "modelo_logit_ord.docx")

library(modelsummary)
modelsummary(logitord, 
             output = "modelo_logit_ord2.docx", 
             stars = c('*' = .05, '**' = .01, '***' = .001))

library(modelsummary)
modelsummary(logitord,
             output = "modelo_logit_ord3.docx",
             statistic = c("p.value" = "p = {p.value}{stars}"),
             stars = TRUE)

library(modelsummary)
modelsummary(logitord,
             output = "modelo_logit_ord4.docx",
             statistic = "{statistic}",
             stars = TRUE)

##Nota sobre RMSE

RMSE significa Root Mean Squared Error (Raíz del Error Cuadrático Medio). Es una métrica común para evaluar la precisión de modelos de predicción, especialmente en regresión. Más bajo = mejor. Un RMSE de 0 indica predicciones perfectas. Este es más usual en casos de variables cuantitativas, para determinar la precisión lo hacemos con la matriz de confusión.

Criterio de información de Akaike (AIC) y criterio de Información Bayesiano (BIC)

##AIC y BIC del Logit

AIC(logitord)

[1] 22123.39

BIC(logitord)

[1] 22185.05

##Ajustar modelo Probit ordenado

library("Rchoice")
probitord <- Rchoice(intensidad_dpp ~ lingrl + anios_esc + edad + t_hijos + etnia + area, 
                  data = data, 
                  na.action = na.omit, 
                  family = ordinal('probit'))

summary(probitord)


Model: ordinal
Model estimated on: sáb. jun. 07 23:47:43 2025 

Call:
Rchoice(formula = intensidad_dpp ~ lingrl + anios_esc + edad + 
    t_hijos + etnia + area, data = data, na.action = na.omit, 
    family = ordinal("probit"), method = "bfgs")


Frequencies of categories:
y
     1      2      3 
0.7798 0.1088 0.1114 
The estimation took: 0h:0m:1s 

Coefficients:
           Estimate Std. Error z-value Pr(>|z|)    
kappa.1    0.454873   0.010062  45.205  < 2e-16 ***
constant  -1.404867   0.057158 -24.579  < 2e-16 ***
lingrl     0.001106   0.004054   0.273   0.7850    
anios_esc -0.006361   0.002774  -2.293   0.0218 *  
edad       0.020295   0.001839  11.037  < 2e-16 ***
t_hijos    0.023193   0.010864   2.135   0.0328 *  
etnia      0.196005   0.034697   5.649 1.61e-08 ***
area       0.071133   0.023881   2.979   0.0029 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Optimization of log-likelihood by BFGS maximization
Log Likelihood: -11050
Number of observations: 16451
Number of iterations: 82
Exit of MLE: successful convergence

library(modelsummary)
modelsummary(probitord,
             output = "modelo_probitord_ord.docx",
             statistic = "{statistic}",
             stars = TRUE)

Comparación AIC y BIC

AIC(probitord)

[1] 22109.98

BIC(probitord)

[1] 22171.64

aic_logit_ord <- AIC(logitord)
aic_probit_ord <- AIC(probitord)

bic_logit_ord <- BIC(logitord)
bic_probit_ord <- BIC(probitord)

##Mostrar resultados

cat("AIC Logit Ordenado:", aic_logit_ord, " | AIC Probit Ordenado:", aic_probit_ord, "\n")

AIC Logit Ordenado: 22123.39  | AIC Probit Ordenado: 22109.98

cat("BIC Logit Ordenado:", bic_logit_ord, " | BIC Probit Ordenado:", bic_probit_ord, "\n")

BIC Logit Ordenado: 22185.05  | BIC Probit Ordenado: 22171.64

Conclusión

El modelo con menor AIC/BIC es el preferido. - Si Probit tiene menor AIC y BIC → es mejor. - Si Logit tiene menor AIC y BIC → es mejor.

###Efectos marginales###### ## Logit

X <- cbind(1, logitord$mf[, -1])

coeficientes<- logitord$coefficients

coeficientes_1<- coeficientes[2:8]

ai <- crossprod(t(X), coeficientes_1) #(producto cruzado por los coeficientes)

deriv1<-mean(dnorm(ai)) #(saco derivada)

##Ingreso

library(msm)
gamma_hatq1_1 <- coeficientes_1["lingrl"]
AME_hat_q1lc_1   <- round(mean(dnorm(ai)) * (gamma_hatq1_1),3);AME_hat_q1lc_1

lingrl 
     0

AME_hat_q1lc_se_1 <-  round(deltamethod(~ deriv1 * x1, coef(logitord), vcov(logitord), ses =  TRUE),3);AME_hat_q1lc_se_1

[1] 0.003

AME_hat_q1lc_cil_1 <-  round(c(AME_hat_q1lc_1 - 1.96 * AME_hat_q1lc_se_1),3);AME_hat_q1lc_cil_1

lingrl 
-0.006

AME_hat_q1lc_ciu_1 <-  round(c(AME_hat_q1lc_1 + 1.96 * AME_hat_q1lc_se_1),3);AME_hat_q1lc_ciu_1

lingrl 
 0.006

AME_hat_q1lc_pr_1 <-  round(2 * pnorm(-abs(AME_hat_q1lc_1/AME_hat_q1lc_se_1)),3);AME_hat_q1lc_pr_1

lingrl 
     1

AME_hat_q1lc_star_1<- " "  
lingrl_q1<-cbind(AME_hat_q1lc_1,AME_hat_q1lc_star_1, AME_hat_q1lc_se_1,AME_hat_q1lc_pr_1,AME_hat_q1lc_cil_1,AME_hat_q1lc_ciu_1)

##Escolaridad

library(msm)
gamma_hatq2_2 <- coeficientes_1["anios_esc"]
AME_hat_q2lc_2   <- round(mean(dnorm(ai)) * (gamma_hatq2_2),3);AME_hat_q2lc_2

anios_esc 
   -0.002

AME_hat_q2lc_se_2 <-  round(deltamethod(~ deriv1 * x2, coef(logitord), vcov(logitord), ses =  TRUE),3);AME_hat_q2lc_se_2

[1] 0.018

AME_hat_q2lc_cil_2 <-  round(c(AME_hat_q2lc_2 - 1.96 * AME_hat_q2lc_se_2),3);AME_hat_q2lc_cil_2

anios_esc 
   -0.037

AME_hat_q2lc_ciu_2 <-  round(c(AME_hat_q2lc_2 + 1.96 * AME_hat_q2lc_se_2),3);AME_hat_q2lc_ciu_2

anios_esc 
    0.033

AME_hat_q2lc_pr_2 <-  round(2 * pnorm(-abs(AME_hat_q2lc_2/AME_hat_q2lc_se_2)),3);AME_hat_q2lc_pr_2

anios_esc 
    0.912

AME_hat_q2lc_star_2<- " "  
anios_esc_q2<-cbind(AME_hat_q2lc_2,AME_hat_q2lc_star_2, AME_hat_q2lc_se_2,AME_hat_q2lc_pr_2,AME_hat_q2lc_cil_2,AME_hat_q2lc_ciu_2)

Edad

library(msm)
gamma_hatq3_3 <- coeficientes_1["edad"]
AME_hat_q3lc_3   <- round(mean(dnorm(ai)) * (gamma_hatq3_3),3);AME_hat_q3lc_3

 edad 
0.006

AME_hat_q3lc_se_3 <-  round(deltamethod(~ deriv1 * x3, coef(logitord), vcov(logitord), ses =  TRUE),3);AME_hat_q3lc_se_3

[1] 0.001

AME_hat_q3lc_cil_3 <-  round(c(AME_hat_q3lc_3 - 1.96 * AME_hat_q3lc_se_3),3);AME_hat_q3lc_cil_3

 edad 
0.004

AME_hat_q3lc_ciu_3 <-  round(c(AME_hat_q3lc_3 + 1.96 * AME_hat_q3lc_se_3),3);AME_hat_q3lc_ciu_3

 edad 
0.008

AME_hat_q3lc_pr_3 <-  round(2 * pnorm(-abs(AME_hat_q3lc_3/AME_hat_q3lc_se_3)),3);AME_hat_q3lc_pr_3

edad 
   0

AME_hat_q3lc_star_3<- "*** "  
edad_q3<-cbind(AME_hat_q3lc_3,AME_hat_q3lc_star_3, AME_hat_q3lc_se_3,AME_hat_q3lc_pr_3,AME_hat_q3lc_cil_3,AME_hat_q3lc_ciu_3)

Numero de Hijos

library(msm)
gamma_hatq4_4 <- coeficientes_1["t_hijos"]
AME_hat_q4lc_4   <- round(mean(dnorm(ai)) * (gamma_hatq4_4),3);AME_hat_q4lc_4

t_hijos 
  0.007

AME_hat_q4lc_se_4 <-  round(deltamethod(~ deriv1 * x4, coef(logitord), vcov(logitord), ses =  TRUE),3);AME_hat_q4lc_se_4

[1] 0.001

AME_hat_q4lc_cil_4 <-  round(c(AME_hat_q4lc_4 - 1.96 * AME_hat_q4lc_se_4),3);AME_hat_q4lc_cil_4

t_hijos 
  0.005

AME_hat_q4lc_ciu_4 <-  round(c(AME_hat_q4lc_4 + 1.96 * AME_hat_q4lc_se_4),3);AME_hat_q4lc_ciu_4

t_hijos 
  0.009

AME_hat_q4lc_pr_4 <-  round(2 * pnorm(-abs(AME_hat_q4lc_4/AME_hat_q4lc_se_4)),3);AME_hat_q4lc_pr_4

t_hijos 
      0

AME_hat_q4lc_star_4<- "*** "  
t_hijos_q4<-cbind(AME_hat_q4lc_4,AME_hat_q4lc_star_4, AME_hat_q4lc_se_4,AME_hat_q4lc_pr_4,AME_hat_q4lc_cil_4,AME_hat_q4lc_ciu_4)

Etnia

library(msm)
gamma_hatq5_5 <- coeficientes_1["etnia"]
AME_hat_q5lc_5   <- round(mean(dnorm(ai)) * (gamma_hatq5_5),3);AME_hat_q5lc_5

etnia 
0.061

AME_hat_q5lc_se_5 <-  round(deltamethod(~ deriv1 * x5, coef(logitord), vcov(logitord), ses =  TRUE),3);AME_hat_q5lc_se_5

[1] 0.001

AME_hat_q5lc_cil_5 <-  round(c(AME_hat_q5lc_5 - 1.96 * AME_hat_q5lc_se_5),3);AME_hat_q5lc_cil_5

etnia 
0.059

AME_hat_q5lc_ciu_5 <-  round(c(AME_hat_q5lc_5 + 1.96 * AME_hat_q5lc_se_5),3);AME_hat_q5lc_ciu_5

etnia 
0.063

AME_hat_q5lc_pr_5 <-  round(2 * pnorm(-abs(AME_hat_q5lc_5/AME_hat_q5lc_se_5)),3);AME_hat_q5lc_pr_5

etnia 
    0

AME_hat_q5lc_star_5<- "*** "  
etnia_q5<-cbind(AME_hat_q5lc_5,AME_hat_q5lc_star_5, AME_hat_q5lc_se_5,AME_hat_q5lc_pr_5,AME_hat_q5lc_cil_5,AME_hat_q5lc_ciu_5)

Area

library(msm)
gamma_hatq6_6 <- coeficientes_1["area"]
AME_hat_q6lc_6 <- round(mean(dnorm(ai)) * gamma_hatq6_6, 3)
AME_hat_q6lc_se_6 <- round(deltamethod(~ deriv1 * x6, coef(logitord), vcov(logitord), ses = TRUE), 3)
AME_hat_q6lc_cil_6 <- round(AME_hat_q6lc_6 - 1.96 * AME_hat_q6lc_se_6, 3)
AME_hat_q6lc_ciu_6 <- round(AME_hat_q6lc_6 + 1.96 * AME_hat_q6lc_se_6, 3)
AME_hat_q6lc_pr_6 <- round(2 * pnorm(-abs(AME_hat_q6lc_6 / AME_hat_q6lc_se_6)), 3)
area_q6 <- cbind(AME_hat_q6lc_6, "***", AME_hat_q6lc_se_6, AME_hat_q6lc_pr_6, AME_hat_q6lc_cil_6, AME_hat_q6lc_ciu_6)

##Unir y exportar efectos marginales a excel

resultados_ame <- rbind(lingrl_q1, anios_esc_q2, edad_q3, t_hijos_q4, etnia_q5, area_q6)
write.csv2(resultados_ame, file = "ame_logitord.csv")

####Matriz de Confusión y Exactitud del Modelo # Obtener predicciones del modelo logit ordenado

pred_logit_ord <- crossprod(t(X), coeficientes_1)

Convertir las probabilidades en categorías predichas (tomando la categoría con mayor probabilidad)

pred_clases <- apply(pred_logit_ord, 1, which.max)

Crear matriz de confusión

conf_matrix <- table(Predicho = pred_clases, Real = as.numeric(data$intensidad_dpp))

Mostrar matriz de confusión

print(conf_matrix)

        Real
Predicho     1     2     3
       1 12828  1790  1833

Calcular exactitud del modelo

exactitud <- sum(diag(conf_matrix)) / sum(conf_matrix)
cat("Exactitud del modelo Logit Ordenado:", exactitud, "\n")

Exactitud del modelo Logit Ordenado: 0.7797702

###Efectos marginales###### ## Probit

X <- cbind(1, probitord$mf[, -1])

coeficientes<- probitord$coefficients

coeficientes_1<- coeficientes[2:8]

ai <- crossprod(t(X), coeficientes_1) #(producto cruzado por los coeficientes)

deriv1<-mean(dnorm(ai)) #(saco derivada)