---
title: "Getting started with birtRcpp"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Getting started with birtRcpp}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

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

`birtRcpp` fits joint models of response accuracy and response time. The main estimator is
`ecm()`: maximum marginal likelihood by the ECM algorithm (Meng & Rubin, 1993), with standard
errors, AIC and BIC, usually in about a second. `gibbs()` samples the posterior of the same
model, starting from the ECM solution. Both are written in C++. The package ports the Julia
package ExtendedRtIrtModeling.jl, with corrected full conditionals. Notation: j indexes
persons, i items.

The model is named by `model =` in `ecm()` and `gibbs()`:

| `model` | model |
|---|---|
| `"irt"` | 2PL (or 1PL, `itemtype = "1pl"`) with a latent regression of ability |
| `"null"` | 2PL + log-normal response times (van der Linden, 2007) |
| `"latreg"` | ... with latent regressions of ability and speed |
| `"latent"` | speed regressed on ability and covariates; `quantile =` (Gibbs) for quantile regression |
| `"cross"` | cross-relations of ability with each log response time; `quantile =` (Gibbs) |

## Data

`input_data()` takes responses, response times (in seconds) and centred covariates, either
as matrices or as column names of a data frame:

```{r}
library(birtRcpp)
set.seed(1)
cond <- set_cond(n_subj = 300, n_item = 8, n_feat = 2)
sim <- sim_data(cond, sim_para(cond, "latent"), "latent")
raw <- data.frame(sim$Y, sim$T, sim$X)
names(raw) <- c(paste0("y", 1:8), paste0("t", 1:8), "x1", "x2")
dat <- input_data(resp = paste0("y", 1:8), time = paste0("t", 1:8), cov = c("x1", "x2"), data = raw)
dat
```

Data from an experiment program often come in long format (one row per person and item),
with times in milliseconds, text responses and categorical covariates. `input_data()` refuses
long data and text responses with the code to fix them; the steps are:

```r
wide <- reshape(long[c("pid", "item", "correct", "rt_ms")], idvar = "pid", timevar = "item",
                direction = "wide")                       # one row per person
wide <- merge(wide, unique(long[c("pid", "group")]), by = "pid")   # person-level columns back
items <- paste0("correct.q", 1:8); times <- paste0("rt_ms.q", 1:8)   # names made by reshape()
wide[times] <- wide[times] / 1000                         # milliseconds -> seconds (only lambda changes)
wide <- wide[stats::complete.cases(wide[c(items, times)]), ]   # complete data only
groups <- stats::model.matrix(~ group, wide)[, -1]        # k - 1 dummies for k groups
wide[colnames(groups)] <- scale(groups, scale = FALSE)    # centred
dat <- input_data(resp = items, time = times, cov = colnames(groups), id = "pid", data = wide)
```

Section 11 of the intro tutorial runs this from a CSV file to a table of the regression
coefficients and person scores merged back into the data.

## Maximum likelihood with ECM

ECM integrates ability over Gauss-Hermite nodes and speed in closed form. Each M step is
split into conditional maximizations, one per block of parameters, and each has a closed
form: the logistic item block through the Polya-Gamma bound (Polson, Scott & Windle, 2013),
the response-time and structural blocks as normal models. The log-likelihood never
decreases; SQUAREM accelerates it, and the standard errors come from the observed
information.

```{r}
fit <- ecm(dat, model = "latent")
fit
estimates(fit, pars = "beta|b_ability") # est, se, z, Wald 95% interval
head(scores(fit)) # EAP and posterior SD of ability and speed
fit_indices(fit) # logLik, npar, N, AIC, BIC (also logLik(fit), AIC(fit), BIC(fit))
```

Tables print with 3 decimals (counts without decimals, p-values below 0.001 as `<0.001`);
the values keep full precision: `estimates(fit)$est`, `as.data.frame(estimates(fit))`.

`anova()` compares fits of the same data: AIC and BIC, and for nested models a
likelihood-ratio test; `model = "null"` is `model = "latreg"` without covariates (it leaves
them out):

```{r}
f0 <- ecm(dat, model = "null")
f1 <- ecm(dat, model = "latreg")
anova(f0, f1)
```

Other settings: `se`, and `prior = "map"` with `priors = rtirt_priors(...)` for the posterior
mode in the item parameters (small samples); `control = ecm_control()` sets `nodes`
(Gauss-Hermite nodes, default 41), `max_iter`, `tol` and `accelerate`.

## Gibbs from the ECM solution

`init = fit` starts every chain near the ECM solution:

```{r}
quiet <- gibbs_control(verbose = FALSE)
g <- gibbs(dat, model = "latent", init = fit, n_iter = 1500, n_chain = 2, seed = 1, control = quiet)
estimates(g, pars = "beta|b_ability") # posterior mean, SD, 95% interval, R-hat, ESS
convergence(g)
```

## Quantile models

Quantile models use the asymmetric Laplace (ALD) working likelihood: the estimates describe a
quantile of the log times (`"cross"`) or of speed (`"latent"`), not a model of the whole
distribution. `ecm()` maximizes the working likelihood and gives sandwich standard errors:

```{r}
fit_q <- ecm(dat, model = "latent", quantile = 0.25)
estimates(fit_q, pars = "beta|b_ability")
```

`gibbs()` samples the same model. For several quantile levels, give a vector; the result holds one
model per level:

```{r}
qs <- gibbs(dat,
  model = "latent", quantile = c(0.25, 0.5, 0.75), n_iter = 1000, n_chain = 2,
  seed = 1, control = quiet
)
qs
```

Use longer chains in practice (the defaults are `n_iter = 5000`, `n_chain = 4`) and check
`convergence()`. The [intro tutorial](https://github.com/jiewenTsai/birtExamples/blob/main/birtRcpp/tutorial/birtRcpp_intro.md)
(in Traditional Chinese; read it in the browser) explains every line and output; its source is
`system.file("tutorial", "birtRcpp_intro.qmd", package = "birtRcpp")`.

## References

Meng, X.-L., & Rubin, D. B. (1993). Maximum likelihood estimation via the ECM algorithm: A
general framework. *Biometrika, 80*, 267-278.

Polson, N. G., Scott, J. G., & Windle, J. (2013). Bayesian inference for logistic models using
Polya-Gamma latent variables. *Journal of the American Statistical Association, 108*, 1339-1349.

van der Linden, W. J. (2007). A hierarchical framework for modeling speed and accuracy on test
items. *Psychometrika, 72*, 287-308.