1
$\begingroup$

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

radiation

$\endgroup$

Your Answer

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

Browse other questions tagged or ask your own question.