simmodelfitshelp: Simulate multiple datasets, fit and save model estimates
In chaddlewick/spr: Simulation-based Parameter Recovery

1	simmodelfits(variablespecification, effectspecifications, outcomespecification, analyticmodel, ntimes)

`variablespecification`	A list of 'observed' variables to create a design specification. Must include 'participant' as one variable.
`effectspecifications`	A list of effects. Effect names must be variable names from variablespecification prefixed with 'b'. Must specify level at which to apply as second argument.
`outcomespecification`	A list containing the outcome distribution and any parameters to be computed based on the columns contained in the simulated dataset. Note that all such parameters should be prefixed with 'dataset$' and then call the variable name
`analyticmodel`	A character object of the model specification to be fit to the simulated data. Must contain variable names as columns in the simulated data.
`ntimes`	Number of times to simulate the data and fit the model

Please note these parameter values are provided for illustrative purposes only

observedvariables = as.list(c(participant = "rep(1:20, each = 40)",
                              qriscore = "rnorm(participant, 10, 2)",
                              hlvascore = "rnorm(participant, 8, 0.5)",
                              texts = "rep(1:10, times = 20, each = 4)",
                              question = "rep(1:800)"))

effectvariables = as.list(c(intercept = "0.15",
                            bparticipant = "rnorm(participant, mean=0, sd=0.04)",
                            bqriscore = "rnorm(participant, 0.025, 0.001)",
                            bhlvascore = "rnorm(participant, 0.02, 0.001)",
                            btexts = "rnorm(texts, 0, 0.02)",
                            bquestion = "rnorm(question, 0, 0.015)"))

outcomegeneration = as.list(c(outcome= "rbinom(observation, 1, dataset$py)",
                              py = "dataset$intercept + dataset$bparticipant + dataset$bqriscore*dataset$qriscore + dataset$bhlvascore*dataset$hlvascore + dataset$btexts + dataset$bquestion"))

analyticmodel = "brm(outcome ~ (1|participant) + (1|texts) + qriscore + hlvascore, data=dataset, family = bernoulli(), cores = 2)"

simmodelfits(observedvariables, effectvariables, outcomegeneration, analyticmodel, 5)