Predicting response from new curves using fda package in R

Basically all I want to do is predict a scalar response using some curves. I've got as far as doing a regression (using fRegress from the fda package) but have no idea how to apply the results to a NEW set of curves (for prediction).

I have N=536 curves, and 536 scalar responses. Here's what I've done so far:

• I've created a basis for the curves.
• I've created a fdPar object to introduce a penalty
• I've created the fd object using smooth.basis to smooth the curves with the chosen penalty on the specified basis.
• I've ran a regression using fRegress(), regressing the curves on the scalar response.

Now, all I'd like to do, is use that regression to produce predictions for a new set of data that I have. I can't seem to find an easy way to do this.

Cheers

-
Even a description of how to calculate the predictions manually from the basis, smoothed (fd) objects and the regression estimates from fRegress() would be very helpful. – dcl Nov 18 '11 at 5:52
Just checking: have you tried using predict.fRegress using the newdata option (from the fda manual here)? – Mike Wierzbicki Nov 18 '11 at 6:33
I have, it's just that I'm not exactly sure what class or format the 'newdata' has to be. It won't accept an fd or an fdSmooth object which are the smoothed curves I wish to predict from. And it won't let me input raw arguments and covariate values. – dcl Nov 18 '11 at 9:08
I remember having a similar problem about a year ago when I played around with the fda package. I was writing up a response that involved getting predictions manually, but a big chunk of it was lost due to not saving it. If someone else doesn't beat me to it, I should have a solution for you in a couple days. – Mike Wierzbicki Nov 20 '11 at 1:56
I would be most grateful Mike. Thank you :) – dcl Nov 20 '11 at 13:42

I don't care for fda's use of Inception-like list-within-list-within-list object structures, but my response will abide by the system the package writers have created.

I think it's instructive to first think about what we're doing exactly. Based on your description of what you've done so far, this is what I believe you're doing (let me know if I have misinterpreted something). I'll be continuing using notation and, due to lack of real data, an example from Ramsay and Silverman's Functional Data Analysis and Ramsay, Hooker, and Graves's Functional Data Analysis with R and MATLAB (Some of the following equations and code are lifted directly from these books).

We are modelling a scalar response via a functional linear model, i.e.

$y_i = \beta_0 + \int_0^T X_i(s)\beta(s)ds + \epsilon_i$

We expand the $\beta$ in some basis. We use, say, $K$ basis functions. So,

$\beta(s) = \sum\limits_{k=1}^K b_k\theta_k (s)$

In matrix notation, this is $\beta (s)=\boldsymbol{\theta}'(s)\mathbf{b}$.

We also expand the covariate functions in some basis, as well (say $L$ basis functions). So,

$X_i(s) = \sum\limits_{k=1}^L c_{ik}\psi_k (s)$

Again, in matrix notation, this is $X(s)=\mathbf{C}\boldsymbol{\psi}(s)$.

And thus, if we let $\mathbf{J}=\int\boldsymbol{\psi}(s)\boldsymbol{\theta}'(s)ds$, our model can be expressed as

$y = \beta_0 + \mathbf{CJb}$.

And if we let $\mathbf{Z} = [\mathbf{1}\quad \mathbf{CJ}]$ and $\boldsymbol{\xi} = [\beta_0\quad \mathbf{b}']'$, our model is

$\mathbf{y} = \mathbf{Z}\boldsymbol{\xi}$

And this looks much more familiar to us.

Now I see you're adding some sort of regularization. The fda package works with roughness penalties of the form

$P = \lambda\int[L\beta(s)]^2 ds$

for some linear differential operator $L$. Now it can be shown (details left out here -- it's really not hard to show this) that if we define the penalty matrix $\mathbf{R}$ as

$\mathbf{R}= \lambda\left(\begin{array}{cccc} 0 & 0 & \cdots & 0 \\ 0 & R_1 & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & R_K \end{array}\right)$

where $R_i$ is in terms of the basis expansion of $\beta_i$, then we minimize the penalized sum of squares:

$\left(\mathbf{y}-\mathbf{Z}\boldsymbol{\xi}\right)'\left(\mathbf{y}-\mathbf{Z}\boldsymbol{\xi}\right) + \lambda\boldsymbol{\xi}'\mathbf{R}\boldsymbol{\xi}$,

and so our problem is merely a ridge regression with solution:

$\hat{\boldsymbol{\xi}}=(\mathbf{Z}'\mathbf{Z}+\lambda\mathbf{R})^{-1}\mathbf{Z}'\mathbf{y}$.

I walked through the above because, (1) I think it's important we understand what we're doing, and (2) some of the above is necessary to understand some of the code I'll use later on. On to the code...

Here's a data example with R code. I'm using the Canadian weather dataset provided in the fda package. We'll model the log annual precipitation for a number of weather stations via a functional linear model and we'll use temperature profiles (temperatures were recorded once a day for 365 days) from each station as functional covariates. We'll proceed similarly to the way you describe in your situation. Data was recorded at 35 stations. I'll break the dataset up into 34 stations, which will be used as my data, and the last station, which will be my "new" dataset.

I continue via R code and comments (I assume you are familiar enough with the fda package such that nothing in the following is too surprising -- if this isn't the case, please let me know):

# pick out data and 'new data'
dailydat <- daily$precav[,2:35] dailytemp <- daily$tempav[,2:35]
dailydatNew <- daily$precav[,1] dailytempNew <- daily$tempav[,1]

# set up response variable
annualprec <- log10(apply(dailydat,2,sum))

# create basis objects for and smooth covariate functions
tempbasis <- create.fourier.basis(c(0,365),65)
tempSmooth <- smooth.basis(day.5,dailytemp,tempbasis)
tempfd <- tempSmooth$fd # create design matrix object templist <- vector("list",2) templist[[1]] <- rep(1,34) templist[[2]] <- tempfd # create constant basis (for intercept) and # fourier basis objects for remaining betas conbasis <- create.constant.basis(c(0,365)) betabasis <- create.fourier.basis(c(0,365),35) betalist <- vector("list",2) betalist[[1]] <- conbasis betalist[[2]] <- betabasis # set roughness penalty for betas Lcoef <- c(0,(2*pi/365)^2,0) harmaccelLfd <- vec2Lfd(Lcoef, c(0,365)) lambda <- 10^12.5 betafdPar <- fdPar(betabasis, harmaccelLfd, lambda) betalist[[2]] <- betafdPar # regress annPrecTemp <- fRegress(annualprec, templist, betalist)  Now when I was first taught about functional data a year or so ago, I played around with this package. I was also unable to get predict.fRegress to give me what I wanted. Looking back at it now, I still don't know how to make it behave. So, we'll just have to get the predictions semi-manually. I'll be using pieces that I pulled straight out of the code for fRegress(). Again, I continue via code and comments. First, the set-up: # create basis objects for and smooth covariate functions for new data tempSmoothNew <- smooth.basis(day.5,dailytempNew,tempbasis) tempfdNew <- tempSmoothNew$fd

# create design matrix object for new data
templistNew <- vector("list",2)
templistNew[[1]] <- rep(1,1)
templistNew[[2]] <- tempfdNew

# convert the intercept into an fd object
onebasis <- create.constant.basis(c(0,365))
templistNew[[1]] <- fd(matrix(templistNew[[1]],1,1), onebasis)


Now to get the predictions

$\hat{\mathbf{y}}_{\mathrm{new}}=\mathbf{Z}_{\mathrm{new}}\hat{\boldsymbol{\xi}}$

I just take the code that fRegress uses to calculate yhatfdobj and edit it slightly. fRegress calculates yhatfdobj by estimating the integral $\int_0^T X_i (s) \beta(s)$ via the trapezoid rule (with $X_i$ and $\beta$ expanded in their respective bases).

Normally, fRegress calculates the fitted values by looping through the covariates stored in annPrecTemp$xfdlist. So for our problem, we replace this covariate list with the corresponding one in our new covariate list, i.e., templistNew. Here's the code (identical to the code found in fRegress with two edits, some deletions of unneeded code, and a couple comments added): # set up yhat matrix (in our case it's 1x1) yhatmat <- matrix(0,1,1) # loop through covariates p <- length(templistNew) for(j in 1:p){ xfdj <- templistNew[[j]] xbasis <- xfdj$basis
xnbasis    <- xbasis$nbasis xrng <- xbasis$rangeval
nfine      <- max(501,10*xnbasis+1)
tfine      <- seq(xrng[1], xrng[2], len=nfine)
deltat     <- tfine[2]-tfine[1]
xmat       <- eval.fd(tfine, xfdj)
betafdParj <- annPrecTemp$betaestlist[[j]] betafdj <- betafdParj$fd
betamat    <- eval.fd(tfine, betafdj)
# estimate int(x*beta) via trapezoid rule
fitj       <- deltat*(crossprod(xmat,betamat) -
0.5*(outer(xmat[1,],betamat[1,]) +
outer(xmat[nfine,],betamat[nfine,])))
yhatmat    <- yhatmat + fitj
}


(note: if you look at this chunk and surrounding code in fRegress, you'll see the steps I outlined above).

I tested the code by re-running the weather example using all 35 stations as our data and compared the output from the above loop to annPrecTemp$yhatfdobj and everything matches up. I also ran it a couple times using different stations as my "new" data and everything seems reasonable. Let me know if any of the above is unclear or if anything is not working correctly. Sorry for the overly detailed response. I couldn't help myself :) And if you don't already own them, check out the two books I used to write up this response. They are really good books. - This looks like it's exactly what I need. Thank you. I assume I won't have to play around with the nfine/tine/deltat stuff right? Should I assume the integration is being doing accurately enough? – dcl Nov 23 '11 at 7:08 Also, I notice you didn't penalise the 'new' covariate or the 'old' covariates directly. It's all done with penalizing beta (and the number of basis functions I guess). The penalty lambda is applied to the beta. Do you achieve the same effect by penalizing the smoothes before regression? (with the same value of lambda) – dcl Nov 23 '11 at 12:29 The grid used to approximate the integral is pretty fine, so the approximation should be pretty good. You could always increase nfine and see how much the integral changes but I'm guessing it won't do much. As far as the penalization, yes we're directly penalizing the$\xi$'s instead of$\beta$'s in this case. Ramsay and Silverman discuss another penalty method that estimates$\hat{\beta}$without basis functions where we apply the penalty directly to$\beta$. Both ways are inducing a smoothness constraint on the$\beta$functions, but I'm unsure whether you'll get the 'same effect.' – Mike Wierzbicki Nov 23 '11 at 14:47 I've tried manipulating the code to produce predictions for multiple curves, but I'm not sure I've done it correctly. For starters, yhatmat is not constant for all curvers after the first iteration of the loop... Is this this meant to be equivalent to$\beta_0$? – dcl Nov 23 '11 at 15:50 @dcl In the loop, when$j=1$, it is adding$\hat{\beta_0}$to$\hat{y}\$ (assuming the first list in your Xlist corresponds to the intercept term). Can you add the snippet of code you're using to your question so I can look at it? – Mike Wierzbicki Nov 23 '11 at 16:49
{"extraction_info": {"found_math": true, "script_math_tex": 0, "script_math_asciimath": 0, "math_annotations": 0, "math_alttext": 0, "mathml": 0, "mathjax_tag": 0, "mathjax_inline_tex": 1, "mathjax_display_tex": 0, "mathjax_asciimath": 1, "img_math": 0, "codecogs_latex": 0, "wp_latex": 0, "mimetex.cgi": 0, "/images/math/codecogs": 0, "mathtex.cgi": 0, "katex": 0, "math-container": 0, "wp-katex-eq": 0, "align": 0, "equation": 0, "x-ck12": 0, "texerror": 0, "math_score": 0.5892781019210815, "perplexity": 1208.7475873425417}, "config": {"markdown_headings": false, "markdown_code": true, "boilerplate_config": {"ratio_threshold": 0.18, "absolute_threshold": 10, "end_threshold": 15, "enable": true}, "remove_buttons": true, "remove_image_figures": true, "remove_link_clusters": true, "table_config": {"min_rows": 2, "min_cols": 3, "format": "plain"}, "remove_chinese": true, "remove_edit_buttons": true, "extract_latex": true}, "warc_path": "s3://commoncrawl/crawl-data/CC-MAIN-2015-48/segments/1448398463315.2/warc/CC-MAIN-20151124205423-00129-ip-10-71-132-137.ec2.internal.warc.gz"}