Skip to content

How to get the penalized log likelihood from smooth.spline()?

2 messages · Henrik Bengtsson, Martin Maechler

Henrik Bengtsson
# 
I use smooth.spline(x, y) in package modreg and I would like to get
value of penalized log likelihood and preferable also its two parts. To
make clear what I am asking for (and make sure that I am asking for the
right thing) I clarify my problem trying to use the same notation as in
help(smooth.spline):

I want to find the natural cubic spline f(x) such that

   L(f) = \sum_{k=1}{n} w[k](y[k] - f(x[k])^2 + \lambda J(f)

is minimized, where J(f) := \int f''(t)^2 dt is the quadratic roughness
functional use. Since J(f) is quadratic one can find a matrix \Sigma such
that J(g) = c^T{\Sigma}c where c is the vector of spline coefficients.
With J(f) defined as above the elements of \Sigma becomes

  \Sigma_{ij} = \int \beta_i''(t)\beta_j''(t) dt

where \beta(t) is the vector of B-spline base functions. Finally, writing
the matrix W as W := diag(\sqrt{w}) one can write L(f) as

  L(f) = (y - f)^T W^2 (y - f) + \lambda c^T{\Sigma}c

which is the form used in help(smooth.spline). So back to my question,
using smooth.spline(), how can I get 

  1) the value of L(f(x)) for a my fitted (x,y) data,
  2) the value of roughness penalty c^T{\Sigma}c,
  3) the value of (y - f)^T W^2 (y - f)?

Does smooth.spline() return any of these directly or indirectly?

Thanks

Henrik Bengtsson

Dept. of Mathematical Statistics @ Centre for Mathematical Sciences 
Lund Institute of Technology/Lund University, Sweden (+2h UTC)
Office: P316, +46 46 222 9611 (phone), +46 46 222 4623 (fax)
h b @ m a t h s . l t h . s e
http://www.maths.lth.se/bioinformatics/





-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !)  To: r-help-request at stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
6 days later
Martin Maechler
# 
HenrikB> I use smooth.spline(x, y) in package modreg and I
    HenrikB> would like to get value of penalized log likelihood
    HenrikB> and preferable also its two parts. To make clear
    HenrikB> what I am asking for (and make sure that I am
    HenrikB> asking for the right thing) I clarify my problem
    HenrikB> trying to use the same notation as in
    HenrikB> help(smooth.spline):


    HenrikB> I want to find the natural cubic spline f(x) such that
    HenrikB>      L(f) = \sum_{k=1}{n} w[k](y[k] - f(x[k])^2 + \lambda J(f)

    HenrikB> is minimized, where J(f) := \int f''(t)^2 dt is the
    HenrikB> quadratic roughness functional use. Since J(f) is
    HenrikB> quadratic one can find a matrix \Sigma such that
    HenrikB> J(g) = c^T{\Sigma}c where c is the vector of spline
    HenrikB> coefficients.  With J(f) defined as above the
    HenrikB> elements of \Sigma becomes

    HenrikB>   \Sigma_{ij} = \int \beta_i''(t)\beta_j''(t) dt

    HenrikB> where \beta(t) is the vector of B-spline base
    HenrikB> functions. Finally, writing the matrix W as W :=
    HenrikB> diag(\sqrt{w}) one can write L(f) as

    HenrikB>     L(f) = (y - f)^T W^2 (y - f) + \lambda c^T{\Sigma}c

    HenrikB> which is the form used in help(smooth.spline). So
    HenrikB> back to my question, using smooth.spline(), how can
    HenrikB> I get

    HenrikB> 1) the value of L(f(x)) for a my fitted (x,y) data,
    HenrikB> 2) the value of roughness penalty c^T{\Sigma}c,
    HenrikB> 3) the value of (y - f)^T W^2 (y - f)?

    HenrikB> Does smooth.spline() return any of these directly or indirectly?

it does so, partly as the `$pen.crit' and `$crit' components of
the result. Since the result also contains lambda, the residuals
and the leverages, it should be pretty straightforward to
compute all three quantities from one of them; I think
res$pen.crit == "1)".

smooth.spline is more extensively documented since R version 1.4
almost exactly using your notation above.
Are you using R version 1.4.1 and have looked at   help(smooth.spline)  ?

{try the result of  help(smooth.spline, offline = TRUE)
 since the Math will be nicer}.

By looking at the source of smooth.spline() you may learn even
more.
I hope this helps; I don't have time at the moment work at the
details.

Martin Maechler <maechler at stat.math.ethz.ch>	http://stat.ethz.ch/~maechler/
Seminar fuer Statistik, ETH-Zentrum  LEO C16	Leonhardstr. 27
ETH (Federal Inst. Technology)	8092 Zurich	SWITZERLAND
phone: x-41-1-632-3408		fax: ...-1228			<><
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !)  To: r-help-request at stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._