fda (version 2.4.0)

fRegress: Functional Regression Analysis

Description

This function carries out a functional regression analysis, where either the dependent variable or one or more independent variables are functional. Non-functional variables may be used on either side of the equation. In a simple problem where there is a single scalar independent covariate with values $z_i, i=1,\ldots,N$ and a single functional covariate with values $x_i(t)$, the two versions of the model fit by fRegress are the scalar dependent variable model $$y_i = \beta_1 z_i + \int x_i(t) \beta_2(t) \, dt + e_i$$ and the concurrent functional dependent variable model $$y_i(t) = \beta_1(t) z_i + \beta_2(t) x_i(t) + e_i(t).$$ In these models, the final term $e_i$ or $e_i(t)$ is a residual, lack of fit or error term. In the concurrent functional linear model for a functional dependent variable, all functional variables are all evaluated at a common time or argument value $t$. That is, the fit is defined in terms of the behavior of all variables at a fixed time, or in terms of "now" behavior. All regression coefficient functions $\beta_j(t)$ are considered to be functional. In the case of a scalar dependent variable, the regression coefficient for a scalar covariate is converted to a functional variable with a constant basis. All regression coefficient functions can be forced to be smooth through the use of roughness penalties, and consequently are specified in the argument list as functional parameter objects.

Usage

fRegress(y, ...)
## S3 method for class 'formula':
fRegress(y, data=NULL, betalist=NULL, wt=NULL,
                 y2cMap=NULL, SigmaE=NULL,
                 method=c('fRegress', 'model'), sep='.', ...)
## S3 method for class 'character':
fRegress(y, data=NULL, betalist=NULL, wt=NULL,
                 y2cMap=NULL, SigmaE=NULL,
                 method=c('fRegress', 'model'), sep='.', ...)
## S3 method for class 'fd':
fRegress(y, xfdlist, betalist, wt=NULL,
                     y2cMap=NULL, SigmaE=NULL, returnMatrix=FALSE, ...)
## S3 method for class 'fdPar':
fRegress(y, xfdlist, betalist, wt=NULL,
                     y2cMap=NULL, SigmaE=NULL, returnMatrix=FALSE, ...)
## S3 method for class 'numeric':
fRegress(y, xfdlist, betalist, wt=NULL,
                     y2cMap=NULL, SigmaE=NULL, returnMatrix=FALSE, ...)

Arguments

y

the dependent variable object. It may be an object of five possible classes:

• character or formula

{ a formula object or a character object that can be coerced into a formula providin

Value

• These functions return either a standard fRegress fit object or or a model specification:
• fRegress fit
  a list of class fRegress with the following components:

  • y

  { the first argument in the call to fRegress (coerced to class fdPar) } xfdlist{ the second argument in the call to fRegress. } betalist{ the third argument in the call to fRegress. } betaestlist{ a list of length equal to the number of independent variables and with members having the same functional parameter structure as the corresponding members of betalist. These are the estimated regression coefficient functions. } yhatfdobj{ a functional parameter object (class fdPar) if the dependent variable is functional or a vector if the dependent variable is scalar. This is the set of predicted by the functional regression model for the dependent variable. } Cmatinv{ a matrix containing the inverse of the coefficient matrix for the linear equations that define the solution to the regression problem. This matrix is required for function fRegress.stderr that estimates confidence regions for the regression coefficient function estimates. } wt{ the vector of weights input or inferred } If class(y) is numeric, the fRegress object also includes: df{ equivalent degrees of freedom for the fit. } OCV{ the leave-one-out cross validation score for the model. } gcv{ the generalized cross validation score. } If class(y) is either fd or fdPar, the fRegress object returned also includes 5 other components: y2cMap{ an input y2cMap } SigmaE{ an input SigmaE } betastderrlist{ an fd object estimating the standard errors of betaestlist } bvar{ a covariance matrix } c2bMap{ a map }

item

• data
• xfdlist
• fd
• betalist
• wt
• y2cMap
• SigmaE
• method
• sep
• returnMatrix
• ...
• model specification
• type
• nbasis
• xVars

code

fd

itemize

• xfdlist0

Details

Alternative forms of functional regression can be categorized with traditional least squares using the following 2 x 2 table: lcccc{ explanatory variable response | scalar | function | | scalar | lm | fRegress.numeric | | function | fRegress.fd or | fRegress.fd or | fRegress.fdPar | fRegress.fdPar or linmod } For fRegress.numeric, the numeric response is assumed to be the sum of integrals of xfd * beta for all functional xfd terms. fRegress.fd or .fdPar produces a concurrent regression with each beta being also a (univariate) function. linmod predicts a functional response from a convolution integral, estimating a bivariate regression function. In the computation of regression function estimates in fRegress, all independent variables are treated as if they are functional. If argument xfdlist contains one or more vectors, these are converted to functional data objects having the constant basis with coefficients equal to the elements of the vector. Needless to say, if all the variables in the model are scalar, do NOT use this function. Instead, use either lm or lsfit. These functions provide a partial implementation of Ramsay and Silverman (2005, chapters 12-20).

References

Ramsay, James O., Hooker, Giles, and Graves, Spencer (2009) Functional Data Analysis in R and Matlab, Springer, New York. Ramsay, James O., and Silverman, Bernard W. (2005), Functional Data Analysis, 2nd ed., Springer, New York.

See Also

fRegress.formula, fRegress.stderr, fRegress.CV, linmod

Examples

###
###
### scalar response and explanatory variable
###   ... to compare fRegress and lm
###
###
# example from help('lm')
     ctl <- c(4.17,5.58,5.18,6.11,4.50,4.61,5.17,4.53,5.33,5.14)
     trt <- c(4.81,4.17,4.41,3.59,5.87,3.83,6.03,4.89,4.32,4.69)
     group <- gl(2,10,20, labels=c("Ctl","Trt"))
     weight <- c(ctl, trt)
lm.D9 <- lm(weight ~ group)

fRegress.D9 <- fRegress(weight ~ group)

(lm.D9.coef <- coef(lm.D9))

(fRegress.D9.coef <- sapply(fRegress.D9$betaestlist, coef))

stopifnot(
all.equal(as.numeric(lm.D9.coef), as.numeric(fRegress.D9.coef))
)

###
###
### vector response with functional explanatory variable
###
###

##
## set up
##
annualprec <- log10(apply(CanadianWeather$dailyAv[,,"Precipitation.mm"],
                          2,sum))
# The simplest 'fRegress' call is singular with more bases
# than observations, so we use a small basis for this example
smallbasis  <- create.fourier.basis(c(0, 365), 25)
# There are other ways to handle this,
# but we will not discuss them here
tempfd <- smooth.basis(day.5,
          CanadianWeather$dailyAv[,,"Temperature.C"], smallbasis)$fd

##
## formula interface
##

precip.Temp.f <- fRegress(annualprec ~ tempfd)

##
## Get the default setup and modify it
##

precip.Temp.mdl <- fRegress(annualprec ~ tempfd, method='m')
# First confirm we get the same answer as above:
precip.Temp.m <- do.call('fRegress', precip.Temp.mdl)
stopifnot(
all.equal(precip.Temp.m, precip.Temp.f)
)

#  set up a smaller basis than for temperature
nbetabasis  <- 21
betabasis2.  <- create.fourier.basis(c(0, 365), nbetabasis)
betafd2.     <- fd(rep(0, nbetabasis), betabasis2.)
# add smoothing
betafdPar2.  <- fdPar(betafd2., lambda=10)

precip.Temp.mdl2 <- precip.Temp.mdl
precip.Temp.mdl2[['betalist']][['tempfd']] <- betafdPar2.

# Now do it.
precip.Temp.m2 <- do.call('fRegress', precip.Temp.mdl2)

# Compare the two fits
precip.Temp.f[['df']] # 26
precip.Temp.m2[['df']]# 22 = saved 4 degrees of freedom

(var.e.f <- mean(with(precip.Temp.f, (yhatfdobj-yfdPar)^2)))
(var.e.m2 <- mean(with(precip.Temp.m2, (yhatfdobj-yfdPar)^2)))
# with a modest increase in lack of fit.

##
## Manual construction of xfdlist and betalist
##

xfdlist <- list(const=rep(1, 35), tempfd=tempfd)

# The intercept must be constant for a scalar response
betabasis1 <- create.constant.basis(c(0, 365))
betafd1    <- fd(0, betabasis1)
betafdPar1 <- fdPar(betafd1)

betafd2     <- with(tempfd, fd(basisobj=basis, fdnames=fdnames))
# convert to an fdPar object
betafdPar2  <- fdPar(betafd2)

betalist <- list(const=betafdPar1, tempfd=betafdPar2)

precip.Temp <- fRegress(annualprec, xfdlist, betalist)
stopifnot(
all.equal(precip.Temp, precip.Temp.f)
)

###
###
### functional response with vector explanatory variables
###
###

##
## simplest:  formula interface
##
daybasis65 <- create.fourier.basis(rangeval=c(0, 365), nbasis=65,
                  axes=list('axesIntervals'))
Temp.fd <- with(CanadianWeather, smooth.basisPar(day.5,
                dailyAv[,,'Temperature.C'], daybasis65)$fd)
TempRgn.f <- fRegress(Temp.fd ~ region, CanadianWeather)

##
## Get the default setup and possibly modify it
##
TempRgn.mdl <- fRegress(Temp.fd ~ region, CanadianWeather, method='m')# make desired modifications here
# then run
TempRgn.m <- do.call('fRegress', TempRgn.mdl)

# no change, so match the first run
stopifnot(
all.equal(TempRgn.m, TempRgn.f)
)

##
## More detailed set up
##
region.contrasts <- model.matrix(~factor(CanadianWeather$region))
rgnContr3 <- region.contrasts
dim(rgnContr3) <- c(1, 35, 4)
dimnames(rgnContr3) <- list('', CanadianWeather$place, c('const',
   paste('region', c('Atlantic', 'Continental', 'Pacific'), sep='.')) )

const365 <- create.constant.basis(c(0, 365))
region.fd.Atlantic <- fd(matrix(rgnContr3[,,2], 1), const365)region.fd.Continental <- fd(matrix(rgnContr3[,,3], 1), const365)
region.fd.Pacific <- fd(matrix(rgnContr3[,,4], 1), const365)
region.fdlist <- list(const=rep(1, 35),
     region.Atlantic=region.fd.Atlantic,
     region.Continental=region.fd.Continental,
     region.Pacific=region.fd.Pacific)
beta1 <- with(Temp.fd, fd(basisobj=basis, fdnames=fdnames))
beta0 <- fdPar(beta1)
betalist <- list(const=beta0, region.Atlantic=beta0,
             region.Continental=beta0, region.Pacific=beta0)

TempRgn <- fRegress(Temp.fd, region.fdlist, betalist)

stopifnot(
all.equal(TempRgn, TempRgn.f)
)

###
###
### functional response with
###            (concurrent) functional explanatory variable
###
###

##
##  predict knee angle from hip angle;  from demo('gait', package='fda')
##
## formula interface
##
(gaittime <- as.numeric(dimnames(gait)[[1]])*20)
gaitrange <- c(0,20)
gaitbasis <- create.fourier.basis(gaitrange, nbasis=21)
harmaccelLfd <- vec2Lfd(c(0, (2*pi/20)^2, 0), rangeval=gaitrange)
gaitfd <- smooth.basisPar(gaittime, gait,
       gaitbasis, Lfdobj=harmaccelLfd, lambda=1e-2)$fd
hipfd  <- gaitfd[,1]
kneefd <- gaitfd[,2]

knee.hip.f <- fRegress(kneefd ~ hipfd)
##
## manual set-up
##
#  set up the list of covariate objects
const  <- rep(1, dim(kneefd$coef)[2])
xfdlist  <- list(const=const, hipfd=hipfd)

beta0 <- with(kneefd, fd(basisobj=basis, fdnames=fdnames))
beta1 <- with(hipfd, fd(basisobj=basis, fdnames=fdnames))

betalist  <- list(const=fdPar(beta0), hipfd=fdPar(beta1))

fRegressout <- fRegress(kneefd, xfdlist, betalist)

stopifnot(
all.equal(fRegressout, knee.hip.f)
)

#See also the following demos:

#demo('canadian-weather', package='fda')
#demo('gait', package='fda')
#demo('refinery', package='fda')
#demo('weatherANOVA', package='fda')
#demo('weatherlm', package='fda')