Why is my design matrix rank deficient? (modelling seasonal data with a cyclic spline)

Question

I am trying to fit a hurdle model to seasonal data for which the day of the year is known.

To this end, I tried to construct a spline with cyclic restrictions (ends meet beginnings). Here is a MWE using lm instead of a hurdle model to simplify the problem:

# 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)

However, if I try to use the resulting spline basis functions as predictors for my model, it ends up being rank deficient:

> lm(y ~ cyclicSpline)

Call:
lm(formula = y ~ cyclicSpline)

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

I tried orthogonalizing the spline basis functions, but this still results in rank deficiency:

> 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 

I have also tried this with different numbers of knots, but it seems to always result in rank deficiency. I was under the impression that once orthogonalized, I end up with independent basis functions. If so, why does that result in a rank deficient design matrix?

Answer 1

cyclicSpline already contains the constant vector in its span so if you additionally add an intercept it'll be rank deficient.

Matrix::rankMatrix(cyclicSpline) shows that cyclicSpline is full rank by itself and doing lm(y ~ cyclicSpline - 1) will fix the issue.

To confirm that this really is the case, we can compute the hat matrix explicitly as U = svd(cyclicSpline)$u; H = U %*% t(U) and then check that H %*% rep(1, n) is within numerical rounding of rep(1,n) (i.e. H acts as the identity on the vector of all 1s so that vector is entirely within the column space of cyclicSpline).

Answer 2

As jld writes, your splines contain the constant 1 in their span, so you should fit a model without an intercept.

As a matter of fact, the splines sum rowwise to 1: rowSums(cyclicSpline) gives you a constant vector of 1s.

Here is a plot of x against its spline transform, which shows very nicely how the splines add to 1:

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

Answer 3

I'm sure this lack of identification is due to your model also including a constant (intercept). mgcv is certainly applying identifiability constraints on cyclic splines it creates using cSplineDes() but that happens in other functions as you want to apply those constraints through the smooth.construct.xxxx functions while cSplineDes is a lower-level function that just creates the basis.

Remove the intercept and things will start working:

> # 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