The smoothing spline models are considered as penalized model, and hence, hypothesis testing of the regression coefficients becomes more complicated due to the shrinkage imposed on the coefficients. Similarly, when jointly testing more than one coefficients, for example testing the penalized part of a spline function, the degree of freedom is unclear. Defaulting to the number of coefficients will overestimate the degree of freedom, and hence introduce under-powered test. Using the well established concept “effective degree of freedom” can alleviate the problem, but introduce another one. When using Wald type of test, the test statistics follows a \(\chi^2\) distribution with degree of freedom \(k\), where \(k\) is a positive integer. However, this does not applied to the effective degree of freedom which is not necessarily an integer. Hence, establishing the limiting distribution of the test statistics becomes a new challenge.
Overall, using hypothesis testing for function selection has two-fold challenges: first of all, establish the effective degree of freedom; secondly, find the sampling distribution of the test statistics.
Estimating effective degree of freedom (EDF) is a common problem encountered in complex models, such as penalized models and Bayesian hierarchical models. It is used to establish the complexity of a model. The problem of estimating EDF is well studied. Previous work on estimating EDF include. In our GAM model, we follow Wood (2017, PP. 252) to estimate EDF.
The estimated EDF of a GAM model is \[ \tau = \text{tr}(2A-AA), \] where A is \(X(X^TX+S_\lambda)^{-1}X^T\) and \(S_\lambda = \sum_j \lambda_j S_j\).
As our Bayesian model is not a penalized model, the penalty parameter \(\lambda_j\) is not directly available, and hence needed to be estimate. Following our mixture normal model specification, we can write the SS prior as a \(l_2\) penalty following Equation X
\[ log(p(\beta_j|\gamma_j)) = \]
n <- 100
p_1 <- 5
p_2 <- 10
x_1 <- MASS::mvrnorm(n,rep(0, p_1), diag(p_1))
x_2 <- MASS::mvrnorm(n,rep(0, p_2), diag(p_2))
x <- cbind(x_1, x_2)
A <- x%*%solve(crossprod(x) + diag(p_1+p_2))%*%t(x)
A_1 <- x_1%*%solve(crossprod(x_1) + diag(p_1))%*%t(x_1)
A_2 <- x_2%*%solve(crossprod(x_2) + diag(p_2))%*%t(x_2)
tau <- sum(diag(2*A))-sum(diag(A%*%A))
tau_1 <- sum(diag(2*A_1))-sum(diag(A_1%*%A_1))
tau_2 <- sum(diag(2*A_2))-sum(diag(A_2%*%A_2))