4.1. Introducción

Los modelos estadísticos comunes permiten medir el efecto de un conjunto de covariables X₁,…,X_p en una variable respuesta Y. Un supuesto fundamental en estos modelos es que la variable respuesta y las covariables se miden sin error. Sin embargo, esto no siempre es posible, especial-mente cuando la variable respuesta mide atributos sensibles.

Una variable sensible es la medición de un atributo sensible que se refiere a información personal del entrevistado que no es bien visto por la sociedad, o sobre la práctica de actividades ilícitas como la evasión de impuestos, o la opinión acerca de la práctica del aborto ilegal.

La técnica de respuesta aleatorizada (TRA) es una herramienta para realizar entrevistas diseñada para proteger la privacidad del encuestado y evitar

el estigma social o el miedo a las represalias y al mismo tiempo reduce las respuestas deshonestas y el sesgo de no respuesta. En la TRA la respuesta a una pregunta sensible depende de 1) el verdadero estado del entrevistado respecto de la pregunta sensible y 2) el resultado de un dispositivo de aleatorización (Cruyff et al., 2008). El primero en proponer una TRA fue Stanley Warner en 1965. La TRA de Warner tiene a) un mecanismo para obtener información sobre una variable sensible binaria, y b) un estimador de la prevalencia del atributo sensible binario (Warner, 1965). La TRA de Warner consta de dos afirmaciones complementarias: Tengo el atributo sensible y No tengo el atributo sensible. El entrevistado responde solo a una de las dos afirmaciones y la elección de cuál debe responder se realiza con un dispositivo aleatorio, como un

dado o una ruleta. La elección de la pregunta a contestar se oculta al entrevistador, que solo recibe una respuesta de Sí o No sin saber qué pregunta fue respondida.

Las variables sensibles pueden ser dicotómicas, politómicas o cuantitativas según el atributo sensible estudiado. En el primer y segundo caso, el parámetro de interés es la proporción de individuos (prevalencia) en cada categoría del atributo sensible. En el último caso el parámetro de interés es la media poblacional.

Varias modificaciones se han propuesto a la TRA de Warner. Una familia de modificaciones cambia la pregunta complementaria por una pregunta incorrelacionada (Greenberg et al., 1969) o por una respuesta forzada (Boruch, 1971). Otra familia de modificaciones cambia la respuesta del entrevistado por una respuesta inocua como Rojo (Verde) en lugar de No (Sí) (Kuk, 1990) o por un valor codificado o mezclado (Eichhorn y Hayre, 1983; Bar-Lev et al. 2004).

Al usar una TRA se introduce un ruido aleatorio adicional a la respuesta del entrevistado lo que agrega más incertidumbre a los estimadores. Aun así, es posible obtener estimadores insesgados de parámetros unidimensionales, prevalencias o medias, pero estos estimadores tienen mayores varianzas estimadas que en el caso no sensible (Ardah y Oral, 2017).

La mayoría de artículos se enfocan en la estimación de la prevalencia de una variable aleatorizada binaria sensible. Un reducido número de artículos estiman la asociación entre una variable aleatorizada binaria y otra del mismo tipo, o con una variable binaria no sensible (Barabesi et al., 2012; Drane, 1976; Ewemooje y Amahia, 2015; Lee et al., 2013; Tamhane, 1981). Muy pocos artículos estiman el efecto de covariables no sensibles y factores en una variable respuesta binaria sensible. Este escenario es más interesante y desafiante porque el ruido aleatorio añadido mediante el uso de una TRA en la variable respuesta sensible provoca varianzas estimadas muy grandes, además porque es necesario estimar parámetros adicionales.

La regresión logística es el modelo más frecuentemente usado para modelar el efecto de covariables no sensibles en una variable respuesta aleatorizada binaria. Este modelo se ha utilizado

cuando la variable aleatorizada binaria se obtuvo bajo los diseños de Warner, de preguntas incorrelacionadas de Greenberg, de respuesta forzada de Boruch y de respuesta enmascarada de Kuk (Scheers y Dayton, 1988; Blair et al., 2015, van der Heijden y van Gils, 1996; Lensvelt-Mulders et al., 2006; van den Hout et al., 2007).

Abul-Ela et al. (1967) y Eriksson (1973) fueron los primeros en estimar las prevalencias de las categorías de variables aleatorizadas politómicas. Algunas preguntas que se han utilizado para obtener una variable aleatorizada ordinal son:

¿Cuántos días usó drogas ilegales la semana pasada?, donde las posibles respuestas son 0, 1, 2 o ≥ 3 (Kim y Warde, 2005), ¿Cuántos abortos tuvo una mujer?, ¿En cuántas ocasiones se tomaron drogas?, y ¿Cuánto dinero no reportaste en tu declaración de impuestos? (Eichhorn y Hayre, 1983). Otras preguntas producen natural-mente variables ordinales aleatorizadas, como

¿Qué tan de acuerdo estás con la interrupción legal del embarazo antes de las doce semanas?, donde las respuestas pueden ser muy en desacuerdo, en desacuerdo, de acuerdo y muy de acuerdo.

El número de artículos en que se estima la prevalencia de las categorías de variables ordinales aleatorizadas es mucho menor que para el caso binario, incluso hay menos artículos que estiman asociaciones entre una variable ordinal aleatorizada y otra del mismo tipo, o con covariables no sensibles (Kim y Warde, 2005).

El modelo de regresión ordinal se utiliza cuando la variable respuesta es ordinal y no sensible. La estimación de los parámetros de este modelo cuando la variable respuesta es no sensible se ha realizado con ambos paradigmas estadísticos, frecuentista y Bayesiano. Sin embargo, no se encontraron artículos donde la variable respuesta ordinal sea sensible.

Este trabajo tiene como objetivo proponer un modelo para estudiar el efecto de covariables no sensibles en una variable respuesta ordinal aleatorizada bajo el diseño de respuesta forzada, usar la inferencia Bayesiana para estimar los parámetros y evaluar el desempeño de los estimadores propuestos.

El objetivo de este trabajo es usar el enfoque Bayesiano para medir el efecto de covariables no

sensibles en una variable respuesta ordinal aleatorizada obtenida bajo el diseño de respuesta forzada y evaluar el desempeño de los estimadores propuestos.

4.2. Material y métodos

Considere que hay un atributo ordinal sensible que se mide por una variable ordinal Z con categorías 1, 2,…,J, es decir, tiene un total de J categorías. Si al entrevistado se le hace una pregunta directa sobre su pertenencia a una categoría de Z, él o ella puede no responder honestamente o puede no contestar. La información sobre la pertenencia a una de las J categorías de Z se puede obtener mediante una TRA que generará una variable ordinal aleatorizada Y estrechamente relacionada con Z. El uso de la TRA permite motivar la participación del entrevistado al darle confianza para responder su verdadera pertenencia a la categoría de la variable ordinal aleatorizada Y y al mismo tiempo proteger su privacidad.

Para cada individuo i, sea 𝑍_𝑖 = 𝑗; j=1,…,J, la categoría verdadera pero desconocida para la variable ordinal sensible y 𝑌_𝑖 = 𝑗 la categoría de la variable ordinal aleatorizada observada obtenida utilizando el dispositivo de aleatorización de Hopkins que se describe a continuación.

Dispositivo de aleatorización de Hopkins

Liu y Chow (1976) propusieron un modelo de respuesta aleatorizado cuantitativo discreto utilizando el dispositivo de aleatorización de Hopkins. El dispositivo consiste en una jarra que contiene bolas de dos colores diferentes, por ejemplo, verde y azul, en proporciones 𝑔 y 1 − 𝑔 , respectivamente. Las bolas azules están marcadas con los números 1,…,J en proporciones q₁ ,…, q_J (Figura 1). Se le pide al entrevistado que dé la vuelta a la jarra, la agite bien y permita que

Figura 1: Dispositivo de aleatorización de Hopkins para el diseño de respuesta forzada para variables ordinales.

aparezca una bola en el cuello. Si la bola es verde, el entrevistado debe reportar su verdadera categoría del atributo ordinal sensible, pero si la bola es azul debe informar el número marcado en la bola. q_j tiene el efecto de aumentar las respuestas para la j-ésima categoría. Este proceso se lleva a cabo sin la observación del entrevistador lo que protege la privacidad del entrevistado y motiva su participación. El método descrito implementa un diseño de respuesta forzada para obtener respuestas sobre un atributo ordinal sensible.

El dispositivo de aleatorización de Hopkins relaciona la probabilidad de la j-ésima categoría para la variable ordinal sensible para el i-ésimo

donde 𝜀_𝑖~𝑁(0,1) ; 𝑖 = 1, … , 𝑛, ; 𝒙 ^𝑇 = (𝑥_𝑖1, … , 𝑥 _𝑖𝑝) es el vector de observaciones correspondientes a p covariables no sensibles para el i-ésimo individuo. Sea 𝛾₀, … , 𝛾 _𝐽 un conjunto de umbrales, la j-ésima categoría de Z se obtiene cuando 𝛾_𝑗−1 < 𝑍^𝐿 ≤ 𝛾_𝑗 para j=1,…,J. Para evitar problemas de estimación se supone que 𝛾₀ = −∞ y 𝛾_𝐽 = ∞ , por lo que la probabilidad de la primera y la última categoría ordinal están bien definidas.

�unable to handle picture here, no embed or link𝑖El modelo de regresión probit ordinal está dado por (Kruschke, 2014):

�unable to handle picture here, no embed or link∗ 𝑖 𝛾 𝑗−𝒙^𝑇𝜷

𝜋_𝑖𝑗 = 𝑃(𝑍_𝑖 = 𝑗|𝒙_𝑖) = Φ ( _𝜎 ) −

individuo, ^∗ (

), a la probabilidad de

𝛾𝑗−1−𝒙^𝑇𝜷

(2)

𝜋_𝑖𝑗 = 𝑃 𝑍_𝑖 = 𝑗

la j-ésima categoría de la variable ordinal aleatorizada observada, 𝜋_𝑖𝑗 = 𝑃(𝑌_𝑖 = 𝑗) mediante la ecuación:

Φ ( ^𝑖 )

�unable to handle picture here, no embed or link𝜎

�unable to handle picture here, no embed or link𝑖La Ecuación 2 usa el hecho de que Z^L tiene distribución normal con media 𝒙^𝑇𝜷, y desviación

�unable to handle picture here, no embed or link𝑖𝑗𝜋_𝑖𝑗 = 𝑔𝜋^∗ + (1 − 𝑔)𝑞_𝑗 (1) donde 𝑖 = 1, … , 𝑛 𝑗 = 1, … , 𝐽

El término 𝑔 está relacionado con la probabilidad de responder a una pregunta directa y no debe ser cero ni uno. Si 𝑔 = 0 , solo se obtienen respuestas forzadas y no se obtiene información del atributo ordinal sensible, por otro lado, si 𝑔 = 1 solo se utilizan preguntas directas sobre el atributo ordinal sensible y el entrevistado

estándar σ.

Modelo de regresión probit ordinal aleatorizado

�unable to handle picture here, no embed or link𝑖𝑗El objetivo de este trabajo es medir el efecto de las covariables no sensibles en la verdadera pero desconocida variable ordinal sensible Z, usando la variable ordinal aleatorizada observada Y y la relación entre sus probabilidades, 𝜋^∗ y 𝜋_𝑖𝑗, dadas por las Ecuaciones 1 y 2.

El modelo de regresión probit ordinal aleatorizado propuesto es:

sentirá amenazada su privacidad y se negará a

𝑃(𝑌 = 𝑗|𝒙 ) = 𝑔𝑃(𝑍

= 𝑗|𝒙 ) + (1 − 𝑔)𝑞

contestar, o lo hará de forma deshonesta.

𝑖 𝑖

𝑖 𝑖 𝑗

∗

Modelo de regresión probit ordinal

𝜋_𝑖𝑗 = 𝑔𝜋_𝑖𝑗 + (1 − 𝑔)𝑞_𝑗

𝛾𝑗−𝒙^𝑇𝜷 𝛾𝑗−1−𝒙^𝑇𝜷

El modelo de regresión probit ordinal permite medir el efecto de las covariables no sensibles en una variable respuesta ordinal no sensible. La motivación del modelo es simple. Se supone que

�unable to handle picture here, no embed or link𝑖=1una variable ordinal Z es una discretización de una variable latente continua Z^L que es una

𝜋 = 𝑔 [Φ ( ^𝑖 ) − Φ ( ^𝑖 )] +

�unable to handle picture here, no embed or link𝑖𝑗𝜎 𝜎

(1 − 𝑔 )𝑞_𝑗 (3)

La función de verosimilitud para el Modelo 3 es:

combinación lineal de covariables no sensibles,

𝑋₁, … , 𝑋 _𝑝 a la que se agrega un error aleatorio con

𝐿(𝜷, 𝜎, 𝜸|𝑿, 𝒚) = ∏^𝑛

𝐽

�unable to handle picture here, no embed or link∑𝑗=1

{𝐼(𝑌_𝑖 = 𝑗) × 𝜋_𝑖𝑗}

(4)

distribución normal estándar:

�unable to handle picture here, no embed or link𝑖 𝑖𝑍^𝐿 = 𝒙^𝑇𝜷 + 𝜀_𝑖,

donde I(·) es la función indicadora y 𝜸 = (𝛾₀, 𝛾₁, … , 𝛾_𝐽−1, 𝛾_𝐽) ^𝑇

En la revisión de la literatura no se encontraron artículos que aborden la estimación frecuentista o

Bayesiana de los parámetros del modelo de regresión probit ordinal aleatorizado.

Modelación Bayesiana

El teorema de Bayes establece que la distribución posterior del parámetro 𝜽 dados los datos D, π(𝜽|D), es proporcional al producto de la verosimilitud de los datos, L( 𝜽 |D), y la distribución a priori del parámetro 𝜽, π( 𝜽), es decir:

𝜋(𝜽|𝑫) ∝ 𝐿(𝜽|𝑫)𝜋(𝜽)

Un estimador Bayesiano de θ es el que minimiza la esperanza bajo la distribución posterior de la función de pérdida. Cuando la función de pérdida es el error cuadrado o lineal absoluta, entonces el estimador Bayesiano de 𝜽 es la media o la mediana de la distribución posterior, según corresponda.

La distribución a priori conjunta se da como un producto de distribuciones marginales a prioris independientes:

�unable to handle picture here, no embed or link𝑘=1𝜋(𝜶, 𝜸, 𝜎) = 𝜋(𝛼₀) ∏ ^𝑝 𝜋(𝛼_𝑘) ×

�unable to handle picture here, no embed or link𝑗=1∏^𝐽−1 𝜋(𝛾_𝑗) × 𝜋 (𝜎) (5)

Las estimaciones de parámetros unidimensionales como son la prevalencia o media de variables sensibles después de usar alguna de las diferentes TRA tienen varianzas estimadas grandes. Un efecto similar ocurrirá al estimar vectores de parámetros como los efectos de covariables no sensibles en una variable respuesta ordinal aleatoria. Esto motiva a usar a prioris no informativas o a prioris que reduzcan las estimaciones de los parámetros α_j.

Se estudian cuatro a prioris conjuntas, todas

ellas usan 𝛼₀~𝑁(^1+𝐽, 𝐽 ²) , 𝜎~𝑁( ^𝐽 , 10𝐽) , y

2 1000

En el modelo de regresión probit ordinal aleatorizado 𝜽 = (𝜷^𝑇, 𝜸 ^𝑇, 𝜎 )^𝑇 y 𝑫 = {𝑿, 𝒚 } , donde X es la matriz diseño que contiene la información de las covariables no sensibles. La distribución posterior es:

𝜋(𝜷, 𝜎, 𝜸|𝑿, 𝒚) ∝ 𝐿(𝜷, 𝜎, 𝜸|𝑿, 𝒚) 𝜋(𝜷, 𝜎, 𝜸)

La verosimilitud de los datos bajo el diseño de respuesta forzada se da en las Ecuaciones 3 y 4. Obtener la media de la distribución posterior no es analíticamente posible. En su lugar, la distribución posterior se aproxima por Cadenas de Markov Monte Carlo (MCMC, por sus siglas en inglés). Para mejorar la convergencia de las MCMC se estandariza X. Sea X^* la matriz estandarizada, es decir, la k-ésima columna es

𝛾_𝑗~𝑁(𝑗 + 0.5, 𝐽 ²) ; j = 2,…,J-2. La diferencia entre las cuatro a prioris conjuntas es la distribución de los 𝛼_𝑘 , k = 1,…,p. Estas distribuciones a prioris son Normal (N), Doble exponencial (DE), t y Cauchy. La lista de las

distribuciones a prioris utilizadas y los parámetros adicionales necesarios son:

dada por 𝑥^∗ = (𝑥

− 𝑥̅

)/𝑠

, donde 𝑥̅

y 𝑠

�unable to handle picture here, no embed or link𝛼_𝑘 ~ 𝑁 (0, 𝐽 ²)(6)𝛼_𝑘 ~ 𝐷𝐸 (0, 𝜆 ), 𝜆 ~ 𝑈(0.001,10)(7)𝛼 _𝑘 ~ 𝑡 (0, 𝜎 ², 𝑣 ^∗ + 1 ),𝑡𝑣^∗ ~ 𝐸𝑥𝑝 (1/29),𝜎_𝑡 ~ 𝑈(0.001, 1000)(8)𝛼 _𝑘 ~ 𝑡 (0, 𝜎 ², 1 ), 𝜎 ~ 𝑈(0.001, 10 𝑡 𝑡(9)𝑖𝑘

𝑖𝑘

𝑘 𝑥𝑘

𝑘 𝑥𝑘

En la literatura, las a prioris para los 𝛼

que se

son la media muestral y la desviación estándar para la j-ésima covariable. Para evitar problemas de identificabilidad en la estimación de los umbrales, dos de ellos se fijan: γ₁ = 1.5 y γ_J-1 = J-

0.5 (Kruschke, 2014). Asociado con las covariables estandarizadas, se introduce un nuevo

𝑗

usan cuando se modela una variable respuesta ordinal no sensible tienen una distribución normal con media cero y varianzas 100, 0.5 o 0.1 (Van, 2017; Xie et al., 2009). La distribución a priori

𝛼_𝑘 ~ 𝑁 (0, 𝐽 ²) es vaga en la escala estandarizada

vector de parámetros

𝑇

�unable to handle picture here, no embed or link.𝜶 = (𝛼₀, 𝛼₁, … , 𝛼_𝑝)

para el modelo de regresión probit ordinal

aleatorizado y la varianza del hiperparámetro es

La distribución posterior se convierte en:

𝜋(𝜶, 𝜸, 𝜎|𝑿^∗, 𝒚) ∝ 𝐿(𝜶, 𝜸, 𝜎|𝑿^∗, 𝒚) × 𝜋(𝜶, 𝜸, 𝜎)

totalmente definido por el número de categorías. La distribución a priori 7 se usa comúnmente en estimación Lasso Bayesiana para reducir las

�unable to handle picture here, no embed or link𝑡pendientes estimadas porque su masa se concentra cerca de cero y el hiperparámetro λ toma cualquier valor en el intervalo (0.001, 10). La distribución a priori 9 es un caso especial de las distribuciones a priori 8 cuando los grados de libertad son 1. Ambas a prioris tienen media cero, de modo que 𝜎² está en un rango amplio de 10^-6 a 10⁶.

Después de obtener los 𝛼′𝑠 en las MCMC, los parámetros en la escala original se obtuvieron con las ecuaciones:

𝑝

distribución posterior y cuando el tamaño de la muestra es grande, lo hace la verosimilitud de los datos. Un tamaño de la muestra pequeño determinó de acuerdo con el número de parámetros en el modelo, y el tamaño de la muestra grande como el triple del tamaño de muestra pequeño.

Número de covariables no sensibles, p. Se estudiaron una y dos covariables no sensibles, p=1,2.

Por lo tanto, cuando Y tiene J=3 categorías ordinales y el modelo p=1 covariable no sensible,

𝛼_𝑘

𝛽₀ = 𝛼₀ − ∑ _𝑠

𝑥̅_𝑘,

los tamaños de muestra fueron n=125, 375 y si p=2 los tamaños de muestra fueron n=150, 450.

𝛽_𝑘

𝛼_𝑘

=

𝑠_𝑥

𝑘=1

𝑥𝑘

𝑘 = 1, … , 𝑝

Cuando Y tiene J=4 categorías ordinales y el modelo tenía p=1 covariable no sensible, los

𝑘

Estudio de simulación

El desempeño de los estimadores Bayesianos bajo las distribuciones a prioris estudiadas se evaluó mediante simulación Monte Carlo. El dispositivo de aleatorización de Hopkins que induce el diseño de respuesta forzada para un atributo ordinal sensible se usó con 𝑔 = 2 ⁄3 y

tamaños de muestra eran n=150, 450, y si p=2, los tamaños de muestra eran n=175, 525 (Tabla 1).

Generación de datos

Las p covariables no sensibles X_k se generaron de forma independiente de una distribución uniforme en (0.9, 1.9). Los parámetros verdaderos se eligieron de acuerdo con el número de covariables no sensibles en el modelo, si p=1

𝑞₁

1

_{�unable to handle picture here, no embed or link}= ⋯ = 𝑞_𝐽 = _𝐽

lo que induce respuestas

directas con probabilidad 2/3.

Factores estudiados. Muchos factores pueden afectar el rendimiento de los estimadores Bayesianos del modelo de regresión probit ordinal aleatorizado, pero en esta investigación solo se consideran tres. Los factores estudiados en la simulación Monte Carlo fueron el número de categorías de la variable dependiente, el tamaño de la muestra y el número de covariables no sensibles.

Número de categorías de la variable respuesta ordinal aleatorizada, J. Se consideraron tres y cuatro categorías ordinales, J=3,4.

Tamaño de la muestra, n. En el enfoque frecuentista cuando el tamaño de la muestra es grande, el estimador de máxima verosimilitud tiene propiedades óptimas. Además, en el enfoque Bayesiano cuando el tamaño de la muestra es pequeño, la distribución a priori domina la

Tabla 1: Tamaños de muestra, parámetros y umbrales utilizados en la simulación.

entonces β = (-10,10)^T y si p=2 entonces β = (-10,10,5)^T. De esta manera, 𝛽₁ = 10 y 𝛽₂ = 5 están lejos de cero y se espera que los estimadores

Bayesianos no sean afectados por el efecto de contracción producido por las a prioris. Los umbrales se definen apropiadamente para que todas las categorías tengan observaciones (Tabla 1).

�unable to handle picture here, no embed or link𝑖La variable latente continua se obtuvo con 𝑍^𝐿 =

�unable to handle picture here, no embed or link𝑖𝒙^𝑇𝜷 + 𝜀_𝑖, 𝜀_𝑖~𝑁(0,1/2²) . Después la variable ordinal sensible Z se obtuvo utilizando los umbrales.

Para simular el proceso de asignar al entrevistado a responder la pregunta directa con probabilidad 𝑔, o a dar una respuesta forzada con probabilidad 1 − 𝑔 , se generó una variable W con

Desempeño de los estimadores

Un criterio para comparar estimadores es el error cuadrático medio (ECM) porque mide la dispersión de las estimaciones respecto al valor verdadero del parámetro y tiene la propiedad de que es la suma del cuadrado del sesgo y la varianza del estimador. Otros criterios usados son dos propiedades de los intervalos de credibilidad que son la longitud (L) y cobertura (C). Es deseable que los intervalos de credibilidad sean de la menor longitud para que sean más informativos y que la proporción de intervalos que contienen el parámetro sea aproximadamente del valor nominal, por ejemplo 95%. Estos criterios se definen como

distribución Bernoulli con parámetro 𝑔. Es decir, para cada Z_i = j, se generó una variable auxiliar W distribuida Bernoulli con parámetro 𝑔. Si W=1

𝑅

^{�unable to handle picture here, no embed or link�unable to handle picture here, no embed or link}𝐸^̂𝐶𝑀(𝛽̃) = ¹ ∑ ¹

𝑅 𝑝

𝑟=1

𝑝

2

∑(𝛽̃_𝑘,𝑟 − 𝛽_𝑘)

𝑘=1

entonces Y_i = j, de lo contrario Y_i es una realización de una distribución multinomial con

1

𝐿̂(𝛽̃) =

𝑅

1 1

_{�unable to handle picture here, no embed or link ∑}

𝑅 𝑝

𝑝

∑{𝐿𝑆 (𝛽̃_𝑘,𝑟

) − 𝐿𝐼 (𝛽̃_𝑘,𝑟)}

�unable to handle picture here, no embed or linkparámetros 𝜋₁ = ⋯ = 𝜋 _𝐽 = _𝐽.

𝑟=1

𝑘=1

𝑅 𝑝

Se estudiaron cuatro estimadores Bayesianos,

1 1

𝐶̂(𝛽̃) = ∑ ∑[𝐿𝐼 (𝛽̃

) ≤ 𝛽

definidos por la distribución a priori usada para α_j. El estimador Bayesiano EB.N se obtiene al usar

𝑅 𝑝

𝑟=1

𝑘=1

𝑘,𝑟 𝑘

̃

la a priori Normal (6), el EB.DE se obtiene con la a priori Doble Exponencial (7), el EB.T se obtiene con la a priori t (8) y el EB.C con la a priori Cauchy (9).

Las distribuciones marginales posteriores se obtuvieron con el paquete R2jags (Su y Yajima, 2015) implementado en el lenguaje R (R Core Team, 2016). Las MCMC se obtuvieron con tres cadenas de 10000 iteraciones, un quemado de 5000, y un adelgazamiento de uno de cada cinco. Aunque no se muestran los resultados, la convergencia de MCMC se verificó mediante el factor de reducción de escala potencial multivariado de Gelman y Brooks implementado en el paquete coda (Plummer et al., 2006). Los intervalos de credibilidad del 95% para 𝛽_𝑘,

k=1,…,p

𝐼𝐶𝑟 [𝛽̃_𝑘] = [𝐼𝐶𝑟𝐼 (𝛽̃_𝑘), 𝐼𝐶𝑟𝑆 (𝛽̃_𝑘)],

�unable to handle picture here, no embed or link𝑘donde 𝐼𝐶𝑟𝐼 (𝛽^̃ ) es el cuantil 0.025 de la distribución posterior marginal de 𝛽_𝑘 y

𝐼𝐶𝑟𝑆 (𝛽̃_𝑘) el cuantil 0.975. El proceso se repitió R = 1000 veces.

≤ 𝐿𝑆 (𝛽_𝑘,𝑟)]

donde 𝛽̃_𝑘,𝑟 es el r-ésimo valor del estimador 𝛽̃_𝑘

de 𝛽_𝑘.

Para comparar un estadístico 𝑇̂ bajo dos condiciones, A y B, se utilizó el porcentaje de ganancia relativa (PGR) del estadístico cuando se pasa de la condición A a la condición B:

𝑇̂ (𝐴) − 𝑇̂ (𝐵)

�unable to handle picture here, no embed or link𝑃𝐺𝑅. 𝑇̂ (𝐴, 𝐵) = × 100

𝑇̂ (𝐵)

donde 𝑇̂ corresponde a los criterios 𝐸^̂𝐶𝑀 , 𝐿̂ y 𝐶̂.

4.3. Resultados

Los resultados de la estimación del intercepto no se informan porque interesa solo el efecto de las covariables no sensibles en la variable respuesta ordinal aleatorizada. La Tabla 2 muestra el ECM,

la longitud y la cobertura de los estimadores Bayesianos y sus intervalos de credibilidad.

Efecto del número de categorías

Todos los estimadores Bayesianos tienen ECM y longitud de intervalo creíble más pequeño cuando la variable respuesta ordinal aleatorizada tiene cuatro categorías (Tabla 2).

La Tabla 3 muestra el porcentaje de ganancia relativa para los tres criterios comparando J=3 y J=4 categorías de la variable respuesta ordinal aleatorizada. El efecto del número de categorías es muy fuerte en el ECM de los estimadores Bayesianos porque el porcentaje de ganancia relativa es mayor, entre 94% y 619%, cuando Y tiene tres categorías. También es fuerte el efecto en la longitud de los intervalos de credibilidad en los estimadores Bayesianos pues el porcentaje de ganancia relativa es mayor, entre 42% y 142%, cuando Y tiene tres categorías. Y finalmente, es débil en la cobertura de intervalos de credibilidad ya que el porcentaje de ganancia relativa suele ser mayor cuando Y tiene cuatro categorías.

Efecto del tamaño de muestra

Todos los estimadores Bayesianos tienen menor ECM, longitud y cobertura de intervalos de credibilidad cuando el tamaño de la muestra es grande (Tabla 2).

La Tabla 4 muestra el PGR para los tres criterios comparando el tamaño de muestra pequeño (TMP) con el tamaño de muestra grande (TMG). De acuerdo con la teoría asintótica, los estimadores Bayesianos muestran un mejor desempeño con tamaño de muestra grande. El efecto del tamaño de la muestra es muy fuerte en el ECM de los estimadores Bayesianos porque el porcentaje de ganancia relativa es entre 226% y 883% mayor que con el tamaño de muestra pequeño. El efecto del tamaño de la muestra es fuerte en la longitud de los intervalos de credibilidad de los estimadores Bayesianos porque el porcentaje de ganancia relativa es entre 87% y 154% mayor cuando el tamaño de la muestra es pequeño. Finalmente, el efecto del tamaño de la muestra es débil en la cobertura de los intervalos de credibilidad de los estimadores Bayesianos lo que

Tabla 2: Valores de los criterios utilizados para comparar los estimadores Bayesianos bajo los factores estudiados.

Tabla 3: Porcentaje de ganancia relativa en ECM, L y C al comparar tres contra cuatro categorías sensibles.

PGR.𝑇̂ (3,4) = 100 (𝑇̂ (3) − 𝑇̂ (4))⁄𝑇̂ (4)

se evidencia en que la cobertura es mayor más frecuentemente con el tamaño de muestra grande, y el porcentaje de ganancia relativa varía entre 1 y 4%.

Efecto del número de covariables no sensibles

Todos los estimadores Bayesianos tienen mejor ECM, longitud y cobertura de intervalos de credibilidad cuando hay dos covariables no sensibles (Tabla 2).

La Tabla 5 muestra el porcentaje de ganancia relativa para los tres criterios comparando p=1 con p=2. El efecto del número de covariables no sensibles en el ECM de los estimadores Bayesianos es muy fuerte, pues es mayor cuando hay una covariable que cuando hay dos covariables, porque el porcentaje de ganancia relativa varía entre un 25 % y un 157 %. Este efecto del número de covariables no sensibles es fuerte en la longitud de los intervalos de credibilidad de los estimadores Bayesianos, pues

es mayor cuando hay una covariable, ya que el porcentaje de ganancia relativa varía entre 14% y 63%. Y finalmente, es débil en la cobertura de los intervalos de credibilidad de los estimadores Bayesianos, puesto que es más frecuente que el porcentaje de ganancia relativa sea mayor con dos covariables no sensibles.

Comparación del estimador de máxima verosimilitud con los estimadores Bayesianos estudiados

Existen dos enfoques importantes para realizar inferencia, el frecuentista y el Bayesiano. El estimador de máxima verosimiltud (EMV) es uno de los métodos más importantes en el enfoque frecuentista. En el contexto del modelo de regresión probit ordinal aleatorizado es natural comparar el desempeño de los estimadores, de máxima verosimiltud y Bayesianos, en función de sus esperanzas, sus intervalos de confianza (IC) o de credibilidad (ICr), su longitud y su cobertura.

Tabla 4: Porcentaje de ganancia relativa en ECM, L y C cuando se compara el tamaño de muestra pequeño con el grande.

PGR.𝑇̂ (𝑇𝑀𝑃, 𝑇𝑀𝐺) = 100 (𝑇̂ (𝑇𝑀𝑃) − 𝑇̂ (𝑇𝑀𝐺))⁄𝑇̂ (𝑇𝑀𝐺)

Se realizó una simulación con R = 1000 repeticiones en el escenario en que la variable respuesta ordinal aleatorizada tiene 4 categorías que se explican con dos variables independientes

donde 𝑍_1−𝛼/2 es el cuantil 1 − 𝛼/2 de la distribución normal estándar y 𝐸𝐸(𝛽̂_𝑘) es el elemento jj de la inversa de la matriz de Fisher.

Los límites de confianza inferior y superior de 𝛽_𝑘

y se usa un tamaño de muestra de 525 que es el

son 𝐿𝐶𝐼(𝛽̂ ) = 𝛽̂

− 𝑍

𝐸𝐸(𝛽̂ ) y

mayor de los tamaños de muestra considerados en

𝑘

𝐿𝐶𝑆(𝛽̂ ) = 𝛽̂ + 𝑍

𝑘 1−𝛼/2 𝑘

𝐸𝐸(𝛽̂ ) .

este trabajo. Este tamaño de muestra se puede

considerar grande y de acuerdo con la teoría

𝑘 𝑘

1−𝛼/2 𝑘

asintótica, el EMV tiene distribución normal con media el vector de parámetros verdadero y matriz de covarianzas la inversa de la matriz de información de Fisher, esto es, se espera que el EMV tenga un excelente desempeño.

La esperanza e intervalos de confianza (IC) del

100(1 − 𝛼 )% para el EMV de 𝛽_𝑘 son:

Aunque no se muestran, los valores del ECM de los estimadores Bayesianos estudiados son similares a los reportados en la Tabla 2. Los datos se generaron con 𝛽₁ = 10 y 𝛽₂ = 5 por lo que la

esperanza de los EMV, 13.97 y 6.38,

sobreestiman a los valores verdaderos, y muestra que los EMV de los parámetros del modelo de regresión probit ordinal aleatorizado no son

𝐸̂ (𝛽̂_𝑘) =

𝑅

1

�unable to handle picture here, no embed or link_𝑅 ^{∑ 𝛽̂}𝑘,𝑟

𝑟=1

insesgados. Por el contrario, los estimadores Bayesianos son más cercanos a los valores verdaderos y el mejor de ellos es EB.DE, que usa

𝛽̂_𝑘

^{± 𝑍}1−𝛼/2

𝐸𝐸(𝛽̂_𝑘)

la distribución a priori doble exponencial.

Tabla 5: Porcentaje de ganancia relativa en ECM, L y C al comparar una contra dos covariables sensibles.

PGR.𝑇̂ (1,2) = 100 (𝑇̂ (1) − 𝑇̂ (2))⁄𝑇̂ (2)

La sobreestimación de los EMV causa que la media del 𝐿𝐶𝐼(𝛽̂₁) sea 12.10 que es mayor que el valor verdadero 𝛽₁ = 10 y como consecuencia la cobertura del IC para 𝛽₁ es apenas de 0.04, muy lejos del valor nominal de 0.95. Por el contrario, los estimadores Bayesianos tienen ICr con medias

del 𝐿𝐶𝑟𝐼(𝛽̃ ) y 𝐿𝐶𝑟𝑆(𝛽̃ ) que permiten que el

El otro criterio usado para comparar los estimadores fue la longitud de los IC e ICr. Los EMV tienen longitud aproximadamente el doble de la longitud de los ICr para 𝛽₁ y de aproximadamente 2.6 veces mayores para 𝛽₂.

Aún en el escenario más favorable para el EMV por tener un tamaño de muestra grande, se

^{1 1} observa que los estimadores Bayesianos son

valor verdadero 𝛽₁ = 10 esté en el ICr lo que se refleja en que su cobertura sea de 0.93 o 0.94 que son más cercanos al valor nominal de 0.95.

Las medias de 𝐿𝐶𝐼(𝛽̂ ) y 𝐿𝐶𝑆(𝛽̂ ) son 4.74 y

mejores pues son aproximadamente insesgados, con menores longitudes y cobertura de los ICr.

^{2 2} Conclusiones

8.03 que permiten que el valor verdadero 𝛽₂ = 5 y esté contenido en el IC. Estos IC tienen cobertura de 0.62 que sigue siendo menor al valor nominal de 0.95. Por el contrario, los estimadores

Bayesianos tienen ICr con medias de 𝐿𝐶𝑟𝐼(𝛽̃ ) y

Los estimadores Bayesianos propuestos estiman correctamente los parámetros del modelo de regresión probit ordinal aleatorizado aun cuando se usa el diseño de respuesta forzada inducido por

_̃ ² el dispositivo de aleatorización de Hopkins para

𝐿𝐶𝑟𝑆(𝛽₂) que permiten que el valor verdadero

𝛽₂ = 5 esté en el ICr lo que se refleja en que su cobertura sea de 0.94 o 0.95 que coinciden con el valor nominal de 0.95.

producir una variable respuesta ordinal aleatorizada. El valor nominal de cobertura del 95% de los intervalos de credibilidad se alcanzó

Tabla 6. Medias de las esperanzas, límites inferiores y superiores de los intervalos de confianza y de credibilidad, longitud y cubrimiento de los estimadores de máxima verosimilitud y Bayesianos del modelo de regresión probit ordinal aleatorizado.

generalmente con los tamaños de muestra más grandes, lo que es consistente con la teoría asintótica. Un resultado adicional es que todas las covariables no sensibles incluidas en el modelo de regresión probit ordinal aleatorizado rechazan que

𝛽_𝑘 ≠ 0, j=1,…, p porque las longitudes de los intervalos de credibilidad indican que las estimaciones 𝛽̃_𝑘 están lejos de cero.

Los estimadores Bayesianos de los 𝛽_𝑘 de las covariables no sensibles en el modelo de

regresión ordinal aleatorizado tienen óptimo error cuadrático medio, longitud y cobertura de los intervalos de credibilidad cuando se utilizan más categorías de la variable respuesta ordinal aleatorizada, tamaño de muestra más grande y mayor número de covariables no sensibles.

Sorprendentemente, los estimadores Bayesianos tienen menor error cuadrático medio usando la a priori doble exponencial para cada coeficiente de las covariables no sensibles. Por lo tanto, el mejor estimador Bayesiano en el modelo de regresión probit ordinal aleatorizado se obtiene utilizando la a priori doble exponencial.

4.4. Referencias

Abul-Ela, A.L.A., Greenberg, G.G., Horvitz, D.G. (1967). A Multiproportions Randomized Response Model. Journal of the American Statistical Association, 62, 990-1008. https://doi.org/10.2307/2283687

Ardah, H.I. Oral, E. (2017). Model Selection in Randomized Response Techniques for Binary Responses. Communications in Statistics - Theory and Methods, 47, 3305-3323. https://doi.org/10.1080/03610926.2017.135362 6

Bar-Lev, S.K., Bobovitch, E., Boukai, B. (2004). A Note on Randomized Response Models for Quantitative Data. Metrika, 60, 225-250. https://doi.org/10.1007/s001840300308

Barabesi, L., Franceschi, S., Marcheselli, M. (2012). A Randomized Response Procedure for Multiple Sensitive Questions. Statistical Papers, 53, 703-718.

Blair, G., Imai, K., Zhou, Y.-Y. (2015). Design and Analysis of the Randomized Response Technique. Journal of the American Statistical Association, 110, 1304-1319. https://doi.org/10.1080/01621459.2015.105002 8

Boruch, R.F. (1971). Assuring Confidentiality of Responses in Social Research: a Note on Strategies. The American Sociologist, 6, 308-

311.

Cruyff, M.J.L.F., van den Hout, A., van der Heijden, P.G.M. (2008). The Analysis of Randomized Response Sum Score Variables. Journal of the Royal Statistical Society, Series B, 71, 21-30.

http://dx.doi.org/10.1111/j.1467-9868.2007.00624.x

Drane, W. (1976). N the Theory of Randomized Responses to Two Sensitive Questions. Communications in Statistics - Theory and Methods, 5, 565-574. https://doi.org/10.1080/03610927608827375

Eichhorn, B.H., Hayre, L.S. (1983). Scrambled Randomized Response Methods for Obtaining Sensitive Quantitative Data. Journal of Statistical planning and Inference, 7, 307–316.

Eriksson, S.A. (1973). A New Model for Randomized Response. International Statistical Review, 41, 101-113.

https://doi.org/10.2307/1402791

Ewemooje, O.S., Amahia, G.N. (2015). Improved Randomized Response Technique for Two Sensitive Attributes. Afrika Statistika, 10, 839-

853.

https://doi.org/10.16929/as/2015.639.78 Greenberg, B., Abul-Ela, A., Simmons, W.,

Horvitz, D.G. (1969). The Unrelated Question Randomized Response Model: Theoretical Framework. Journal of the American Statistical Association, 64, 520-539. https://doi.org/10.2307/2283636

Kim, J.-M., Warde, W.D. (2005). Some New Results on the Multinomial Randomized Response Model. Communications in Statistics. Theory and Methods, 847-856. https://doi.org/10.1081/STA-200054378

Kruschke, J.K. (2014). Doing Bayesian Data Analysis: A Tutorial with R, JAGS, and Stan. Academic Press.

Kuk, A.Y.C. (1990). Asking Sensitive Question Indirectly. Biometrika, 77, 436-438. https://doi.org/10.1093/biomet/77.2.436

Lee, C.S., Sedory, S.A., Singh, S. (2013). Estimating at Least Seven Measures of Qualitative Variables from a Single Sample using Randomized Response Technique. Statistics and Probability Letters, 83, 399-409. https://doi.org/10.1016/j.spl.2012.10.004

Lensvelt-Mulders, G.J., van der Heijden, P.G., Laudy, O., van Gils, G. (2006). A Validation of a Computer-Assisted Randomized Response

Survey to Estimate the Prevalence of Fraud in Social Security. Journal of the Royal Statistical Society A., 169, 305-318.

A Validation of a Computer-Assisted Randomized Response Survey to Estimate the Prevalence of Fraud in Social Security on JSTOR

Liu, P.T., Chow, L.P. (1976). A new discrete quantitative randomized response model. ACM SIGSIM Simulation Digest, 7, 30-31. https://doi.org/10.1080/01621459.1976.104814 79

Plummer, M., Best, N., Cowles, K., Vines, K. (2006). CODA: Convergence diagnosis and output analysis for MCMC. R News, 6, 7-11.

CODA: Convergence diagnosis and output analysis for MCMC

R Core Team (2016). R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing, Vienna, Austria.

https://www.R-project.org/

Scheers, N., Dayton, C. (1988). Covariate Randomized Response Models. Journal of the American Statistical Association, 83, 969-974. https://doi.org/10.1080/01621459.1988.104786 86

Su, Y.-S., Yajima, M. (2015). R2jags: Using R to run JAGS.

Tamhane, A. (1981). Randomized Response Techniques for Multiple Sensitive Attributes. Journal of the American Statistical Association, 76, 916-923.

https://doi.org/10.1080/01621459.1981.104777 41

van den Hout, A., van der Heijden, P.G.M., Gilchrist, R. (2007). The Logistic Regression Model with Response Variables Subject to Randomized Response. Computational Statistics & Data Analysis, 51, 6060-6069. https://doi.org/10.1016/j.csda.2006.12.002

van der Heijden, P., van Gils, G. (1996). Some logistic regression models for randomized response data. In Proceedings of the 11th International Workshop on Statistical Modeling, 341-348.

Van, M.H. (2017). Determinants of the Levels of Development Based on the Human Development Index: Bayesian Ordered Probit Model. International Journal of Economics and Financial Issues, 7, 425.

Warner, S. (1965). Randomized Response: A Survey Technique for Eliminating Evasive

Answer Bias. Journal of the American Statistical Association, 60, 63-69. https://doi.org/10.1080/01621459.1965.104807 75

Xie, Y., Zhang, Y., Liang, F. (2009). Crash Injury Severity Analysis Using Bayesian Ordered Probit Models. Journal of Transportation Engineering, 135, 18-25.