Deber_07_junio

Author

LuisR_LeninJ_07_junio

1 REGRESIÓN MÚLTIPLE CON VARIABLE DISCRETA ORDENADA

# Cargar librerías
library(MASS)
library(VGAM)

Cargando paquete requerido: stats4

Cargando paquete requerido: splines

library(margins)
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(haven)
library(modelsummary)


Adjuntando el paquete: 'modelsummary'

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

    Max

library(officer)
library(flextable)
library(broom)
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/

library(msm)
library(car)

Cargando paquete requerido: carData


Adjuntando el paquete: 'car'

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

    logit

# Leer los datos (ajustar la ruta del archivo)
data <- read_dta("Data1_R.dta")

# Ver distribución de la variable dependiente
table(data$intensidad_dpp)  # variable ordinal


    1     2     3 
12828  1790  1833

2 Ejemplo 1: Modelo Probit Ordenado

# Estimación del 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 11:23:36 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

# Exportar resumen del modelo
modelsummary(probitord,
             output = "modelo_probit_ord.docx",
             statistic = "{statistic}",
             stars = TRUE)

# Criterios de información
AIC(probitord)

[1] 22109.98

BIC(probitord)

[1] 22171.64

El modelo probitord tiene un AIC de 22109.98 y un BIC de 22171.64, lo que sirve como referencia para comparar con modelos alternativos. El modelo más adecuado será aquel con los valores de AIC y BIC más bajos entre todos los modelos considerados.

3 Efectos marginales promedio (AME)

# Preparar X y coeficientes
X <- model.matrix(probitord)[, -1]
coeficientes <- coef(probitord)
coeficientes_1 <- coeficientes[1:ncol(X)]
ai <- as.vector(X %*% coeficientes_1)
deriv1 <- mean(dnorm(ai))

# Función de AME
calc_marginal <- function(pos) {
  gamma <- coeficientes_1[pos]
  ame <- round(deriv1 * gamma, 3)
  formula_str <- paste0("~", deriv1, "*x", pos)
  se <- tryCatch(
    round(deltamethod(as.formula(formula_str), coeficientes_1, vcov(probitord)[1:length(coeficientes_1), 1:length(coeficientes_1)]), 3),
    error = function(e) NA
  )
  cil <- round(ame - 1.96 * se, 3)
  ciu <- round(ame + 1.96 * se, 3)
  pval <- round(2 * pnorm(-abs(ame / se)), 3)
  sig <- ifelse(pval < 0.001, "***", ifelse(pval < 0.01, "**", ifelse(pval < 0.05, "*", " ")))
  cbind(AME = ame, Sig = sig, SE = se, p = pval, CI_Low = cil, CI_Up = ciu)
}

# Calcular AME
nvars <- length(coeficientes_1)
resultados_ame <- do.call(rbind, lapply(1:nvars, calc_marginal))
rownames(resultados_ame) <- colnames(X)

# Mostrar y guardar
print(resultados_ame)

          AME      Sig   SE  p     CI_Low   CI_Up   
lingrl    "0.002"  "***" "0" "0"   "0.002"  "0.002" 
anios_esc "-0.005" "***" "0" "0"   "-0.005" "-0.005"
edad      "0"      NA    "0" "NaN" "0"      "0"     
t_hijos   "0"      NA    "0" "NaN" "0"      "0"     
etnia     "0"      NA    "0" "NaN" "0"      "0"     
area      "0"      NA    "0" "NaN" "0"      "0"

write.csv2(resultados_ame, file = "ame_probitord.csv")

El modelo probit ordenado muestra que solo lingrl (AME = 0.002) y anios_esc (AME = -0.005) tienen efectos marginales significativos (p < 0.001), mientras que las demás variables no presentan efectos detectables.

4 Predicción y matriz de confusión

# Definir pred_probit_ord antes de usarlo
pred_probit_ord <- as.vector(X %*% coeficientes_1)

# Extraer umbrales
umbrales <- coef(probitord)[grepl("^tau", names(coef(probitord)))]
umbrales <- sort(umbrales)

# Clasificar predicciones
pred_clases <- cut(pred_probit_ord,
                   breaks = c(-Inf, umbrales, Inf),
                   labels = 1:(length(umbrales) + 1),
                   right = TRUE)

# Real
real_clases <- as.numeric(as.character(data$intensidad_dpp))

# Matriz de confusión y exactitud
conf_matrix <- table(Predicho = pred_clases, Real = real_clases)
print(conf_matrix)

        Real
Predicho     1     2     3
       1 12828  1790  1833

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

Exactitud del modelo Probit Ordenado: 0.78

La matriz de confusión muestra que el modelo predice correctamente los casos de menor intensidad de depresión postparto, pero tiende a confundir los de categorías superiores. La exactitud general del modelo es del 78%, aunque se recomienda revisar otras métricas como precisión o sensibilidad debido al posible desbalance en las clases.