样条曲线用于回归建模以对可能复杂的非线性函数形式进行建模。样条平滑趋势由分段连续多项式组成，其前导系数在每个断点或节点处发生变化。可以根据趋势的多项式次数以及断点来指定样条。协变量的样条表示将观察值的单个向量扩展为一个矩阵，该矩阵的维数是多项式次数加上节点数。

样条曲线的周期性版本仅仅是任何回归的周期性版本：数据被切割成周期长度的重复。因此，例如，在大鼠的多日实验中模拟昼夜趋势需要将实验时间重新编码为 24 小时增量，因此第 154 小时将是 10 的模 24 值（154 = 6*24 + 10）。如果您对切割数据进行线性回归，它将估计趋势的锯齿波形。如果您在周期的某处拟合阶跃函数，它将是适合该系列的方波。样条能够表达更复杂的小波。对于它的价值，在splines包中，有一个功能periodicSpline正是这样做的。

我没有发现 R 的默认样条“bs”实现对解释有用。所以我在下面写了自己的脚本。对于具有次样条曲线，此表示为前列提供标准多项式表示，第列 ( ) 简单地计算为其中是节点的实际向量。$p$$n_k$$p$$p+i$$i \le n_k$$S_{p+i} = (X - k_i)^p\mathcal{I}(X<k_i)$$k$

myspline <- function(x, degree, knots) {
  knots <- sort(knots)
  val <- cbind(x, outer(x, knots, `-`))
  val[val < 0] <- 0
  val <- val^degree
  if(degree > 1)
    val <- cbind(outer(x, 1:{degree-1}, `^`), val)
  colnames(val) <- c(
    paste0('spline', 1:{degree-1}, '.1'),
    paste0('spline', degree, '.', seq(length(knots)+1))
  )
  val
}

（或 ）的域上插入正弦趋势，如下所示：$2\pi$$\tau$

x <- seq(0, 2*pi, by=pi/2^8)
y <- sin(x)
plot(x,y, type='l')
s <- myspline(x, 2, pi)
fit <- lm(y ~ s)
yhat <- predict(fit)
lines(x,yhat)

你会发现他们很默契。此外，命名约定可以解释。在回归输出中，您会看到：

> summary(fit)

Call:
lm(formula = y ~ s)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.04564 -0.02050  0.00000  0.02050  0.04564 

Coefficients:
             Estimate Std. Error  t value Pr(>|t|)    
(Intercept) -0.033116   0.003978   -8.326 7.78e-16 ***
sspline1.1   1.268812   0.004456  284.721  < 2e-16 ***
sspline2.1  -0.400520   0.001031 -388.463  < 2e-16 ***
sspline2.2   0.801040   0.001931  414.878  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.02422 on 509 degrees of freedom
Multiple R-squared:  0.9988,    Adjusted R-squared:  0.9988 
F-statistic: 1.453e+05 on 3 and 509 DF,  p-value: < 2.2e-16

我的 spline1.1 度的第一组协变量是第一个断点后面的第一个域的多项式趋势。线性项是原点处切线的斜率，X=0。这几乎是 1，这将由正弦曲线的导数 (cos(0) = 1) 表示，但我们必须记住，这些是近似值，并且将二次趋势外推到很容易出现错误错误。二次项表示负的凹形。spline2.2 项表示与第一个二次斜率的差异，导致 0.4 的正前导系数表示向上的凸形。所以我们现在可以解释样条输出，并可以相应地判断推断和估计。$\pi/2$

我将假设您知道手头数据的周期性。如果数据缺少增长或移动平均组件，您可以将较长的时间序列转换为持续时间为 1 个周期的较短序列的复制品。您现在有了重复，可以使用数据分析来估计重复趋势。

假设我生成了以下有点嘈杂、很长的时间序列：

x <- seq(1, 100, by=0.01)
y <- sin(x) + rnorm(length(x), 0, 10)
xp <- x %% (2*pi)
s <- myspline(xp, degree=2, knots=pi)
lm(y ~ s)

结果输出显示了合理的性能。

> summary(fit)

Call:
lm(formula = y ~ s)

Residuals:
    Min      1Q  Median      3Q     Max 
-39.585  -6.736   0.013   6.750  37.389 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) -0.48266    0.38155  -1.265 0.205894    
sspline1.1   1.52798    0.42237   3.618 0.000299 ***
sspline2.1  -0.44380    0.09725  -4.564 5.09e-06 ***
sspline2.2   0.76553    0.18198   4.207 2.61e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 9.949 on 9897 degrees of freedom
Multiple R-squared:  0.006406,  Adjusted R-squared:  0.006105 
F-statistic: 21.27 on 3 and 9897 DF,  p-value: 9.959e-14

我最近正在寻找这个问题的答案，并使用最近的包找到了以下解决方案splines2。有一个函数可以计算周期性的 m 样条（m 样条是归一化的 b 样条）。用法与函数非常相似bs。假设我们有一个 24 小时的噪声固定信号，在 2 天内以固定间隔测量：

library(ggplot2)
library(splines2)

t <- seq(0, 48, length.out = 500)
y <- sin(time/2*pi/6) + rnorm(500, sd = 0.5)

df <- data.frame(t = t, y = y)

ggplot(df, aes(x = t, y = y)) + geom_point() + theme_minimal()

现在我们可以在这些数据上拟合一个周期样条，并为我们的定期间隔创建预测：

# (boundary knots determine the period)
pspline_fit <- lm(y ~ mSpline(x = t, 
                              df = 4, 
                              periodic = TRUE, 
                              Boundary.knots = c(0, 24)), data = df)
df <- cbind(df, as.data.frame(predict(pspline_fit, interval = "prediction")))
pred_plot <- 
  ggplot(df, aes(x = t, y = y)) + 
  geom_ribbon(aes(ymin = lwr, ymax = upr), alpha = 0.4) + 
  geom_line(aes(y = fit), size = 1, colour = "blue") + 
  geom_point() + 
  theme_minimal()
pred_plot

周期性样条曲线的优点在于 24 小时标记处没有不连续性，您可以使用极坐标对其进行可视化：

pred_plot + xlim(0, 24) + coord_polar()