学习函数回归模型，关注r语言代码：This is the first of two chapters on the functional linear model. Here we have a dependent or response variable whose value is to be predicted or approximated on the basis of a set of independent or covariate variables, and at least one of these is functional in nature. The focus here is on linear models, or functional analogues of linear regression analysis. This chapter is confined to considering the prediction of a scalar response on the basis of one or more functional covariates, as well as possible scalar covariates.
Confidence intervals are developed for estimated regression functions in order to permit conclusions about where along the t axis a covariate plays a strong role in predicting a functional responses. The chapter also offers some permutation tests of hypotheses.
More broadly, we begin here the study of input/output systems. This and the next chapter lead in to Chapter 11, where the response is the derivative of the output from the system.
9.1	Functional Linear Regression with a Scalar Response
We have so far focused on representing a finite number of observations as a func-tional data object, which can theoretically be evaluated at an infinite number of points, and on graphically exploring the variation and covariation of populations of functions. However, we often need to model predictive relationships between differ-ent types of data, and here we expect that some of these data will be functional.
In classic linear regression, predictive models are often of the form
p
yi = XxijBj + a, i= 1,…,N	(9.1)
j=0
that model the relationship between a response yi, and covariates xij as a linear structure. The dummy covariate with j = 0, which has the value one for all i, is
J.O. Ramsay et al., Functional Data Analysis with R and MATLAB, Use R,
DOI: 10.1007/978-0-387-98185-7_9,
© Springer Science + Business Media, LLC 2009 usually included because origin of the response variable and/or one or more of the independent variables can be arbitrary, and 00 codes the constant needed to allow for this. It is often called the intercept term.
The term & allows for sources of variation considered extraneous, such as mea-surement error, unimportant additional causal factors, sources of nonlinearity and so forth, all swept under the statistical rug called error. The a are assumed to add to the prediction of the response, and are usually considered to be independently and identically distributed.
In this chapter we replace at least one of the p observed scalar covariate vari¬ables on the right side of the classic equation by a functional covariate. To simplify the exposition, though, we will describe a model consisting of a single functional independent variable plus an intercept term.
9.2	A Scalar Response Model for Log Annual Precipitation
Our test-bed problem in this section is to predict the logarithm of annual precipita-tion for 35 Canadian weather stations from their temperature profiles. The response in this case is, in terms of the fda package in R,
annualprec = log10(apply(daily$precav,2,sum))
We want to use as the predictor variable the complete temperature profile as well as a constant intercept. These two covariates can be stored in a list of length 2, or in Matlab as a cell array. Here we set up a functional data object for the 35 temperature profiles, called tempfd. To keep things simple and the computation rapid, we will use 65 basis functions without a roughness penalty. This number of basis functions has been found to be adequate for most purposes, and can, for example, capture the ripples observed in early spring in many weather stations.
tempbasis =create.fourier.basis(c(0,365),65)
tempSmooth=smooth.basis(day.5,daily$tempav,tempbasis) tempfd =tempSmooth$fd
9.3	Setting Up the Functional Linear Model
But what can we do when the vector of covariate observations x i = (x〃,..., Xip) in (9.1) is replaced by a function Xi(t)? A first idea might be to discretize each of the N functional covariates Xi (t) by choosing a set of times 11,..., tq and consider fitting the model
q
yi = «O + f Xi(tj )0j + &i.
j=1 
But which times tj are important, given that we must have q < N?
If we choose a finer and finer mesh of times, the summation approaches an inte-gral equation:

We now have a finite number N of observations with which to determine the infinite-dimensional B(t). This is an impossible problem: it is almost always possible to find a B(t) so that (9.2) is satisfied with Ei = 0. More importantly, there are always an infinite number of possible regression coefficient functions B(t) that will produce exactly the same predictions yi. Even if we expand each functional covariate in terms of a limited number of basis functions, it is entirely possible that the total number of basis functions will exceed or at least approach N.
9.4	Three Estimates of the Regression Coefficient Predicting Annual Precipitation
Three strategies have been developed to deal with this underdetermination issue. The first two redefine the problem using a basis coefficient expansion of 0:
K
0(t) = £ Ck帆 t) = c0 Mt).	(9.3)
k
The third replaces the potentially high-dimensional covariate functions by a low-dimensional approximation using principal components analysis. The first two ap-proaches will be illustrated using function fRegress. Function fRegress in R and Matlab requires at least three arguments:
yfdPar This object contains the response variable. It can be a functional param¬eter object, a functional data object, or simply a vector of N scalar responses. In this chapter we restrict ourselves to the scalar response situation.
xfdlist This object contains all the functional and scalar covariate functions used to predict the response using the linear model. Each covariate is an element or component in a list object in R or an entry in a cell array in Matlab.
betalist This is a list object in R or a cell array in Matlab of the same length as the second argument; it specifies the functional regression coefficient objects. Because it is possible that any or all of them can be subject to a roughness penalty, they are all in principle functional parameter objects, although fRegress will by default convert both functional data objects and basis objects into functional pa-rameter objects.
Here we store the two functional data covariates required for predicting log an-nual precipitation in a list of length two, which we here call templist, to be used for the argument xfdlist.
templist	=	vector("list",2)
templist[[1]] = rep(1,35)
templist[[2]] = tempfd
Notice that the intercept term is here set up as a constant function with 35 repli-cations.
9.4.1	Low-Dimensional Regression Coefficient Function p
The simplest strategy for estimating p is just to keep the dimensionality K of p in (9.3) small relative to N. In our test-bed expansion, we will work with five Fourier basis functions for the regression coefficient p multiplying the temperature profiles; we will also use a constant function for a, the multiplier of the constant intercept covariate set up above.
conbasis	= create.constant.basis(c(0,365))
betabasis = create.fourier.basis(c(0,365),5) betalist	= vector("list",2)
betalist[[1]] = conbasis
betalist[[2]] = betabasis
Now we can call function fRegress, which returns various results in a list object that we call fRegressList.
fRegressList = fRegress(annualprec,templist,betalist)
The command names(fRegressList) reveals a component betaestlist containing the estimated regression coefficient functions. Each of these is a func¬tional parameter object. We can plot the estimate of the regression function for the temperature profiles with the commands
betaestlist = fRegressList$betaestlist tempbetafd = betaestlist[[2]]$fd plot(tempbetafd, xlab="Day", ylab="Beta for temperature")
Figure 9.1 shows the result. The intercept term can be obtained from coef( betaestlist[[1]]); its value in this case is 0.0095. We will defer commenting on these estimates until we consider the next more sophisticated strategy.
We need to assess the quality of this fit. First, let us extract the fitted values defined by this model and compute the residuals. We will also compute error sums of squares associated with the fit as well as for the fit using only a constant or intercept.
annualprechat1 = fRegressList$yhatfdobj annualprecres1 = annualprec - annualprechatl SSE1.1 = sum(annualprecres1"2)
SSE0 = sum((annualprec - mean(annualprec))"2)
 
0	100	200	300
Day
Fig. 9.1 Estimated p(t) for predicting log annual precipitation from average daily temperature using five Fourier basis functions.

We can now compute the squared multiple correlation and the usual F-ratio for comparing these two fits.
RSQ1 = (SSE0-SSE1.1)/SSE0
Fratiol = ((SSE0-SSE1)/5)/(SSE1/29)
The squared multiple correlation is 0.80, and the corresponding F-ratio with 5 and 29 degrees of freedom is 22.6, suggesting a fit to the data that is far better than we would expect by chance.
9.4.2 Coefficient p Estimate Using a Roughness Penalty
There are two ways to obtain a smooth fit. The simplest is to use a low-dimensional basis for P(t). However, we can get more direct control over what we mean by “smooth” by using a roughness penalty. The combination of a high-dimensional basis with a roughness penalty reduces the possibilities that either (a) important features are missed or (b) extraneous features are forced into the image by using a basis set that is too small for the application.
Suppose, for example, that we fit (9.2) by minimizing the penalized sum of squares
PENSSEa(an,p) = £加一a - ZXi(t)p(t)dt]2 + Z Z[Lp(t)]2dt.	(9.4) This allows us to shrink variation in 0 as close as we wish to the solution of the dif-ferential equation L0 = 0. Suppose, for example, that we are working with periodic data with a known period. As noted with expression (5.11), the use of a harmonic acceleration operator,
L0 = (rn 2) D0 + D3 0,
places no penalty on a simple sine wave and increases the penalty on higher-order harmonics in a Fourier approximation approximately in proportion to the sixth power of the order of the harmonic. (In this expression,① is determined by the pe-riod, which is assumed to be known.) Thus, increasing the penalty 入 in (9.4) forces 0 to look more and more like 0 (t) = c 1 + c2 sin(①t) + c3 cos(①t).
More than one functional covariate can be incorporated into this model and scalar covariates may also be included. Let us suppose that in addition to yi we have measured p scalar covariates zi = (zi 1,...,zip) and q functional covariates Xi 1 (t),...,xiq(t). We can put these into a linear model as follows
q ,
a0 + 4a + 工J Xij(t)0j(t)dt + a,
where zi is the p - vector of scalar covariates. A separate smoothing penalty may be employed for each of the 0j (t), j = 1,..., q.
Using (9.5), we can define a least-squares estimate as follows. We define Z by:
—JX 11(t画(t)dt ••• XX1 q(t鸣(t)-
Z =.	：	.
.	.
_Z'nRXn 1(t)01(t)dt …RXnq(t陶(t)_
Here 0k is the basis expansion used to represent 0k(t). We also define a penalty matrix
0 	
0 ZIR1 •••	0
R(A)= .	.	.	.	(9.6)
.	.	.	.
一 0…3
where Rk is the penalty matrix associated with the smoothing penalty for 0k and Zk is the corresponding smoothing parameter. With these objects, we can define
b = °Z0Z + R(Z))-1 Z0y
to hold the vector of estimated coefficients a along with the coefficients defining each estimated coefficient function 0k (t) estimated by penalized least squares. These are then extracted to form the appropriate functional data objects.
Now let us apply this approach to predicting the log annual precipitations. First, we set up a harmonic acceleration operator, as we did already in Chapter 5.
Lcoef = c(0,(2*pi/365)-2,0)
harmaccelLfd = vec2Lfd(Lcoef, c(0,365))
Now we replace our previous choice of basis for defining the 0 estimate by a func-tional parameter object that incorporates both this roughness penalty and a level of smoothing:
betabasis = create.fourier.basis(c(0, 365), 35) lambda = 10"12.5
betafdPar = fdPar(betabasis, harmaccelLfd, lambda) betalist[[2]] = betafdPar
These now allow us to invoke fRegress to return the estimated functional coeffi-cients and predicted values:
annPrecTemp = fRegress(annualprec, templist, betalist)
betaestlist2 = annPrecTemp$betaestlist annualprechat2 = annPrecTemp$yhatfdobj
The command print(fRegressList$df) indicates the degrees of freedom for this model, including the intercept, is 4.7, somewhat below the value of 6 that we used for the simple model above.
Now we compute the usual R2 and F-ratio statistics to assess the improvement in fit achieved by including temperature as a covariate.
SSE1.2	= sum((annualprec-annualprechat2)"2)
RSQ2	= (SSE0 - SSE1.2)/SSE0
Fratio2 = ((SSE0-SSE1.2)/3.7)/(SSE1/30.3)
The squared multiple correlation is now 0.75, a small drop from the value for the simple model, due partly to using fewer degrees of freedom. The F-ratio is 25.1 with 3.7 and 30.3 degrees of freedom, and is even more significant than for the simple model. The reader should note that because a smoothing penalty has been used, the F-distribution only represents an approximation to the null distribution for this model. Figure 9.2 compares predicted and observed values of log annual pre-cipitation. Figure 9.3 plots the coefficient 0(t) along with the confidence intervals derived below. Comparing this version with that in Figure 9.1 shows why the rough-ness penalty approach is to be preferred over the fixed low-dimension strategy; now we see that only the autumn months really matter in defining the relationship, and that the substantial oscillations over other parts of the year in Figure 9.1 are actually extraneous.
To complete the picture, we should ask whether we could do just as well with a constant value for 0. Here we use the constant basis, run fRegress, and redo the comparison using this fit as a benchmark. The degrees of freedom for this model is now 2.
betalist[[2]] = fdPar(conbasis)
fRegressList = fRegress(annualprec, templist, betalist)
betaestlist = fRegressList$betaestlist