I'm no expert in this kind of analysis, but I will sketch out a sort of answer. My feeling is that doing five models isn't the most efficient way of doing this because you're not using the information that the plausible values come from the same distribution (for the same individual).
As I see it you have some pseudo-response values that are drawn from a distribution centered at the true response value. You know that this distribution is a normal distribution. It is fairly straight-forward to use this information in a Bayesian model.
The main idea is to treat the true response variable that you're interested in as an unknown parameter in the model you're estimating, and the plausible values are treated as data to inform on this parameter. I provide an example in R and Stan. You might be able to do this sort of thing with some of the friendlier packages like brms but I wouldn't know how.
First we generate some fake data
set.seed(2023-11-17)
# True model is y = 1 + x + noise
x = rnorm(50)
y = 1 + x + rnorm(50, sd=.05)
We observe some true $y$ values generated as above. For some reason (privacy perhaps?) we choose to provide some randomly drawn plausible values instead of the true value so that the plausible values are centered on the true value:
y_plausible = t(sapply(y, function(yy) { rnorm(5, mean=yy, sd=.025) }))
head(y_plausible)
#> [,1] [,2] [,3] [,4] [,5]
#> [1,] 2.1016310 2.0761254 2.1356810 2.1256719 2.0619211
#> [2,] 0.7671114 0.7383365 0.7434291 0.7278975 0.7367159
#> [3,] 1.0430924 1.0673116 1.1229851 1.0795531 1.0617507
#> [4,] 1.4613296 1.4723853 1.4703300 1.4365836 1.4870355
#> [5,] 2.7055201 2.6943164 2.7165702 2.7290942 2.6345115
#> [6,] 0.6336314 0.6393465 0.6268552 0.6292515 0.6270424
This is all that the analyst has to work with, the true $y$ being hidden.
Then formulate a Bayesian model for the data generating process
Below is the code from the file model.stan that I use further down. Basically it says that true $y$, an unknown parameter, is distributed as normal in a typical linear regression way. "Observed" data, $y_{plausible}$, is distributed as normal, centered on true $y$.
data {
int<lower=0> N;
vector[N] x;
vector[5] y_plausible[N];
}
// The parameters accepted by the model. Notice that true obscured y becomes
// a parameter to be estimated
parameters {
// regression parameters
real beta_0;
real beta_1;
real<lower=0> sigma;
// parameters of the distributions from which plausible values came.
// you can give them individual standard deviations if you want. I just do a
// common one to all of them, which probably isn't very realistic
vector[N] y;
real<lower=0> sigma_plausible;
}
// The model to be estimated.
model {
y ~ normal(beta_0 + beta_1*x, sigma);
for (i in 1:N) {
y_plausible[i] ~ normal(y[i], sigma_plausible);
}
// I just let stan set some default priors, which is probably not advisable in
// general but for this simple model it seems to work fine
}
Now we fit the model and look at some estimates
Our posterior distribution for the regression slope is more or less centered on the true value of 1:
library(rstan)
dataset = list(N=50, x=x, y_plausible=y_plausible)
stan_fit = stan("model.stan", data=dataset)
draws = extract(stan_fit)
# the intercept, truth marked as a vertical bar
hist(draws$beta_1, breaks = 50, prob=T, border = "lightgrey",
main="Draws from beta_1 posterior")
abline(v=1, lwd=1.5)
The posterior for the first true $y$ is also pretty close:
# Let's look at the first y value, truth marked as vertical bar
estimated_y = draws$y
hist(estimated_y[, 1], breaks = 50, prob=T, border = "lightgrey",
main="Draws from y[1] posterior")
abline(v=y[1], lwd=1.5)
Created on 2023-11-17 with reprex v2.0.2
