¿Cómo incorporar la incertidumbre de los coeficientes del modelo en el intervalo de predicción de una regresión lineal múltiple?

Estoy lidiando con el modelado de pequeños conjuntos de datos de física experimental (específicamente la pegajosidad de los compuestos de pegamento). Como la mayoría del trabajo experimental no genera miles de muestras, sino un puñado, necesito ser inventivo en cómo manejar esta pequeña cantidad de conjuntos de datos (digamos 10-20). En este punto, tengo un modelo-marco (ver regresión a continuación en PSS) que puede manejar esto bastante bien.

Sin embargo, para tener una mejor idea de la precisión de mis predicciones, quiero tener una barra de error en mis valores predichos, esto para verificar qué tan bien mis predicciones predicen nuevos experimentos. Como este trabajo es de naturaleza numérica, la barra de error se originará en el modelo teórico subyacente, ¿cómo se propagan estos errores (es decir, el análisis de errores como se acostumbra en la física experimental)

En aras de la simplicidad, supongamos que estoy tratando con un modelo de regresión lineal múltiple, digamos (en realidad habrá muchos más términos):

$$y = \beta_0 + \beta_1 x_1 + \beta_2 x_2 \tag{1}$$ Lo que estoy buscando es una forma algebraica de calcular (numéricamente) las barras de error (en realidad, es el intervalo de predicción (PI) o el intervalo de confianza (CI), ya que ambos están relacionados). En la literatura estadística, hay referencias a tal problema y ejemplos de cómo se pueden calcular el PI y el IC. Sin embargo, estos sólo consideran la variabilidad de la$x$'s. El PI y el CI se relacionan entonces con (cf. pregunta 147242 ): $$\hat{V}_f=s^2\cdot\mathbf{x_0}\cdot\mathbf{(X^TX)^{-1}}\cdot\mathbf{x_0^T} + s^2 \tag{2}$$

En contraste con estos, cada uno de los coeficientes de mi modelo [ver: PSS a continuación] ($\beta_0, \beta_1$y$\beta_2$) en este caso tienen una barra de error (extraído mediante bootstrapping de una distribución, siendo las distribuciones de naturaleza numérica no analítica, y las distribuciones son específicas para cada uno de los tres coeficientes). ¿Hay alguna manera de incorporar la incertidumbre de la$\beta_i$'s (cq las "barras de error") en el cálculo del PI (y CI).

En pocas palabras, ¿cómo puede la ecuación

$$\hat{V}_f=s^2\cdot\mathbf{x_0}\cdot\mathbf{(X^TX)^{-1}}\cdot\mathbf{x_0^T} + s^2 \tag{3}$$ modificarse para incorporar también el hecho de que los coeficientes mismos son una media de una distribución.

(PD: Uno podría crear un conjunto de varias instancias modelo con el$\beta_i$extraídas de sus respectivas distribuciones, y en base a la distribución de las obtenidas$y_0$calcular el IC de la$y_0$, pero esto no es realmente eficiente desde el punto de vista computacional y trae muchos otros problemas que me gustaría evitar).

(PPS: El modelo de regresión presentado no es el resultado de una regresión directa hacia un solo conjunto de datos, sino que se construye de la siguiente manera:

1. Cree un conjunto de N conjuntos de datos.
2. En cada conjunto de datos, una regresión da lugar a un modelo lineal como se indica en la publicación anterior. Esto da lugar a N valores para cada uno de los coeficientes$\beta$.
3. Se calcula la media de cada uno de los tres conjuntos.
4. Estos tres coeficientes medios son los coeficientes del modelo presentado anteriormente.
5. El objetivo aquí: encontrar el intervalo de predicción para el modelo promediado anterior teniendo en cuenta el hecho de que los coeficientes$\beta$se calculan a partir de distribuciones numéricas).