simmr posterior_predictive

[Craigdux]

@andrewcparnell --another question! I finally have my complete dataset (18N and 15O isotopes for groundwater samples). My output for the posterior_predictive for a first trial with a subset of wells results in this:

print(post_pred) $table interval.1 interval.2 data outside 1 1.366803 4.299307 1.3 TRUE 2 5.259295 8.197375 6.0 FALSE

$prop_outside [1] 0.5

I need to dig into the model output, but as a first cut, is this telling me half (50%) of the data fall outside the model results?

Thank you.

[andrewcparnell]

Hi @Craigdux. Thanks for the query. The idea of this table is that, if the model is doing a good job capturing the uncertainty in the data, then approximately 50% of the observations should lie outside the 50% prediction intervals. For your data this is exactly what you're getting so this would give you confidence that the model is fitting well.

However, you do only have two observations so it's not a particularly robust statistic. I would recommend only really using this tool if you have at least 10 observations, and preferably more.

[Craigdux]

Hi @andrewcparnell --thanks again. for the clarification. I must be setting my model up incorrectly--I have eight wells with, on average, about 20 observations (two isotopes) per well. I even put "group = 8" in my code?:

Where am I going wrong? thanks

Setup the model

well_simmr <- simmr_load(mixtures = as.matrix(targets[, 1:2]), source_names = sources$Sources, source_means = as.matrix(sources[,2:3]), source_sds = as.matrix(sources[,4:5]),

correction_means = as.matrix(tefs[,2:3]),

                    #correction_sds = as.matrix(tefs[,4:5]),
                    #concentration_means = as.matrix(concep[,2:3]),
                    group = targets$well)

##################################

Run the simmr model:

well_simmr_out <- simmr_mcmc(well_simmr)

summary(well_simmr_out, type = "quantiles", group = 8)

############################

Check model:

posterior_predictive(well_simmr_out, group = 8)

post_pred <- posterior_predictive(well_simmr_out) print(post_pred)

[andrewcparnell]

Hmm. That's weird. What's the output of:

print(well_simmr)

and what does the plot look like when you run posterior_predictive? If it only shows 2 observations then something is missing when you load in the data.

[Craigdux]

@andrewcparnell --I am sorry, I had just copied the first and last few lines. And I found every well ("group!"). And the Gelman diagnostics looked favorable. And my results match our hypotheses! Thank you again!

[andrewcparnell]

Not worries. I'm glad it's all working.