Prédire la réponse de nouvelles courbes à l'aide du package fda dans R

10

Fondamentalement, tout ce que je veux faire, c'est prédire une réponse scalaire à l'aide de certaines courbes. Je suis allé jusqu'à faire une régression (en utilisant fRegress du package fda) mais je n'ai aucune idée de comment appliquer les résultats à un NOUVEAU jeu de courbes (pour la prédiction).

J'ai N = 536 courbes et 536 réponses scalaires. Voici ce que j'ai fait jusqu'à présent:

  • J'ai créé une base pour les courbes.
  • J'ai créé un objet fdPar pour introduire une pénalité
  • J'ai créé l'objet fd en utilisant smooth.basis pour lisser les courbes avec la pénalité choisie sur la base spécifiée.
  • J'ai effectué une régression en utilisant fRegress (), régressant les courbes sur la réponse scalaire.

Maintenant, tout ce que j'aimerais faire, c'est utiliser cette régression pour produire des prédictions pour un nouvel ensemble de données que j'ai. Je n'arrive pas à trouver un moyen facile de le faire.

À votre santé

dcl
la source
Même une description de la façon de calculer manuellement les prévisions à partir de la base, des objets lissés (fd) et des estimations de régression de fRegress () serait très utile.
DCL
Juste vérification: avez-vous essayé d'utiliser predict.fRegressl' newdataoption (du manuel de la FDA ici )?
J'ai, c'est juste que je ne suis pas exactement sûr de la classe ou du format des 'newdata'. Il n'acceptera pas un objet fd ou fdSmooth qui sont les courbes lissées à partir desquelles je souhaite prédire. Et cela ne me permettra pas de saisir des arguments bruts et des valeurs de covariables.
DCL
1
Je me souviens d'avoir eu un problème similaire il y a environ un an lorsque j'ai joué avec le fdapackage. J'écrivais une réponse qui impliquait d'obtenir des prédictions manuellement, mais une grande partie de celle-ci a été perdue en raison de ne pas l'avoir enregistrée. Si quelqu'un d'autre ne me bat pas, je devrais avoir une solution pour vous dans quelques jours.

Réponses:

14

Je ne me soucie pas de fdal'utilisation des structures d'objets list-within-list-within-list de type Inception , mais ma réponse respectera le système que les rédacteurs de packages ont créé.

Je pense qu'il est instructif de penser d'abord à ce que nous faisons exactement. Sur la base de votre description de ce que vous avez fait jusqu'à présent, c'est ce que je crois que vous faites (faites-moi savoir si j'ai mal interprété quelque chose). Je continuerai à utiliser la notation et, en raison du manque de données réelles, un exemple de l' analyse fonctionnelle des données de Ramsay et Silverman et l' analyse fonctionnelle des données de Ramsay, Hooker et Graves avec R et MATLAB (certaines des équations et du code suivants sont directement extraits de ces livres).

Nous modélisons une réponse scalaire via un modèle linéaire fonctionnel, c'est-à-dire

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

Nous développons la dans une certaine base. Nous utilisons, disons, fonctions de base. Alors,KβK

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

En notation matricielle, il s'agit de .β(s)=θ(s)b

Nous développons également les fonctions de covariation sur une certaine base également (disons les fonctions de base ). Alors,L

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

Encore une fois, en notation matricielle, il s'agit de .X(s)=Cψ(s)

Et donc, si nous laissons , notre modèle peut être exprimé commeJ=ψ(s)θ(s)ds

y=β0+CJb .

Et si nous laissons et , notre modèle estZ=[1CJ]ξ=[β0b]

y=Zξ

Et cela nous semble beaucoup plus familier.

Maintenant, je vois que vous ajoutez une sorte de régularisation. Le fdapackage fonctionne avec des pénalités de rugosité du formulaire

P=λ[Lβ(s)]2ds

pour un certain opérateur différentiel linéaire . Maintenant, il peut être montré (détails laissés ici - ce n'est vraiment pas difficile à montrer) que si nous définissons la matrice de pénalité commeLR

R=λ(0000R1000RK)

où est en termes d'expansion de base de , alors nous minimisons la somme pénalisée des carrés:Riβi

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

et donc notre problème est simplement une régression de crête avec solution:

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

J'ai parcouru ce qui précède parce que, (1) je pense qu'il est important que nous comprenions ce que nous faisons, et (2) une partie de ce qui précède est nécessaire pour comprendre une partie du code que j'utiliserai plus tard. Passons au code ...

Voici un exemple de données avec le code R. J'utilise l'ensemble de données météorologiques canadiennes fourni dans le fdapackage. Nous modéliserons les précipitations annuelles logarithmiques pour un certain nombre de stations météorologiques via un modèle linéaire fonctionnel et nous utiliserons les profils de température (les températures ont été enregistrées une fois par jour pendant 365 jours) de chaque station comme covariables fonctionnelles. Nous procéderons de la même manière que vous décrivez dans votre situation. Les données ont été enregistrées dans 35 stations. Je vais diviser l'ensemble de données en 34 stations, qui seront utilisées comme mes données, et la dernière station, qui sera mon "nouveau" ensemble de données.

Je continue via le code R et les commentaires (je suppose que vous êtes assez familier avec le fdapackage pour que rien de ce qui suit n'est trop surprenant - si ce n'est pas le cas, faites-le moi savoir):

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

Maintenant, quand j'ai appris pour la première fois les données fonctionnelles il y a environ un an, j'ai joué avec ce package. Je n'ai pas non plus pu me predict.fRegressdonner ce que je voulais. En y repensant maintenant, je ne sais toujours pas comment le faire se comporter. Donc, nous devrons simplement obtenir les prédictions semi-manuellement. Je vais utiliser des morceaux que j'ai extraits du code fRegress(). Encore une fois, je continue via le code et les commentaires.

Tout d'abord, la configuration:

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

Maintenant pour obtenir les prédictions

y^new=Znewξ^

Je prends juste le code qui l' fRegressutilise pour le calculer yhatfdobjet le modifier légèrement. fRegresscalcule yhatfdobjen estimant l'intégrale via la règle trapézoïdale (avec et étendus dans leurs bases respectives). 0TXi(s)β(s)Xiβ

Normalement, fRegresscalcule les valeurs ajustées en parcourant les covariables stockées dans annPrecTemp$xfdlist. Donc, pour notre problème, nous remplaçons cette liste de covariables par celle correspondante dans notre nouvelle liste de covariables, c'est-à-dire templistNew. Voici le code (identique au code trouvé fRegressavec deux modifications, quelques suppressions de code inutile et quelques commentaires ajoutés):

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

(Remarque: si vous regardez ce morceau et le code qui l'entoure fRegress, vous verrez les étapes que j'ai décrites ci-dessus).

J'ai testé le code en réexécutant l'exemple météo en utilisant les 35 stations comme données et j'ai comparé la sortie de la boucle ci-dessus à annPrecTemp$yhatfdobjet tout correspond. Je l'ai également exécuté plusieurs fois en utilisant différentes stations comme mes "nouvelles" données et tout semble raisonnable.

Faites-moi savoir si l'un des éléments ci-dessus n'est pas clair ou si quelque chose ne fonctionne pas correctement. Désolé pour la réponse trop détaillée. Je ne pouvais pas m'en empêcher :) Et si vous ne les possédez pas déjà, consultez les deux livres que j'ai utilisés pour rédiger cette réponse. Ce sont de très bons livres.


la source
On dirait que c'est exactement ce dont j'ai besoin. Je vous remercie. Je suppose que je n'aurai pas à jouer avec les trucs nfine / tine / deltat non? Dois-je supposer que l'intégration se déroule avec suffisamment de précision?
dcl
De plus, je remarque que vous n'avez pas pénalisé directement la «nouvelle» covariable ou les «anciennes» covariables. Tout est fait avec la bêta pénalisante (et le nombre de fonctions de base je suppose). La pénalité lambda est appliquée à la bêta. Réalisez-vous le même effet en pénalisant les lissages avant régression? (avec la même valeur de lambda)
dcl
1
La grille utilisée pour approximer l'intégrale est assez fine, donc l'approximation devrait être assez bonne. Vous pouvez toujours augmenter nfineet voir à quel point les modifications intégrales changent, mais je suppose que cela ne fera pas grand-chose. En ce qui concerne la pénalisation, oui, nous pénalisons directement les au lieu des dans ce cas. Ramsay et Silverman discutent d'une autre méthode de pénalité qui estime sans fonctions de base où nous appliquons la pénalité directement à . Les deux façons induisent une contrainte de fluidité sur les fonctions , mais je ne sais pas si vous obtiendrez le «même effet». ξββ^ββ
J'ai essayé de manipuler le code pour produire des prédictions pour plusieurs courbes, mais je ne suis pas sûr de l'avoir fait correctement. Pour commencer, yhatmat n'est pas constant pour tous les courbes après la première itération de la boucle ... Est-ce que cela est censé être équivalent à ? β0
dcl
1
@dcl Dans la boucle, lorsque , il ajoute à (en supposant que la première liste de votre Xlist correspond au terme d'interception). Pouvez-vous ajouter l'extrait de code que vous utilisez à votre question afin que je puisse l'examiner? j=1β0^y^