Prevendo a resposta de novas curvas usando o pacote fda em R

Basicamente, tudo o que quero fazer é prever uma resposta escalar usando algumas curvas. Eu tenho tanto quanto fazer uma regressão (usando fRegress do pacote fda), mas não tenho idéia de como aplicar os resultados a um novo conjunto de curvas (para previsão).

Eu tenho N = 536 curvas e 536 respostas escalares. Aqui está o que eu fiz até agora:

• Eu criei uma base para as curvas.
• Eu criei um objeto fdPar para introduzir uma penalidade
• Eu criei o objeto fd usando smooth.basis para suavizar as curvas com a penalidade escolhida na base especificada.
• Fiz uma regressão usando fRegress (), regredindo as curvas na resposta escalar.

Agora, tudo que eu gostaria de fazer é usar essa regressão para produzir previsões para um novo conjunto de dados que eu tenho. Não consigo encontrar uma maneira fácil de fazer isso.

Felicidades

Eu não ligo para o fdauso 's de Inception -como estruturas de objetos lista-dentro-de-lista-dentro-lista, mas a minha resposta vai respeitar o sistema os escritores pacote criaram.

Acho instrutivo pensar primeiro no que estamos fazendo exatamente. Com base na sua descrição do que você fez até agora, é isso que acredito que esteja fazendo (deixe-me saber se interpretei algo errado). Continuarei usando a notação e, devido à falta de dados reais, um exemplo da Análise de Dados Funcionais de Ramsay e Silverman e Análise de Dados Funcionais de Ramsay, Hooker e Graves com R e MATLAB (Algumas das equações e códigos a seguir são levantados diretamente desses livros).

Estamos modelando uma resposta escalar por meio de um modelo linear funcional, ou seja,

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

Expandimos o de alguma forma. Usamos, digamos, funções básicas. Então,$\beta$$K$

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

Na notação matricial, isso é .$\beta (s)=\boldsymbol{\theta}'(s)\mathbf{b}$

Também expandimos as funções covariáveis ​​de alguma maneira (digamos , funções básicas de ). Então,$L$

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

Novamente, na notação matricial, esse é .$X(s)=\mathbf{C}\boldsymbol{\psi}(s)$

E assim, se permitirmos , nosso modelo poderá ser expresso como$\mathbf{J}=\int\boldsymbol{\psi}(s)\boldsymbol{\theta}'(s)ds$

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

E se permitirmos e , nosso modelo é$\mathbf{Z} = [\mathbf{1}\quad \mathbf{CJ}]$$\boldsymbol{\xi} = [\beta_0\quad \mathbf{b}']'$

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

E isso parece muito mais familiar para nós.

Agora vejo que você está adicionando algum tipo de regularização. O fdapacote funciona com penalidades de rugosidade do formulário

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

para alguns operador diferencial linear . Agora pode ser mostrado (detalhes deixados aqui - não é realmente difícil mostrar isso) que se definirmos a matriz de penalidade 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)$

onde é em termos da expansão de , minimizamos a soma penalizada dos quadrados:$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}$ ,

e, portanto, nosso problema é apenas uma regressão de crista com solução:

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

Eu caminhei pelo exposto porque, (1) acho que é importante entendermos o que estamos fazendo e (2) alguns dos itens acima são necessários para entender parte do código que usarei mais tarde. Para o código ...

Aqui está um exemplo de dados com código R. Estou usando o conjunto de dados climáticos canadense fornecido no fdapacote. Modelaremos a precipitação anual de log para várias estações meteorológicas por meio de um modelo linear funcional e usaremos perfis de temperatura (as temperaturas foram registradas uma vez por dia durante 365 dias) de cada estação como covariáveis ​​funcionais. Vamos proceder da mesma forma que você descreve em sua situação. Os dados foram gravados em 35 estações. Dividirei o conjunto de dados em 34 estações, que serão usadas como meus dados, e a última estação, que será meu "novo" conjunto de dados.

Continuo via código R e comentários (presumo que você esteja familiarizado o suficiente com o fdapacote, de modo que nada a seguir seja muito surpreendente - se esse não for o caso, informe-me):

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

Agora, quando eu fui ensinado sobre dados funcionais há mais ou menos um ano, brinquei com este pacote. Também não consegui predict.fRegressme dar o que queria. Olhando para trás agora, ainda não sei como fazê-lo se comportar. Portanto, teremos que obter as previsões semi-manualmente. Vou usar peças que extraí direto do código fRegress(). Novamente, continuo via código e comentários.

Primeiro, a configuração:

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

Agora, para obter as previsões

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

Eu apenas pego o código que fRegressusa para calcular yhatfdobje editar um pouco. fRegresscalcula yhatfdobjestimando a integral através da regra trapezoidal (com e expandidos em suas respectivas bases). $\int_0^T X_i (s) \beta(s)$$X_i$$\beta$

Normalmente, fRegresscalcula os valores ajustados percorrendo as covariáveis ​​armazenadas em annPrecTemp$xfdlist. Portanto, para o nosso problema, substituímos essa lista covariável pela correspondente em nossa nova lista covariável, ou seja templistNew,. Aqui está o código (idêntico ao encontrado em fRegressduas edições, algumas exclusões de código desnecessário e alguns comentários adicionados):

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

(observação: se você olhar para este bloco e o código ao redor fRegress, verá as etapas descritas acima).

Testei o código reexecutando o exemplo climático usando todas as 35 estações como nossos dados e comparando a saída do loop acima annPrecTemp$yhatfdobje tudo corresponde. Também executei algumas vezes usando estações diferentes como meus "novos" dados e tudo parece razoável.

Deixe-me saber se alguma das opções acima não está clara ou se algo não está funcionando corretamente. Desculpe pela resposta excessivamente detalhada. Eu não pude evitar :) E se você ainda não os possui, confira os dois livros que eu escrevi para esta resposta. São livros realmente bons.