Within-group residual correlation (R-side autocorrelation)

ar(), ma(), arma(), cosy() and unstr() are written as terms of the model formula, next to the fixed and random effects. They come in two forms. The covariance form (cov = TRUE, and always for cosy() and unstr()) makes the residuals of one group a single correlated draw instead of independent ones:

Value

ar(), ma(), arma(), cosy() and unstr() are formula terms, not free-standing functions: bf() reads them at parse time, and the value they contribute is the fitted autocorrelation block of the model, reachable through autocor_matrix(), VarCorr() and confint_varcorr(). This page itself documents the term grammar and returns nothing.

Details

$$y_g \sim N(\mu_g,\; D_g R D_g), \qquad D_g = \mathrm{diag}(\sigma_i, i \in g).$$

R is a unit-diagonal correlation matrix over the time points and $\sigma_i$ is the marginal residual standard deviation of row $i$. Under cosy() and unstr() that is sigma. Under ar(), ma() and arma() sigma is the standard deviation of the INNOVATIONS, as in brms, and the marginal one is larger: for an AR(1), $\sigma_i = \sigma / \sqrt{1 - \phi^2}$; see "The scale of sigma". A sigma ~ ... distributional model enters through the diagonal. Nothing is added to the linear predictor, so fitted(), predict.frmtmb_fit() and se.fit are unchanged; what changes is the likelihood. This is the model nlme::gls(correlation = corAR1()) fits, and the one brms fits under cov = TRUE.

The residual-regression form is brms's DEFAULT for ar(), ma() and arma(), so it is what those terms mean without cov = TRUE. It is described in its own section below.

Structures

ar(time, gr, p = 1, cov = TRUE): Stationary AR(p) correlation. p = 1 gives $R_{ij} = \rho^{|i-j|}$.
ma(time, gr, q = 1, cov = TRUE): Invertible MA(q): correlation dies after lag q.
arma(time, gr, p = 1, q = 1, cov = TRUE): Stationary, invertible ARMA(p, q).
cosy(time, gr): Compound symmetry - one correlation shared by every pair. Equivalent in fit to a random intercept per group, but the correlation may also be negative (down to -1 / (d - 1)), which a variance component cannot express.
unstr(time, gr): One free correlation per pair of time levels, through the same Cholesky parameterization the us() random-effect structure uses.
fcor(M): A FIXED residual matrix over all the observations, with no parameter; see "A fixed residual matrix".

brms's sar(), the spatial simultaneous autoregressive term, is in the extension package frmtmb.spatial. The AR and MA coefficients are estimated through partial autocorrelations (the Monahan/Jones transform nlme::corARMA also uses), so every parameter value is a stationary and invertible process and the optimizer cannot leave the parameter space. The ARMA autocorrelation function is exact, not a truncated MA($\infty$) expansion.

Arguments

The argument order is brms's, so the FIRST positional argument is time and the second is gr: write cosy(gr = subj), not cosy(subj).

time: ONE variable name, whose levels index the correlation. Omit it (brms's time = NA) to use each row's position within its group. An expression is refused rather than evaluated, as it is in brms: the distinct values of the index are the time points and the distance between them is the lag, so ar(x + t, g) would be indexed by the sum. Compute the column first and name it.
gr: The grouping variable the residual factorizes over, or several crossed with :. Omit it to treat the whole data set as one series. gr = a:b crosses the two, on numeric codes as well as on factors; every other operator is refused rather than evaluated, again as in brms ("It may contain only variable names combined by the symbol ':'"), because gr = a/b would otherwise group by the quotient. gr = factor(g) and gr = interaction(a, b) are refused for the same reason: the grouping variable is made a factor here anyway, so write gr = g and gr = a:b. The crossing PASTES the levels with an underscore, as brms's combine_groups() does, so it is not injective: levels ("x_1", "2") and ("x", "1_2") both become "x_1_2" and share one series, where interaction() would keep them apart. brms merges them too; make the column yourself if the labels can collide.
p, q: Autoregressive and moving-average orders.
cov: FALSE, the default as in brms, gives the residual-regression form; TRUE gives the covariance form. They are different likelihoods, not two implementations of one.

The residual-regression form (`cov = FALSE`)

This is the model brms 2.23.0 fits for ar(), ma() and arma() when cov is not set. Within each group, in time order, let $r_t = y_t - \mu_t$ and $$e_t = r_t - \sum_{i=1}^{q} \theta_i e_{t-i},$$ with $e_s = 0$ before the group's first row. The mean of row $t$ becomes $$\mu^*_t = \mu_t + \sum_{i=1}^{q} \theta_i e_{t-i} + \sum_{i=1}^{p} \phi_i e_{t-i},$$ and every row keeps the family's own density at $\mu^*_t$. The MA part enters $e$ and the AR part does not, exactly as in brms's Stan code. So:

It is a CONDITIONAL likelihood. A group starts with no residuals behind it, so lag $i$ first reaches the row at within-group position $i + 1$: the FIRST row of a group gets no lagged term, a row at position $k$ gets only the lags up to $k - 1$, and from position $\max(p, q) + 1$ onward a row gets all of them. None of the opening rows gets the stationary distribution the covariance form gives them. For one series this is the conditional sum of squares of stats::arima(method = "CSS").
The coefficients ar[i] and ma[i] are unconstrained reals, as in brms, because nothing in this likelihood needs a stationary or invertible process. A prior or a bound with set_prior(class = "ar") acts on them directly, at any order.
sigma is the sd of each row around $\mu^*_t$, which for a correct model is the innovation sd.
The lag is counted in ROWS of the group, sorted by time, as brms counts it: a group that skips a time point treats the next row as one step later. The covariance form counts levels instead (see below).
Families: gaussian() and student(), the two families brms gives a residual y - mu. A predicted nu ~ ... works. Every other family has the model of "Latent residuals" instead. With a non-identity link the term is added to mu on the response scale, after the inverse link, as brms's Stan code adds it.
It combines with weights(), cens(), trunc(), mi() on the response (the residual is then taken against the imputed value), random effects and rescor = TRUE, as in brms. se() is refused, as brms refuses it.
After the fit, fitted(), frm_linpred() and predict(type = "response") give brms's one-step mean $\mu^*_t$, which reads the OBSERVED earlier residuals; so residuals() are against it. newdata must then carry the response, and its rows form their own groups in their own time order. predict.frmtmb_fit() draws each row around its one-step mean, as brms's posterior_predict() does. simulate() instead runs the recursion over its own draws, so a replicate carries the fitted serial dependence; on newdata it starts each group from an empty past and needs no response, which makes it the forecasting route. conditional_effects() and emmeans() drop the term, as brms's do. autocor_matrix() and residuals(type = "osa") are refused: this form defines no correlation matrix, and its tape reads the response as data. The pointwise log-density of this form is each row's own density given its group's observed past, which is what brms's log_lik() returns. arma_cond_resp() and arma_cond_dpars() (frmtmb-sampling-api) are what a sampling extension reads to build it.

The pointwise log-density

A residual covariance term (the covariance form of ar(), ma() and arma(), cosy(), unstr(), fcor() and the terms of frmtmb_register_autocor()) has no density per row. brms's log_lik() gives the leave-one-out conditional density instead, the density of one row given the other rows of its group, after the loo package's vignette on non-factorized models, and frmtmb.sample::log_lik() gives the same, through autocor_row_loglik() (frmtmb-sampling-api). With $Q$ the precision matrix of a group's response and $e$ its residual, row $i$ is normal with mean $y_i - (Qe)_i / Q_{ii}$ and variance $1 / Q_{ii}$, or for student() a t with $\nu + k - 1$ degrees of freedom, $k$ the size of the group. The values are brms's, so the row sum is NOT the log-likelihood of the fit, and loo() on them leaves out one row with the rest of its group kept. A group missing a time point gives the conditional of the model this page describes (lags from the time levels), which differs from brms's for ar(), ma() and arma() (lags by position) exactly as the two likelihoods differ.

A fixed residual matrix, `fcor(M)`

fcor(M) makes the residuals of all the rows one draw with covariance $\sigma^2 M$, or $D M D$ with $D = \mathrm{diag}(\sigma_i)$ under a predicted sigma ~ ..., as brms 2.23.0 fits it. student() gets the multivariate-t analog, with a constant nu. M is a known matrix in data2 over every row of the data, and the rows the model drops (missing values) are dropped from it, as brms drops them. M may be a correlation or a covariance; for a covariance, sigma is a scale factor, as in brms. A vector is the diagonal. M must be symmetric and positive definite. brms 2.23.0 checks the second with min(eigen(M)$values <= 0), which refuses only a matrix whose eigenvalues are all non-positive; frmtmb refuses every matrix that is not positive definite.

The term has no parameter, and the mean of the response is mu. fitted(), predict() and simulate() take newdata with one row per row of M, read in their order, as brms takes it, so a counterfactual such as x + 1 everywhere works; when the fit dropped rows for missing values, newdata covers them too, and a row with a missing predictor gives NA. Other newdata is refused; brms refuses it too unless a matrix over the new rows comes in its newdata2. frm_linpred() gives mu on any new data. Pearson residuals divide by the marginal sd, $\sigma \sqrt{M_{ii}}$. autocor_matrix() returns M over the fitted rows.

Families

This section and the two after it describe the covariance form on gaussian() and student(), the two families with a residual y - mu, and exactly the two brms treats this way. student() gets the multivariate-t analog, with one nu per group, so a predicted nu ~ ... is refused. Every other family has the model of "Latent residuals".

Time points, gaps and ragged groups

The lag between two rows is the distance between their positions in the GLOBAL set of time levels, so a group missing week 3 gets cor(week2, week4) = rho^2. That is nlme's reading (corAR1(form = ~ week | subj)) and what the agreement tests pin down; brms instead indexes ar()/ma()/arma()/cosy() by position within the group and treats a missing row as no gap (only its unstr() carries a time index). On complete balanced groups the two coincide. Time levels that are whole numbers but not consecutive warn, because the lag is then not the one the labels suggest.

Groups of different sizes are handled by construction: the density is evaluated one PATTERN at a time (a pattern being a distinct set of present time levels), so a balanced design costs one Cholesky per gradient evaluation and ragged data costs one per distinct pattern. Duplicate (gr, time) pairs are refused, as they are in brms.

What cannot be combined with it

The likelihood is a joint density over each group, so it no longer factorizes into per-row contributions. weights(), cens(), trunc(), se() and mi() on the same response are therefore refused, as are rescor = TRUE, mixture families, quadrature = TRUE, frm_simulate() and residuals(type = "osa"). brms refuses weights(), cens(), trunc() and rescor = TRUE too; it fits se() here (see "Divergence from brms"). Random effects ARE allowed and are the point of the feature: the marginal likelihood is a Laplace approximation over the modes with the correlated residual density inside, which reproduces nlme::lme(random = ~ 1 | subj, correlation = corAR1()).

Latent residuals

On a family with no residual y - mu, these terms add a latent residual to the linear predictor of mu, one per row, as brms 2.23.0 does: $$\eta_i = \dots + err_i.$$ sderr is the standard deviation of the latent residuals, with brms's prior class "sderr". The terms give err its law:

ar(), ma(), arma() with cov = TRUE: Order one only. Each group is one draw of the stationary process, and sderr is the sd of its innovations, as sigma is on gaussian(). brms refuses a higher order, and so does frm().
ar() with cov = FALSE: Any order. The latents z are independent, and $err_i = sderr (z_i + \sum_j ar_j z_{i-j})$. Each ar[j] is in (-1, 1). This is a moving average of the latents, by brms's own program. On a family whose density brms writes on the response scale (Gamma(), Beta(), a non-canonical link, a cens() or trunc() response), brms adds the sum to mu after the inverse link, and only frm_sample() fits it; see "Where the Laplace approximation fails". ma() and arma() are refused, in brms's words, because AR and MA structures cannot be distinguished on one vector of latents.
cosy(), unstr(): sderr is the marginal sd. cosy is in (0, 1), as brms bounds it. unstr() reads the time levels of each row.

fcor() and sar() are refused, in brms's words, and so are the multi-category families and the mixtures, where the term would not sit on mu, and a nonlinear mu. A nonlinear shape or zi beside a linear mu fits, as brms builds it. The lag counts rows of a group, as in brms, and not time levels as under the covariance form on gaussian(). If every group has one row, frm() refuses the term: brms builds it, but the correlation is not seen in the data, so its chains walk a ridge.

The latents are integrated by the Laplace approximation, on sparse coordinates that leave the approximation unchanged, so a long series costs a banded inner Hessian and not a dense one. They are not random effects of a group: ranef() and VarCorr() do not list them, and re_formula = NA keeps them, as brms keeps err. summary() and confint_varcorr() report sderr and the coefficients, and variables() names them as brms does. Under frm_sample() the draws of the latents are named latent_state[k]. These are the chain's own states and not err: ma(), arma() and cosy() have one extra state per group, so 60 rows in 12 groups have 72 states.

After the fit, fitted(), frm_linpred(), residuals(), predict() and simulate() read the conditional mode of err on the fitted rows. On new rows, under cov = TRUE:

predict() and simulate() draw new latents per group of the new data, as brms does, at the estimates of sderr and the coefficients.
fitted() is the mean over the latent, as brms's posterior_epred() is, by Gauss-Hermite quadrature at each row's marginal sd. With a log link and a marginal sd of 1.15 this is 1.95 times the value at err = 0. marginaleffects reads the same mean. Its delta method moves beta and holds sderr and the coefficients at their estimates, so its standard error leaves the latent's draw out, as on any fit with a group-level term: predictions() conditions on what get_coef() does not hold. Est.Error of fitted() holds the draw beside the coefficients', as on a new me() row, and its Q columns read that draw: on poisson Est.Error 0.133 became the draws' 1.005.
frm_linpred() stays at err = 0, because it is the predictor.

Under cov = FALSE new rows are refused, in brms's words. conditional_effects() and emmeans() drop the term, as brms's do.

The joint covariance of the estimates and the latents is held as the sparse Cholesky factor of their joint precision, which is banded along a chain, so a standard error reads it at a cost that grows with the rows, not with their square. fitted(), residuals() and frm_linpred(se.fit = TRUE) on the fitted rows read each row's variance through the selected inverse of the factor: on one poisson series, 0.4 s at 8000 rows and 1.7 s at 32000. A hurdle or zero-inflated mean reads the latent of mu and a level of hu ~ (1 | id) or zi ~ (1 | id) on one row, a pair that no row's likelihood couples; that pair is read from a column of the inverse solved once per level, 3.3 s at 32000 rows. Rows in newdata read the estimates and the levels of any other group-level block alone, frm_curve() and conditional_effects() too.

A dense block of the joint covariance is formed only where a function returns one or must form one: frm_joint_cov(), frm_lp_basis() on the fitted rows, and the category probabilities of an ordinal family's fitted rows. A block of more than getOption("frmtmb.latres_dense_rows", 20000L)^2 numbers is refused, with its size in the message. The same option bounds a wider block: on any fit with group-level effects, latent or not, a vcov = that marginaleffects carries into those effects inverts the inner block, one row per effect, densely, and past that many inner rows the effects are held at their estimates with a warning (frmtmb-marginaleffects). Above 3000 rows the category probabilities' standard errors are NA, with a message of class frmtmb_latres_se_skipped; options(frmtmb.latres_se_rows = ) moves that limit. A fit whose covariance of the estimates is held as a factor, which an ill-conditioned design gives, keeps the dense joint covariance, and both limits apply to every standard error of its fitted rows.

A coefficient near an end of its range, or an sderr near 0 against the response's scale, is at a limit only if the fit's numbers say that end is not determined: the standard error of log sderr is above 1, or the row's own is not finite. ar[1] at 0.9992 on 2000 rows has a small standard error, and reports it. At a limit the fit says so with a message of class frmtmb_boundary_fit, summary() gives NA for Est.Error and the interval, and diagnose() lists the limit. At ar[1] = 1 with sderr = 0 the data read as cosy() or (1 | group), and only the marginal sd is known. At ar[1] = -1 they read as one size per group whose sign alternates by row, which (1 | group) does not fit. When ar[1] and ma[1] cancel, the data show no correlation and neither coefficient is known.

One decision that differs from brms 2.23.0, under cov = FALSE only: brms's posterior_epred(), log_lik() and predict() read err and leave out the AR sum of its own Stan program, so their row densities do not add up to the model. Summed over the rows, brms's log_lik() misses the program's log_prob by -6.78 nats on poisson and -113.48 on Gamma(log). The methods here read the model, whose row densities add up to log_prob. At ar 0.6 and sderr 0.5, the two means differ by up to a factor of 1.62 on poisson and 5.18 on Gamma(log).

Where the Laplace approximation fails

The inner problem of a latent residual reads what one row tells its latent. The share of a row is $d v / (1 + d v)$, with d the curvature of the row's log density in its linear predictor and v the latent's variance. It has no units. Where it is small, or where a dispersion trades off against the latent, the approximation is biased. Median over 40 data sets, against the exact maximum likelihood estimate by forward recursion:

Design	Share	`sderr`	Exact
bernoulli, 20 groups of 10	0.11	0.46	1.71
cumulative, 5 categories	0.13	0.43	0.73
poisson, mean about 0.25	0.11	0.55	0.33
binomial, 3 trials	0.34	0.67	0.81
`Beta()`, phi 4	0.18	0.26	0.42
`beta_binomial()`, 10 trials, phi 10	0.16	0.26	0.65
`zero_one_inflated_beta()`, phi 10	0.24	0.28	0.59
negbinomial, mean 1, shape 2	0.19	0.46	0.57
negbinomial, mean 5, shape 2	0.28	0.43	0.55
`Gamma(link = "log")`, shape 5, `sderr` 0.05	0.12	0.12	0.05
binomial, 5 trials	0.49	0.74	0.80
poisson, mean about 1, `sderr` 0.5	0.25	0.47	0.44
negbinomial, mean 20, shape 4	0.48	0.47	0.52
`Gamma()`, shape 50	0.66	0.31	0.31
`poisson("identity")`, mean 50	0.27	3.79	3.73
lognormal, sigma 2, where Laplace is exact	0.28	0.95	0.95

frm() warns, with class frmtmb_latres_laplace, when the rows are ordinal, whose error is not local, or when the first-order Laplace correction moves log sderr or ar[1] by more than 0.15, or cannot be read. The correction is C = -1/8 sum(h4 v^2) + 1/8 sum(h3 h3 v v v) + 1/12 sum(h3 h3 v^3), with h3 and h4 the third and fourth derivatives of a row's negative log density in its predictor and v the covariance of the latents at the mode. The exact value is about Laplace plus C: the median error falls from 0.340 to 0.012 on designs the rule was not read on, and the fit moved by that correction falls from a median of 0.106 to 0.010 in log sderr. Over 84 designs and 3130 fits, with a fit wrong when sderr is off the exact estimate by more than a factor of 1.25 or ar by more than 0.3:

Set of designs	Rule	Warns on a wrong fit	Warns on a right fit
designs no side had seen	this rule	0.816	0.169
designs the rule was read on	this rule	0.960	0.147
all	this rule	0.865	0.157

Above frmtmb.latres_se_rows rows, where the correction's dense v is the cost that limit avoids, the warning goes by the share above instead. It reads the rows' share of the latent, a row's trials and the dispersion's trade-off. suppressWarnings(expr, classes = "frmtmb_latres_laplace") silences it, and check_se does not reach it. frmtmb.sample::frm_sample() samples the same model and agrees with brms. A silent skew_normal() is right to be silent: at a coverage of 0.725, the exact estimate with this fit's standard error covers at 0.725 too, so the fault is the Wald interval and not the approximation.

frm() refuses cov = FALSE after the inverse link before the fit, with its reason: no number from a fit says whether it reads true. frm_sample() samples it. On benign designs, 25 groups of 6, 20 data sets each:

`ar()` on	Margin	`sderr`	Exact	`ar[1]`	Exact	Read true
`Gamma("log")`, mean 20	23.6	0.41	0.50	-0.86	0.47	2
`inverse.gaussian("log")`	33.7	0.28	0.40	-0.55	0.52	3
`weibull()`	27.4	0.35	0.40	-1.00	-0.04	6
`exponential()`	52.2	0.29	0.25	-0.55	0.37	13
`poisson("sqrt")`, mean 49	72.5	0.50	0.50	-0.34	0.40	14
`Beta()`, phi 50	6.6	0.11	0.12	0.49	0.48	16

beta_binomial() at 10 trials and phi 10 stops at a nonzero code on 20 of 20, with a gap of 9.0 nats to the exact value, and zero_one_inflated_beta() at phi 10 reads sderr 0.28 for 0.38 with code 0.

On a link that does not hold the mean under cov = TRUE (identity, inverse, 1/mu^2), the margin is the smallest distance of eta from the link's edge, in marginal sds of the latent. frm() refuses after the fit below a margin of 2. On 680 fits against the exact marginal:

Margin	Fits	Read true	Coverage of log `sderr`	A start reads NaN
below 1	62	0.00	0.00	0.89
1 to 2	88	0.16	0.13	0.09
2 to 3	82	0.50	0.61	0
3 to 5	182	0.62	0.81	0
5 and above	266	0.55	0.98	0

A dispersion the rows do not ask for can run off beside a latent per row, and past a log of about 25 the objective is rounding: logLik read +6.5e17 on beta_binomial(). A constant dispersion is held at a log of 15, where the value agrees with the fit without the box to 1e-3, and a dispersion with a predictor reads NaN above 25. A fit at the box warns with class frmtmb_latres_dispersion.

frm() also refuses mi() on the response beside a latent residual, whose approximation has not been measured.

Where the parameters appear

On the internal (unconstrained) scale they are the thetaac_* rows of confint.frmtmb_fit() and of vcov(full = TRUE). On the natural scale, summary() prints them under "Within-group residual correlation" and confint_varcorr() reports one row per parameter under brms's names - ar[1], ma[1], cosy, cortime__<t1>__<t2> - with a delta-method interval (Fisher-z for the bounded ones, the identity for the coefficients of a higher-order AR/MA process). hypothesis() sees them under the same names, ar[1] included, which brms's renaming lets an expression carry. autocor_matrix() returns the fitted R. They are NOT part of VarCorr(), which reports random-effect blocks. set_prior() targets them by brms's classes, "ar", "ma", "cosy" and "cortime", one prior slot per parameter as in brms; see its section Residual correlation.

The scale of `sigma`

Under the covariance form of ar(), ma() and arma(), sigma is the standard deviation of the innovations, as in brms 2.23.0, whose cholesky_cor_ar1(), cholesky_cor_ma1() and cholesky_cor_arma1() are the autocovariance of the process at unit innovation variance. The marginal residual standard deviation is sigma times $1 / \sqrt{1 - \phi^2}$ for an AR(1), $\sqrt{1 + \theta^2}$ for an MA(1), and $\sqrt{(1 + \theta^2 + 2\phi\theta) / (1 - \phi^2)}$ for an ARMA(1, 1); at a higher order the factor is the square root of the process variance at unit innovation variance. So sigma(), VarCorr(), summary(), the standard errors, frm_linpred(dpar = "sigma"), a prior on class "sigma" and the coefficients of a sigma ~ x model are about the innovation sd, as brms's are, and pearson residuals, simulate() and the pointwise log-density read the marginal one. Under sigma ~ x the marginal sd of row $i$ is $\sigma_i$ times the same factor, as in brms's normal_time_het_lpdf(). Under cosy(), unstr() and the cov = FALSE form, sigma is what it is in brms too: the marginal sd, and the sd of a row around its one-step mean.

A fit saved by frmtmb 0.76.0 or earlier with one of these terms held the marginal sd as sigma. A method whose answer reads sigma or its coefficients refuses it; one that reads only the mean, such as fitted() or logLik(), answers, since the model is the same. update() fits it again.

conditional_effects() with method = "posterior_predict" draws its band at sigma, so here at the innovation sd, as brms draws it: brms's conditional_effects() drops every autocorrelation term.

Divergence from brms

brms limits cov = TRUE to order one ("Covariance formulation of ARMA structures is only possible for effects of maximal order one"); higher p and q are supported here, with sigma the innovation sd there too. brms 2.23.0 fits se() beside a cov = TRUE term, adding $se_i^2$ to the diagonal of each group's covariance (normal_time_hom_se_lpdf() and its relatives); frmtmb refuses that combination.

On a nonlinear body with cov = FALSE, brms 2.23.0's post-processing does not follow its own model. Its Stan program carries the ARMA term in mu, but predictor.bprepnl() has no predictor_ac(), so its posterior_epred() is the body alone and its log_lik() is the density at an ARMA coefficient of 0: on one model with a real AR part, -185.87 against -150.99 from its own log_prob(). fitted(), predict() and logLik() here follow the model, as they do for a linear mu, where brms and frmtmb agree.

Examples

set.seed(1)
d <- expand.grid(week = 1:5, subj = factor(1:30))
d$x <- rnorm(150)
e <- as.vector(vapply(1:30, function(i) {
  as.vector(stats::filter(rnorm(5), 0.6, "recursive"))
}, numeric(5)))
d$y <- 1 + 0.5 * d$x + e

fit <- frm(bf(y ~ x + ar(week, subj, cov = TRUE)) + gaussian(),
           data = d)
summary(fit)
#>  Family: gaussian 
#>  Links: mu = identity; sigma = log
#> 
#> Formula: y ~ x + ar(week, subj, cov = TRUE) 
#>    Data: d (Number of observations: 150) 
#>  Method: ML   logLik: -214.257   AIC: 436.514   BIC: 448.556 
#> 
#> Correlation Structures:
#>       Estimate Est.Error l-95% CI u-95% CI
#> ar[1]     0.52      0.07     0.37     0.65
#> 
#> Regression Coefficients:
#>           Estimate Est.Error l-95% CI u-95% CI z value Pr(>|z|)
#> Intercept     1.11      0.14     0.83     1.38    7.93  2.1e-15
#> x             0.45      0.08     0.30     0.61    5.63  1.8e-08
#> 
#> Further Distributional Parameters:
#>       Estimate Est.Error l-95% CI u-95% CI
#> sigma     0.98      0.06     0.87      1.1
autocor_matrix(fit)
#>            1         2         3         4          5
#> 1 1.00000000 0.5226004 0.2731112 0.1427280 0.07458974
#> 2 0.52260044 1.0000000 0.5226004 0.2731112 0.14272804
#> 3 0.27311122 0.5226004 1.0000000 0.5226004 0.27311122
#> 4 0.14272804 0.2731112 0.5226004 1.0000000 0.52260044
#> 5 0.07458974 0.1427280 0.2731112 0.5226004 1.00000000

# brms's default residual-regression form: mu gains ar * (y - mu) of
# the previous row of the same subject
fit0 <- frm(bf(y ~ x + ar(week, subj)) + gaussian(), data = d)
summary(fit0)
#>  Family: gaussian 
#>  Links: mu = identity; sigma = log
#> 
#> Formula: y ~ x + ar(week, subj) 
#>    Data: d (Number of observations: 150) 
#>  Method: ML   logLik: -214.929   AIC: 437.857   BIC: 449.9 
#> 
#> Correlation Structures:
#>       Estimate Est.Error l-95% CI u-95% CI
#> ar[1]     0.52      0.08     0.37     0.68
#> 
#> Regression Coefficients:
#>           Estimate Est.Error l-95% CI u-95% CI z value Pr(>|z|)
#> Intercept     1.12      0.13     0.86     1.38    8.35  < 2e-16
#> x             0.47      0.08     0.31     0.63    5.83  5.6e-09
#> 
#> Further Distributional Parameters:
#>       Estimate Est.Error l-95% CI u-95% CI
#> sigma     1.01      0.06     0.91     1.14
head(fitted(fit0))
#>       Estimate  Est.Error      Q2.5     Q97.5
#> [1,] 0.8236406 0.14500849 0.5394291 1.1078520
#> [2,] 1.3714204 0.08119577 1.2122796 1.5305611
#> [3,] 0.7952423 0.10168233 0.5959486 0.9945360
#> [4,] 1.7156570 0.17627991 1.3701548 2.0611593
#> [5,] 0.6977275 0.11448821 0.4733348 0.9221203
#> [6,] 0.7316354 0.15168752 0.4343334 1.0289375

# compound symmetry, and the unstructured correlation over the five
# weeks
frm(bf(y ~ x + cosy(week, subj)) + gaussian(), data = d)
#>  Family: gaussian 
#>  Links: mu = identity; sigma = log
#> 
#> Formula: y ~ x + cosy(week, subj) 
#>    Data: d (Number of observations: 150) 
#>  Method: ML   logLik: -220.426   AIC: 448.852   BIC: 460.894 
#> 
#> Correlation Structures:
#>      Estimate Est.Error l-95% CI u-95% CI
#> cosy     0.36      0.09     0.17     0.52
#> 
#> Regression Coefficients:
#>           Estimate Est.Error l-95% CI u-95% CI z value Pr(>|z|)
#> Intercept     1.09      0.15     0.80     1.38    7.43  1.1e-13
#> x             0.41      0.09     0.23     0.59    4.52  6.2e-06
#> 
#> Further Distributional Parameters:
#>       Estimate Est.Error l-95% CI u-95% CI
#> sigma     1.15      0.08        1     1.32
frm(bf(y ~ x + unstr(week, subj)) + gaussian(), data = d)
#>  Family: gaussian 
#>  Links: mu = identity; sigma = log
#> 
#> Formula: y ~ x + unstr(week, subj) 
#>    Data: d (Number of observations: 150) 
#>  Method: ML   logLik: -213.317   AIC: 452.634   BIC: 491.773 
#> 
#> Correlation Structures:
#>               Estimate Est.Error l-95% CI u-95% CI
#> cortime__1__2     0.44      0.14     0.14     0.67
#> cortime__1__3     0.28      0.16    -0.05     0.56
#> cortime__1__4     0.24      0.18    -0.12     0.55
#> cortime__1__5     0.10      0.19    -0.26     0.44
#> cortime__2__3     0.61      0.10     0.36     0.77
#> cortime__2__4     0.38      0.16     0.04     0.64
#> cortime__2__5     0.14      0.18    -0.21     0.46
#> cortime__3__4     0.55      0.12     0.28     0.75
#> cortime__3__5     0.31      0.16    -0.02     0.58
#> cortime__4__5     0.49      0.14     0.17     0.71
#> 
#> Regression Coefficients:
#>           Estimate Est.Error l-95% CI u-95% CI z value Pr(>|z|)
#> Intercept     1.12      0.14     0.84     1.40    7.82  5.2e-15
#> x             0.44      0.09     0.26     0.61    4.94  7.8e-07
#> 
#> Further Distributional Parameters:
#>       Estimate Est.Error l-95% CI u-95% CI
#> sigma     1.15      0.08        1     1.32

# a random intercept alongside the correlated residual is allowed
frm(bf(y ~ x + (1 | subj) + ar(week, subj, cov = TRUE)) + gaussian(),
    data = d)
#>  Family: gaussian 
#>  Links: mu = identity; sigma = log
#> 
#> Formula: y ~ x + (1 | subj) + ar(week, subj, cov = TRUE) 
#>    Data: d (Number of observations: 150) 
#>  Method: ML   logLik: -214.232   AIC: 438.465   BIC: 453.518 
#> 
#> Multilevel Hyperparameters:
#> ~subj (Number of levels: 30) 
#>               Estimate Est.Error l-95% CI u-95% CI
#> sd(Intercept)     0.27      0.58        0    18.55
#> 
#> Correlation Structures:
#>       Estimate Est.Error l-95% CI u-95% CI
#> ar[1]     0.49      0.14     0.17     0.73
#> 
#> Regression Coefficients:
#>           Estimate Est.Error l-95% CI u-95% CI z value Pr(>|z|)
#> Intercept     1.11      0.14     0.83     1.38    7.81  5.7e-15
#> x             0.45      0.08     0.29     0.61    5.54  3.0e-08
#> 
#> Further Distributional Parameters:
#>       Estimate Est.Error l-95% CI u-95% CI
#> sigma     0.97      0.07     0.84     1.11

# a family with no residual: brms's latent residual, with its SD sderr
d$cnt <- rpois(150, exp(0.5 + 0.3 * d$x + e))
summary(frm(bf(cnt ~ x + ar(week, subj, cov = TRUE)), family = poisson(),
            data = d))
#>  Family: poisson 
#>  Links: mu = log
#> 
#> Formula: cnt ~ x + ar(week, subj, cov = TRUE) 
#>    Data: d (Number of observations: 150) 
#>  Method: ML   logLik: -343.263   AIC: 694.526   BIC: 706.568 
#> 
#> Correlation Structures:
#>       Estimate Est.Error l-95% CI u-95% CI
#> ar[1]      0.6       0.1     0.37     0.76
#> sderr      0.9       0.1     0.72     1.13
#> 
#> Regression Coefficients:
#>           Estimate Est.Error l-95% CI u-95% CI z value Pr(>|z|)
#> Intercept     0.65      0.17     0.32     0.97    3.87  0.00011
#> x             0.27      0.10     0.07     0.48    2.62  0.00867