Run the longitudinal hierarchical model with MCMC.

```
hbl_mcmc_hierarchical(
data,
response = "response",
study = "study",
study_reference = max(data[[study]]),
group = "group",
group_reference = min(data[[group]]),
patient = "patient",
rep = "rep",
rep_reference = min(data[[rep]]),
covariates = grep("^covariate", colnames(data), value = TRUE),
constraint = FALSE,
s_delta = 30,
s_beta = 30,
s_sigma = 30,
s_lambda = 1,
s_mu = 30,
s_tau = 30,
covariance_current = "unstructured",
covariance_historical = "unstructured",
control = list(max_treedepth = 17, adapt_delta = 0.99),
...
)
```

- data
Tidy data frame with one row per patient per rep, indicator columns for the response variable, study, group, patient, rep, and covariates. All columns must be atomic vectors (e.g. not lists).

- response
Character of length 1, name of the column in

`data`

with the response/outcome variable.`data[[response]]`

must be a continuous variable, and it*should*be the change from baseline of a clinical endpoint of interest, as opposed to just the raw response. Treatment differences are computed directly from this scale, please supply change from baseline unless you are absolutely certain that treatment differences computed directly from this quantity are clinically meaningful.- study
Character of length 1, name of the column in

`data`

with the study ID.- study_reference
Atomic of length 1, element of the

`study`

column that indicates the current study. (The other studies are historical studies.)- group
Character of length 1, name of the column in

`data`

with the group ID.- group_reference
Atomic of length 1, element of the

`group`

column that indicates the control group. (The other groups may be treatment groups.)- patient
Character of length 1, name of the column in

`data`

with the patient ID.- rep
Character of length 1, name of the column in

`data`

with the rep ID.- rep_reference
Atomic of length 1, element of the

`rep`

column that indicates baseline, i.e. the first rep chronologically. (The other reps may be post-baseline study visits or time points.)- covariates
Character vector of column names in

`data`

with the columns with baseline covariates. These can be continuous, categorical, or binary. Regardless,`historicalborrowlong`

derives the appropriate model matrix.Each baseline covariate column must truly be a

*baseline*covariate: elements must be equal for all time points within each patient (after the steps in the "Data processing" section). In other words, covariates must not be time-varying.A large number of covariates, or a large number of levels in a categorical covariate, can severely slow down the computation. Please consider carefully if you really need to include such complicated baseline covariates.

- constraint
Logical of length 1, whether to pool all study arms at baseline (first rep). Appropriate when the response is the raw response (as opposed to change from baseline) and the first rep (i.e. time point) is prior to treatment.

- s_delta
Numeric of length 1, prior standard deviation of the study-by-group effect parameters

`delta`

.- s_beta
Numeric of length 1, prior standard deviation of the fixed effects

`beta`

.- s_sigma
Numeric of length 1, prior upper bound of the residual standard deviations.

- s_lambda
shape parameter of the LKJ priors on the unstructured correlation matrices.

- s_mu
Numeric of length 1, prior standard deviation of

`mu`

.- s_tau
Numeric of length 1, Upper bound on

`tau`

.- covariance_current
Character of length 1, covariance structure of the current study. Possible values are

`"unstructured"`

for fully parameterized covariance matrices,`"ar1"`

for AR(1) covariance matrices, and`"diagonal"`

for residuals independent across time within each patient. In MCMC (e.g.`hbl_mcmc_hierarchical()`

), the covariance structure affects computational speed. Unstructured covariance is slower than AR(1), and AR(1) is slower than diagonal. This is particularly true for`covariance_historical`

if there are many historical studies in the data.- covariance_historical
Same as

`covariance_current`

, but for the covariance structure of each separate historical study. Each historical study has its own separate covariance matrix.- control
A named

`list`

of parameters to control the sampler's behavior. It defaults to`NULL`

so all the default values are used. First, the following are adaptation parameters for sampling algorithms. These are parameters used in Stan with similar names here.`adapt_engaged`

(`logical`

)`adapt_gamma`

(`double`

, positive, defaults to 0.05)`adapt_delta`

(`double`

, between 0 and 1, defaults to 0.8)`adapt_kappa`

(`double`

, positive, defaults to 0.75)`adapt_t0`

(`double`

, positive, defaults to 10)`adapt_init_buffer`

(`integer`

, positive, defaults to 75)`adapt_term_buffer`

(`integer`

, positive, defaults to 50)`adapt_window`

(`integer`

, positive, defaults to 25)

In addition, algorithm HMC (called 'static HMC' in Stan) and NUTS share the following parameters:

`stepsize`

(`double`

, positive)`stepsize_jitter`

(`double`

, [0,1])`metric`

(`string`

, one of "unit_e", "diag_e", "dense_e")

For algorithm NUTS, we can also set:

`max_treedepth`

(`integer`

, positive)

For algorithm HMC, we can also set:

`int_time`

(`double`

, positive)

For

`test_grad`

mode, the following parameters can be set:`epsilon`

(`double`

, defaults to 1e-6)`error`

(`double`

, defaults to 1e-6)

- ...
Other optional parameters:

`chain_id`

(`integer`

)`init_r`

(`double`

, positive)`test_grad`

(`logical`

)`append_samples`

(`logical`

)`refresh`

(`integer`

)`save_warmup`

(`logical`

)deprecated:

`enable_random_init`

(`logical`

)

`chain_id`

can be a vector to specify the chain_id for all chains or an integer. For the former case, they should be unique. For the latter, the sequence of integers starting from the given`chain_id`

are used for all chains.`init_r`

is used only for generating random initial values, specifically when`init="random"`

or not all parameters are initialized in the user-supplied list or function. If specified, the initial values are simulated uniformly from interval [-`init_r`

,`init_r`

] rather than using the default interval (see the manual of (cmd)Stan).`test_grad`

(`logical`

). If`test_grad=TRUE`

, Stan will not do any sampling. Instead, the gradient calculation is tested and printed out and the fitted`stanfit`

object is in test gradient mode. By default, it is`FALSE`

.`append_samples`

(`logical`

). Only relevant if`sample_file`

is specified*and*is an existing file. In that case, setting`append_samples=TRUE`

will append the samples to the existing file rather than overwriting the contents of the file.`refresh`

(`integer`

) can be used to control how often the progress of the sampling is reported (i.e. show the progress every`refresh`

iterations). By default,`refresh = max(iter/10, 1)`

. The progress indicator is turned off if`refresh <= 0`

.Deprecated:

`enable_random_init`

(`logical`

) being`TRUE`

enables specifying initial values randomly when the initial values are not fully specified from the user.`save_warmup`

(`logical`

) indicates whether to save draws during the warmup phase and defaults to`TRUE`

. Some memory related problems can be avoided by setting it to`FALSE`

, but some diagnostics are more limited if the warmup draws are not stored.

A tidy data frame of parameter samples from the
posterior distribution. Columns `.chain`

, `.iteration`

,
and `.draw`

have the meanings documented in the
`posterior`

package.

Before running the MCMC, dataset is pre-processed. This includes expanding the rows of the data so every rep of every patient gets an explicit row. So if your original data has irregular rep IDs, e.g. unscheduled visits in a clinical trial that few patients attend, please remove them before the analysis. Only the most common rep IDs should be added.

After expanding the rows, the function fills in missing values for every column except the response. That includes covariates. Missing covariate values are filled in, first with last observation carried forward, then with last observation carried backward. If there are still missing values after this process, the program throws an informative error.

Other mcmc:
`hbl_convergence()`

,
`hbl_mcmc_independent()`

,
`hbl_mcmc_pool()`

,
`hbl_mcmc_sge()`

```
if (!identical(Sys.getenv("HBL_TEST", unset = ""), "")) {
set.seed(0)
data <- hbl_sim_hierarchical(
n_study = 2,
n_group = 2,
n_patient = 5,
n_rep = 3
)$data
tmp <- utils::capture.output(
suppressWarnings(
mcmc <- hbl_mcmc_hierarchical(
data,
chains = 1,
warmup = 10,
iter = 20,
seed = 0
)
)
)
mcmc
}
```