I am following the discussion on chapter 5 of Elements of statistical learning which discusses basis expansion using splines. The data set I used is the Ozone data which can be found at :http://web.stanford.edu/~hastie/ElemStatLearn/
It has a response variable "ozone",and three predictor variables "radiation","temperature" and "wind".
Essentially what I am trying to do is get the basis matrix for each one of the predictor variables and combine them to train a generalized linear additive model to predict the ozone levels.
To be more specific the underlying model is : $Y = \theta_{0} + \boldsymbol{h}_{1}^{T}\boldsymbol{\theta}_{1} + \boldsymbol{h}_{2}^{T}\boldsymbol{\theta}_{2}+ \boldsymbol{h}_{3}^{T}\boldsymbol{\theta}_{3} = \theta_{0} + f_{1}(X_{1}) +f_{2}(X_{2})+f_{3}(X_{3})$ , where Where $\theta_{0}$ is the intercept term. $X_{i}$ is the $i^{th}$ feature of data X.
And $\boldsymbol{h}_{1}(X_{1})$ is of form $\begin{bmatrix}
h_{11}(X_{1})\\
h_{12}(X_{1})\\
h_{13}(X_{1})\\
h_{14}(X_{1})\\
\vdots
\end{bmatrix} , $. Each entry of $\boldsymbol{h}_{1}(X_{1})$ is a spline basis function. That is each variable is expanded by the spline basis (natural cubic spline or smoothing spline etc.) and then do a linear regression on the whole model.
The smoothing spline is found by minimizing $RSS(f,\lambda) = \sum_{i}^{N}\{y_{i}-f(x_{i})\}^{2} + \lambda\int\{f''(t)\}dt$, and the solution is given by $\hat{\theta} = (N^{T}N + \lambda\Omega_{N})^{-1}N^{T}\boldsymbol{y}$, where $\{\Omega_{N}\}_{jk} = \int N''_{j}(t)N''_{k}(t) dt$
I am really stuck on how to calculate the $\Omega$ matrix here. Currently I am using R and with the built-in function "ns(x,knots = )" from which I can get a basis matrix with specified knot locations. Since smoothing splines place knots at each distinctive x value , the command I use is ns(Ozone$radiation,knots = sort(Ozone$radiation)[2:110]). The data matrix is of size (111,3), since ns function does not allow placing knots at boundary points I only select the values from index 2 to index 110. Not surprisingly the resulting basis matrix is very sparse thus it is not viable to take the inverse and need to add the $\Omega$
matrix in order to get the $\hat{\theta}$.But I am not sure how to calculate that , what is the range of integration? And is there libraries that I can use to calculate the derivatives of the basis functions?
Aside
I have tried a simpler model with natural cubic splines with degree of freedom = 4. And below is the code I used. The figure plots "radiation v.s. f(radiation)" where f(radiation) is the $f_{1}(X_{1})$ in my model and it is equal to $\boldsymbol{h}_{1}^{T}\boldsymbol{\theta}_{1}$ and the red lines are the confidence band obtained following the procedures in Section 5.2.2 in ESL.
require(stats); require(graphics)
library(splines)
library(ggplot2)
OZ=read.table("ozone.data",sep="", header = TRUE)
train = OZ[,c(2,3,4)]
yTrain=OZ$ozone
yTrain = scale(yTrain,center = TRUE, scale = TRUE)
X=cbind(ns(train$radiation,4),ns(train$temperature,4),ns(train$wind,4))
X = scale(X,center = TRUE, scale = TRUE)
X
df_2=data.frame(yTrain,X)
regression = lm(yTrain~.,data = df_2)
B = cbind(1,X)
N = length(yTrain)
p = 4*3+1
#########################
# Yhat, theta_hat #
yhat = B%*%(solve(t(B)%*%B))%*%(t(B))%*%(yTrain)
theta_hat = (solve(t(B)%*%B))%*%(t(B))%*%(yTrain)
#########################
# J #
J = (solve(t(B)%*%B))%*%(t(B))
#########################
sigma_val = (((t(yTrain - yhat))%*%(yTrain - yhat))/(N-p-1))[1]
V_y = sigma_val * diag(N)
#########################
# Sigma_hat #
Sigma_hat = (J%*%V_y)%*%(t(J))
#########################
#########################
# Radiation #
h_1 = B[,2:5]
f_1=h_1%*%theta_hat[2:5]
sigma_hat = Sigma_hat[2:5,2:5]
ycov=h_1%*%sigma_hat%*%t(h_1)
yvar=diag(ycov)
low=f_1-2*sqrt(yvar)
high=f_1+2*sqrt(yvar)
plot_df = as.data.frame(cbind(train$radiation,
f_1,high,low,yTrain))
ggplot() + geom_line(aes(V1,V2),plot_df)+ geom_line(aes(V1,V3),plot_df,color = "red")+
geom_line(aes(V1,V4),plot_df,color = "red")+
labs(y= "f(radiation)", x = "radiation")
