Regresión logística bayesiana

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           12

NOS 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: 4

Nota. 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-27

hoslem.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.01593

Nota. 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.05917289

Nota. 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.024774

Nota. 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.7

Nota. 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.5278626

Nota. 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 file

library(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 3.1e-05 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.31 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: 0.231 seconds (Warm-up)
## Chain 1:                0.114 seconds (Sampling)
## Chain 1:                0.345 seconds (Total)
## Chain 1: 
## 
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 2).
## Chain 2: 
## Chain 2: Gradient evaluation took 1.1e-05 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.11 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.22 seconds (Warm-up)
## Chain 2:                0.111 seconds (Sampling)
## Chain 2:                0.331 seconds (Total)
## Chain 2: 
## 
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 3).
## Chain 3: 
## Chain 3: Gradient evaluation took 1.2e-05 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.12 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.267 seconds (Warm-up)
## Chain 3:                0.113 seconds (Sampling)
## Chain 3:                0.38 seconds (Total)
## Chain 3: 
## 
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 4).
## Chain 4: 
## Chain 4: Gradient evaluation took 1.1e-05 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.11 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.289 seconds (Warm-up)
## Chain 4:                0.108 seconds (Sampling)
## Chain 4:                0.397 seconds (Total)
## Chain 4:

# Resumen de los resultados
summary(ajuste)

## $summary
##               mean     se_mean        sd        2.5%          25%          50%
## alpha   1.94519751 0.048098555 1.4249642  -0.7647737   0.99619400   1.94748044
## beta1  -0.03628245 0.001104015 0.0328699  -0.1049694  -0.05709114  -0.03605348
## lp__  -20.60457146 0.035251438 1.0711631 -23.5719868 -20.98951172 -20.25911884
##                75%        97.5%    n_eff     Rhat
## alpha   2.83973954   5.02417426 877.6955 1.002323
## beta1  -0.01531548   0.02740069 886.4354 1.001657
## lp__  -19.85695788 -19.57674122 923.3313 1.003887
## 
## $c_summary
## , , chains = chain:1
## 
##          stats
## parameter         mean         sd        2.5%          25%          50%
##     alpha   2.02748840 1.49048336  -0.6725260   1.08582884   1.94422865
##     beta1  -0.03790701 0.03413455  -0.1179062  -0.05634573  -0.03553694
##     lp__  -20.61975876 1.15135823 -23.9100236 -21.01593986 -20.19828395
##          stats
## parameter          75%        97.5%
##     alpha   2.78874712   5.49915534
##     beta1  -0.01829783   0.02595669
##     lp__  -19.80094941 -19.56877506
## 
## , , chains = chain:2
## 
##          stats
## parameter         mean        sd        2.5%          25%          50%
##     alpha   1.87931764 1.4490929  -0.9320417   0.89135207   1.90998401
##     beta1  -0.03494794 0.0336431  -0.1032937  -0.05748087  -0.03503116
##     lp__  -20.65696966 1.0981207 -23.7506596 -21.06413173 -20.31942431
##          stats
## parameter          75%        97.5%
##     alpha   2.86196782   4.81827891
##     beta1  -0.01218442   0.03059206
##     lp__  -19.87223068 -19.58448170
## 
## , , chains = chain:3
## 
##          stats
## parameter         mean         sd        2.5%         25%          50%
##     alpha   1.94487181 1.34625860  -0.5255522   0.9772843   1.95742929
##     beta1  -0.03638578 0.03129487  -0.0999430  -0.0588545  -0.03674804
##     lp__  -20.59371251 0.93662440 -23.0827274 -21.0083788 -20.30666669
##          stats
## parameter          75%        97.5%
##     alpha   2.90071264   4.69999871
##     beta1  -0.01444731   0.02279752
##     lp__  -19.91162673 -19.58138087
## 
## , , chains = chain:4
## 
##          stats
## parameter         mean        sd        2.5%        25%          50%
##     alpha   1.92911217 1.4081848  -0.8041675   1.024572   1.96614259
##     beta1  -0.03588909 0.0323107  -0.1008598  -0.055792  -0.03708014
##     lp__  -20.54784489 1.0853243 -23.5621502 -20.875716 -20.19871855
##          stats
## parameter          75%        97.5%
##     alpha   2.78338973   4.83427174
##     beta1  -0.01523633   0.03022772
##     lp__  -19.84106022 -19.56460927

Nota. 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: 4

Nota. 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 = NA

Nota. 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.2894509

Nota. 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.9839981

Nota. 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.6333333

Nota. 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 rows

pred

##         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.6481549

Nota. 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.181

Nota. 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.

Regresión logística bayesiana

Yesica CASTILLO oLIVARES

2023-09-23

ASUMIENDO NORMALIDAD

CONSIDERANDO LOS SIGUIENTES DATOS

NOS PIDEN

ASUMIENDO UNIFORMIDAD

CONSIDERANDO LOS DATOS ANTERIORES

NOS PIDEN