Simulacion Modelo Skew Probit, Logit y Probit bajo un enfoque Bayesiano
Manuel Rain manuelrain@hotmail.com
7 de agosto de 2019
Introducción
Los modelos lineales generalizados permiten el tratamiento de problemas de regresión, en los que se puede elegir la distribucion de la variable respuesta (McCullagh, 1984). El enlace Logit (Bliss, 1935) y Probit (Berkson, 1944), (Nagler, 1994), (Bazan et al., 2006), (Chen et al., 1999) y (Chen, 2004), entre otros. mostraron la importancia de elegir apropiadamente la funcion de enlace y que tan sensible es la inferencia si se utiliza incorrectamente una funcion de enlace asimetrico en lugar de un enlace asimetrico.
En forma general este estudio incorpora un parámetro que controla la asimetría de la función de enlace y se basa en la Skew Normal propuesta por Azzalini (1985). En consecuencia, el metodo se basa en un enfoque de variables latentes de Albert and Chib (1993).
Por otro lado, mediante la introducción de variables latentes, los algoritmos Monte Carlo de cadena de Markov (MCMC) pueden ser facilmente implementados para muestrear eficientemente la distribución a posteriori de los parametros.
Objetivos
Objetivo General
Implementar y estimar el modelo Skew Probit bajo un enfoque Bayesiano.
Objetivos Especificos
Estudiar la distribución Skew-Normal.
Construir el modelo Skew-Probit.
Implementar computacionalmente los modelos Skew Probit, Logit y Probit bajo una perspectiva Bayesiana en el software R Studio (2012).
Comparar los modelos Skew Probit, Logit y Probit.
Distribución Skew-Normal
Esta fue introducida por (Azzalini, 1985), la familia incluye la distribución normal como un caso particular y permite modelar la asimetría presente en un conjunto de observaciones.
Una variable aleatoria \(Z\) tiene una distribución Skew-Normal si su función de densidad de probabilidad (fdp en adelante) está dada por
\[ f_z(z)=2 \phi(z)\Phi(\lambda z), \quad z \in \mathbb{R} \quad \lambda \in \mathbb{R}, \]
donde \(\phi\) y \(\Phi\) son la fdp y la función de distribución acumulada (fda en adelante) de una normal estándar, respectivamente, \(\lambda\) es el parámetro de asimetria. Consideremos la notación para denotar una variable aleatoria \(Z\) con distribución Skew Normal estándar \(Z \sim \mathcal{SN} (0, 1,\lambda)\).
Una variable más flexible se puede construir añadiendo parámetros de localización y escala, como \(Y = \xi + \tau Z\). Entonces, \(Y\) tiene un fdp dada por
\[f_Y(y)=\frac{1}{\tau} \phi \left(\frac{y-\xi}{\tau} \right) \Phi \left(\lambda \frac{y-\xi}{\tau}\right), \quad \xi \in \mathbb{R}, \tau >0 \]
Consideramos la notación para denotar una variable aleatoria \(Y\) con distribución Skew Normal \(Y \sim \mathcal{SN} (\xi, \tau^2,\lambda)\), llamado el modelo de tres parámetros.
Si \(Z\sim \mathcal{SN}(0,1,\lambda)\), entonces las siguientes propiedades son válidas:
\(\mathcal{SN}(0,1,0)= \mathcal{N}(0,1)\);
\(\mathcal{SN}(0,1,\infty)= \mathcal{HN}(0,1)\);
\(|Z| \sim \mathcal{HN}(0,1)\);
\(Z^{2} \sim \chi^2_{1}\)
El siguiente resultado ha sido presentado repetidamente en la literatura, con una prueba que incluyen algunos autores (Ellison, 1964) y (Zacks et al., 1981).
Sea \(\mathcal{U} \sim \mathcal{N}(0,1)\) entonces
\[\mathbb{E}\left\lbrace\Phi(h \mathcal{U} + k)\right\rbrace=\Phi\left(\frac{k}{\sqrt{1+h^2}}\right), \quad h,k \in \mathbb{R} \]
A partir de este resultado, la función generadora de momentos de \(Y\) se obtiene, es decir
\[ M_{Y}(t)= \mathbb{E} exp(tY)\] \[ M_{Y}(t)= \mathbb{E} exp(t(\xi+\tau Z) \] \[ M_{Y}(t)= \mathbb{E} exp(\xi t+\tau Zt) \] \[ M_{Y}(t)=2 \exp(\xi t+\frac{1}{2}\tau^2t^2)\int_{\mathbb{R}}\phi(z-\tau t)\Phi(\lambda z)dz \] \[ M_{Y}(t)=2 \exp(\xi t+\frac{1}{2}\tau^2t^2) \Phi(\delta \tau t) \] donde \(\delta=\delta(\lambda)=\frac{\lambda}{\sqrt{1+\lambda^2}}, \quad \lambda \in \mathbb{R} \Rightarrow \delta \in (-1,1)\)
Simulaciones Skew-Normal
library("sn")## Warning: package 'sn' was built under R version 3.5.3## Loading required package: stats4##
## Attaching package: 'sn'## The following object is masked from 'package:stats':
##
## sdx<-seq(-5,5,length=1000)
par(mfrow=c(1,2))
plot(x,dsn(x,0,1,1/2),ylab =" ",type = "l", lwd = 2, col="red", main="", ylim=c(0,1.3))
lines(x,dsn(x,0,1,2),type = "l", lwd = 2, col="blue")
lines(x,dsn(x,0,1,14),type = "l", lwd = 2, col="orange")
legend("topleft", inset=c(0,0), c( expression(lambda==1/2),expression(lambda==2),expression(lambda==14)),bty="n",fill=c("red","blue","orange"))
plot(x,dsn(x,0,1,-1/2),ylab =" ",type = "l", lwd = 2, col="green", main="", ylim=c(0,1.3))
lines(x,dsn(x,0,1,-2),type = "l", lwd = 2, col="purple")
lines(x,dsn(x,0,1,-14),type = "l", lwd = 2, col="brown")
legend("topleft", inset=c(0,0), c( expression(lambda==-1/2),expression(lambda==-2),expression(lambda==-14)),bty="n",fill=c("green","purple","brown"))
Para valores positivos de \(\lambda\), la \(\mathcal{SN}(\xi=0,\tau=1,\lambda)\) es sesgada a la derecha, caso contrario es sesgada a la izquierda.
Modelo
Sea \(\mathbf{y}=(y_{1},y_{2},...,y_{n})^{'}\) un vector \(n \times 1\) de \(n\) variables aleatorias independientes dicotómicas. Además, \(\mathbf{x}_{i}=(x_{i1},...,x_{ik})^{'}\) es un vector \(k \times 1\) de covariables, donde \(x_{i1}\) puede ser \(1\) (que corresponda a un intercepto), \(\mathrm{X}\) denotamos la matriz \(n\times k\) con filas \(\mathbf{x}_{i}^{'}\) , y $ = ( {1},…,{k}) ^{’}$ es un vector de \(k\times1\) coeficientes de regresión.
Supongamos que \(y_{i}=1\) con probabilidad \(p_{i}\) y \(y_{i}=0\) con \(1-p_{i}\).
Es un modelo de respuesta dicotómica. En el caso general, asumimos que
\[ \begin{align} p_{i}= F(\mathbf{x}_{i}^{'} \boldsymbol{ \beta}), \end{align} \]
donde \(F(\cdot)\) es una función de distribución acumulada \((cdf)\) y \(F^{-1}\) es llamada función de enlace en un modelo lineal generalizado.
Sea \(\mathbf{w}=(w_{1}.w_{2},...,w_{n})'\) un vector de variables latentes independiente. Entonces nuestro modelo de enlace propuesto por se formula como
\[ \begin{align} y_{i}= \left \{ \begin{matrix} \quad 0 \quad , si & w_{i} <0, \\ \\\quad 1 \quad, si & w_{i} \geq 0, \end{matrix}\right. \end{align}\]
donde \[w_{i}=\mathbf{x} _{i}^{'}\boldsymbol{\beta} + \delta z_{i}+ \varepsilon_{i} \] \(z_{i} \sim \mathcal{HN}(0,1)\) y \(\varepsilon_{i} \sim \mathcal{N}(0,1)\) con función de densidad \(g(z)= \frac{2}{ \sqrt{2\pi}}e^{-(\frac{z^{2}}{2})}, \quad z>0\).
\[\begin{align} p_{i}= \int_{0}^{\infty} \Phi (\mathbf{x} _{i}^{'}\boldsymbol{\beta}+\delta z_{i})g(z_{i})dz_{i} \end{align}\]
La Función de Verosimilitud
Sea \(D_{obs}=(n,\mathbf{y},\mathbf{X})\) denota los datos observados. Entonces, la función de verosimilitud para el modelo de enlace asimétrico está dado por,
\[ \begin{align} L(\boldsymbol{\beta},\delta|D_{obs})= \prod_{i=1}^{n} \int _{-\infty}^{\infty} [F(\mathbf{x} _{i}^{'}\boldsymbol{\beta}+ \delta z_{i})]^{y_{i}}\times [1-F(\mathbf{x} _{i}^{'}\boldsymbol{\beta} + \delta z_{i})]^{1-y_{i}}g(z_{i})dz_{i}, \end{align}\]
Sea \(\mathbf{z}=(z_{1},...,z_{n})^{'}\), y por \(D=(n,\mathbf{y},\mathbf{X},\mathbf{z})\) denotamos los datos completos. Entonces, la función de verosimilitud de los datos completos para los parámetros \((\beta, \delta)\) se pueden escribir como
\[\begin{align} L(\beta, \delta | D)= \prod_{i=1}^{n}[F(x_{i}^{'}\beta +\delta z_{i})]^{y_{i}}\times [1-F(x_{i}^{'}\beta +\delta z_{i})]^{1-y_{i}}g(z_{i}), \end{align} \] es equivalente a la representación de un modelo de regresión binomial estándar con \(\theta=(\beta. \delta)\) un vector de coeficiente de regresión \((x_{i}^{'},z_{i})^{'}\) y vector de covariables.
Distribución a priori
La distribuciones de a priori no informativas para \((\boldsymbol{\beta},\delta)\).
Supongamos que consideramos una a priori impropia para \((\boldsymbol{\beta},\delta)\) de la forma
\[ \begin{align} \pi(\boldsymbol{\beta},\delta)\propto \pi(\delta), \end{align} \]
donde \(\pi(\delta)\) es una a priori, es decir, \(\pi(\delta)\propto 1\). La distribución a posteriori de \((\boldsymbol{\beta},\delta)\) se basa en la observación de los datos \(D_{obs}\) y es dada por,
\[\begin{align*} \begin{split} p (\boldsymbol{\beta},\delta | D_{obs}) & \propto L( \boldsymbol{\beta},\delta| D_{obs})\pi(\delta) \\ & =\left\lbrace \prod_{i=1}^{n} \int_{-\infty}^{\infty} [F(\mathbf{x} _{i}^{'}\boldsymbol{\beta}+\delta z_{i})]^{y_{i}}\times [1- F(\mathbf{x} _{i}^{'}\boldsymbol{\beta}+\delta z_{i})]^{1-y_{i}}g(z_{i})dz_{i} \right\rbrace \pi (\delta). \end{split} \end{align*}\]
A priori no informativa: es relativa a la función de verosimilitud, tiene un impacto mínimo a la distribución a posteriori.
MCMC para inferencia Bayesiana
Para facilitar la explicación, consideramos el caso para el cual se usa una a priori no informativa. Para la distribución a posterior \(p(\boldsymbol{\beta},\delta|D_{obs})\), se introducen las variables latentes \(\mathbf{z}=(z_{1},...,z_{n})\). Entonces la distribución a posteriori conjunta de \((\boldsymbol{\beta},\delta,\mathbf{z})\) está dada por
\[ \begin{align} p(\boldsymbol{\beta},\delta,\mathbf{z} |D_{obs})\propto \left\lbrace \prod_{i=1}^{n} [F(\mathbf{x} _{i}^{'}\boldsymbol{\beta}+ \delta z_{i})]^{y_{i}}\times [1-F( \mathbf{x} _{i}^{'}\boldsymbol{\beta}+ \delta z_{i})]^{1-y_{i}}g(z_{i})\right\rbrace \pi (\delta). \end{align}\] Para ejecutar muestreo de Gibss, necesitamos muestrear las siguientes distribuciones condicionales completas de: \([z_{i}|\boldsymbol{\beta},\delta,D_{obs}]\) por \(i=1,2,...,n\), \([\boldsymbol{\beta}|\delta,z,D_{obs}]\) y \([\delta|\boldsymbol{\beta}.z,D_{obs}]\).
Simulación
Se realiza un estudio de simulación para evaluar el desempeño relativo del procedimiento de estimación en términos de recuperación de parámetros del modelo considerando: dos tamaños de muestra \(n = (200)\) y el uso de media o mediana a posterior como medida de resumen para \(\lambda\).
De modo que para cada uno, se obtienen \(n\) variables de respuesta de Bernoulli independientes \(y_i\) para el modelo Skew Probit, con una covariable \(x_i\) y el valor de \(\beta_1= 1\), \(\beta_2= 2\) y \(\lambda = 0,5\).
Para cada conjunto de datos generado, nos adaptamos al modelo verdadero correspondiente. Para cada muestra simulada. Se quemó \(5000\) en los \(100.000\) valores de cadena.
Estimaciones de las simulaciones
Modelo Skew Probit
library("BRugs")## Warning: package 'BRugs' was built under R version 3.5.3## Welcome to BRugs connected to OpenBUGS version 3.2.3library("boot")
library("coda")
modelFile <- "C:/Users/Manuel/Desktop/200/campeon.txt"
bodega <- read.table("C:/Users/Manuel/Desktop/200/datos.txt",
header=TRUE, sep="", na.strings="NA",
dec=".", strip.white=TRUE)
data <- c(as.list(bodega), list(n=200,k=2))
inits<-function(){list(beta=c(0.0,0.0),delta=0.5)}
parametersToSave <- c("beta","delta")
BRugsFit(modelFile, data, inits, numChains=3,
parametersToSave, nBurnin=5000, nIter=
100000, nThin=50, coda=FALSE, DIC=TRUE,
working.directory=NULL, digits=3,
seed=NULL, BRugsVerbose=getOption("BRugsVerbose"))## model is syntactically correct## data loaded## model compiled## [1] "C:\\Users\\Manuel\\AppData\\Local\\Temp\\Rtmp0EQ6i4/inits1.txt"
## [2] "C:\\Users\\Manuel\\AppData\\Local\\Temp\\Rtmp0EQ6i4/inits2.txt"
## [3] "C:\\Users\\Manuel\\AppData\\Local\\Temp\\Rtmp0EQ6i4/inits3.txt"## Initializing chain 1:## initial values loaded but chain contain uninitialized variables## Initializing chain 2:## initial values loaded but chain contain uninitialized variables## Initializing chain 3:## initial values loaded but chain contain uninitialized variables## initial values generated, model initialized## 5000 updates took 8 s## deviance set## monitor set for variable 'beta'## monitor set for variable 'delta'## 100000 updates took 253 s## $Stats
## mean sd MC_error val2.5pc median val97.5pc start sample
## beta[1] 1.00800 0.5748 0.01846 -0.00869 0.9840 2.1070 5001 6000
## beta[2] 2.07700 0.4031 0.01860 1.37900 2.0460 2.9600 5001 6000
## delta 0.09549 0.5939 0.01583 -0.94420 0.1411 0.9675 5001 6000
##
## $DIC
## Dbar Dhat DIC pD
## er 277.30 277.30 277.30 7.907e-10
## z 51.26 52.61 49.91 -1.349e+00
## total 328.50 329.90 327.20 -1.349e+00mychain <- buildMCMC('*')
summary(mychain)##
## Iterations = 5001:104951
## Thinning interval = 50
## Number of chains = 3
## Sample size per chain = 2000
##
## 1. Empirical mean and standard deviation for each variable,
## plus standard error of the mean:
##
## Mean SD Naive SE Time-series SE
## beta[1] 1.00822 0.5749 0.007421 0.01842
## beta[2] 2.07717 0.4031 0.005204 0.01929
## delta 0.09549 0.5940 0.007668 0.01648
##
## 2. Quantiles for each variable:
##
## 2.5% 25% 50% 75% 97.5%
## beta[1] -0.008704 0.5554 0.984 1.4456 2.106
## beta[2] 1.378688 1.7974 2.046 2.3151 2.959
## delta -0.944260 -0.4156 0.141 0.6354 0.967plot(mychain, mfrow = c(1, 3))
Modelo Logit
library("BRugs")
library("boot")
library("coda")
modelFile <- "C:/Users/Manuel/Desktop/200/logit.txt"
bodega <- read.table("C:/Users/Manuel/Desktop/200/datos.txt",
header=TRUE, sep="", na.strings="NA",
dec=".", strip.white=TRUE)
data <- c(as.list(bodega), list(n=200,k=2))
inits<-function(){list(beta=c(0,0))}
parametersToSave <- "beta"
BRugsFit(modelFile, data, inits, numChains=3, parametersToSave, nBurnin=5000, nIter=100000, nThin=50, coda=FALSE, DIC=TRUE, working.directory=NULL, digits=3,
seed=NULL, BRugsVerbose=getOption("BRugsVerbose"))## model is syntactically correct## data loaded## model compiled## [1] "C:\\Users\\Manuel\\AppData\\Local\\Temp\\Rtmp0EQ6i4/inits1.txt"
## [2] "C:\\Users\\Manuel\\AppData\\Local\\Temp\\Rtmp0EQ6i4/inits2.txt"
## [3] "C:\\Users\\Manuel\\AppData\\Local\\Temp\\Rtmp0EQ6i4/inits3.txt"## Initializing chain 1:## initial values loaded and chain initialized but another chain contain uninitialized variables## Initializing chain 2:## initial values loaded and chain initialized but another chain contain uninitialized variables## Initializing chain 3:## model is initialized## model is already initialized## 5000 updates took 2 s## deviance set## monitor set for variable 'beta'## 100000 updates took 52 s## $Stats
## mean sd MC_error val2.5pc median val97.5pc start sample
## beta[1] 1.87 0.5793 0.008326 0.8724 1.814 3.123 5001 6000
## beta[2] 3.89 0.8543 0.011260 2.4490 3.815 5.805 5001 6000
##
## $DIC
## Dbar Dhat DIC pD
## y 53.77 51.84 55.69 1.925
## total 53.77 51.84 55.69 1.925mychain <- buildMCMC('*')
summary(mychain)##
## Iterations = 5001:104951
## Thinning interval = 50
## Number of chains = 3
## Sample size per chain = 2000
##
## 1. Empirical mean and standard deviation for each variable,
## plus standard error of the mean:
##
## Mean SD Naive SE Time-series SE
## beta[1] 1.87 0.5793 0.007479 0.007685
## beta[2] 3.89 0.8544 0.011030 0.011030
##
## 2. Quantiles for each variable:
##
## 2.5% 25% 50% 75% 97.5%
## beta[1] 0.8724 1.459 1.814 2.243 3.119
## beta[2] 2.4489 3.278 3.815 4.421 5.802plot(mychain)
Probit
library("BRugs")
library("boot")
library("coda")
modelFile <- "C:/Users/Manuel/Desktop/200/probit.txt"
bodega <- read.table("C:/Users/Manuel/Desktop/200/datos.txt",
header=TRUE, sep="", na.strings="NA",
dec=".", strip.white=TRUE)
data <- c(as.list(bodega), list(n=200,k=2))
inits<-function(){list(beta=c(0,0))}
parametersToSave <- "beta"
BRugsFit(modelFile, data, inits, numChains=3, parametersToSave, nBurnin=5000, nIter=
100000, nThin=50, coda=FALSE, DIC=TRUE, working.directory=NULL, digits=3,
seed=NULL, BRugsVerbose=getOption("BRugsVerbose"))## model is syntactically correct## data loaded## model compiled## [1] "C:\\Users\\Manuel\\AppData\\Local\\Temp\\Rtmp0EQ6i4/inits1.txt"
## [2] "C:\\Users\\Manuel\\AppData\\Local\\Temp\\Rtmp0EQ6i4/inits2.txt"
## [3] "C:\\Users\\Manuel\\AppData\\Local\\Temp\\Rtmp0EQ6i4/inits3.txt"## Initializing chain 1:## initial values loaded and chain initialized but another chain contain uninitialized variables## Initializing chain 2:## initial values loaded and chain initialized but another chain contain uninitialized variables## Initializing chain 3:## model is initialized## model is already initialized## 5000 updates took 2 s## deviance set## monitor set for variable 'beta'## 100000 updates took 59 s## $Stats
## mean sd MC_error val2.5pc median val97.5pc start sample
## beta[1] 0.882 0.2570 0.003125 0.3987 0.8729 1.423 5001 6000
## beta[2] 1.966 0.3594 0.004098 1.3250 1.9370 2.727 5001 6000
##
## $DIC
## Dbar Dhat DIC pD
## y 54.54 52.6 56.47 1.934
## total 54.54 52.6 56.47 1.934mychain <- buildMCMC('*')
summary(mychain)##
## Iterations = 5001:104951
## Thinning interval = 50
## Number of chains = 3
## Sample size per chain = 2000
##
## 1. Empirical mean and standard deviation for each variable,
## plus standard error of the mean:
##
## Mean SD Naive SE Time-series SE
## beta[1] 0.882 0.2571 0.003319 0.003233
## beta[2] 1.966 0.3594 0.004640 0.004654
##
## 2. Quantiles for each variable:
##
## 2.5% 25% 50% 75% 97.5%
## beta[1] 0.3986 0.7043 0.8728 1.043 1.422
## beta[2] 1.3247 1.7169 1.9374 2.194 2.724plot(mychain)
Conclusión
Esta clase de modelo introduce un parámetro asociado con la asimetría que controla la tasa de aumento (o disminución) de la probabilidad de éxito (o fracaso) de la respuesta binaria.
Debo destacar, que la probabilidad de éxito se obtiene considerando la función de distribución acumulada del predictor lineal y el parámetro de asimetría esta asociado con la distribución elegida y es independiente del predictor lineal.
La especificación de la distribución a priori en la regresión binaria.
En el primer ejemplo, los datos se generaron a partir del modelo Skew Probit concluimos igual que Prentice (1976),Czado and Santner (1992), que el modelo con enlace asimétrico es más conveniente y tienen mejor desempeño que los enlaces Logit y Probit.
Como lo demuestra (Chen, 2004), la elección de los enlaces es importante ya que un enlace mal especificado puede conducir a una predicción deficiente.
Referencias
Albert, J. H. and Chib, S. (1993). Bayesian analysis of binary and polychotomous response data. Journal of the American statistical Association, 88(422):669-679.
Azzalini, A. (1985). A class of distributions which includes the normal ones. Scandinavian journal of statistics, pages 171-178.
Bazán, J. L., Branco, M. D., Bolfarine, H., et al. (2006). A skew item response model. Bayesian analysis, 1(4):861-892.
Berkson, J. (1944). Application of the logistic function to bio-assay. Journal of the American Statistical Association, 39(227):357-365.
Bliss, C. I. (1935). The calculation of the dosage-mortality curve. Annals of Applied Biology, 22(1):134-167.
Buxton, O. M. and Marcelli, E. (2010). Short and long sleep are positively associated with obesity, diabetes, hypertension, and cardiovascular disease among adults in the united states. Social science & medicine, 71(5):1027-1036.
Chen, M. (2004). Skewed link models for categorical response data. Skew-elliptical distributions and their applications: a journey beyond normality, pages 131-151.
Chen, M.-H. and Dey, D. K. (1998). Bayesian modeling of correlated binary responses via scale mixture of multivariate normal link functions. Sankhyiaa: The Indian Journal of Statistics, Series A, pages 322-343.
McCullagh, P. (1984). Generalized linear models. European Journal of Operational Research, 16(3):285-292.
Ogden, C. L., Carroll, M. D., Curtin, L. R., McDowell, M. A., Tabak, C. J., and Flegal, K. M. (2006). Prevalence of overweight and obesity in the united states, 1999-2004. Jama, 295(13):1549-1555.
Studio, R. (2012). Rstudio: integrated development environment for r. RStudio Inc, Boston, Massachusetts.
Stukel, T. A. (1988). Generalized logistic models. Journal of the American Statistical Association, 83(402):426-431.
Zacks, S., Lakshmikantham, V., and Tsokos, C. (1981). Parametric statistical inference.