Overview

The Romeb package fits robust median-based Bayesian linear growth curve models for longitudinal outcomes. The main function, Romeb(), supports complete data as well as data with missing values under MCAR, MAR, and MNAR assumptions.

The current implementation fits a linear growth model without covariates in the growth model. Auxiliary variables, when supplied, are used only in the MNAR missingness model. This vignette introduces the basic data structure, the main function arguments, and typical model calls.

Model fitting is performed with rjags, so JAGS must be installed on the system before fitting models.

Installation

After installing the package, load it with:

library(Romeb)

If the package is installed from CRAN, use:

install.packages("Romeb")

Data format

Romeb() expects a matrix or data frame containing repeated outcome measurements and, optionally, auxiliary variables. Measurement times are supplied separately through the time argument.

In earlier versions, users had to place auxiliary variables in the first K columns and outcome variables immediately afterward. The updated interface allows outcome and auxiliary variables to be selected directly by column names or column indices.

The following simulated data set has three repeated outcome variables and one auxiliary variable:

set.seed(123)

n <- 100

dat <- data.frame(
  y0 = rnorm(n, mean = 0.0, sd = 1.0),
  y1 = rnorm(n, mean = 0.5, sd = 1.0),
  y2 = rnorm(n, mean = 1.0, sd = 1.0),
  aux1 = rnorm(n)
)

time <- c(0, 1, 2)

head(dat)
#>            y0         y1        y2       aux1
#> 1 -0.56047565 -0.2104066 3.1988103 -0.7152422
#> 2 -0.23017749  0.7568837 2.3124130 -0.7526890
#> 3  1.55870831  0.2533081 0.7348549 -0.9385387
#> 4  0.07050839  0.1524574 1.5431941 -1.0525133
#> 5  0.12928774 -0.4516186 0.5856601 -0.4371595
#> 6  1.71506499  0.4549723 0.5237531  0.3311792

The outcome columns are y0, y1, and y2. The time scores are 0, 1, and 2, so the latent slope is interpreted as the change in the fitted median per one-unit increase in time.

Complete-data model

For complete data, set Missing_Type = "no missing". The selected outcome variables must contain no missing values.

The following code shows a complete-data model call. The chunk is not evaluated when building the vignette, because it runs MCMC sampling through JAGS.

fit_complete <- Romeb(
  Missing_Type = "no missing",
  data = dat,
  time = time,
  seed = 123,
  outcome_vars = c("y0", "y1", "y2"),
  chain = 1,
  Niter = 6000,
  burnIn = 3000,
  n_adapt = 1000
)

print(fit_complete)

The returned object has class RomebResult and contains posterior summaries, diagnostic information, MCMC samples, a parameter map, and model metadata.

names(fit_complete)

fit_complete$parameter_map
fit_complete$model_info

MCAR and MAR models

For incomplete outcome data under MCAR or MAR assumptions, use Missing_Type = "MCAR" or Missing_Type = "MAR". In the current implementation, MCAR and MAR use the observed-data outcome model. No explicit missingness model is fitted unless Missing_Type = "MNAR".

The following example creates a copy of the data with missing values:

dat_miss <- dat

set.seed(456)
missing_id <- sample(seq_len(nrow(dat_miss)), size = 20)
dat_miss$y1[missing_id] <- NA

head(dat_miss)
#>            y0         y1        y2       aux1
#> 1 -0.56047565 -0.2104066 3.1988103 -0.7152422
#> 2 -0.23017749  0.7568837 2.3124130 -0.7526890
#> 3  1.55870831  0.2533081 0.7348549 -0.9385387
#> 4  0.07050839  0.1524574 1.5431941 -1.0525133
#> 5  0.12928774 -0.4516186 0.5856601 -0.4371595
#> 6  1.71506499  0.4549723 0.5237531  0.3311792

A MAR model call can be written as follows:

fit_mar <- Romeb(
  Missing_Type = "MAR",
  data = dat_miss,
  time = time,
  seed = 123,
  outcome_vars = c("y0", "y1", "y2"),
  chain = 1,
  Niter = 6000,
  burnIn = 3000,
  n_adapt = 1000
)

print(fit_mar)

Similarly, an MCAR model can be fitted with:

fit_mcar <- Romeb(
  Missing_Type = "MCAR",
  data = dat_miss,
  time = time,
  seed = 123,
  outcome_vars = c("y0", "y1", "y2"),
  chain = 1,
  Niter = 6000,
  burnIn = 3000,
  n_adapt = 1000
)

print(fit_mcar)

MNAR model without auxiliary variables

For MNAR missingness, set Missing_Type = "MNAR". Without auxiliary variables, the missingness model includes an intercept, a coefficient for the previous outcome, and a coefficient for the current outcome.

fit_mnar <- Romeb(
  Missing_Type = "MNAR",
  data = dat_miss,
  time = time,
  seed = 123,
  outcome_vars = c("y0", "y1", "y2"),
  chain = 1,
  Niter = 6000,
  burnIn = 3000,
  n_adapt = 1000
)

print(fit_mnar)

MNAR model with auxiliary variables

Auxiliary variables can be included in the MNAR missingness model through the auxiliary_vars argument. Auxiliary variables are currently supported only for Missing_Type = "MNAR".

fit_mnar_aux <- Romeb(
  Missing_Type = "MNAR",
  data = dat_miss,
  time = time,
  seed = 123,
  outcome_vars = c("y0", "y1", "y2"),
  auxiliary_vars = "aux1",
  chain = 1,
  Niter = 6000,
  burnIn = 3000,
  n_adapt = 1000
)

print(fit_mnar_aux)

The same variables can also be selected by column indices:

fit_mnar_aux_idx <- Romeb(
  Missing_Type = "MNAR",
  data = dat_miss,
  time = time,
  seed = 123,
  outcome_vars = 1:3,
  auxiliary_vars = 4,
  chain = 1,
  Niter = 6000,
  burnIn = 3000,
  n_adapt = 1000
)

Prior specification

The priors argument allows users to override the default weakly informative prior settings. It should be a named list. Unknown names trigger an error so that misspecified prior inputs are not silently ignored.

Supported prior elements include:

muLS_mean
muLS_prec
sigma_shape
sigma_rate
wishart_R
wishart_df
missing_coef_mean
missing_coef_prec
aux_coef_mean
aux_coef_prec

In JAGS, normal priors are parameterized by precision rather than variance.

The following example modifies the prior for the latent intercept and slope means and the residual precision parameter:

custom_priors <- list(
  muLS_mean = c(0, 0),
  muLS_prec = c(0.01, 0.01),
  sigma_shape = 0.01,
  sigma_rate = 0.01
)

fit_custom_priors <- Romeb(
  Missing_Type = "no missing",
  data = dat,
  time = time,
  seed = 123,
  outcome_vars = c("y0", "y1", "y2"),
  priors = custom_priors,
  chain = 1,
  Niter = 6000,
  burnIn = 3000,
  n_adapt = 1000
)

print(fit_custom_priors)

For an MNAR model with one auxiliary variable, auxiliary coefficient priors can be specified as follows:

custom_mnar_priors <- list(
  missing_coef_mean = c(0, 0, 0),
  missing_coef_prec = c(0.01, 0.01, 0.01),
  aux_coef_mean = c(0),
  aux_coef_prec = c(0.01)
)

fit_custom_mnar <- Romeb(
  Missing_Type = "MNAR",
  data = dat_miss,
  time = time,
  seed = 123,
  outcome_vars = c("y0", "y1", "y2"),
  auxiliary_vars = "aux1",
  priors = custom_mnar_priors,
  chain = 1,
  Niter = 6000,
  burnIn = 3000,
  n_adapt = 1000
)

Initial values

The inits argument allows users to supply initial values for JAGS. If inits is omitted, chain-specific RNG initial values are generated internally.

For a single chain, inits can be a named list:

fit_with_inits <- Romeb(
  Missing_Type = "no missing",
  data = dat,
  time = time,
  seed = 123,
  outcome_vars = c("y0", "y1", "y2"),
  inits = list(
    muLS = c(0, 0)
  ),
  chain = 1,
  Niter = 6000,
  burnIn = 3000,
  n_adapt = 1000
)

For multiple chains, inits can be a list of named lists with length equal to the number of chains:

fit_two_chains <- Romeb(
  Missing_Type = "no missing",
  data = dat,
  time = time,
  seed = 123,
  outcome_vars = c("y0", "y1", "y2"),
  inits = list(
    list(muLS = c(0, 0)),
    list(muLS = c(0.1, 0.1))
  ),
  chain = 2,
  Niter = 6000,
  burnIn = 3000,
  n_adapt = 1000
)

MCMC settings

The main MCMC-related arguments are:

chain: number of MCMC chains.
Niter: number of saved posterior sampling iterations per chain.
burnIn: number of saved iterations discarded before posterior summaries are computed.
n_adapt: number of JAGS adaptation iterations before posterior sampling.

The adaptation iterations controlled by n_adapt are not saved. They are separate from burnIn, which discards saved posterior samples during post-processing.

For real analyses, users should choose sufficiently large Niter, burnIn, and chain values and inspect convergence diagnostics.

Output structure

A fitted Romeb object contains:

names(fit_complete)

The main components are:

quantiles: posterior means, standard deviations, and quantiles.
geweke: Geweke diagnostic z-scores.
credible_intervals: equal-tail 95 percent credible intervals.
hpd_intervals: 95 percent highest posterior density intervals.
samps_full: full coda::mcmc.list object.
parameter_map: mapping between output names and model parameters.
model_info: model settings, selected variables, priors, and MCMC settings.

Users can access the full MCMC samples through samps_full:

coda::traceplot(fit_complete$samps_full[, "latent_intercept_mean"])
coda::summary(fit_complete$samps_full)

Legacy data layout

The legacy K argument is still supported. If outcome_vars is omitted, the function uses the first K columns as auxiliary variables and the following length(time) columns as outcome variables.

For new analyses, it is recommended to use outcome_vars and auxiliary_vars because they are clearer and do not require a restrictive column order.

fit_legacy <- Romeb(
  Missing_Type = "MNAR",
  data = dat_miss[, c("aux1", "y0", "y1", "y2")],
  time = time,
  seed = 123,
  K = 1,
  chain = 1,
  Niter = 6000,
  burnIn = 3000,
  n_adapt = 1000
)

Practical recommendations

For applied analyses:

Use outcome_vars to specify longitudinal outcomes explicitly.
Use auxiliary_vars only for MNAR models.
Check that length(time) equals the number of selected outcome variables.
Increase Niter, burnIn, and chain for final analyses.
Inspect trace plots, Geweke diagnostics, credible intervals, and HPD intervals.
Report the missingness assumption and prior settings clearly.

Session information

sessionInfo()
#> R version 4.5.1 (2025-06-13)
#> Platform: aarch64-apple-darwin20
#> Running under: macOS Sonoma 14.1
#> 
#> Matrix products: default
#> BLAS:   /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/lib/libRblas.0.dylib 
#> LAPACK: /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/lib/libRlapack.dylib;  LAPACK version 3.12.1
#> 
#> locale:
#> [1] C/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
#> 
#> time zone: America/New_York
#> tzcode source: internal
#> 
#> attached base packages:
#> [1] stats     graphics  grDevices utils     datasets  methods   base     
#> 
#> loaded via a namespace (and not attached):
#>  [1] digest_0.6.37     R6_2.6.1          fastmap_1.2.0     xfun_0.53        
#>  [5] cachem_1.1.0      knitr_1.50        htmltools_0.5.8.1 rmarkdown_2.29   
#>  [9] lifecycle_1.0.4   cli_3.6.5         sass_0.4.10       jquerylib_0.1.4  
#> [13] compiler_4.5.1    rstudioapi_0.17.1 tools_4.5.1       evaluate_1.0.5   
#> [17] bslib_0.9.0       yaml_2.3.10       rlang_1.1.6       jsonlite_2.0.0

Romeb Package: An Introduction

Dandan Tang

Xin Tong

2026-07-02