[R] how lambda is computed in smoot.spline given _df_

Martin Maechler maechler at stat.math.ethz.ch
Mon Nov 8 14:51:39 CET 2004

>>>>> "Eryk" == Eryk Wolski <wolski at molgen.mpg.de>
>>>>>     on Mon, 08 Nov 2004 10:39:26 +0100 writes:

    Eryk> Hi,

    Eryk> I posted some days ago a question concerning the
    Eryk> computation of lambda in the smooth.spline function
    Eryk> (which I repreat at the bottom of the mail) given _df_.

I had thought one should be able to figure it out.
The help page cannot give all the theory behind smoothing splines!
The 'References' of  help(smooth.spline) do explain how "df" is
defined, and even the help page mentions that 
df = trace( <smoother matrix> ). 

    Eryk> Unfortunately the documentation is not clear to me.
    Eryk> Maybee someone can help to answer in my view the
    Eryk> basic question:

    Eryk> If the penalized log likelihood is

    Eryk>  L = (y - f)' W (y - f) + lambda c' Sigma c 
    Eryk>  how the _lambda_ in the above equation is computed if _df_ is
    Eryk>  given and _spar_ not?

Well, we have a computer for solving such problems:
spar (or lambda, equivalently) is determined such that the
resulting df matches the desired  df  by a traditional zero-finder
(Brent).  Basically the same algorithm as used by uniroot().

    Eryk> And, Is there a way to define lambda directly?

Almost {note: this is high school math}:
Since the help page gives the linear relation
between log(lambda) and spar, you call smooth.spline() once with
e.g. lambda = 0.5; and then know the formula how to compute spar
from lambda:

Using the 'cars' example from  help(smooth.spline):
assume you want to set  

my.lambda <- 1e-4

data(cars)# needed in R versions before 2.0.0
## Call it with an arbitrary fixed 'spar':
(sspl <- smooth.spline(speed, dist, spar= .5))

## now solve the linear relationship for the only unknown  s0 :

s0 <- with(sspl, spar - 0.0601*log(lambda))

ssp. <- smooth.spline(speed, dist, spar= s0 + 0.0601*log(my.lambda))

str(ssp.)## confirms that lambda *is* 1e-4 = my.lambda


Hoping this helped:
Martin Maechler Seminar fuer Statistik, ETH Zurich

More information about the R-help mailing list