library(haven)
data <- read_dta("Data1_R.dta")
View(data)clase 01 de junio
REGRESIÓN MÚLTIPLE CON VARIABLE DISCRETA ORDENADA
install.packages(“MASS”)
install.packages(“VGAM”)
install.packages(“margins”)
install.packages(“pROC”)
install.packages(“ggplot2”)
install.packages(“nnet”)
install.packages(“modelsummary”)
install.packages(“pandoc”)
install.packages(“officer”)
install.packages(“flextable”)
install.packages(“broom”)
install.packages(c(“officer”, “flextable”, “broom”))
install.packages(“Rchoice”)
Cargar librerías
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 (ajustar la ruta del archivo)
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 Analisis
indican que 12828 observaciones no han tenido depresion; mientras que 3623 si han presenato depresion
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: ju. jun. 05 G-InnoVa 21:22:51 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 no tiene ninguna influencia sobre la depresión post parto, no hay evidencia estadística.
Los años de escolaridad cuando aumentan, reduce la probabilidad de sufrir de depresión posparto, es estadisticamente significativo.
La edad aumenta la probabilidad de sufrir depresión post parto
A medida que tienen un hijo mas, es mayormente probable que sufran depresion postparto.
Las mujeres de la etnia indigena tienen mayor probabilidad de padecer de depresion postparto.
Las mujeres que habitan en zonas rurales tienen mayor probabilidad de sufrir depresion postparto.
Exportan 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)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'))Mostrar resumen de los modelos
summary(probitord)
Model: ordinal
Model estimated on: ju. jun. 05 G-InnoVa 21:22:56 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_ord4.docx",
statistic = "{statistic}",
stars = TRUE)Comando para unir los modelos
library(modelsummary)
modelsummary(list("Logit" = logitord, "Probit" = probitord),
output = "modelos_logit_probit.docx",
stars = TRUE)Comparar 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 Analisis
El modelo con menor AIC/BIC es el preferido, Ya que tiene mejor ajuste, por lo tanto, el modelo probit ordenado es el mejor.
Efectos marginales
PROBIT
X <- cbind(1, probitord$mf[, -1])coeficientes<- probitord$coefficientscoeficientes_1<- coeficientes[2:8]ai <- crossprod(t(X), coeficientes_1)deriv1<-mean(dnorm(ai))INGRESO
library(msm)
gamma_hatq1_1 <- coeficientes_1["lingrl"]
AME_hat_q1lc_1 <- round(mean(dnorm(ai)) * (gamma_hatq1_1),3);AME_hat_q1lc_1lingrl
0 AME_hat_q1lc_se_1 <- round(deltamethod(~ deriv1 * x1, coef(logitord), vcov(logitord), ses = TRUE),3);AME_hat_q1lc_se_1 [1] 0.005AME_hat_q1lc_cil_1 <- round(c(AME_hat_q1lc_1 - 1.96 * AME_hat_q1lc_se_1),3);AME_hat_q1lc_cil_1lingrl
-0.01 AME_hat_q1lc_ciu_1 <- round(c(AME_hat_q1lc_1 + 1.96 * AME_hat_q1lc_se_1),3);AME_hat_q1lc_ciu_1lingrl
0.01 AME_hat_q1lc_pr_1 <- round(2 * pnorm(-abs(AME_hat_q1lc_1/AME_hat_q1lc_se_1)),3);AME_hat_q1lc_pr_1lingrl
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_2anios_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.029AME_hat_q2lc_cil_2 <- round(c(AME_hat_q2lc_2 - 1.96 * AME_hat_q2lc_se_2),3);AME_hat_q2lc_cil_2anios_esc
-0.059 AME_hat_q2lc_ciu_2 <- round(c(AME_hat_q2lc_2 + 1.96 * AME_hat_q2lc_se_2),3);AME_hat_q2lc_ciu_2anios_esc
0.055 AME_hat_q2lc_pr_2 <- round(2 * pnorm(-abs(AME_hat_q2lc_2/AME_hat_q2lc_se_2)),3);AME_hat_q2lc_pr_2anios_esc
0.945 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.002AME_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.002 AME_hat_q3lc_ciu_3 <- round(c(AME_hat_q3lc_3 + 1.96 * AME_hat_q3lc_se_3),3);AME_hat_q3lc_ciu_3edad
0.01 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.003 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_4t_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.001AME_hat_q4lc_cil_4 <- round(c(AME_hat_q4lc_4 - 1.96 * AME_hat_q4lc_se_4),3);AME_hat_q4lc_cil_4t_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_4t_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_4t_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_5etnia
0.057 AME_hat_q5lc_se_5 <- round(deltamethod(~ deriv1 * x5, coef(logitord), vcov(logitord), ses = TRUE),3);AME_hat_q5lc_se_5 [1] 0.001AME_hat_q5lc_cil_5 <- round(c(AME_hat_q5lc_5 - 1.96 * AME_hat_q5lc_se_5),3);AME_hat_q5lc_cil_5etnia
0.055 AME_hat_q5lc_ciu_5 <- round(c(AME_hat_q5lc_5 + 1.96 * AME_hat_q5lc_se_5),3);AME_hat_q5lc_ciu_5etnia
0.059 AME_hat_q5lc_pr_5 <- round(2 * pnorm(-abs(AME_hat_q5lc_5/AME_hat_q5lc_se_5)),3);AME_hat_q5lc_pr_5etnia
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_6 area
0.021 AME_hat_q6lc_se_6 <- round(deltamethod(~ deriv1 * x6, coef(logitord), vcov(logitord), ses = TRUE),3);AME_hat_q6lc_se_6 [1] 0.005AME_hat_q6lc_cil_6 <- round(c(AME_hat_q6lc_6 - 1.96 * AME_hat_q6lc_se_6),3);AME_hat_q6lc_cil_6 area
0.011 AME_hat_q6lc_ciu_6 <- round(c(AME_hat_q6lc_6 + 1.96 * AME_hat_q6lc_se_6),3);AME_hat_q6lc_ciu_6 area
0.031 AME_hat_q6lc_pr_6 <- round(2 * pnorm(-abs(AME_hat_q6lc_6/AME_hat_q6lc_se_6)),3);AME_hat_q6lc_pr_6area
0 AME_hat_q6lc_star_6<- "*** "
area_q6<-cbind(AME_hat_q6lc_6,AME_hat_q6lc_star_6, AME_hat_q6lc_se_6,AME_hat_q6lc_pr_6,AME_hat_q6lc_cil_6,AME_hat_q6lc_ciu_6)
resultados_ame<-rbind(lingrl_q1,anios_esc_q2, edad_q3, t_hijos_q4,etnia_q5,area_q6);resultados_ame AME_hat_q1lc_1 AME_hat_q1lc_star_1 AME_hat_q1lc_se_1
lingrl "0" " " "0.005"
anios_esc "-0.002" " " "0.029"
edad "0.006" "*** " "0.002"
t_hijos "0.007" "*** " "0.001"
etnia "0.057" "*** " "0.001"
area "0.021" "*** " "0.005"
AME_hat_q1lc_pr_1 AME_hat_q1lc_cil_1 AME_hat_q1lc_ciu_1
lingrl "1" "-0.01" "0.01"
anios_esc "0.945" "-0.059" "0.055"
edad "0.003" "0.002" "0.01"
t_hijos "0" "0.005" "0.009"
etnia "0" "0.055" "0.059"
area "0" "0.011" "0.031" Exportar efectos marginales a excel
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 1833Valores predichos
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 Analisis
El 77.97%, lo cual significa que el modelo predijo en un 77.97% las observaciones.
MODELO PROBIT
Efectos marginales Probit
X <- cbind(1, probitord$mf[, -1])coeficientes<- probitord$coefficientscoeficientes_1<- coeficientes[2:8]ai <- crossprod(t(X), coeficientes_1)deriv1<-mean(dnorm(ai))INGRESO
library(msm)
gamma_hatq1_1 <- coeficientes_1["lingrl"]
AME_hat_q1lc_1 <- round(mean(dnorm(ai)) * (gamma_hatq1_1),3);AME_hat_q1lc_1lingrl
0 AME_hat_q1lc_se_1 <- round(deltamethod(~ deriv1 * x1, coef(probitord), vcov(probitord), ses = TRUE),3);AME_hat_q1lc_se_1 [1] 0.003AME_hat_q1lc_cil_1 <- round(c(AME_hat_q1lc_1 - 1.96 * AME_hat_q1lc_se_1),3);AME_hat_q1lc_cil_1lingrl
-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_1lingrl
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_1lingrl
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_2anios_esc
-0.002 AME_hat_q2lc_se_2 <- round(deltamethod(~ deriv1 * x2, coef(probitord), vcov(probitord), ses = TRUE),3);AME_hat_q2lc_se_2 [1] 0.017AME_hat_q2lc_cil_2 <- round(c(AME_hat_q2lc_2 - 1.96 * AME_hat_q2lc_se_2),3);AME_hat_q2lc_cil_2anios_esc
-0.035 AME_hat_q2lc_ciu_2 <- round(c(AME_hat_q2lc_2 + 1.96 * AME_hat_q2lc_se_2),3);AME_hat_q2lc_ciu_2anios_esc
0.031 AME_hat_q2lc_pr_2 <- round(2 * pnorm(-abs(AME_hat_q2lc_2/AME_hat_q2lc_se_2)),3);AME_hat_q2lc_pr_2anios_esc
0.906 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(probitord), vcov(probitord), ses = TRUE),3);AME_hat_q3lc_se_3 [1] 0.001AME_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_3edad
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_4t_hijos
0.007 AME_hat_q4lc_se_4 <- round(deltamethod(~ deriv1 * x4, coef(probitord), vcov(probitord), ses = TRUE),3);AME_hat_q4lc_se_4 [1] 0.001AME_hat_q4lc_cil_4 <- round(c(AME_hat_q4lc_4 - 1.96 * AME_hat_q4lc_se_4),3);AME_hat_q4lc_cil_4t_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_4t_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_4t_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_5etnia
0.057 AME_hat_q5lc_se_5 <- round(deltamethod(~ deriv1 * x5, coef(probitord), vcov(probitord), ses = TRUE),3);AME_hat_q5lc_se_5 [1] 0.001AME_hat_q5lc_cil_5 <- round(c(AME_hat_q5lc_5 - 1.96 * AME_hat_q5lc_se_5),3);AME_hat_q5lc_cil_5etnia
0.055 AME_hat_q5lc_ciu_5 <- round(c(AME_hat_q5lc_5 + 1.96 * AME_hat_q5lc_se_5),3);AME_hat_q5lc_ciu_5etnia
0.059 AME_hat_q5lc_pr_5 <- round(2 * pnorm(-abs(AME_hat_q5lc_5/AME_hat_q5lc_se_5)),3);AME_hat_q5lc_pr_5etnia
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_6 area
0.021 AME_hat_q6lc_se_6 <- round(deltamethod(~ deriv1 * x6, coef(probitord), vcov(probitord), ses = TRUE),3);AME_hat_q6lc_se_6 [1] 0.003AME_hat_q6lc_cil_6 <- round(c(AME_hat_q6lc_6 - 1.96 * AME_hat_q6lc_se_6),3);AME_hat_q6lc_cil_6 area
0.015 AME_hat_q6lc_ciu_6 <- round(c(AME_hat_q6lc_6 + 1.96 * AME_hat_q6lc_se_6),3);AME_hat_q6lc_ciu_6 area
0.027 AME_hat_q6lc_pr_6 <- round(2 * pnorm(-abs(AME_hat_q6lc_6/AME_hat_q6lc_se_6)),3);AME_hat_q6lc_pr_6area
0 AME_hat_q6lc_star_6<- "*** "
area_q6<-cbind(AME_hat_q6lc_6,AME_hat_q6lc_star_6, AME_hat_q6lc_se_6,AME_hat_q6lc_pr_6,AME_hat_q6lc_cil_6,AME_hat_q6lc_ciu_6)Modelo probit
data$intensidad_dpp <- ordered(data$intensidad_dpp)library(MASS)
probitord <- polr(intensidad_dpp ~ lingrl + anios_esc + edad + t_hijos + etnia + area,
data = data,
method = "probit",
Hess = TRUE)Mostrar resumen de los modelos
AIC(probitord)[1] 22109.98BIC(probitord)[1] 22171.64Efectos marginales
pred_clases <- predict(probitord, type = "class")conf_matrix <- table(Predicho = pred_clases, Real = data$intensidad_dpp)
print(conf_matrix) Real
Predicho 1 2 3
1 12828 1790 1833
2 0 0 0
3 0 0 0exactitud <- sum(diag(conf_matrix)) / sum(conf_matrix)
cat("Exactitud del modelo Probit Ordenado:", exactitud, "\n")Exactitud del modelo Probit Ordenado: 0.7797702