7
$\begingroup$

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?

$\endgroup$
0
8
$\begingroup$

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

$\endgroup$
2
  • $\begingroup$ I see, I noticed that it was not rank deficient without an intercept, but I was unsure that was a correct approach. But essentially, including all basis functions is equivalent to including an intercept? $\endgroup$ – Frans Rodenburg Feb 23 at 7:18
  • 1
    $\begingroup$ @FransRodenburg yeah, these basis functions have the constant function in their span so it's just like how if we have $$X = \begin{bmatrix}1&0\\1&0\\0&1\\0&1\end{bmatrix}$$ we don't need to include an intercept because it'd be collinear and anything that the intercept could do we can also do with a combination of the other columns $\endgroup$ – jld Feb 23 at 17:28
6
$\begingroup$

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:

splines

plot(x,cyclicSpline[,1],pch=19)
points(x,cyclicSpline[,2],pch=19,col=2)
points(x,cyclicSpline[,3],pch=19,col=3)
$\endgroup$
1
  • $\begingroup$ The plot is indeed helpful, as is rowSums, thanks! $\endgroup$ – Frans Rodenburg Feb 23 at 7:15
5
$\begingroup$

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
$\endgroup$
2
  • $\begingroup$ Thanks Gavin, your answers have been enormously helpful in better understanding GAMs and splines. Do you think excluding the intercept is the identifiability constraint imposed by mgcv, or something more complex? $\endgroup$ – Frans Rodenburg Feb 23 at 7:20
  • 1
    $\begingroup$ No, the intercept is left in the model, you really want an intercept when dealing with parametric factor terms for example. Simon applies a sum-to-zero constraint over the range of the covariate(s) to the basis for each spline (which is why the plots of the smooths are always centred about zero) for most of the spline types (some special ones get handled differently IIRC). By imposing that sum-to-zero constraint you remove the unidentifiability due to the intercept and the constant term in the basis. It's also why for example, the max EDF of the smooth is always less than the k you set. $\endgroup$ – Gavin Simpson Feb 26 at 16:53

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.