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} = β_{1} z_{i} + \int x_{i} (t) β_{2} (t) d t + e_{i}

and the concurrent functional dependent variable model

y_{i} (t) = β_{1} (t) z_{i} + β_{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, ...)
## S3 method for class 'fdPar':
fRegress(y, xfdlist, betalist, wt=NULL,
                     y2cMap=NULL, SigmaE=NULL, ...)
## S3 method for class 'numeric':
fRegress(y, xfdlist, betalist, wt=NULL,
                     y2cMap=NULL, SigmaE=NULL, ...)

Arguments

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
...
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.

Examples

Run this code

###
###
### 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')

#  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)

# Now do it.
precip.Temp.m <- do.call('fRegress', precip.Temp.mdl)

# With the change in betalist, the answers may be different
# Without that change, the answes will be the same.
stopifnot(
all.equal(precip.Temp.m, precip.Temp.f)
)

##
## 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')

Run the code above in your browser using DataLab

Last chance! 50% off unlimited learning