Modélisation d'une spline dans le temps - matrice de conception et étude des approches

8

Une variable de réponse y est une fonction non linéaire d'un certain nombre de variables prédictives X (dans mes données réelles, la réponse est distribuée de façon binomiale, mais ici j'utilise une valeur normalement distribuée pour plus de simplicité). Je peux modéliser les relations entre les prédicteurs et la réponse à l'aide de splines / lissages (par exemple, les modèles GAM dans le mgcvpackage dans R).

Jusqu'ici tout va bien. Cependant, chaque réponse est le résultat de processus qui évoluent avec le temps. C'est-à-dire que la relation entre les prédicteurs X et la réponse y change avec le temps. Pour chaque réponse, j'ai des données pour les prédicteurs sur un certain nombre de points de temps autour de la réponse. Autrement dit, il y a une réponse par groupe de points dans le temps (pas que la réponse évolue avec le temps).

Certaines illustrations sont probablement utiles à ce stade. Voici quelques données avec des paramètres connus (code ci-dessous), puis tracées à l'aide de ggplot2 (spécifiant la méthode GAM et un lisseur approprié), avec du temps dans les facettes. À titre d'illustration, y est une fonction quadratique de x1, et le signe et l'ampleur de cette relation changent en fonction du temps.

EDIT: depuis que je suis nouveau je ne peux pas poster d'image;  exécutez le code pour le voir.

La relation entre x2 et y est circulaire, correspondant à une augmentation de y avec une certaine direction de x2. L'amplitude de cette relation module au fil du temps. (Modélisé dans ggplot en utilisant un gam spécifiant un lisseur cubique circulaire "cc").

EDIT: depuis que je suis nouveau je ne peux pas poster d'image;  exécutez le code pour le voir.

Je voudrais modéliser le changement (non linéaire) de chaque prédicteur en fonction du temps en utilisant quelque chose comme une spline bidimensionnelle.

J'ai envisagé d'utiliser un lissage bidimensionnel dans le package mgcv (quelque chose comme te(x1,t)), sauf que cela nécessiterait les données sous une forme longue (c'est-à-dire une seule colonne de points temporels). Je pense que cela est inapproprié, car une réponse est associée à tous les points dans le temps - donc organiser les données sous une forme longue (dupliquant ainsi la même réponse sur plusieurs lignes de la matrice de conception) violerait l'indépendance des observations. Mes données sont actuellement organisées en colonnes (y, x1.t1, x1.t2, x1.t3, ..., x2.t1, x2.t2, ...)et je pense que c'est le format le plus approprié.

J'aimerais savoir:

  1. existe-t-il une meilleure façon de modéliser ces données
  2. dans l'affirmative, à quoi ressemblerait la matrice / formule de conception du modèle. En fin de compte, je voudrais estimer les coefficients du modèle en utilisant l'inférence bayésienne dans un package mcmc comme JAGS, donc je voudrais savoir comment écrire une spline bidimensionnelle.

Code R pour reproduire mon exemple:

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
la source
1
Les données des séries chronologiques ne sont jamais des données indépendantes, car l'évolution dans le temps corrèle les observations de la même série temporelle (auto-corrélation temporelle), sauf si, bien sûr, vous avez affaire à du bruit. Mais cette corrélation temporelle n'a rien à voir avec le fait que les 2 séries temporelles sont indépendantes les unes des autres. Le temps est le même pour x1 et x2, donc le ième score ou l'entrée de la variable x1 doit être dans la même ligne de la matrice du modèle que le ième score de la variable x2. Je ne vois pas pourquoi cela viole l'hypothèse que les séries chronologiques sont indépendantes. En cas de doute, calculez la corrélation croisée.
nukimov

Réponses:

1

Je suis d'accord avec vous que vous devrez peut-être tenir compte des termes d'erreur des répondants individuels dans le temps, en particulier si vous n'avez pas de résultat pour toute la période pour chaque répondant.

Une façon de le faire est avec le BayesX . Il permet des effets spatiaux avec des splines où vous pouvez avoir du temps dans une dimension et la valeur de covariable dans l'autre. De plus, vous pouvez ajouter un effet aléatoire pour chaque observation. Potentiellement, jetez un œil à ce document .

Cependant, je suis sûr que vous devrez mettre votre modèle en format long. De plus, vous devrez ajouter une idcolonne pour le répondant ou l'effet aléatoire.

Benjamin Christoffersen
la source