Vorhersage der Reaktion neuer Kurven mit dem FDA-Paket in R.

10

Grundsätzlich möchte ich nur anhand einiger Kurven eine skalare Reaktion vorhersagen. Ich bin so weit gekommen, eine Regression durchzuführen (mit fRegress aus dem FDA-Paket), habe aber keine Ahnung, wie ich die Ergebnisse auf einen NEUEN Satz von Kurven anwenden soll (zur Vorhersage).

Ich habe N = 536 Kurven und 536 skalare Antworten. Folgendes habe ich bisher getan:

  • Ich habe eine Basis für die Kurven erstellt.
  • Ich habe ein fdPar-Objekt erstellt, um eine Strafe einzuführen
  • Ich habe das fd-Objekt mit glatt.basis erstellt, um die Kurven mit der gewählten Strafe auf der angegebenen Basis zu glätten.
  • Ich habe mit fRegress () eine Regression ausgeführt und dabei die Kurven der Skalarantwort zurückgeführt.

Jetzt möchte ich nur noch diese Regression verwenden, um Vorhersagen für einen neuen Datensatz zu erstellen, den ich habe. Ich kann anscheinend keinen einfachen Weg finden, dies zu tun.

Prost

dcl
quelle
Selbst eine Beschreibung der manuellen Berechnung der Vorhersagen anhand der Basis, der geglätteten (fd) Objekte und der Regressionsschätzungen von fRegress () wäre sehr hilfreich.
dcl
Nur überprüfen: Haben Sie versucht, predict.fRegressdie newdataOption zu verwenden (aus dem FDA-Handbuch hier )?
Ich bin mir nur nicht sicher, welche Klasse oder welches Format die 'neuen Daten' haben müssen. Es wird kein fd- oder fdSmooth-Objekt akzeptiert, das die geglätteten Kurven sind, aus denen ich vorhersagen möchte. Und ich kann keine rohen Argumente und kovariaten Werte eingeben.
dcl
1
Ich erinnere mich, dass ich vor ungefähr einem Jahr ein ähnliches Problem hatte, als ich mit dem fdaPaket herumgespielt habe. Ich habe eine Antwort geschrieben, bei der Vorhersagen manuell abgerufen wurden, aber ein großer Teil davon ging verloren, weil sie nicht gespeichert wurden. Wenn mich jemand anderes nicht schlägt, sollte ich in ein paar Tagen eine Lösung für Sie haben.

Antworten:

14

Die fdaVerwendung von Inception- ähnlichen Objektstrukturen für Listen innerhalb von Listen innerhalb von Listen ist mir egal , aber meine Antwort wird sich an das System halten, das die Paketschreiber erstellt haben.

Ich finde es lehrreich, zuerst darüber nachzudenken, was wir genau tun. Basierend auf Ihrer Beschreibung dessen, was Sie bisher getan haben, glaube ich, dass Sie dies tun (lassen Sie mich wissen, wenn ich etwas falsch interpretiert habe). Ich werde weiterhin die Notation verwenden und aufgrund fehlender realer Daten ein Beispiel aus der Funktionsdatenanalyse von Ramsay und Silverman sowie die Funktionsdatenanalyse von Ramsay, Hooker und Graves mit R und MATLAB (Einige der folgenden Gleichungen und Codes werden direkt aufgehoben aus diesen Büchern).

Wir modellieren eine skalare Antwort über ein funktionales lineares Modell, dh

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

Wir erweitern die in gewisser Weise. Wir verwenden beispielsweise Basisfunktionen. Damit,βK

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

In der Matrixnotation ist dies .β(s)=θ(s)b

Wir erweitern auch die kovariaten Funktionen in einigen Basen (sagen wir Basisfunktionen). Damit,L

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

In der Matrixnotation ist dies wiederum .X(s)=Cψ(s)

Wenn wir also , kann unser Modell ausgedrückt werden alsJ=ψ(s)θ(s)ds

y=β0+CJb .

Und wenn wir und , ist unser ModellZ=[1CJ]ξ=[β0b]

y=Zξ

Und das kommt uns viel vertrauter vor.

Jetzt sehe ich, dass Sie eine Art Regularisierung hinzufügen. Das fdaPaket arbeitet mit Rauheitsstrafen des Formulars

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

für einige linearen Differentialoperator . Jetzt kann gezeigt werden (Details hier weggelassen - es ist wirklich nicht schwer, dies zu zeigen), dass, wenn wir die Strafmatrix als definierenLR

R=λ(0000R1000RK)

wo in Bezug auf die Basis Ausbau der ist , dann minimieren wir die bestraft Summe der Quadrate:Riβi

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

und so ist unser Problem nur eine Gratregression mit Lösung:

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

Ich habe das oben Gesagte durchgearbeitet, weil (1) ich denke, es ist wichtig, dass wir verstehen, was wir tun, und (2) einige der oben genannten Dinge notwendig sind, um einen Teil des Codes zu verstehen, den ich später verwenden werde. Weiter zum Code ...

Hier ist ein Datenbeispiel mit R-Code. Ich verwende den im fdaPaket enthaltenen kanadischen Wetterdatensatz . Wir modellieren den logarithmischen Jahresniederschlag für eine Reihe von Wetterstationen über ein funktionelles lineares Modell und verwenden Temperaturprofile (die Temperaturen wurden 365 Tage lang einmal täglich aufgezeichnet) von jeder Station als funktionale Kovariaten. Wir werden ähnlich vorgehen, wie Sie es in Ihrer Situation beschreiben. Die Daten wurden an 35 Stationen aufgezeichnet. Ich werde den Datensatz in 34 Stationen aufteilen, die als meine Daten verwendet werden, und die letzte Station, die mein "neuer" Datensatz sein wird.

Ich fahre mit R-Code und Kommentaren fort (ich gehe davon aus, dass Sie mit dem fdaPaket so vertraut sind, dass nichts im Folgenden zu überraschend ist - wenn dies nicht der Fall ist, lassen Sie es mich bitte wissen):

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

Als ich vor ungefähr einem Jahr zum ersten Mal über funktionale Daten unterrichtet wurde, habe ich mit diesem Paket herumgespielt. Ich konnte predict.fRegressmir auch nicht geben, was ich wollte. Wenn ich jetzt zurückblicke, weiß ich immer noch nicht, wie ich es dazu bringen soll, sich zu verhalten. Wir müssen die Vorhersagen also nur halb manuell abrufen. Ich werde Teile verwenden, für die ich direkt aus dem Code gezogen habe fRegress(). Wieder gehe ich über Code und Kommentare weiter.

Zunächst das Setup:

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

Nun zu den Vorhersagen

y^new=Znewξ^

Ich nehme nur den Code, fRegressmit dem er berechnet yhatfdobjund leicht bearbeitet wird. fRegressberechnet yhatfdobjdurch Schätzen des Integrals über die Trapezregel (wobei und in ihren jeweiligen Basen erweitert sind). 0TXi(s)β(s)Xiβ

Normalerweise fRegresswerden die angepassten Werte berechnet, indem die in gespeicherten Kovariaten durchlaufen werden annPrecTemp$xfdlist. Für unser Problem ersetzen wir diese Kovariatenliste durch die entsprechende in unserer neuen Kovariatenliste, d templistNew. H. Hier ist der Code (identisch mit dem Code, der fRegressmit zwei Änderungen, einigen Löschungen von nicht benötigtem Code und einigen Kommentaren hinzugefügt wurde):

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

(Hinweis: Wenn Sie sich diesen Block und den umgebenden Code in fRegressansehen, werden die oben beschriebenen Schritte angezeigt.)

Ich habe den Code getestet, indem ich das Wetterbeispiel mit allen 35 Stationen als Daten erneut ausgeführt und die Ausgabe der obigen Schleife mit verglichen habe annPrecTemp$yhatfdobjund alles übereinstimmt. Ich habe es auch ein paar Mal mit verschiedenen Stationen als "neue" Daten ausgeführt und alles scheint vernünftig.

Lassen Sie mich wissen, wenn einer der oben genannten Punkte unklar ist oder wenn etwas nicht richtig funktioniert. Entschuldigen Sie die zu detaillierte Antwort. Ich konnte mir nicht helfen :) Und wenn Sie sie noch nicht besitzen, lesen Sie die beiden Bücher, mit denen ich diese Antwort geschrieben habe. Das sind wirklich gute Bücher.


quelle
Das sieht so aus, als wäre es genau das, was ich brauche. Vielen Dank. Ich nehme an, ich muss nicht mit dem nfine / tine / deltat-Zeug herumspielen, oder? Sollte ich davon ausgehen, dass die Integration genau genug erfolgt?
dcl
Außerdem stelle ich fest, dass Sie die "neue" Kovariate oder die "alten" Kovariaten nicht direkt bestraft haben. Es ist alles mit der Bestrafung der Beta (und der Anzahl der Basisfunktionen, denke ich) erledigt. Das Straf-Lambda wird auf die Beta angewendet. Erreichen Sie den gleichen Effekt, indem Sie die Glättungen vor der Regression bestrafen? (mit dem gleichen Wert von Lambda)
dcl
1
Das Gitter, das zur Approximation des Integrals verwendet wird, ist ziemlich gut, daher sollte die Approximation ziemlich gut sein. Sie könnten immer zunehmen nfineund sehen, wie sehr sich das Integral ändert, aber ich vermute, es wird nicht viel bewirken. Was die Bestrafung angeht, ja, wir bestrafen in diesem Fall direkt die 's anstelle von ' s. Ramsay und Silverman diskutieren eine andere Strafmethode, die ohne Basisfunktionen schätzt, wobei wir die Strafe direkt auf anwenden . Beide Methoden führen zu einer Einschränkung der Glätte der Funktionen, aber ich bin mir nicht sicher, ob Sie den gleichen Effekt erzielen. ξββ^ββ
Ich habe versucht, den Code zu manipulieren, um Vorhersagen für mehrere Kurven zu erstellen, bin mir aber nicht sicher, ob ich es richtig gemacht habe. Für den Anfang ist Yhatmat nach der ersten Iteration der Schleife nicht für alle Kurven konstant ... Soll dies ? β0
dcl
1
@dcl Wenn in der Schleife , wird zu hinzugefügt (vorausgesetzt, die erste Liste in Ihrer X-Liste entspricht dem Intercept-Term). Können Sie Ihrer Frage den Codeausschnitt hinzufügen, den Sie verwenden, damit ich ihn mir ansehen kann? j=1β0^y^