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

10

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

dcl
fonte
Mesmo uma descrição de como calcular as previsões manualmente a partir da base, objetos suavizados (fd) e estimativas de regressão de fRegress () seria muito útil.
dcl
Apenas verificando: você já tentou usar predict.fRegressa newdataopção (no manual da fda aqui )?
Eu tenho, é que não tenho muita certeza de qual classe ou formato os 'novos dados' devem ser. Ele não aceita um objeto fd ou fdSmooth que são as curvas suavizadas das quais desejo prever. E isso não me permite inserir argumentos brutos e valores covariáveis.
dcl
11
Lembro-me de ter um problema semelhante há cerca de um ano atrás, quando brinquei com o fdapacote. Eu estava escrevendo uma resposta que envolvia obter previsões manualmente, mas uma grande parte dela foi perdida por não salvá-la. Se alguém não me derrotar, eu devo ter uma solução para você em alguns dias.

Respostas:

14

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,

yi=β0+0TXi(s)β(s)ds+ϵi

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

β(s)=k=1Kbkθk(s)

Na notação matricial, isso é .β(s)=θ(s)b

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

Xi(s)=k=1Lcikψk(s)

Novamente, na notação matricial, esse é .X(s)=Cψ(s)

E assim, se permitirmos , nosso modelo poderá ser expresso comoJ=ψ(s)θ(s)ds

y=β0+CJb .

E se permitirmos e , nosso modelo éZ=[1CJ]ξ=[β0b]

y=Zξ

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=λ[Lβ(s)]2ds

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 comoLR

R=λ(0000R1000RK)

onde é em termos da expansão de , minimizamos a soma penalizada dos quadrados:Riβi

(yZξ)(yZξ)+λξRξ ,

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

ξ^=(ZZ+λR)1Zy .

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

y^new=Znewξ^

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). 0TXi(s)β(s)Xiβ

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.


fonte
Parece que é exatamente o que eu preciso. Obrigado. Suponho que não precisarei brincar com as coisas nfine / tine / deltat, certo? Devo assumir que a integração está sendo feita com precisão suficiente?
dcl 23/11
Além disso, notei que você não penalizou a covariável 'nova' ou as covariáveis ​​'velhas' diretamente. Tudo é feito com a penalização da versão beta (e o número de funções básicas, eu acho). A penalidade lambda é aplicada à versão beta. Você consegue o mesmo efeito penalizando as suavizações antes da regressão? (com o mesmo valor de lambda)
dcl 23/11
11
A grade usada para aproximar a integral é bastante fina, portanto a aproximação deve ser muito boa. Você sempre pode aumentar nfinee ver o quanto a integral muda, mas acho que não fará muito. Quanto à penalização, sim, estamos penalizando diretamente os em vez de neste caso. Ramsay e Silverman discutem outro método de penalidade que estima sem funções básicas em que aplicamos a penalidade diretamente em . Ambas as formas estão induzindo uma restrição de suavidade nas funções , mas não tenho certeza se você obterá o 'mesmo efeito'. ξββ^ββ
Tentei manipular o código para produzir previsões para várias curvas, mas não tenho certeza se o fiz corretamente. Para iniciantes, o yhatmat não é constante para todas as curvas após a primeira iteração do loop ... Isso significa que é equivalente a ? β0
dcl 23/11
11
@dcl No loop, quando , ele adiciona a (assumindo que a primeira lista na sua lista X corresponda ao termo de interceptação). Você pode adicionar o trecho de código que está usando à sua pergunta para que eu possa ver? j=1β0^y^