Models with Uncertainty in Y

Question

Let's say I have data with predictors x1, x2, x3...xn for a variable y. I have essentially imputed y using a Bayesian analysis, which means I have a posterior distribution for each value of y. To calculate a follow up model predicting y from x1 to xn while keeping uncertainty in y, I'd normally use something like the brm_multiple function from the brms package, which allows for us to calculate the same model on different values of y from the posterior distribution and then combine them into the same model at the end. Unfortunately, I have so much data that a Bayesian approach is just not feasible.

How would I do something similar with frequentist statistics? Is there a way to run a bunch of models and then combine the models at the end, or possibly have each y as a distribution in the original model, etc.?

Here's R code to generate a toy dataset that takes the form of what I'm talking about:

# Set seed for reproducibility
set.seed(42)

# Create toy data with predictors
n <- 100
x1 <- rnorm(n)
x2 <- rnorm(n, mean = 3)
x3 <- rnorm(n, mean = -2)

# Simulate posterior samples for y (e.g., 5 samples from a posterior distribution)
y_post_1 <- 2 + 1.5*x1 - 0.7*x2 + 0.5*x3 + rnorm(n)
y_post_2 <- 2 + 1.5*x1 - 0.7*x2 + 0.5*x3 + rnorm(n)
y_post_3 <- 2 + 1.5*x1 - 0.7*x2 + 0.5*x3 + rnorm(n)
y_post_4 <- 2 + 1.5*x1 - 0.7*x2 + 0.5*x3 + rnorm(n)
y_post_5 <- 2 + 1.5*x1 - 0.7*x2 + 0.5*x3 + rnorm(n)

# Combine into a data frame
toy_data <- data.frame(x1, x2, x3, y_post_1, y_post_2, y_post_3, 
                       y_post_4, y_post_5)

# Display the first few rows of the dataset
head(toy_data)

Answer 1

If you just have an uncertainty estimate, one approach could be to weight linear regression by the inverse of standard deviation - calculate mean and sd for each imputed value, regress on the mean, and use the weights argument in lm() to weight by 1/sd.

I'd caution you though - summarizing posterior preditions this way loses some information, particularly any covariation within each set of imputed data. In your example, you are using the same parameter values to sample the posterior, but if you actually using a posterior from a bayesian regression, every imputed data set will be generated from a different posterior draw. So you might have a wide range of possible values for any y among multiple imputations, but stronger relationships of y with predictor variables within each different imputation. So really keeping it bayesian is the best thing if you can, especially if your earlier workflow is bayesian regression.