Simulating Observations from the Prior Distribution?

Can you directly simulate data from Bayesian Priors and "add them as new rows" to your existing dataset - and then fit a statistical model to this new dataset?

In the context of statistical modelling, I could see this as having two advantages:

1) Simplify computations and avoid the need of MCMC sampling

2) Better exploit real world knowledge about your variables (Bayesian Regression requires you to specify prior distributions on the model parameters (i.e. regression beta coefficients) - when in real life, you are more likely to have knowledge about the prior distributions of the actual variables, and not the model parameters).

For example, suppose you have the following data about the height, weight and age of giraffes. Suppose you are interested in creating a regression model that predicts "age" as a function of "height and weight". Suppose I have the following measurements:

     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


Now, suppose a giraffe expert tells me the following information that I think can be used as priors:

• The average age of a giraffe is 50 years, with a "bell curve shape"
• The average height of a giraffe is 16 feet, with a "bell curve shape"
• The average weight of a giraffe is 2500 lbs, with a "bell curve shape"

Could I simulate data from this (e.g. suppose I assume that each of these variables has a normal distribution and (for simplicity sake) is not correlated to the other variables), and add them to the original data? For example (using the R programming language):

age = rnorm(10, 50, 1)
height = rnorm(10, 16,1)
weight = rnorm(10, 2500, 100)

bayesian_data = data.frame(age, height, weight)
bayesian_data$source = "simulated" old_data$source = "observed"

combined_data = rbind(old_data, bayesian_data)


The new data ("combined_data") would look as follows:

      weight   height      age    source
2998.958 15.26611 53.00000  observed
3002.208 18.08711 52.00000  observed
3008.171 16.70896 49.00000  observed
3002.374 17.37032 55.00000  observed
3000.658 18.04860 50.00000  observed
3002.688 17.24797 45.00000  observed
3004.923 16.45360 47.00000  observed
2987.264 16.71712 47.00000  observed
3011.332 17.76626 50.00000  observed
2983.783 18.10337 42.00000  observed
3007.167 18.18355 50.00000  observed
3007.049 18.11375 53.00000  observed
3002.656 15.49990 42.00000  observed
2986.710 16.73089 47.00000  observed
2998.286 17.12075 52.00000  observed
2563.850 15.93615 49.55841 simulated
2646.080 15.85238 49.71639 simulated
2336.807 16.24114 49.35653 simulated
2467.514 15.67480 49.67832 simulated
2501.552 16.55515 48.07730 simulated
2582.426 16.33497 50.51050 simulated
2549.077 16.20062 48.55126 simulated
2595.276 16.80354 50.72527 simulated
2429.517 16.90839 49.48397 simulated
2497.516 14.25892 49.33416 simulated


Then, I could directly fit a regression model to this data that would "Bayesian in Spirit":

model <- lm(age ~ height + weight, data = combined_data)

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

Residuals:
Min     1Q Median     3Q    Max
-7.570 -1.602  0.416  1.091  5.879

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 45.512105  10.789639   4.218 0.000354 ***
height       0.560421   0.742336   0.755 0.458290
weight      -0.002040   0.003024  -0.675 0.506924
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 3.132 on 22 degrees of freedom
Multiple R-squared:  0.02989,   Adjusted R-squared:  -0.0583
F-statistic: 0.3389 on 2 and 22 DF,  p-value: 0.7162


My Question: Does the approach I used make sense?

• I do not think what I have done is mathematically correct - I am not sure if you can just simulate data from your priors and add it to your original data.

• The approach I listed would require you to assume priors on all your variables, even if you only believed that you had reasonable prior information on only some of your variables. In a traditional Bayesian setting, you can just place "uninformative priors" on variables you did not have reliable information about, and place "informative priors" on the variables that you did have reliable information about (I suppose for any given row, you could just randomly sample your existing data and re-use existing values for some variables and simulate values for the other variables).

• The approach I listed partly defeats the purpose of Bayesian Modelling - my approach does not provide the "much desired" posterior distributions and credible intervals of the parameters, as is traditionally associated with Bayesian Models. Although I suppose you re-run the simulation 1000 times, create 1000 new datasets where the results of each simulation is appended to the original data - then you can fit a linear regression model to each of the 1000 datasets, and plot a histogram for the estimates of each parameter (i.e. the regression beta coefficients) ... but this will probably become computationally expensive and defeat the original purpose of this approach (i.e. the original purpose of this approach was to reduce computation time by avoiding MCMC and directly simulating from the priors).

• And finally, in the above example, I had 15 observations and simulated 10 observations from the priors. But why 10? Would 11 have been a better choice - why not 8? In the traditional Bayesian setting, you are not faced with this decision.

I am highly skeptical myself if what I have done is correct (as in the reasons I have just listed) - but can someone please provide some additional comments?

Thanks!