crossbasis: Generate a cross-basis matrix for a DLNM

Description

Generate the basis functions for the two spaces of predictor and lags, choosing among a set of possible bases. Then, these functions are combined in order to create the related cross-basis matrix, which can be included in a model formula to fit a distributed lag non-linear model (DLNM).

Usage

crossbasis(x, vartype="ns", vardf=1, vardegree=1, varknots=NULL,
	varbound=range(x), varint=FALSE, cen=TRUE, cenvalue=mean(x),
	maxlag=0, lagtype="ns", lagdf=1, lagdegree=1,lagknots=NULL,
	lagbound=c(0,maxlag), lagint=TRUE, group=NULL)

## S3 method for class 'crossbasis':
summary(object, ...)

Arguments

the predictor variable, defined as a numeric vector representing a complete series of ordered observations.

vartype, lagtype

type of basis. See Details below for the list of possible choices.

vardf, lagdf

dimension of the basis, equivalent to number of degrees of freedom spent to specify the relationship in each space. They depend on knots if provided, or on degree for type="poly".

vardegree, lagdegree

degree of polynomial. Used only for type equal to "bs" (degree of the piecewise polynomial for the B-spline) or "poly" (degree of the polynomial).

varknots, lagknots

knots location for the basis. They specify the position of the internal knots for "ns" and "bs", the cut-off points for "strata" (defining right-open intervals) and the threshold(s)/cut-off points for "lthr"

varbound, lagbound

boundary knots (sometimes called external knots). Used only for type equal to "ns" and "bs".

varint, lagint

logical. If TRUE, an 'intercept' is included in the basis, with different meanings depending on type above: see Details below. The default values should not be changed: see Warnings below.

cen

logical. If TRUE, the basis functions for the space of predictor are centered. See Note below.

cenvalue

centering value, used as a reference point for the predicted effects.

maxlag

a positive value defining the maximum lag.

group

a factor defining groups of observations, representing multiple series. Each series must be consecutive, complete and ordered.

object

a object of class "crossbasis".

...

additional arguments to be passed to summary.

Value

A matrix object of class "crossbasis" which can be included in a model formula in order to fit a DLNM. It contains the attributes crossdf (global number of degrees of freedom) and range (range of the original vector of observations). Additional attributes are returned that correspond to the arguments to crossbasis, and explicitly give type, df, degree, knots, bound, cen, cenvalue and maxlag related to the corresponding basis ( with prefix var- or lag-) for use of crosspred. The function summary.crossbasis returns a summary of the cross-basis matrix and the related attributes, and can be used to check the options for the bases chosen for the two dimensions.

Warnings

Meaningless combinations of arguments (for example the inclusion of knots lying outside the range for type equal to "strata" or thr-type) could lead to collinear variables, with identifiability problems in the model and the exclusion of some of them. It is strongly recommended to avoid the inclusion of an intercept in the basis for x (varint should be FALSE), otherwise a rank-deficient cross-basis matrix will be specified, causing some of the cross-variables to be excluded in the regression model. Conversely, an intercept should always be included in the basis for the space of lags when lagtype is equal to "ns", "bs", "strata" or "poly".

Details

The value in type defines the basis for each space (predictor and lags). It must be one of: "ns": natural cubic B-splines (constrained to be linear beyond the boundary knots). Specified by knots (internal knots) and bound (boundary or external knots). If knots is provided, the dimension df is set to length(knots)+1+int. An intercept is included if int=T (corresponding to a vector of 1's if df=1, involving a more complex parameterization for df>1). See the functions ns for additional information. The transformed variables can be centered at cenvalue. "bs": B-splines characterized by degree (degree of the piecewise polynomial). Specified by knots (internal knots) and bound (boundary or external knots). If knots is provided, the dimension df is set to length(knots)+degree+int; if not, df must be higher than degree+int. An intercept is included if int=T (corresponding to a vector of 1's if df=1, involving a more complex parameterization for df>1). See the functions bs for additional information. The transformed variables can be centered at cenvalue. "strata": strata variables (dummy parameterization) determined by internal cut-off values specified in knots, which represent the lower boundaries for the right-open intervals. Intervals containing no observation are automatically discarded. If knots is provided, the dimension df is set to length(knots)+int. A dummy variable for the reference stratum (the first one by default) is included if int=T, generating a full rank basis. Never centered. "poly": polynomial with power specified by degree. The dimension df is set to to degree+int. An intercept, corresponding to a vector of 1's (the power 0 of the polynomial) is included if int=T. The transformed variables can be centered at cenvalue. "integer": strata variables (dummy parameterization) for each integer values, expressly created to specify an unconstrained function in the space of lags. df is set automatically to the number of integer values minus 1 plus int. A dummy variable for the reference stratum (the first one by default) is included if int=T, generating a full rank basis. Never centered. "hthr", "lthr": high and low threshold parameterization, with a linear relationship above or below the threshold, respectively, and flat otherwise. The threshold is chosen by knots: if more than one is provided, a piecewise linear relationship is applied above the first knot or below the last one, respectively, with the slope changing at each further knot. df is automatically set to length(knots)+int. An intercept (corresponding to a vector of 1's) is included if int=T. Never centered. "dthr": double threshold parameterization (2 independent linear relationships above the second and below the first threshold, flat between them). The thresholds are chosen by knots. If only one is provided, the threshold is unique (V-model). If more than 2 are provided, the first and the last ones are chosen. df is automatically set to 2+int. An intercept (corresponding to a vector of 1's) is included if int=T. Never centered. "lin": linear relationship (untransformed apart from optional centering). df is automatically set to 1+int. An intercept (corresponding to a vector of 1's) is included if int=T. It can be centered at cenvalue. Some arguments can be automatically changed for not sensible combinations, or set to NULL if not required. Use summary.crossbasis to check the result. The argument group defines groups of observations representing independent series. Each series must be consecutive, complete and ordered. crossbasis is run on each of them applying the same cross-basis functions: default choices (knots position, range, etc.) are taken considering the pooled distribution. For a detailed illustration of the use of the functions, see: vignette("dlnmOverview")

References

Gasparrini A. Distributed lag linear and non-linear models in R: the package dlnm. Journal of Statistical Software. 2011; 43(8):1-20. URL http://www.jstatsoft.org/v43/i08/. Gasparrini A., Armstrong, B.,Kenward M. G. Distributed lag non-linear models. Statistics in Medicine. 2010; 29(21):2224-2234.

Examples

Run this code

# Example 1. See crosspred and plot.crosspred for other examples

### simple DLM
### space of predictor: linear effect for PM10
### space of predictor: 5df natural cubic spline for temperature
### lag function: 4th degree polynomial for PM10 up to lag15
### lag function: strata intervals at lag 0 and 1-3 for temperature

# CREATE THE CROSS-BASIS FOR EACH PREDICTOR AND CHECK WITH SUMMARY
basis.pm <- crossbasis(chicagoNMMAPS$pm10, vartype="lin", lagtype="poly",
	lagdegree=4, cen=FALSE, maxlag=15)
basis.temp <- crossbasis(chicagoNMMAPS$temp, vardf=5, lagtype="strata",
	lagknots=1, cenvalue=21, maxlag=3)
summary(basis.pm)
summary(basis.temp)

# RUN THE MODEL AND GET THE PREDICTION FOR PM10
library(splines)
model <- glm(death ~ basis.pm + basis.temp + ns(time, 7*14) + dow,
	family=quasipoisson(), chicagoNMMAPS)
pred.pm <- crosspred(basis.pm, model, at=0:20, cumul=TRUE)

# PLOT THE LINEAR EFFECT OF PM10 ALONG LAGS
plot(pred.pm, "slices", var=10, col=3, ylab="RR", ci.arg=list(density=15,lwd=2),
	main="Effect of a 10-unit increase in PM10 along lags")
plot(pred.pm, "slices", var=10, cumul=TRUE, ylab="Cumulative RR",
	main="Cumulative effect of a 10-unit increase in PM10 along lags")

# GET THE FIGURES FOR THE OVERALL EFFECT, WITH CI
pred.pm$allRRfit["10"]
cbind(pred.pm$allRRlow, pred.pm$allRRhigh)["10",]

### See the vignette 'dlnmOverview' for a detailed explanation of this example

Run the code above in your browser using DataLab