¿Cómo se deciden los priores bayesianos en la vida real?

Question

Siempre he tenido la siguiente pregunta: ¿Cómo se deciden los priores bayesianos en la vida real?

He creado el siguiente escenario para plantear mi pregunta: Suponga que es usted investigador y está interesado en estudiar si se puede predecir la edad de una jirafa por su peso y altura (por ejemplo, modelo de regresión lineal: edad = b_o + b_1 altura + b_2 peso). Llegas a un parque nacional para registrar las mediciones de las jirafas, pero después de tomar las medidas de unas pocas jirafas, se produce una terrible tormenta y tienes que detener tu estudio. Sólo has tenido tiempo de medir 15 jirafas:

     weight   height age
1  2998.958 15.26611  53
2  3002.208 18.08711  52
3  3008.171 16.70896  49
4  3002.374 17.37032  55
5  3000.658 18.04860  50
6  3002.688 17.24797  45
7  3004.923 16.45360  47
8  2987.264 16.71712  47
9  3011.332 17.76626  50
10 2983.783 18.10337  42
11 3007.167 18.18355  50
12 3007.049 18.11375  53
13 3002.656 15.49990  42
14 2986.710 16.73089  47
15 2998.286 17.12075  52

Lamentablemente, esta información no es suficiente para completar su estudio. Sin embargo, investigas un poco y descubres que este tipo de mediciones se han realizado en jirafas en el pasado. Por ejemplo:

Estudio 1 : En 1800 se realizó un estudio en el que se midieron 1000 jirafas y se descubrió que la altura media de esas jirafas era de 17 pies, el peso medio era de 2800 libras y la edad media era de 35 años. Sin embargo, esto se hizo en el siglo XIX y es dudoso que las mediciones no hayan sido tan precisas en ese entonces, y los problemas en el medio ambiente (por ejemplo, la caza furtiva) podrían haber causado que las jirafas cambiaran de tamaño.

Estudio 2 : Un estudio se hizo en 2010 fueron 50 jirafas en los zoológicos de todo el mundo y su altura era de 16 pies, el peso era de 300 libras y la edad era de 50 años. Este estudio es más reciente, pero es escéptico que las jirafas en los zoológicos podrían ser diferentes de las jirafas en la naturaleza.

Estudio 3 : Un experto en jirafas cree firmemente que la edad, la altura y el peso de las jirafas tienen distribuciones en forma de campana. El experto también cree que las jirafas siguen creciendo durante toda su vida, es decir, que a medida que aumenta la edad, también lo hacen el peso y la altura. No tiene cifras concretas, pero se le considera el principal experto.

etc.

Pregunta: En este problema, ¿es posible complementar las limitadas mediciones que se tienen, junto con los conocimientos previos disponibles sobre las jirafas (teniendo en cuenta su fiabilidad)? ¿Es este problema un ejemplo de cómo se pueden utilizar los modelos bayesianos (por ejemplo, la regresión bayesiana) en la vida real, o este problema carece fundamentalmente de suficientes datos con los que trabajar?

Supongamos que se consultan varios estudios en los que se registran las alturas y se evalúa manualmente la credibilidad de estos estudios (asignando "pesos bajos" a los estudios considerados no creíbles, por ejemplo altura_ajustada = puntuación_de_credibilidad * altura_media_registrada_en_el_estudio ):

head(my_data)

 average_recorded_height_in_study credibility_score study_number adjusted_height
1                         13.74253         1.0000000            1       13.742525
2                         20.08053         0.3222523            2        6.470999
3                         13.25037         0.5132335            3        6.800532
4                         15.74946         0.2625349            4        4.134783
5                         11.68657         0.5966327            5        6.972592
6                         17.27276         1.0000000            6       17.272759

Hay muchas herramientas/paquetes (por ejemplo, utilizando el lenguaje de programación R) que pueden intentar explorar esta "información previa" y ajustar la distribución

library(fitdistrplus)
library(patchwork)
library(ggplot2)

 fg <- fitdist(my_data$adjusted_height, "gamma")
 fln <- fitdist(my_data$adjusted_height, "lnorm")
fg <- fitdist(my_data$adjusted_height, "gamma")
 fw <- fitdist(my_data$adjusted_height, "weibull")

 par(mfrow = c(2, 2))
 plot.legend <- c("Weibull", "lognormal", "gamma")

a <- denscomp(list(fw, fln, fg), legendtext = plot.legend, plotstyle = "ggplot")
b <- qqcomp(list(fw, fln, fg), legendtext = plot.legend, plotstyle = "ggplot")
c <- cdfcomp(list(fw, fln, fg), legendtext = plot.legend, plotstyle = "ggplot")
d <- ppcomp(list(fw, fln, fg), legendtext = plot.legend, plotstyle = "ggplot")

a+b+c+d

El análisis anterior podría repetirse para las demás variables del estudio. En este caso, podríamos ver qué distribución se ajusta mejor a los datos (por ejemplo, utilizando la - verosimilitud), y registrar las estimaciones de los parámetros de esta distribución.

¿Es ésta la idea correcta de cómo se incorporan los priores a los modelos bayesianos en el mundo real? En este ejemplo que he creado, ¿se puede analizar la información de estudios anteriores y utilizarla para crear priores para una Regresión Lineal Bayesiana?

Gracias

Nota: Supongamos que las 15 jirafas que ha medido resultan ser jirafas enfermas y que sus medidas de altura/peso no son representativas de la población general de jirafas, pero quizás la información codificada en los datos a priori representa una amplia gama de jirafas. Por lo tanto, la combinación de sus mediciones con la información a priori podría dar lugar a un modelo más realista que podría generalizar a una población más grande de jirafas (este hecho es desconocido para usted en este momento).

Answer 1

Hay dos grandes enfoques para este problema. En primer lugar, utilizar los datos relevantes del pasado para crear de alguna manera "automática" una prioridad (o incluir de alguna manera estos datos relevantes en un modelo único con nuestros nuevos datos). Esta opción suele considerarse atractiva, porque "tiene cierta objetividad". En segundo lugar, preguntar a los expertos (después de mostrarles los datos relevantes que puedan tener en cuenta). Por último, pero quizás menos relevante, está la opción de utilizar priores débilmente informativos (o priores que intentan ser poco informativos).

En la primera clase de enfoques, el (robusta) de predicción meta-analítica (MAP) de Schmidli et al. ya se ha mencionado y se utiliza con bastante frecuencia - especialmente en la versión robusta con un componente extra de mezcla débil/no informativa añadido-, pero hay varios variantes , alternativas como priores de potencia adaptativos , ideas para ajustar un único modelo a los datos antiguos y a los nuevos de una manera robusta a los conflictos de datos anteriores y otras ideas similares.

En la segunda clase de enfoques, hay muchas formas de obtener opiniones previas de los expertos de manera que se minimicen los sesgos a los que están sujetas las personas (incluidos los expertos) (= "elicitación de expertos" ). Uno de estos marcos es SHELF sobre el que puedes encontrar un curso completo en su página web y para el que también hay un Paquete R . Menciono ese específicamente, porque lo uso en la práctica, pero hay otros con diferentes sabores/filosofías.

A continuación se presentan algunos ejemplos de fijación de valores a priori en la práctica, en su mayoría extraídos de ensayos clínicos/desarrollo de fármacos (simplemente porque soy el que está más familiarizado con ello; para más ejemplos, véase por ejemplo este libro ): para un estudio de prueba de concepto en EPOC para un prueba de concepto en artritis reumatoide (y otro también para la AR ), para un riesgo exponencial a partir de datos históricos , para efectos del tratamiento en los ensayos clínicos y para predecir tasas de eventos y parámetro de dispersión para resultados de conteo . En la industria farmacéutica, el uso de la información previa y el conocimiento de los expertos es especialmente común para analizar los estudios en las primeras etapas del desarrollo clínico (por ejemplo, el análisis de los estudios de prueba de concepto y la decisión de seguir adelante) o para la toma de decisiones posterior Mientras que es más raro en el caso de los estudios confirmatorios destinados a respaldar la aprobación reglamentaria (en parte, una prioridad demasiado optimista es más un problema para la empresa cuando se trata de la toma de decisiones interna, mientras que las autoridades reglamentarias someten las prioridades elegidas para los estudios confirmatorios a un escrutinio mucho mayor).

Answer 2

OP aquí, sólo quería añadir algo de material complementario y demostrar lo siguiente: una comparación entre Regresión frecuencial y Regresión bayesiana utilizando R

#cool trick to directly bring this data into R

my_data <- data.frame(read.table(header=TRUE,
row.names = 1,
text="
                         weight   height age
                      1  2998.958 15.26611  53
                      2  3002.208 18.08711  52
                      3  3008.171 16.70896  49
                      4  3002.374 17.37032  55
                      5  3000.658 18.04860  50
                      6  3002.688 17.24797  45
                      7  3004.923 16.45360  47
                      8  2987.264 16.71712  47
                      9  3011.332 17.76626  50
                      10 2983.783 18.10337  42
                      11 3007.167 18.18355  50
                      12 3007.049 18.11375  53
                      13 3002.656 15.49990  42
                      14 2986.710 16.73089  47
                      15 2998.286 17.12075  52
"))

1. Regresión frecuencial : Así es como funciona un modelo de regresión frecuentista (es decir, un modelo de regresión en el que los parámetros se estiman mediante mínimos cuadrados ordinarios (MCO), lo que todos aprendemos en la escuela).

En primer lugar, ajuste el modelo de regresión:

#fit regression model

model_1 <- lm(age ~ weight + height, data = my_data)

A continuación, vea los resultados:

#view results

 summary(model_1)

Call:
lm(formula = age ~ weight + height, data = my_data)

Residuals:
    Min      1Q  Median      3Q     Max 
-6.2369 -1.8688  0.3864  2.1065  5.6170 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)
(Intercept) -525.2843   369.9144  -1.420    0.181
weight         0.1875     0.1238   1.515    0.156
height         0.6871     1.0859   0.633    0.539

Residual standard error: 3.796 on 12 degrees of freedom
Multiple R-squared:  0.1954,    Adjusted R-squared:  0.06135 
F-statistic: 1.457 on 2 and 12 DF,  p-value: 0.2712

Opcional : Visualizar los resultados

library(scatterplot3d)

s3d <- scatterplot3d(my_data$weight, my_data$height,my_data$age, pch = 19, type = c("p"), color = "darkgrey",
                     main = "Regression Plane", grid = TRUE, box = FALSE,  
                     mar = c(2.5, 2.5, 2, 1.5), angle = 55)

# regression plane
s3d$plane3d(model_1, draw_polygon = TRUE, draw_lines = TRUE, 
            polygon_args = list(col = rgb(.1, .2, .7, .5)))

# overlay positive residuals
wh <- resid(model_1) > 0
s3d$points3d(my_data$height, my_data$weight, my_data$age, pch = 19)

2) Regresión bayesiana: Ahora, intentamos ajustar un modelo de regresión bayesiano a los mismos datos:

#load library

library(rstanarm)
library(see)
library(bayestestR)
library(performance)

En primer lugar, especificamos las variables de Altura y Peso (elegí una distribución normal para ambas - en mi pregunta original, habríamos decidido estas priores utilizando la investigación realizada en jirafas por otros biólogos):

#specify priors
my_prior <- normal(location = c(3000, 17), scale = c(1, 2))

A continuación, ejecutamos el modelo de regresión bayesiano

#run bayesian regression model
model_2 <- stan_glm(age~., data=my_data, prior = my_prior,    seed=111)

Ahora, podemos ver los resultados:

 summary(model_2)

Model Info:
 function:     stan_glm
 family:       gaussian [identity]
 formula:      age ~ .
 algorithm:    sampling
 sample:       4000 (posterior sample size)
 priors:       see help('prior_summary')
 observations: 15
 predictors:   3

Estimates:
              mean       sd         10%        50%        90%     
(Intercept) -9000290.7     3116.3 -9004290.9 -9000230.6 -8996293.9
weight          2999.7        1.0     2998.4     2999.7     3001.1
height            17.0        2.0       14.4       17.0       19.6
sigma           3207.5       65.0     3124.2     3207.2     3291.0

Fit Diagnostics:
           mean    sd      10%     50%     90%  
mean_PPD    55.5   824.4 -1002.3    66.1  1107.1

Mira el rendimiento del modelo:

#model performance 

 performance(model_2)

# Indices of model performance

ELPD     | ELPD_SE |    LOOIC | LOOIC_SE |     WAIC |    R2 | R2 (adj.) |      RMSE |    Sigma
----------------------------------------------------------------------------------------------
-574.459 | 154.366 | 1148.918 |  308.733 | 1160.324 | 0.983 |    -1.000 | 23876.735 | 3207.163
> se <- sqrt(diag(vcov(model_2)))
> se
    (Intercept)      weight      height 
3116.342642    1.038384    2.040471 

**Opcional: Visualizar los resultados**

#MCMC Trace

x <- as.array(model_2, pars = c("(Intercept)", "height", "weight"))
bayesplot::mcmc_trace(x, facet_args = list(nrow = 2))

#Posterior Distributions

plot_title <- ggplot2::ggtitle("Posterior Distributions")
plot(model_2, "hist", "weight", "height") + plot_title

#confidence ellipse
bayesplot::color_scheme_set("green")
plot(model_2, "scatter", pars = c("height", "weight"),
     size = 3, alpha = 0.5) +
    ggplot2::stat_ellipse(level = 0.9) 

Referencias :

• https://rpubs.com/Qsheep/BayesianLinearRegression
• https://www.theoj.org/joss-papers/joss.01541/10.21105.joss.01541.pdf
• https://cran.r-project.org/web/packages/rstanarm/vignettes/priors.html#default-priors-and-scale-adjustments
• https://mc-stan.org/rstanarm/reference/plot.stanreg.html

Nota: Todavía estoy aprendiendo sobre la Regresión Bayesiana - por favor, siéntase libre de corregir cualquier error que pueda haber cometido (por ejemplo, parece que el Modelo de Regresión Bayesiana está funcionando mucho peor que el Modelo de Regresión Lineal debido a mi elección de priores). Cuando ejecuto el Modelo de Regresión Bayesiana con las priores por defecto ("priores débilmente informativos"), por ejemplo model_2 <- stan_glm(age~., data=my_data, seed=111) - los resultados de la regresión lineal bayesiana son comparables con el modelo de regresión lineal. Debo estar haciendo algo mal).

Gracias.