2
$\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$

1 Answer 1

1
$\begingroup$

Not sure I understand your question fully, but why are you use nknots = length(2:110)? that seems like a lot of over-fitting potentially. If you choose fewer knots, I think this will still be smooth but wont have the problem with sparesness that you mention.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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