library(haven)
data <- read_dta("Data1_R.dta")
View(data)Clase 01 de junio
REGRESIÓN MÚLTIPLE CON VARIABLE DISCRETA ORDENADA
setwd(“C:/Users/Garciadaniel/Desktop/CLASE 01 DE JUNIO”)
Carga de librerías
install.packages(c(“officer”, “flextable”, “broom”))
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”)
Carga de base de datos
Análisis de Data
Visualización de primeras filas del dataset
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>, …Distribución de variable dependiente
table(data$intensidad_dpp)
1 2 3
12828 1790 1833 table(data$depresion_pp)
0 1
12828 3623 Modelo con variable dependiente odenada
Ajuste del modelo logit ordenado
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'))Resumen modelo logit ordenado
summary(logitord)
Model: ordinal
Model estimated on: dom. jun. 08 12:49:28 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.
Cuando aumentan los años de escolaridad, se reduce la probabilidad de sufrir de depresión, se desconoce la magnitud.
La edad aumenta la probabilidad de sufrir depresión post parto.
El número de hijos aumenta la probabilidad de sufrir depresión post parto.
Las mujeres indígenas tienen mayor probabilidad de presentar depresión post parto.Similar a las mujeres de zonas rurales.
Se exportan los resultados a un documento Word
Resultados generales
library(modelsummary)
modelsummary(logitord, output = "modelo_logit_ord.docx")Se incluye nivel de significancia (*)
library(modelsummary)
modelsummary(logitord,
output = "modelo_logit_ord2.docx",
stars = c('*' = .05, '**' = .01, '***' = .001))Reporte incluyendo z-value (p value)
library(modelsummary)
modelsummary(logitord,
output = "modelo_logit_ord3.docx",
statistic = c("p.value" = "p = {p.value}{stars}"),
stars = TRUE)Reporte académico
library(modelsummary)
modelsummary(logitord,
output = "modelo_logit_ord4.docx",
statistic = "{statistic}",
stars = TRUE)En este último reporte se puede analizar:
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(logitord) [1] 22123.39BIC(logitord)[1] 22185.05Ajuste 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'))Resumen modelo probit ordenado
summary(probitord)
Model: ordinal
Model estimated on: dom. jun. 08 12:49:34 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 Reporte resultados del modelo
library(modelsummary)
modelsummary(probitord,
output = "modelo_probitord_ord4.docx",
statistic = "{statistic}",
stars = TRUE)Comando para reportar resultados en un mismo documento
library(modelsummary)
modelsummary(list("Logit" = logitord, "Probit" = probitord),
output = "modelos_logit_probit.docx",
stars = TRUE)Criterio de información de Akaike (AIC) y criterio de Información Bayesiano (BIC)
AIC(probitord) [1] 22109.98BIC(probitord)[1] 22171.64Comparación AIC - BIC de ambos modelos
aic_logit_ord <- AIC(logitord)
aic_probit_ord <- AIC(probitord)bic_logit_ord <- BIC(logitord)
bic_probit_ord <- BIC(probitord)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, ya que tiene mejor ajuste.En este caso, el modelo probit ordenado es el mejor.
Cálculo de efectos marginales modelo logit
X <- cbind(1, logitord$mf[, -1])coeficientes<- logitord$coefficientscoeficientes_1<- coeficientes[2:8]Producto cruzado por los coeficientes
ai <- crossprod(t(X), coeficientes_1) Derivada
deriv1<-mean(dnorm(ai)) Variable 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.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)Variable 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.018AME_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.037 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.033 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.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)Variable 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.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)Variable 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)Variable 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.061 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.059 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.063 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)Variable área
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.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)Resultados
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.003"
anios_esc "-0.002" " " "0.018"
edad "0.006" "*** " "0.001"
t_hijos "0.007" "*** " "0.001"
etnia "0.061" "*** " "0.001"
area "0.021" "*** " "0.003"
AME_hat_q1lc_pr_1 AME_hat_q1lc_cil_1 AME_hat_q1lc_ciu_1
lingrl "1" "-0.006" "0.006"
anios_esc "0.912" "-0.037" "0.033"
edad "0" "0.004" "0.008"
t_hijos "0" "0.005" "0.009"
etnia "0" "0.059" "0.063"
area "0" "0.015" "0.027" Se exportan los efectos a archivo Excel
write.csv2(resultados_ame,file="ame_logitord.csv")Matriz de Confusión y Exactitud del Modelo
Se obtienen predicciones del modelo logit ordenado
pred_logit_ord <- crossprod(t(X), coeficientes_1)Se convierten las probabilidades en categorías predichas (tomando la categoría con mayor probabilidad)
pred_clases <- apply(pred_logit_ord, 1, which.max)Matriz de confusión
conf_matrix <- table(Predicho = pred_clases, Real = as.numeric(data$intensidad_dpp))Mostramos la matriz de confusión
print(conf_matrix) Real
Predicho 1 2 3
1 12828 1790 1833Calcular 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 Cálculo de efectos marginales modelo probit
Y <- cbind(1, probitord$mf[, -1])coeficientes_probit<- probitord$coefficientscoeficientes_2<- coeficientes_probit[2:8]Producto cruzado por los coeficientes
ai_p <- crossprod(t(Y), coeficientes_2) Derivada
deriv2<-mean(dnorm(ai_p)) Variable ingreso
library(msm)
gamma_hatq1_7 <- coeficientes_2["lingrl"]
AME_hat_q1lc_7 <- round(mean(dnorm(ai_p)) * (gamma_hatq1_7),3);7[1] 7AME_hat_q1lc_se_7 <- round(deltamethod(~ deriv2 * x1, coef(probitord), vcov(probitord), ses = TRUE),3);AME_hat_q1lc_se_7 [1] 0.003AME_hat_q1lc_cil_7 <- round(c(AME_hat_q1lc_7 - 1.96 * AME_hat_q1lc_se_7),3);AME_hat_q1lc_cil_7lingrl
-0.006 AME_hat_q1lc_ciu_7 <- round(c(AME_hat_q1lc_7 + 1.96 * AME_hat_q1lc_se_7),3);AME_hat_q1lc_ciu_7lingrl
0.006 AME_hat_q1lc_pr_7 <- round(2 * pnorm(-abs(AME_hat_q1lc_7/AME_hat_q1lc_se_7)),3);AME_hat_q1lc_pr_7lingrl
1 AME_hat_q1lc_star_7<- " "
lingrl_q7<-cbind(AME_hat_q1lc_7,AME_hat_q1lc_star_7, AME_hat_q1lc_se_7,AME_hat_q1lc_pr_7,AME_hat_q1lc_cil_7,AME_hat_q1lc_ciu_7)Variable escolaridad
library(msm)
gamma_hatq2_7 <- coeficientes_2["anios_esc"]
AME_hat_q2lc_7 <- round(mean(dnorm(ai_p)) * (gamma_hatq2_7),3);AME_hat_q2lc_7anios_esc
-0.002 AME_hat_q2lc_se_7 <- round(deltamethod(~ deriv2 * x2, coef(probitord), vcov(probitord), ses = TRUE),3);AME_hat_q2lc_se_7 [1] 0.017AME_hat_q2lc_cil_7 <- round(c(AME_hat_q2lc_7 - 1.96 * AME_hat_q2lc_se_7),3);AME_hat_q2lc_cil_7anios_esc
-0.035 AME_hat_q2lc_ciu_7 <- round(c(AME_hat_q2lc_7 + 1.96 * AME_hat_q2lc_se_7),3);AME_hat_q2lc_ciu_7anios_esc
0.031 AME_hat_q2lc_pr_7 <- round(2 * pnorm(-abs(AME_hat_q2lc_7/AME_hat_q2lc_se_7)),3);AME_hat_q2lc_pr_7anios_esc
0.906 AME_hat_q2lc_star_7<- " "
anios_esc_q7<-cbind(AME_hat_q2lc_7,AME_hat_q2lc_star_7, AME_hat_q2lc_se_7,AME_hat_q2lc_pr_7,AME_hat_q2lc_cil_7,AME_hat_q2lc_ciu_7)Variable edad
library(msm)
gamma_hatq3_7 <- coeficientes_2["edad"]
AME_hat_q3lc_7 <- round(mean(dnorm(ai_p)) * (gamma_hatq3_7),3);AME_hat_q3lc_7 edad
0.006 AME_hat_q3lc_se_7 <- round(deltamethod(~ deriv2 * x3, coef(probitord), vcov(probitord), ses = TRUE),3);AME_hat_q3lc_se_7[1] 0.001AME_hat_q3lc_cil_7 <- round(c(AME_hat_q3lc_7 - 1.96 * AME_hat_q3lc_se_7),3);AME_hat_q3lc_cil_7 edad
0.004 AME_hat_q3lc_ciu_7 <- round(c(AME_hat_q3lc_7 + 1.96 * AME_hat_q3lc_se_7),3);AME_hat_q3lc_ciu_7 edad
0.008 AME_hat_q3lc_pr_7 <- round(2 * pnorm(-abs(AME_hat_q3lc_7/AME_hat_q3lc_se_7)),3);AME_hat_q3lc_pr_7edad
0 AME_hat_q3lc_star_7<- "*** "
edad_q7<-cbind(AME_hat_q3lc_7,AME_hat_q3lc_star_7, AME_hat_q3lc_se_7,AME_hat_q3lc_pr_7,AME_hat_q3lc_cil_7,AME_hat_q3lc_ciu_7)Variable hijos
library(msm)
gamma_hatq4_7 <- coeficientes_2["t_hijos"]
AME_hat_q4lc_7 <- round(mean(dnorm(ai_p)) * (gamma_hatq4_7),3);AME_hat_q4lc_7t_hijos
0.007 AME_hat_q4lc_se_7 <- round(deltamethod(~ deriv2 * x4, coef(probitord), vcov(probitord), ses = TRUE),3);AME_hat_q4lc_se_7[1] 0.001AME_hat_q4lc_cil_7 <- round(c(AME_hat_q4lc_7 - 1.96 * AME_hat_q4lc_se_7),3);AME_hat_q4lc_cil_7t_hijos
0.005 AME_hat_q4lc_ciu_7 <- round(c(AME_hat_q4lc_7 + 1.96 * AME_hat_q4lc_se_7),3);AME_hat_q4lc_ciu_7t_hijos
0.009 AME_hat_q4lc_pr_7 <- round(2 * pnorm(-abs(AME_hat_q4lc_7/AME_hat_q4lc_se_7)),3);AME_hat_q4lc_pr_7t_hijos
0 AME_hat_q4lc_star_7<- "*** "
t_hijos_q7<-cbind(AME_hat_q4lc_7,AME_hat_q4lc_star_7, AME_hat_q4lc_se_7,AME_hat_q4lc_pr_7,AME_hat_q4lc_cil_7,AME_hat_q4lc_ciu_7)Variable etnia
library(msm)
gamma_hatq5_7 <- coeficientes_2["etnia"]
AME_hat_q5lc_7 <- round(mean(dnorm(ai_p)) * (gamma_hatq5_7),3);AME_hat_q5lc_7etnia
0.057 AME_hat_q5lc_se_7 <- round(deltamethod(~ deriv2 * x5, coef(probitord), vcov(probitord), ses = TRUE),3);AME_hat_q5lc_se_7[1] 0.001AME_hat_q5lc_cil_7 <- round(c(AME_hat_q5lc_7 - 1.96 * AME_hat_q5lc_se_7),3);AME_hat_q5lc_cil_7etnia
0.055 AME_hat_q5lc_ciu_7 <- round(c(AME_hat_q5lc_7 + 1.96 * AME_hat_q5lc_se_7),3);AME_hat_q5lc_ciu_7etnia
0.059 AME_hat_q5lc_pr_7 <- round(2 * pnorm(-abs(AME_hat_q5lc_7/AME_hat_q5lc_se_7)),3);AME_hat_q5lc_pr_7etnia
0 AME_hat_q5lc_star_7<- "*** "
etnia_q7<-cbind(AME_hat_q5lc_7,AME_hat_q5lc_star_7, AME_hat_q5lc_se_7,AME_hat_q5lc_pr_7,AME_hat_q5lc_cil_7,AME_hat_q5lc_ciu_7)Variable área
library(msm)
gamma_hatq6_7 <- coeficientes_2["area"]
AME_hat_q6lc_7 <- round(mean(dnorm(ai_p)) * (gamma_hatq6_7),3);AME_hat_q6lc_7 area
0.021 AME_hat_q6lc_se_7 <- round(deltamethod(~ deriv2 * x6, coef(probitord), vcov(probitord), ses = TRUE),3);AME_hat_q6lc_se_7[1] 0.003AME_hat_q6lc_cil_7 <- round(c(AME_hat_q6lc_7 - 1.96 * AME_hat_q6lc_se_7),3);AME_hat_q6lc_cil_7 area
0.015 AME_hat_q6lc_ciu_7 <- round(c(AME_hat_q6lc_7 + 1.96 * AME_hat_q6lc_se_7),3);AME_hat_q6lc_ciu_7 area
0.027 AME_hat_q6lc_pr_7 <- round(2 * pnorm(-abs(AME_hat_q6lc_7/AME_hat_q6lc_se_7)),3);AME_hat_q6lc_pr_7area
0 AME_hat_q6lc_star_7<- "*** "
area_q7<-cbind(AME_hat_q6lc_7,AME_hat_q6lc_star_7, AME_hat_q6lc_se_7,AME_hat_q6lc_pr_7,AME_hat_q6lc_cil_7,AME_hat_q6lc_ciu_7)Resultados
resultados_ame_probit<-rbind(lingrl_q7,anios_esc_q7, edad_q7, t_hijos_q7,etnia_q7,area_q7);resultados_ame_probit AME_hat_q1lc_7 AME_hat_q1lc_star_7 AME_hat_q1lc_se_7
lingrl "0" " " "0.003"
anios_esc "-0.002" " " "0.017"
edad "0.006" "*** " "0.001"
t_hijos "0.007" "*** " "0.001"
etnia "0.057" "*** " "0.001"
area "0.021" "*** " "0.003"
AME_hat_q1lc_pr_7 AME_hat_q1lc_cil_7 AME_hat_q1lc_ciu_7
lingrl "1" "-0.006" "0.006"
anios_esc "0.906" "-0.035" "0.031"
edad "0" "0.004" "0.008"
t_hijos "0" "0.005" "0.009"
etnia "0" "0.055" "0.059"
area "0" "0.015" "0.027" Se exportan los efectos a archivo Excel
write.csv2(resultados_ame_probit,file="ame_probitord.csv")Matriz de Confusión y Exactitud del Modelo
Se obtienen predicciones del modelo logit ordenado
pred_probit_ord <- crossprod(t(Y), coeficientes_2)Se convierten las probabilidades en categorías predichas (tomando la categoría con mayor probabilidad)
pred_clases_p <- apply(pred_probit_ord, 1, which.max)Matriz de confusión
conf_matrix_p <- table(Predicho = pred_clases_p, Real = as.numeric(data$intensidad_dpp))Mostramos la matriz de confusión
print(conf_matrix_p) Real
Predicho 1 2 3
1 12828 1790 1833Calcular exactitud del modelo
exactitud_probit <- sum(diag(conf_matrix_p)) / sum(conf_matrix_p)
cat("Exactitud del modelo Probit Ordenado:", exactitud_probit, "\n")Exactitud del modelo Probit Ordenado: 0.7797702