5
$\begingroup$

I'm trying to obtain the smoother matrix out of my smooth.spline fit. I use the bone mineral density data from http://statweb.stanford.edu/~tibs/ElemStatLearn/.

bone <-read.table("bone.data", header=TRUE)
bmd_age <- smooth.spline(bone$age, bone$spnbmd, all.knots=TRUE, cv=TRUE)
bmd_fit <- predict(bmd_age, sort(bone$age))
df <- bmd_age$df

To obtain a column of the smoother matrix, I can replace the response vector (bone$spnbmd) by a vector with a single 1 and the rest filled with 0's. It is both what the professor recommended and what I found online https://stat.ethz.ch/pipermail/r-help/2006-June/108471.html.

So I use

smooth.matrix = function(x){
  n = length(x);
  sm = matrix(0, n, n);  
  for(i in 1:n){
    y = rep(0, n); y[i]=1;
    sm_i = predict(smooth.spline(x, y, df=df),x)$y;
    sm[,i]= sm_i;
  }
  return(sm)
}

sm <- smooth.matrix(bone$age)

If the smoother matrix is correct, the following two quantities should be the same (both fitted values from the smoothing spline model).

fromsm <- sm%*%(bone$spnbmd[order(bone$age)])
fromfit <- bmd_fit$y

However, they are not. I think the problem is in the definition of smooth.matrix function, where 

sm_i = predict(smooth.spline(x, y, df=df),x)$y;

is not using the same smoothing fit as in bmd_age. I've tried fixing the degree of freedom, spar, lambda, cv=FALSE, etc. but no luck so far. How to fix it?

Improve this question
$\endgroup$
6
  • 1
    $\begingroup$ I have done considerable theoretical statistical work on bone mineral density, but have not a clue as to what you are doing. Let's start by asking you to tell us, in simple terms, what you are trying to accomplish. Then, maybe we can figure out what you are stuck on. $\endgroup$
    – Carl
    Commented Sep 13, 2016 at 2:38
  • $\begingroup$ This is a homework assignment. It's not about what I want to find about bone mineral density. bmd_age is the smoothing spline model I obtained from smooth.spline(). I just want to get the correct smoother matrix associated with this model. $\endgroup$
    – Badoe
    Commented Sep 13, 2016 at 2:49
  • 1
    $\begingroup$ If you just want the basis matrix for a cubic (or general higher order spline), you can get this directly with the function bs. $\endgroup$
    – Cliff AB
    Commented Sep 13, 2016 at 2:53
  • $\begingroup$ @CliffAB What I want is not the basis matrix, but the smoother matrix. Any smoothing spline model can be written as $\hat y=S_\lambda \cdot y$, where $S_\lambda$ only depends on my input $x$ and the smoothing parameter $\lambda$. I'm trying to get $S_\lambda$ out of the model. $\endgroup$
    – Badoe
    Commented Sep 13, 2016 at 2:55
  • 1
    $\begingroup$ I have a very hard time believing smooth.spline does not depend on y. In fact, simple experiments demonstrate it produces dramatically different output when you fix x and vary y. Your problem is that smooth.spline does too much: it selects a different $\lambda$ for each vector y. You need to suppress that behavior by overriding the selection. See the help page for the spar and cv arguments. $\endgroup$
    – whuber ♦
    Commented Sep 13, 2016 at 16:02

3 Answers 3

Reset to default
3
$\begingroup$

After many hours of exploration, here is what I found:

Because smooth.spline algorithm chooses spar instead of lambda, it is only possible to (sort of) fix spar. However, lambda is a function of spar and another variable matrix. So fixing spar does not fix lambda necessarily. I have not found an easy way to extract the smoother matrix out of smooth.spline. However, for the purpose of computing variance, the algorithm provided in https://stat.ethz.ch/pipermail/r-help/2006-June/108471.html (fix spar instead of df) is a close estimate of the true smoother matrix. The variance computed from $SyS^T$, where $S$ is the estimated smoother matrix, is pretty close to the one computed from the correct smoother matrix.

Another R package "assist" has a function "ssr()" that also does smooth spline regression. It is not as powerful as smooth.spline. But the built-in function "hat.ssr()" gives the true smoother matrix of the model obtained from "ssr()".

Improve this answer
$\endgroup$
2
$\begingroup$

The above answers only approximate the smoothing matrix. Here is a solution that will get you the exact smoothing matrix from the r function smooth.spline(). The key is to recognize that the smoothing matrix is only a function of the values of $x$ and the penalization parameter $\lambda$, allowing us to smooth a vector $\tilde{y} = (0,0, ..., 0, 1, 0,...,0)^{T}$, and therefore get each column of the smoothing matrix.

library(splines)
x = seq(0, 100, by=0.1)
y = x*sin(x) + rnorm(length(x), 0, 0.1)

#use cross-validation to choose best smoothing parameter
spar = seq(0.01, 1, by = 0.01)
cv = rep_len(NA, length(spar))
for(i in 1:length(spar)){
    tempfit = smooth.spline(x, y, spar = spar[i], cv=TRUE, all.knots = TRUE)
    cv[i] = tempfit$cv.crit
}

#use the optimal smoothing parameter to produce a final fit
fit = smooth.spline(x, y, spar = spar[which(cv == min(cv))], cv=TRUE, all.knots = TRUE)

#calculate the smoothing matrix
L = matrix(nrow = length(x), ncol = length(x))
for(j in 1:length(x)){
    yi = rep_len(0, length(x))
    yi[j] = 1
    L[,j] = predict(smooth.spline(x, yi, lambda = fit$lambda, cv=TRUE,
                                  all.knots = TRUE), x)$y
}

The matrix $L$ is the resulting smoothing matrix.

Improve this answer
$\endgroup$
2
$\begingroup$

The accepted answer isn't correct here - smooth.matrix is working just fine.

The only reason fromsm and fromfit aren't the same in the above example is because of misplaced parentheses.

Replace fromsm <- sm%*%(bone$spnbmd[order(bone$age)]) with fromsm <- (sm%*%bone$spnbmd)[order(bone$age)] and they are the same.

Improve this answer
$\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.