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é
predict.fRegress
l'newdata
option (du manuel de la FDA ici )?fda
package. 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:
Je ne me soucie pas de
fda
l'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
Nous développons la dans une certaine base. Nous utilisons, disons, fonctions de base. Alors,Kβ K
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
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
Et si nous laissons et , notre modèle estZ=[1CJ] ξ=[β0b′]′
Et cela nous semble beaucoup plus familier.
Maintenant, je vois que vous ajoutez une sorte de régularisation. Le
fda
package fonctionne avec des pénalités de rugosité du formulairepour 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é commeL R
où est en termes d'expansion de base de , alors nous minimisons la somme pénalisée des carrés:Ri βi
et donc notre problème est simplement une régression de crête avec solution:
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
fda
package. 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
fda
package pour que rien de ce qui suit n'est trop surprenant - si ce n'est pas le cas, faites-le moi savoir):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.fRegress
donner 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 codefRegress()
. Encore une fois, je continue via le code et les commentaires.Tout d'abord, la configuration:
Maintenant pour obtenir les prédictions
Je prends juste le code qui l'∫T0Xi(s)β(s) Xi β
fRegress
utilise pour le calculeryhatfdobj
et le modifier légèrement.fRegress
calculeyhatfdobj
en estimant l'intégrale via la règle trapézoïdale (avec et étendus dans leurs bases respectives).Normalement,
fRegress
calcule les valeurs ajustées en parcourant les covariables stockées dansannPrecTemp$xfdlist
. Donc, pour notre problème, nous remplaçons cette liste de covariables par celle correspondante dans notre nouvelle liste de covariables, c'est-à-diretemplistNew
. Voici le code (identique au code trouvéfRegress
avec deux modifications, quelques suppressions de code inutile et quelques commentaires ajoutés):(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$yhatfdobj
et 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
nfine
et 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».