Prediciendo la respuesta de nuevas curvas usando el paquete fda en R

Básicamente, todo lo que quiero hacer es predecir una respuesta escalar usando algunas curvas. Llegué a hacer una regresión (usando fRegress del paquete fda) pero no tengo idea de cómo aplicar los resultados a un NUEVO conjunto de curvas (para la predicción).

Tengo N = 536 curvas y 536 respuestas escalares. Esto es lo que he hecho hasta ahora:

• He creado una base para las curvas.
• He creado un objeto fdPar para introducir una penalización.
• He creado el objeto fd usando smooth.basis para suavizar las curvas con la penalización elegida sobre la base especificada.
• He ejecutado una regresión usando fRegress (), regresando las curvas en la respuesta escalar.

Ahora, todo lo que me gustaría hacer es usar esa regresión para producir predicciones para un nuevo conjunto de datos que tengo. Parece que no puedo encontrar una manera fácil de hacer esto.

Salud

No me importa para fdasu uso 's de Inception -como las estructuras objeto de lista dentro de otra lista dentro de otra lista, pero mi respuesta a acatar el sistema de los escritores han creado paquetes.

Creo que es instructivo pensar primero en lo que estamos haciendo exactamente. Según su descripción de lo que ha hecho hasta ahora, esto es lo que creo que está haciendo (avíseme si he malinterpretado algo). Continuaré usando la notación y, debido a la falta de datos reales, un ejemplo del Análisis de datos funcionales de Ramsay y Silverman y el Análisis de datos funcionales de Ramsay, Hooker y Graves con R y MATLAB (Algunas de las siguientes ecuaciones y códigos se eliminan directamente de estos libros).

Estamos modelando una respuesta escalar a través de un modelo lineal funcional, es decir

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

Expandimos el de alguna manera. Usamos, digamos, funciones básicas. Entonces,$\beta$$K$

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

En notación matricial, esto es .$\beta (s)=\boldsymbol{\theta}'(s)\mathbf{b}$

También expandimos las funciones covariables en alguna base, también (digamos funciones base ). Entonces,$L$

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

Nuevamente, en notación matricial, esto es .$X(s)=\mathbf{C}\boldsymbol{\psi}(s)$

Y así, si dejamos , nuestro modelo puede expresarse como$\mathbf{J}=\int\boldsymbol{\psi}(s)\boldsymbol{\theta}'(s)ds$

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

Y si dejamos que y , nuestro modelo esξ = [ β $\mathbf{Z} = [\mathbf{1}\quad \mathbf{CJ}]$$\boldsymbol{\xi} = [\beta_0\quad \mathbf{b}']'$

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

Y esto nos parece mucho más familiar.

Ahora veo que estás agregando algún tipo de regularización. El fdapaquete funciona con penalidades de rugosidad de la forma.

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

para algunos lineal operador diferencial . Ahora se puede mostrar (detalles aquí, no es difícil mostrar esto) que si definimos la matriz de penalización como$L$$\mathbf{R}$

$\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)$

donde es en términos de la expansión base de , entonces minimizamos la suma penalizada de cuadrados:$R_i$$\beta_i$

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

y entonces nuestro problema es simplemente una regresión de cresta con solución:

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

Revisé lo anterior porque, (1) creo que es importante que comprendamos lo que estamos haciendo, y (2) parte de lo anterior es necesario para comprender parte del código que usaré más adelante. En el código ...

Aquí hay un ejemplo de datos con código R. Estoy usando el conjunto de datos meteorológicos canadienses que se proporciona en el fdapaquete. Modelaremos el registro de precipitación anual para varias estaciones meteorológicas mediante un modelo lineal funcional y utilizaremos perfiles de temperatura (las temperaturas se registraron una vez al día durante 365 días) de cada estación como covariables funcionales. Vamos a proceder de manera similar a la forma en que usted describe su situación. Los datos se registraron en 35 estaciones. Dividiré el conjunto de datos en 34 estaciones, que se utilizarán como mis datos, y la última estación, que será mi "nuevo" conjunto de datos.

Continúo a través del código R y los comentarios (supongo que está lo suficientemente familiarizado con el fdapaquete de modo que nada de lo siguiente es demasiado sorprendente; si este no es el caso, hágamelo saber):

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

Ahora, cuando me enseñaron por primera vez sobre datos funcionales hace aproximadamente un año, jugué con este paquete. Tampoco pude predict.fRegressdarme lo que quería. Mirando hacia atrás ahora, todavía no sé cómo hacer que se comporte. Entonces, solo tendremos que obtener las predicciones semi-manualmente. Usaré piezas que saqué directamente del código fRegress(). Nuevamente, continúo por código y comentarios.

Primero, la configuración:

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

Ahora para obtener las predicciones

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

Solo tomo el código que fRegressusa para calcularlo yhatfdobjy editarlo ligeramente. fRegresscalcula yhatfdobjestimando la integral través de la regla trapezoide (con y expandidos en sus respectivas bases). X i$\int_0^T X_i (s) \beta(s)$$X_i$$\beta$

Normalmente, fRegresscalcula los valores ajustados recorriendo las covariables almacenadas en annPrecTemp$xfdlist. Así que para nuestro problema, reemplazamos esta lista de covarianza con la que corresponde en nuestra nueva lista de covarianza, es decir, templistNew. Aquí está el código (idéntico al código que se encuentra fRegresscon dos ediciones, algunas eliminaciones de código innecesario y un par de comentarios agregados):

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

(nota: si observa este fragmento y el código circundante fRegress, verá los pasos que describí anteriormente).

Probé el código volviendo a ejecutar el ejemplo del clima usando las 35 estaciones como nuestros datos y comparé la salida del bucle anterior annPrecTemp$yhatfdobjy todo coincide. También lo ejecuté un par de veces usando diferentes estaciones como mis "nuevos" datos y todo parece razonable.

Avíseme si alguno de los puntos anteriores no está claro o si algo no funciona correctamente. Perdón por la respuesta demasiado detallada. No pude evitarlo :) Y si aún no los posee, consulte los dos libros que usé para escribir esta respuesta. Son realmente buenos libros.