Modellierung eines Splines über die Zeit - Entwurfsmatrix und Überblick über Ansätze

8

Eine Antwortvariable y ist eine nichtlineare Funktion einer Anzahl von Prädiktorvariablen X (in meinen realen Daten ist die Antwort binomial verteilt, aber hier verwende ich der Einfachheit halber einen normalverteilten Wert). Ich kann die Beziehungen zwischen den Prädiktoren und der Antwort mithilfe von Splines / Smooths modellieren (z. B. GAM-Modelle im mgcvPaket in R).

So weit, ist es gut. Jede Antwort ist jedoch das Ergebnis von Prozessen, die sich im Laufe der Zeit entwickeln. Das heißt, die Beziehung zwischen den Prädiktoren X und der Antwort y ändert sich mit der Zeit. Für jede Antwort habe ich Daten für die Prädiktoren über eine Reihe von Zeitpunkten um die Antwort herum. Das heißt, es gibt eine Antwort pro Gruppe von Zeitpunkten (nicht, dass sich die Antwort im Laufe der Zeit entwickelt).

Einige Abbildungen sind an dieser Stelle wahrscheinlich hilfreich. Hier sind einige Daten mit bekannten Parametern (Code unten) aufgeführt, die dann mit ggplot2 (unter Angabe der GAM-Methode und eines geeigneten Glätters) mit der Zeit in den Facetten aufgezeichnet wurden. Zur Veranschaulichung ist y eine quadratische Funktion von x1, und das Vorzeichen und die Größe dieser Beziehung ändern sich als Funktion der Zeit.

BEARBEITEN: Da ich neu bin, kann ich kein Bild posten.  Führen Sie den Code aus, um ihn anzuzeigen.

Die Beziehung zwischen x2 und y ist kreisförmig, was einer Zunahme von y mit einer bestimmten Richtung von x2 entspricht. Die Amplitude dieser Beziehung moduliert sich über die Zeit. (Modelliert in ggplot unter Verwendung eines Spiels, das einen "cc" kreisförmigen kubischen Glätter angibt).

BEARBEITEN: Da ich neu bin, kann ich kein Bild posten.  Führen Sie den Code aus, um ihn anzuzeigen.

Ich möchte die (nichtlineare) Änderung in jedem Prädiktor als Funktion der Zeit mit einem zweidimensionalen Spline modellieren.

Ich habe überlegt, eine zweidimensionale Glättung im mgcv-Paket zu verwenden (so etwas wie te(x1,t)), außer dass dies die Daten in einer langen Form (dh einer einzelnen Spalte von Zeitpunkten) erfordern würde. Ich halte dies für unangemessen, da eine Antwort allen Zeitpunkten zugeordnet ist. Wenn Sie also die Daten in langer Form anordnen (und damit dieselbe Antwort über mehrere Zeilen der Entwurfsmatrix duplizieren), würde dies die Unabhängigkeit der Beobachtungen verletzen. Meine Daten sind derzeit mit Spalten angeordnet, (y, x1.t1, x1.t2, x1.t3, ..., x2.t1, x2.t2, ...)und ich denke, dies ist das am besten geeignete Format.

Ich würde gerne wissen:

  1. Gibt es eine bessere Möglichkeit, diese Daten zu modellieren?
  2. Wenn ja, wie würde die Entwurfsmatrix / -formel des Modells aussehen? Letztendlich möchte ich Modellkoeffizienten mithilfe der Bayes'schen Inferenz in einem mcmc-Paket wie JAGS schätzen, also möchte ich wissen, wie man einen zweidimensionalen Spline schreibt .

R-Code zur Reproduktion meines Beispiels:

library(ggplot2)
library(mgcv)
#-------------------
# start by generating some data with known relationships between two variables,
# one periodic, over time.
set.seed(123)

nTimeBins <- 6
nSamples <- 500

# the relationship between x1, x2 and y are not linear.
# y = 0.4*x1^2 -1.2*x1 + 0.4*sin(x2) + 1.2*cos(x2)

# the relationship between x1, x2 and y evolve over time. 
x1.timeMult <- cos(seq(-pi,pi,length=nTimeBins))
x2.timeMult <- cos(seq(-pi/2,pi/2,length=nTimeBins))

qplot(x=1:nTimeBins,y=x1.timeMult,geom="line") + 
      geom_line(aes(x=1:nTimeBins,y=x2.timeMult,colour="red")) + 
      guides(colour=FALSE) + ylab("multiplier")

df <- data.frame(setup=rep(NA,times=nSamples))
for (time in 1 : nTimeBins){
  text <- paste('df$x1.t',time,' <- runif(nSamples,min=-3,max=3)',sep="")
  eval(parse(text=text))

  text <- paste('df$x2.t',time,' <- runif(nSamples,min=-pi,max=pi)',sep="")
  eval(parse(text=text))
}
df$setup <- NULL

# each y is a function of x over time.
text <- 'y <- '
# replicated from above for reference:
# y = 0.4*x1^2 -1.2*x1 + 0.4*sin(x2) + 1.2*cos(x2)
for (time in 1 : nTimeBins){
  text <- paste(text,'(0.4*x1.t',time,'^2-1.2*x1.t',time,') * 
                x1.timeMult[',time,'] + (0.4*sin(x2.t',time,') + 
                1.2*cos(x2.t',time,'))*x2.timeMult[',time,'] + ',sep="")
}

text <- paste(text,'rnorm(nSamples,sd=0.2)')
attach(df)
eval(parse(text=text))
df$y <- y

#-------------------
# transform into long form data for plotting:
df.long <- data.frame(y=rep(df$y,times=nTimeBins))
textX1 <- 'df.long$x1 <- c('
textX2 <- 'df.long$x2 <- c('

for (time in 1:nTimeBins){
    textX1 <- paste(textX1,'x1.t',time,',',sep="")
    textX2 <- paste(textX2,'x2.t',time,',',sep="")
}
textX1 <- paste(textX1,'NULL)',sep="")
textX2 <- paste(textX2,'NULL)',sep="")
eval(parse(text=textX1))
eval(parse(text=textX2))
# time stamp:
df.long$t <- factor(rep(1:nTimeBins,each=nSamples))

#-------------------
# plot relationships over time using GAM fits in ggplot:
p1 <- ggplot(df.long,aes(x=x1,y=y)) + geom_point() + 
             stat_smooth(method="gam",formula=y ~ s(x,bs="cs",k=4)) + 
             facet_wrap(~ t, ncol=3) + opts(title="x1 versus y over time")
p1

p2 <- ggplot(df.long,aes(x=x2,y=y)) + geom_point() + 
             stat_smooth(method="gam",formula=y ~ s(x,bs="cc",k=5)) + 
             facet_wrap(~ t, ncol=3) + opts(title="x2 versus y over time")
p2
Tsawallis
quelle
1
Daten aus Zeitreihen sind niemals unabhängige Daten, da die zeitliche Entwicklung die Beobachtungen derselben Zeitreihe korreliert (Zeit-Autokorrelation), es sei denn, Sie haben es natürlich mit Rauschen zu tun. Diese Zeitkorrelation hat jedoch nichts damit zu tun, dass die beiden Zeitreihen unabhängig voneinander sind. Die Zeit ist für x1 und x2 gleich, daher sollte sich die i-te Punktzahl oder der Eintrag der Variablen x1 in derselben Zeile der Modellmatrix befinden wie die i-te Punktzahl der Variablen x2. Ich verstehe nicht, warum dies gegen die Annahme verstößt, dass die Zeitreihen unabhängig sind. Berechnen Sie im Zweifelsfall die Kreuzkorrelation.
Nukimov

Antworten:

1

Ich stimme Ihnen zu, dass Sie möglicherweise die Fehlerbedingungen einzelner Befragter im Laufe der Zeit berücksichtigen müssen, insbesondere wenn Sie nicht für jeden Befragten das Ergebnis für den gesamten Zeitraum hätten.

Ein Weg, dies zu tun, ist mit dem BayesX . Es ermöglicht räumliche Effekte mit Splines, bei denen Sie Zeit in einer Dimension und den kovariaten Wert in der anderen haben können. Außerdem können Sie für jede Beobachtung einen zufälligen Effekt hinzufügen. Schauen Sie sich möglicherweise dieses Papier an .

Ich bin mir jedoch ziemlich sicher, dass Sie Ihr Modell in ein Langformat bringen müssen. Außerdem müssen Sie eine idSpalte für den Befragten oder den Zufallseffekt hinzufügen .

Benjamin Christoffersen
quelle