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

[Martin Maechler]

>>>>> "HenrikB" == Henrik Bengtsson <hb at maths.lth.se> writes:

    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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._