library(haven)
data <- read_dta("Data1_R.dta")
View(data)Clase 1 de junio
Setear el documento
setwd("C:/Users/JEFERSON CHAMBA/Desktop/Clase 1 de junio")REGRESIÓN MÚLTIPLE CON VARIABLE DISCRETA ORDENADA
Instalar paquetes
install.packages(“MASS”) Modelos Probit y Logit Ordenado
install.packages(“VGAM”) Modelos Logit y Probit Ordenado
install.packages(“margins”) Efectos marginales install.packages(“pROC”) Curva ROC
install.packages(“ggplot2”) Gráficos install.packages(“nnet”) Matriz de confusión
install.packages(“modelsummary”) Exportar a word install.packages(“pandoc”) Exportar a word
install.packages(“officer”)
install.packages(“flextable”)
install.packages(“broom”)
Cargar librerias
library(MASS)
library(VGAM)
library(margins) library(pROC)
library(ggplot2)
library(nnet)
library(haven)
library(modelsummary)
library(VGAM)
library(officer)
library(flextable)
library(broom)
library(Rchoice)
library(“foreign”)
library(“Rchoice”)
library(“msm”)
library(“car”)
Leer los datos
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)Ver distribución de la variable dependiente
table(data$intensidad_dpp)
1 2 3
12828 1790 1833 table(data$depresion_pp)
0 1
12828 3623 Ejemplo 1: Modelo con variable dependiente odenada
library("Rchoice")Cargando paquete requerido: FormulaCargando paquete requerido: maxLikCargando 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/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 20:09:05 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 tanto no permite conocer su efecto sobre la necesidad de depresión.
a medida que aumenta un año de esolcaridad reduce la intensidad de depresión post parto, dado que es estadisticamente significativo
cuando las mujeres aumentan en un año de edad aumenta la depresión post parto
a medida que tiene un hijo mas es probable que tenga depresión post parto
la mujer indigene tiene mas probabilidad de tener depresión post parto
mujeres que viven en areas rurales tienen mas probabilidad de sufrir depresión post parto
Exportar información a word
library(modelsummary)
modelsummary(logitord, output = "modelo_logit_ord.docx")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)Criterio de información de Akaike (AIC) y criterio de Información Bayesiano (BIC)
AIC(logitord)[1] 22123.39BIC(logitord)[1] 22185.05Ajustar modelo Probit Ordenado
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 20:09:13 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:2s
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)Comando para exportar dos resultados en una misma tabla
library(modelsummary)
modelsummary(list("logit" = logitord, "Probit" = probitord),
output = "modelo_logit.docx",
stars = TRUE)AIC(probitord)[1] 22109.98BIC(probitord)[1] 22171.64Comparar AIC y BIC de ambos modelos
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 Análisis: El modelo con menor AIC/BIC es el preferido, ya que tiene mejor ajuste.Por lo tanto, el modelo probit ordenado es el mejor.