tslm: Fit a linear model with optional autoregressive errors

Description

`tslm()` is a teaching-friendly wrapper for fitting linear models with optional AR(p) error structures. Students specify the mean model using an ordinary formula and add an `ar(p)` term to request autoregressive errors.

Usage

tslm(formula, data = parent.frame(), time, method = "REML", ...)

Value

An object of class `tslm`, containing the original formula, the mean formula fitted internally, the AR order, the time variable if supplied, and the underlying fitted model.

Arguments

formula: a model formula. Use `ar(p)` in the right hand side to specify AR(p) errors, for example `y ~ x + ar(1)`.
data: an optional data frame containing the variables in the model. If omitted, variables are taken from the calling environment.
time: optional unquoted or quoted name of the time variable in `data` or in the calling environment. If omitted for an AR model, the row order of the model data is used.
method: fitting method passed to [nlme::gls()] for AR models. Defaults to `"REML"`.
...: additional arguments passed to [stats::lm()] or [nlme::gls()].

Details

When no `ar(p)` term is present, `tslm()` fits an ordinary [stats::lm()] model. When an `ar(p)` term is present, `tslm()` fits a [nlme::gls()] model with an AR(p) correlation structure using [nlme::corARMA()]. The `ar(p)` term changes the error model, not the mean-model terms printed in the formula.

The formula describes the mean model, just as it does for [stats::lm()]. The special term `ar(p)` is removed from the mean model before fitting and is used only to specify the correlation structure for the errors. For example, `log(passengers) ~ t + month + ar(1)` fits a trend and seasonal mean model with AR(1) errors.

For AR-error models, `time` should usually name the variable giving the time order of the observations. If `time` is omitted, `tslm()` fits the model using the row order of the model data and gives a warning so that this assumption is visible.

Diagnostic methods for AR-error models use normalised residuals by default, because these residuals account for the fitted correlation structure. Use `residualType = "response"` when the raw response residuals are required. `"normalised"` and `"normalized"` are both accepted for compatibility.

Examples

Run this code

data(beer.df)
fit = tslm(beer ~ t + ar(1), data = beer.df, time = t)
coef(fit)

data(airpass.df)
fitAr = tslm(log(passengers) ~ t + month + ar(1),
  data = airpass.df,
  time = t
)
summary(fitAr)
anova(fitAr)

plot(fitAr)
plot(fitAr, residualType = "response")