---
title: "Romeb Package: An Introduction"
author:
  - Dandan Tang
  - Xin Tong
date: "`r Sys.Date()`"
output:
  rmarkdown::html_vignette:
    toc: true
    toc_depth: 2
vignette: >
  %\VignetteIndexEntry{Romeb Package: An Introduction}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  message = FALSE,
  warning = FALSE
)
```

# 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:

```{r load-package, eval = FALSE}
library(Romeb)
```

If the package is installed from CRAN, use:

```{r install-cran, eval = FALSE}
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:

```{r simulate-data}
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)
```

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.

```{r complete-data-fit, eval = FALSE}
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.

```{r complete-data-output, eval = FALSE}
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:

```{r simulate-missing-data}
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)
```

A MAR model call can be written as follows:

```{r mar-fit, eval = FALSE}
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:

```{r mcar-fit, eval = FALSE}
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.

```{r mnar-no-aux-fit, eval = FALSE}
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"`.

```{r mnar-aux-fit, eval = FALSE}
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:

```{r mnar-aux-fit-indices, eval = FALSE}
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:

```{r prior-example, eval = FALSE}
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:

```{r aux-prior-example, eval = FALSE}
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:

```{r inits-one-chain, eval = FALSE}
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:

```{r inits-multiple-chains, eval = FALSE}
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:

```{r output-structure, eval = FALSE}
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`:

```{r mcmc-samples, eval = FALSE}
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.

```{r legacy-layout, eval = FALSE}
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:

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

# Session information

```{r session-info}
sessionInfo()
```