为什么我的设计矩阵秩不足？（用循环样条建模季节性数据）

机器算法验证 r 季节性 样条 毫克CV
2022-04-11 22:18:08

我正在尝试将障碍模型拟合到已知一年中的哪一天的季节性数据。

为此，我尝试构建具有循环限制的样条曲线（末端与起点相接）。这是一个 MWE 使用lm而不是障碍模型来简化问题：

# Example data
set.seed(1234)
n <- 1000
x <- runif(n, 0, 2 * pi) # some random time in the year (2 * pi being day 365)
y <- sin(x) + rnorm(n)

# Generate cyclic B-spline bases using mgcv
k            <- 4
knots        <- seq(0, 2 * pi, length.out = k)
cyclicSpline <- mgcv::cSplineDes(x, knots = knots)

但是，如果我尝试使用生成的样条基函数作为模型的预测变量，它最终会出现秩不足：

> lm(y ~ cyclicSpline)

Call:
lm(formula = y ~ cyclicSpline)

Coefficients:
  (Intercept)  cyclicSpline1  cyclicSpline2  cyclicSpline3  cyclicSpline4  
     -0.09211       -1.40510        0.03613        1.68947             NA

我尝试正交化样条基函数，但这仍然导致秩不足：

> Q <- svd(t(cyclicSpline))$v
> lm(y ~ Q) 

Call:
lm(formula = y ~ Q)

Coefficients:
(Intercept)           Q1           Q2           Q3           Q4  
      8.173      258.817       -7.297       21.510           NA

我也尝试过使用不同数量的结，但它似乎总是导致等级不足。我的印象是，一旦正交化，我最终会得到独立的基函数。如果是这样，为什么会导致秩不足的设计矩阵？

3个回答

cyclicSpline已经在其跨度中包含常量向量，因此如果您另外添加截距，它将是秩不足的。

Matrix::rankMatrix(cyclicSpline)表明它cyclicSpline本身就是满级，做lm(y ~ cyclicSpline - 1)会解决这个问题。

为了确认确实如此，我们可以显式地计算帽子矩阵U = svd(cyclicSpline)$u; H = U %*% t(U)，然后检查它H %*% rep(1, n)是否在 的数值舍入内rep(1,n)（即H作为所有 1 的向量的恒等式，因此向量完全在 的列空间内cyclicSpline）。

正如jld 所写，您的样条曲线在其跨度中包含常数 1，因此您应该拟合没有截距的模型。

事实上，样条曲线的总和为 1：rowSums(cyclicSpline)给你一个 1s 的常数向量。

这是一个x针对其样条变换的图，它很好地显示了样条如何添加到 1：

样条

plot(x,cyclicSpline[,1],pch=19)
points(x,cyclicSpline[,2],pch=19,col=2)
points(x,cyclicSpline[,3],pch=19,col=3)

我确信这种缺乏识别是由于您的模型还包括一个常数（截距）。mgcv肯定会在它创建的循环样条上应用可识别性约束，cSplineDes()但这会发生在其他函数中，因为您想通过smooth.construct.xxxx函数应用这些约束，而cSplineDes这是一个仅创建基础的较低级别的函数。

删除拦截，事情将开始工作：

> # Example data
> set.seed(1234)
> n <- 1000
> x <- runif(n, 0, 2 * pi) # some random time in the year (2 * pi being day 365)
> y <- sin(x) + rnorm(n)
> 
> # Generate cyclic B-spline bases using mgcv
> k            <- 4
> knots        <- seq(0, 2 * pi, length.out = k)
> cyclicSpline <- mgcv::cSplineDes(x, knots = knots)
> m1 <- lm(y ~ cyclicSpline)
> m2 <- lm(y ~ cyclicSpline - 1) ## remove intercept
> AIC(m1, m2)
   df      AIC
m1  4 2742.641
m2  4 2742.641
其它你可能感兴趣的问题
上一篇制服的绝对值 下一篇为什么每次我重新进行分析时，引导标准误差和 95% 置信区间都会发生变化