机器算法验证 - 随着时间的推移对样条进行建模——设计矩阵和方法调查 - 吾爱随笔录

响应变量 y 是多个预测变量 X 的非线性函数（在我的真实数据中，响应是二项式分布的，但为了简单起见，这里我使用的是正态分布值）。mgcv我可以使用样条/平滑（例如， R 包中的 GAM 模型）对预测变量和响应之间的关系进行建模。

到目前为止，一切都很好。但是，每个响应都是随时间演变的过程的结果。也就是说，预测变量 X 和响应 y 之间的关系随时间而变化。对于每个响应，我都有响应周围多个时间点的预测变量数据。也就是说，每组时间点有一个响应（不是响应随时间演变）。

在这一点上，一些插图可能会有所帮助。这是一些具有已知参数的数据（下面的代码），然后使用 ggplot2 绘制（指定 GAM 方法和适当的平滑器），并在刻面中显示时间。为了说明，y 是 x1 的二次函数，并且这种关系的符号和大小随时间而变化。

编辑：因为我是新来的，所以我不能发布图片；运行代码来查看它。

x2和y之间的关系是圆形的，对应于y随着x2的某个方向增加。这种关系的幅度随时间而变化。（使用指定“cc”圆形三次平滑器的 gam 在 ggplot 中建模）。

编辑：因为我是新来的，所以我不能发布图片；运行代码来查看它。

我想使用类似二维样条的东西将每个预测变量的（非线性）变化建模为时间的函数。

我已经考虑在 mgcv 包中使用二维平滑（类似于te(x1,t)），除了这需要长格式的数据（即单列时间点）。我认为这是不合适的，因为一个响应与所有时间点相关联——因此以长格式排列数据（因此在设计矩阵的多行中复制相同的响应）将违反观察的独立性。我的数据目前按列排列，(y, x1.t1, x1.t2, x1.t3, ..., x2.t1, x2.t2, ...)我认为这是最合适的格式。

我想知道：

有没有更好的方法来建模这些数据
如果是这样，模型的设计矩阵/公式会是什么样子。最终，我想在像 JAGS 这样的 mcmc 包中使用贝叶斯推理来估计模型系数，所以我想知道如何编写二维样条曲线。

重现我的示例的 R 代码：

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