Análisis estadístico de respuestas ocultas en encuestas

A veces se hacen encuestas sobre temas sobre los que los encuestados son reticentes a revelar la verdad (p.e., ¿es Vd. un zombi?). Un procedimiento conocido para recabar tal tipo de información es el siguiente:

Se le invita al encuestado a tirar al aire una moneda con las caras etiquetadas con sí y no; la moneda no es una moneda porque tiene una probabidad conocida (y distinta del 50%) de caer en sí.
El encuestado responde sí si la respuesta a la pregunta y el resultado de la tirada de la moneda coinciden y no en caso contrario.

A partir de la proporción de respuestas positivas y conocida la probabilidad del sí de la moneda, $q$ , es posible estimar la proporción $\theta$ de respuestas positivas a la pregunta de subyacente de interés en la muestra. Efectivamente, los síes tienen una distribución binomial $B(p) = B(q\theta + (1-q)(1-\theta))$ y, una vez estimado (por máxima verosimilitud) $\hat{p}$ , puede despejarse $\hat{p}$ de $\hat{p} = q\hat{\theta} + (1-q)(1-\hat{\theta})$ para obtener

$\hat{\theta} = \frac{1 - q - \hat{p}}{1 - 2q}.$

Así,

set.seed(12)
 
n <- 10000
unknown.par <- 0.2 
coin.par    <- 0.3
biased.coin <- rbinom(n, 1, coin.par)
 
hidden.res  <- rbinom(n, 1, unknown.par)
results <- as.numeric(hidden.res == biased.coin)
 
 
get.theta <- function(p, q) (1 - p - q) / (1 - 2*q)
get.theta(mean(results), coin.par)
# 0.19

Y todo es estupendo.

¿Y si queremos intervalos de confianza, etc., del parámetro estimado? Podemos muestrear una $B(\hat{p})$ y ver cuál sería la distribución resultante de $\theta$ :

muestras <- replicate(1000, mean(rbinom(n, 1, unknown.par) == biased.coin))
muestras <- sapply(muestras, function(x) get.theta(x, coin.par))
hist(muestras)
# etc.

Pero el procedimiento anterior tiene algunos caveats:

La función get.theta es algo inestable: nada garantiza siquiera que dé valores positivos; de hecho, en algunas pruebas, me ha dado valores negativos en algunas ocasiones.
Ni siquiera yo entiendo demasiado bien por qué, a la hora de hacer simulaciones, muestreo $B(\hat{p})$ : ¿por qué $\hat{p}$ y no otro valor? Porque sé que la estimación de $\hat{p}$ está sujeta a error, etc.
Hay gente que no sabría despejar $\hat{\theta}$

Afortunadamente, existe un procedimiento alternativo:

library(rstan)
 
standat <- list(N = n, pcoin = coin.par, x = sum(results))
 
stanmodelcode <- '
data {
  int<lower=1> N;
  real pcoin;
  int x;
}
 
parameters {
  real<lower=0, upper = 1> theta;
}
 
model {
  // prior (no informativa) 
  theta ~ beta(1,1);
 
  // verosimilitud
  x ~ binomial(N, theta*pcoin + (1-theta)*(1-pcoin));
}
'
 
fit <- stan(model_code = stanmodelcode, 
            data = standat, 
            iter=12000, warmup=2000, 
            chains=4, thin=10)
 
tmp <- as.data.frame(fit)
hist(tmp$theta)

Que genera

sin mayores complicaciones (ni teóricas ni prácticas). Además, como ventaja adicional, uno siempre puede incluir la información previamente conocida sobre la distribución de $\theta$ en el lugar correspondiente en el código anterior.