Cómo simular la distribución de una prueba de diagnóstico

Question

Tengo un parámetro descrito como mean (95% CI low-hi) y quisiera construir una distribución que tenga aproximadamente esa media y los límites inferior/superior. Las restricciones son que normalmente no son normales y todos los valores tienen que estar entre 0 y 1. ¿Cómo lo haría?

Un método de fuerza bruta que he pensado es el siguiente, pero me preguntaba si hay un adecuado o, al menos, más elegante manera de hacerlo.

Nota: Una versión anterior de esta pregunta mencionaba que estos parámetros son características de las pruebas de diagnóstico como la sensibilidad y la especificidad.

bMoms <- function(m, sd) {
  ## the moments of the beta distribution   
  v <- sd^2
  t <- (1 - m)/v
  a <- m^2 * t
  b <- m * t - a
  ## method of moments
  t. <- m * (1 - m)/v - 1
  a. <- m * t.
  b. <- (1 - m) * t.
  return(c(a, b, a., b.))
}
m = (lo + 4 * m + hi)/6
m. <- (m - lo)/(hi - lo)  #normalize m to scale 0-1
sd. <- 1/6
ab <- bMoms(m., sd.)  # the moments of the beta distribution
df. <- data.frame(x = rbeta(n, ab[1], ab[2]))
df. <- within(df., x <- lo + x * (hi - lo))

Answer 1

La simulación es un medio viable para inspeccionar las características de funcionamiento de una prueba de diagnóstico. Aunque los estudios preliminares presentan la sensibilidad y la especificidad con intervalos de confianza del 95%, estos IC no representan en realidad una probabilidad bayesiana en lo que respecta a dichos valores. Por lo tanto, cuando se proponga realizar su estudio de simulación, es mejor considerar las características de la prueba como fijas o dadas, y luego simular un proceso aleatorio del estado de la enfermedad y la evaluación de la prueba en función de esos parámetros. De este modo, podrá resumir adecuadamente el rendimiento de la prueba (y su incertidumbre) en diversos escenarios, como una elevada prevalencia de la enfermedad o en muestras más grandes. Los parámetros básicos son los siguientes:

\begin{array}{c|ccc} & D & \bar{D} & \\ \hline T & TP & FP& n_t\\ \bar{T} & FN & TN & n-n_t\\ & n_d & n-n_d \end{array}

Como ejemplo, supongamos que se ha desarrollado una prueba y se ha estudiado en una cohorte con una sensibilidad del 40% ( $TP/n_d$ ) y una especificidad del 90% ( $TN/(n-n_d)$ ). Supongamos además que consulto la literatura y descubro que la muestra de origen era de 100 pacientes con una prevalencia de la enfermedad del 10% (digamos, en un entorno hospitalario). Entonces, la tabla anterior se reduce a

\begin{array}{c|ccc} & D & \bar{D} & \\ \hline T & 4 & 9 & 13\\ \bar{T} & 6 & 81 & 87\\ & 10 & 90 \end{array}

Pero si queremos simular estos resultados, el $TN$ y el $FN$ son simplemente modelos de probabilidad binomial. Usando R puedo simular esta tabla de contingencia con:

sim <- function(n, pd, se, sp) {
  nd <- rbinom(1, n, pd)
  tp <- rbinom(1, nd, se)
  tn <- rbinom(1, n-nd, sp)
  matrix(c(tp, nd-tp, n-nd-tn, tn), 2, 2)
}
set.seed(123)
sim(n=100, pd=0.1, se=0.4, sp=0.9)

se pone:

> sim(n=100, pd=0.1, se=0.4, sp=0.9)
     [,1] [,2]
[1,]    4    8
[2,]    4   84

que muestra cierta variabilidad estocástica del diseño. Realización de muchas 1000s de réplicas:

OUT <- replicate(1000, sim(n=100, pd=0.1, se=0.4, sp=0.9))
sens <- OUT[1,1,] / colSums(OUT[,1,])
mean(sens) + c('lower 95%'=-1, 'mean'=0, 'upper 95%'=1) %o%  qnorm(0.975) * sd(sens)

Lo que da

> mean(sens) + c('lower 95%'=-1, 'mean'=0, 'upper 95%'=1) %o%  qnorm(0.975) * sd(sens) 
                [,1]
lower 95% 0.06350984
mean      0.40045478
upper 95% 0.73739973

Así que, como era de esperar, este diseño es ineficaz para evaluar la fiabilidad de las pruebas. No generamos variabilidad en la prueba, sino que simulamos el resultado de la prueba como un proceso aleatorio que se resume con precisión en la sensibilidad y la especificidad y que varía bajo la experimentación en función del diseño experimental.