Regresión logística bayesiana
Rojas Rojas Henry
2023-09-23
La base de datos Trabajo contiene la información de 30 personas del distrito del Porvenir, las variables se presentan a continuación:
Y: Condición de trabajo ** 1: Trabaja ** 0: No trabaja X1: Casada: ** 0: No ** 1: Si X2: Idade: Edad de la persona X3: Escolaridad: Años de escolaridad de la persona
ASUMIENDO NORMALIDAD
CONSIDERANDO LOS SIGUIENTES DATOS
datos<-read.table("trabajo.txt",header = T,sep = "\t")
attach(datos)
head(datos,5)## Obs Y casada idade escolaridade
## 1 1 1 0 31 16
## 2 2 1 1 34 14
## 3 3 1 1 41 16
## 4 4 0 0 67 9
## 5 5 1 0 25 12NOS PIDEN
- Utilizar la base de datos Trabajo con la variable dependiente Y= Condición de trabajo y la variable predictora X= edad de la persona.
- A. Estimar los parámetros con la ecuación de regresión logística binaria clásica.
# Ajustar el modelo de regresión logística binaria
modelo_logistico1 <- glm(Y ~ idade, data = datos, family = binomial)
# Resumen del modelo
summary(modelo_logistico1)##
## Call:
## glm(formula = Y ~ idade, family = binomial, data = datos)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.90310 1.38839 1.371 0.170
## idade -0.03583 0.03165 -1.132 0.258
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 40.381 on 29 degrees of freedom
## Residual deviance: 39.038 on 28 degrees of freedom
## AIC: 43.038
##
## Number of Fisher Scoring iterations: 4Nota. Se observa que, la variable idade es no significativa (p-value = 0.258 > 5%), lo que da a entender que no influye o repercute en la condición de trabajo; Sin embargo, por motivos de práctica seguimos con el análisis.
- B. Determinar la bondad de ajuste del modelo.
#Test de Hosmer-Lemeshow
library(ResourceSelection)## ResourceSelection 0.3-6 2023-06-27hoslem.test(Y, fitted(modelo_logistico1), g=10)##
## Hosmer and Lemeshow goodness of fit (GOF) test
##
## data: Y, fitted(modelo_logistico1)
## X-squared = 18.806, df = 8, p-value = 0.01593Nota. Según la prueba de Hosmer-Lemeshow, el modelo es adecuado, para un nivel de signficancia del 10%.
#Pseudo-R2
library(DescTools)
PseudoR2(modelo_logistico1,which = c("McFadden","CoxSnell","Nagelkerke"))## McFadden CoxSnell Nagelkerke
## 0.03325252 0.04377177 0.05917289Nota. observamos que el valor de Rcuadrado de Cox y Snell es igual a 0.044 y el valor de Rcuadrado de Nagelkerke es igual a 0.059, el cual nos indica que el 5.9% de las proporciones de la varianza de la variable dependiente (Condición de trabajo) es explicada por la variable predictora edad; dando a entender que el modelo de regresión logística no es bueno.
- C. Hallar el coeficiente OR e interpretar.
exp(cbind(OR = coef(modelo_logistico1), confint(modelo_logistico1, level=0.95)))## Waiting for profiling to be done...## OR 2.5 % 97.5 %
## (Intercept) 6.7066208 0.4946832 130.440833
## idade 0.9648037 0.9020713 1.024774Nota. Nos indica la posibilidad de que un individuo trabaje, por su edad, es decir, es aproximadamente 0.96 veces más probable que el individuo trabaje por tener menos edad, respecto a los que tienen más edad.
#Precisión del modelo
y.pred<-ifelse(modelo_logistico1$fitted.values > 0.43, yes = 1, no = 0)
tab_clasif <-table(Y, y.pred,dnn = c("observ", "predic"))
porcentaje <- (3+18)/30
porcentaje## [1] 0.7Nota. El modelo presenta una precisión del 70%, lo que indica un modelo bueno relativamente.
- D. Realizar predicciones en base a los siguientes valores: X (18,60 y 50) e interpretar.
#Prediccion
nuevo<-data.frame(idade=c(18, 60, 50))
pred<-predict(modelo_logistico1,newdata=nuevo,type="response")
pred## 1 2 3
## 0.7787066 0.4386261 0.5278626Nota. De acuerdo a estos resultados podemos decir que, si el individuo tiene 18 años, tiene una probabilidad de 0.78 que trabaje, si consideramos la categorización para la bonomial a un 50%, esta probabilidad permitiría inferir que le corresponde una clasificacion (trabaja), mientras que para un individuo de 60 años (no trabaja), y para uno de 50 años (si trabaja).
*E. Encontrar los parámetros posteriores para una regresión logística bayesiana asumiendo que siguen una distribución normal: Intercepto -> N (0,8); B1 -> N (0,9).
Para ello consideramos el siguiente código:
#REGRESION LOGISTICA BAYESIANA ASUMIENDO NORMALIDAD
# Instala y carga el paquete "rstan" si aún no lo has hecho
library(rstan)## Loading required package: StanHeaders##
## rstan version 2.26.23 (Stan version 2.26.1)## For execution on a local, multicore CPU with excess RAM we recommend calling
## options(mc.cores = parallel::detectCores()).
## To avoid recompilation of unchanged Stan programs, we recommend calling
## rstan_options(auto_write = TRUE)
## For within-chain threading using `reduce_sum()` or `map_rect()` Stan functions,
## change `threads_per_chain` option:
## rstan_options(threads_per_chain = 1)## Do not specify '-march=native' in 'LOCAL_CPPFLAGS' or a Makevars filelibrary(StanHeaders)
# Definir el modelo en formato Stan
modelo_stan <- "
data {
int<lower=0> N; // Número de observaciones
int<lower=0, upper=1> y[N]; // Variable de respuesta binaria
vector[N] x1; // Variable predictora x1
}
parameters {
real alpha; // Intercepto
real beta1; // Coeficiente para x1
}
model {
// Prior para los parámetros
alpha ~ normal(0, 8);
beta1 ~ normal(0, 9);
// Modelo logístico
for (i in 1:N) {
y[i] ~ bernoulli_logit(alpha + beta1 * x1[i]);
}
}
"
# Convertir los datos en formato necesario para Stan
datos_stan <- list(
N = nrow(datos),
y = Y,
x1 = idade
)
# Compilar el modelo
modelo_compilado <- stan_model(model_code = modelo_stan)
# Ajustar el modelo utilizando MCMC
ajuste <- sampling(modelo_compilado, data = datos_stan, chains = 4, iter = 2000)##
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 1).
## Chain 1:
## Chain 1: Gradient evaluation took 0.000111 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 1.11 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1:
## Chain 1:
## Chain 1: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 1: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 1: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 1: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 1: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 1: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 1: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 1: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 1: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 1: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 1: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 1: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 1:
## Chain 1: Elapsed Time: 1.6 seconds (Warm-up)
## Chain 1: 0.338 seconds (Sampling)
## Chain 1: 1.938 seconds (Total)
## Chain 1:
##
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 2).
## Chain 2:
## Chain 2: Gradient evaluation took 3.4e-05 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.34 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2:
## Chain 2:
## Chain 2: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 2: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 2: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 2: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 2: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 2: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 2: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 2: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 2: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 2: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 2: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 2: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 2:
## Chain 2: Elapsed Time: 0.523 seconds (Warm-up)
## Chain 2: 0.429 seconds (Sampling)
## Chain 2: 0.952 seconds (Total)
## Chain 2:
##
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 3).
## Chain 3:
## Chain 3: Gradient evaluation took 3.6e-05 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.36 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3:
## Chain 3:
## Chain 3: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 3: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 3: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 3: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 3: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 3: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 3: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 3: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 3: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 3: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 3: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 3: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 3:
## Chain 3: Elapsed Time: 0.498 seconds (Warm-up)
## Chain 3: 0.225 seconds (Sampling)
## Chain 3: 0.723 seconds (Total)
## Chain 3:
##
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 4).
## Chain 4:
## Chain 4: Gradient evaluation took 1.9e-05 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.19 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4:
## Chain 4:
## Chain 4: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 4: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 4: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 4: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 4: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 4: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 4: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 4: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 4: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 4: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 4: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 4: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 4:
## Chain 4: Elapsed Time: 0.49 seconds (Warm-up)
## Chain 4: 0.228 seconds (Sampling)
## Chain 4: 0.718 seconds (Total)
## Chain 4:# Resumen de los resultados
summary(ajuste)## $summary
## mean se_mean sd 2.5% 25% 50%
## alpha 2.04266311 0.054122500 1.43208236 -0.7441941 1.06151743 2.04518265
## beta1 -0.03885569 0.001244524 0.03291062 -0.1029084 -0.06157107 -0.03879121
## lp__ -20.59794533 0.032411825 1.07780191 -23.4898117 -20.96227259 -20.27266632
## 75% 97.5% n_eff Rhat
## alpha 3.03151300 4.91488805 700.1326 1.005009
## beta1 -0.01669773 0.02550375 699.3032 1.005981
## lp__ -19.86229745 -19.57074663 1105.7856 1.002652
##
## $c_summary
## , , chains = chain:1
##
## stats
## parameter mean sd 2.5% 25% 50%
## alpha 2.14631360 1.45227188 -0.9759760 1.23869176 2.15299408
## beta1 -0.04120048 0.03344002 -0.1045095 -0.06337669 -0.04205268
## lp__ -20.64623270 1.12886560 -23.7217496 -21.05680151 -20.28280073
## stats
## parameter 75% 97.5%
## alpha 3.16960478 4.94284632
## beta1 -0.02128332 0.03148293
## lp__ -19.84960997 -19.57694901
##
## , , chains = chain:2
##
## stats
## parameter mean sd 2.5% 25% 50%
## alpha 2.05153307 1.53199984 -0.7273070 1.03553681 1.98359360
## beta1 -0.03925165 0.03478055 -0.1142846 -0.06246942 -0.03771441
## lp__ -20.62997270 1.12550028 -23.4851643 -21.01443096 -20.29644704
## stats
## parameter 75% 97.5%
## alpha 2.96329849 5.29039034
## beta1 -0.01508468 0.02349584
## lp__ -19.85562703 -19.56815358
##
## , , chains = chain:3
##
## stats
## parameter mean sd 2.5% 25% 50%
## alpha 2.05563540 1.35907830 -0.54818450 1.05987828 2.0829133
## beta1 -0.03929276 0.03136479 -0.09840946 -0.06237412 -0.0398435
## lp__ -20.57050462 0.98961190 -23.23289583 -20.93759098 -20.2784275
## stats
## parameter 75% 97.5%
## alpha 3.0486004 4.4921119
## beta1 -0.0165514 0.0197115
## lp__ -19.8858400 -19.5730098
##
## , , chains = chain:4
##
## stats
## parameter mean sd 2.5% 25% 50%
## alpha 1.91717037 1.37099298 -0.76179476 0.95145217 1.95062275
## beta1 -0.03567788 0.03174991 -0.09601339 -0.05697264 -0.03626512
## lp__ -20.54507132 1.05961462 -23.37953891 -20.89927642 -20.22639217
## stats
## parameter 75% 97.5%
## alpha 2.85225745 4.54360243
## beta1 -0.01382523 0.02700674
## lp__ -19.85168078 -19.56531782Nota. Observamos que, el modelo bayesiano, quedaría con los parametros alfa = 2.07 y beta = -0.0396
# Visualización de resultados (por ejemplo, trazar distribuciones posteriores)
plot(ajuste)## ci_level: 0.8 (80% intervals)## outer_level: 0.95 (95% intervals)
Nota. Se observa que el modelo es adecuado.
ASUMIENDO UNIFORMIDAD
CONSIDERANDO LOS DATOS ANTERIORES
datos<-read.table("trabajo.txt",header = T,sep = "\t")NOS PIDEN
- Utilizar la base de datos TRABAJO con la variable dependiente Y= Condición de trabajo y la variable predictora X= años de escolaridad de la persona.
- A. Estimar los parámetros con la ecuación de regresión logística binaria clásica.
# Ajustar el modelo de regresión logística binaria
modelo_logistico2 <- glm(Y ~ escolaridade, data = datos, family = binomial)
# Resumen del modelo
summary(modelo_logistico2)##
## Call:
## glm(formula = Y ~ escolaridade, family = binomial, data = datos)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -5.9038 2.7624 -2.137 0.0326 *
## escolaridade 0.5429 0.2380 2.281 0.0225 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 40.381 on 29 degrees of freedom
## Residual deviance: 33.152 on 28 degrees of freedom
## AIC: 37.152
##
## Number of Fisher Scoring iterations: 4Nota. Se observa que, la variable idade es significativa (p-value = 0.0225 < 5%), lo que da a entender que influye o repercute en la condición de trabajo.
- B. Determinar la bondad de ajuste del modelo.
#Test de Hosmer-Lemeshow
library(ResourceSelection)
hoslem.test(Y, fitted(modelo_logistico2), g=10)## Warning in hoslem.test(Y, fitted(modelo_logistico2), g = 10): The data did not
## allow for the requested number of bins.##
## Hosmer and Lemeshow goodness of fit (GOF) test
##
## data: Y, fitted(modelo_logistico2)
## X-squared = NaN, df = 5, p-value = NANota. Según la prueba de Hosmer-Lemeshow, el modelo es adecuado.
#Pseudo-R2
library(DescTools)
PseudoR2(modelo_logistico2,which = c("McFadden","CoxSnell","Nagelkerke"))## McFadden CoxSnell Nagelkerke
## 0.1790045 0.2141146 0.2894509Nota. observamos que el valor de Rcuadrado de Cox y Snell es igual a 0.214 y el valor de Rcuadrado de Nagelkerke es igual a 0.2895, el cual nos indica que el 28.95% de las proporciones de la varianza de la variable dependiente (Condición de trabajo) es explicada por la variable predictora escolaridad; dando a entender que el modelo de regresión logística no es bueno.
- C. Hallar el coeficiente OR e interpretar.
exp(cbind(OR = coef(modelo_logistico2), confint(modelo_logistico2, level=0.95)))## Waiting for profiling to be done...## OR 2.5 % 97.5 %
## (Intercept) 0.002729138 4.795682e-06 0.3330419
## escolaridade 1.720977337 1.142515e+00 2.9839981Nota. Nos indica la posibilidad de que un individuo trabaje, por su año de escolaridad, es decir, es aproximadamente 1.72 veces más probable que el individuo trabaje por tener mas años de escolaridad, respecto a los que tienen menos años de escolaridad.
#Precisión del modelo
y.pred<-ifelse(modelo_logistico2$fitted.values > 0.43, yes = 1, no = 0)
tab_clasif <-table(Y, y.pred,dnn = c("observ", "predic"))
porcentaje <- (5+14)/30
porcentaje## [1] 0.6333333Nota. El modelo presenta una precisión del 63.3%, lo que indica un modelo regular relativamente.
- D. Realizar predicciones en base a los siguientes valores: X (5, 20 y 22) e interpretar.
#Prediccion
nuevo<-data.frame(idade=c(5, 20, 22))
pred<-predict(modelo_logistico2,newdata=nuevo,type="response")## Warning: 'newdata' had 3 rows but variables found have 30 rowspred## 1 2 3 4 5 6 7 8
## 0.9417232 0.8451063 0.9417232 0.2654682 0.6481549 0.6481549 0.8451063 0.3834699
## 9 10 11 12 13 14 15 16
## 0.6481549 0.1735562 0.5170050 0.8451063 0.6481549 0.7602100 0.2654682 0.3834699
## 17 18 19 20 21 22 23 24
## 0.8451063 0.3834699 0.6481549 0.7602100 0.8451063 0.6481549 0.1087549 0.5170050
## 25 26 27 28 29 30
## 0.6481549 0.3834699 0.9037510 0.3834699 0.5170050 0.6481549Nota. De acuerdo a estos resultados podemos decir que, si el individuo tiene 5 años de escolaridad, tiene una probabilidad de 0.0396 que trabaje, si consideramos la categorización para la bonomial a un 50%, esta probabilidad permitiría inferir que le corresponde una clasificacion (no trabaja), mientras que para un individuo que tiene 20 años de escolaridad (si trabaja), y para uno que tiene 22 años de escolaridad (si trabaja).
*E. Encontrar los parámetros posteriores para una regresión logística bayesiana asumiendo una distribución a priori uniforme: burnin=1000, mcmc=25000 .
Para ello consideramos el siguiente código:
#REGRESION LOGISTICA BAYESIANA ASUMIENDO UNIFORMIDAD
library("MASS"); library(lattice);library("coda");library("MCMCpack")##
## Attaching package: 'coda'## The following object is masked from 'package:rstan':
##
## traceplot## ##
## ## Markov Chain Monte Carlo Package (MCMCpack)## ## Copyright (C) 2003-2023 Andrew D. Martin, Kevin M. Quinn, and Jong Hee Park## ##
## ## Support provided by the U.S. National Science Foundation## ## (Grants SES-0350646 and SES-0350613)
## ##library("splines"); library("survival"); library("leaps")
##Codigo para la simulaci?n, suponiendo que tenemos prior uniforme
posterior<-MCMClogit(Y~escolaridade, data=datos,burnin=1000,mcmc=25000)
summary(posterior)##
## Iterations = 1001:26000
## Thinning interval = 1
## Number of chains = 1
## Sample size per chain = 25000
##
## 1. Empirical mean and standard deviation for each variable,
## plus standard error of the mean:
##
## Mean SD Naive SE Time-series SE
## (Intercept) -6.7876 2.969 0.018778 0.059643
## escolaridade 0.6229 0.256 0.001619 0.005144
##
## 2. Quantiles for each variable:
##
## 2.5% 25% 50% 75% 97.5%
## (Intercept) -13.2453 -8.6175 -6.5299 -4.705 -1.604
## escolaridade 0.1737 0.4427 0.5999 0.780 1.181Nota. Observamos que, el modelo bayesiano, quedaría con los parametros alfa = -6.7876 y beta = 0.6229
# Visualización de resultados Trace
plot(posterior,trace=FALSE)
Nota. Se observa que el modelo es adecuado.
# Visualización de resultados Densidad
plot(posterior, density=FALSE)
Nota. Se observa que la densidad de los datos.
# Visualización de resultados autocorrelaciones
autocorr.plot(posterior)
Nota. Se observa las autocorrelaciones moderadas.
A MODO PRÁCTICA