pffr: Penalized function-on-function regression

Description

Implements additive regression for functional and scalar covariates and functional responses. This function is a wrapper for mgcv's gam and its siblings to fit models of the general form $E(Y_i(t)) = g(\mu(t) + \int X_i(s)\beta(s,t)ds + f(z_{1i}, t) + f(z_{2i}) + z_{3i} \beta_3(t) + \dots$ with a functional (but not necessarily continuous) response $Y(t)$, response function $g$, (optional) smooth intercept $\mu(t)$, (multiple) functional covariates $X(t)$ and scalar covariates $z_1$, $z_2$, etc.

Usage

pffr(formula, yind, data = NULL, algorithm = NA,
    method = "REML", tensortype = c("te", "t2"),
    bs.yindex = list(bs = "ps", k = 5, m = c(2, 1)),
    bs.int = list(bs = "ps", k = 20, m = c(2, 1)), ...)

Arguments

formula

a formula with special terms as for gam, with additional special terms ff() and c()

yind

a vector with length equal to the number of columns of the matrix of functional responses giving the vector of evaluation points $(t_1, \dots ,t_{G})$. If formula contains an ff-term which

algorithm

the name of the function used to estimate the model. Defaults to gam if the matrix of functional responses has less than 2e5 data points and to ba

data

an (optional) data.frame or a named list containing the data. The functional response and functional covariates have to be supplied as n by matrices, i.e. each row is one functional observation. The model is t

method

Defaults to "REML"-estimation, including of unknown scale. See gam for details.

bs.yindex

a named (!) list giving the parameters for spline bases on the index of the functional response. Defaults to list(bs="ps", k=5, m=c(2, 1)), i.e. 5 cubic B-splines bases with first order difference penalty.

bs.int

a named (!) list giving the parameters for the spline basis for the global functional intercept. Defaults to list(bs="ps", k=20, m=c(2, 1)), i.e. 20 cubic B-splines bases with first order difference penalty.

tensortype

which typ of tenor product splines to use. One of "te" or "t2", defaults to te

...

additional arguments that are valid for gam or bam. weights, subset, offset are not yet implemented!

Value

a fitted pffr-object, which is a gam-object with some additional information in an pffr-entry. If algorithm is "gamm" or "gamm4", only the $gam part of the returned list is modified in this way.

Details

The routine can estimate

(nonlinear, and possibly multivariate) effects of (one or multiple) scalar covariates$z$that vary smoothly over the index$t$of$Y(t)$(e.g.$f(z_{1i}, t)$, specified in theformulasimply as~s(z1)),
(nonlinear) effects of scalar covariates that are constant over$t$(e.g.$f(z_{2i})$, specified as~c(s(z2)), or$\beta_2 z_{2i}$, specified as~c(z2)),
linear functional effects of scalar (numeric or factor) covariates that vary smoothly over$t$(e.g.$z_{3i} \beta_3(t)$, specified as~z3),
function-on-function regression terms (e.g.$\int X_i(s)\beta(s,t)ds$, specified as~ff(X, yindex=t, xindex=s), seeff).

Use the c()-notation to denote model terms that are constant over the index of the functional response. Internally, univariate smooth terms without a c()-wrapper are expanded into bivariate smooth terms in the original covariate and the index of the functional response. Bivariate smooth terms (

s(),
  te()

or t2()) without a c()-wrapper are expanded into trivariate smooth terms in the original covariates and the index of the functional response. Linear terms for scalar covariates or categorical covariates are expanded into varying coefficient terms, varying smoothly over the index of the functional response. For factor variables, a separate smooth function with its own smoothing parameter is estimated for each level of the factor. Functional random intercepts $B_{0g(i)}(t)$ for a grouping variable g can be specified via ~s(g, bs="re")), functional random slopes $u_i B_{1g(i)}(t)$ in a numeric variable u via

~s(g, u,
  bs="re")

). The marginal spline basis used for the index of the the functional response is specified via the global argument bs.yindex. If necessary, this can be overriden for any specific term by supplying a bs.yindex-argument to that term in the formula, e.g. ~s(x, bs.yindex=list(bs="tp", k=7)) would yield a tensor product spline for which the marginal basis for the index of the response are 7 cubic thin-plate spline functions overriding the global default for the basis and penalty on the index of the response given by the global bs.yindex-argument . Use ~-1 + c(1) + ... to specify a model with only a constant and no functional intercept. The functional response and functional covariates have to be supplied as n by matrices, i.e. each row is one functional observation. The model is then fitted with the data in long format, i.e., the rows of the matrix of the functional response evaluations $Y_i(t)$ are stacked into one long vector and the covariates are expanded/repeated correspondingly. This means the models get quite big fairly fast, since the effective number of rows in the design matrix is number of observations times number of evaluations of $Y(t)$ per observation. Note that pffr overrides the default identifiability constraints ($\sum_{i,t} \hat f(z_i, x_i, t) = 0$) implemented in mgcv for tensor product terms whose marginal terms include the index of the functional response $t$. Instead, $\sum_i \hat f(z_i, x_i, t) = 0$ for all $t$ is enforced, so that effects varying over $t$ can be interpreted as deviations from the global functional intercept. I recommend using centered scalar covariates for terms like $z \beta(t)$ (~z) and centered functional covariates with $\sum_i X_i(t) = 0$ for all $t$ in ff-terms so that the global functional intercept can be interpreted as the global mean function.

References

Ivanescu, A., Staicu, A.-M., Scheipl, F. and Greven, S. (2012). Penalized function-on-function regression. (submitted) http://biostats.bepress.com/jhubiostat/paper240/ Scheipl, F., Staicu, A.-M. and Greven, S. (2012). Additive Mixed Models for Correlated Functional Data. (submitted) http://arxiv.org/abs/1207.5947

Examples

Run this code

###############################################################################
# univariate model:
# Y(t) = f(t)  + \int X1(s)\beta(s,t)ds + eps
data1 <- pffrSim(scenario="1")
t <- attr(data1, "yindex")
s <- attr(data1, "xindex")
m1 <- pffr(Y ~ ff(X1, yind=t, xind=s), data=data1)
summary(m1)
plot(m1, pers=TRUE)

###############################################################################
# multivariate model:
# Y(t) = f0(t)  + \int X1(s)\beta1(s,t)ds + \int X2(s)\beta2(s,t)ds +
#		xlin \beta3(t) + f1(xte1, xte2) + f2(xsmoo, t) + beta4 xconst + eps
data2 <- pffrSim(scenario="2", n=200)
t <- attr(data2, "yindex")
s <- attr(data2, "xindex")
m2 <- pffr(Y ~  ff(X1, xind=s, yind=t) + #linear function-on-function
                ff(X2, xind=s, yind=t) + #linear function-on-function
                xlin  +  #varying coefficient term
                c(te(xte1, xte2)) + #bivariate smooth term in xte1 & xte2, const. over Y-index
                s(xsmoo) + #smooth effect of xsmoo varying over Y-index
                c(xconst), # linear effect of xconst constant over Y-index
        data=data2)
summary(m2)
plot(m2, pers=TRUE)

Run the code above in your browser using DataLab