A couple days of detective work have gone along way to reveal how one  
may
recover the actual knots and coefficients of the splines used in ppr.
It's not finished yet, but I may be due some mid-course correction.

First, ... I was asked, ... why?

On Jun 27, 2006, at 11:50 PM, Prof Brian Ripley wrote:

> It is normal to report smooth curves via plots of smooth curves.  
> There are examples for ppr in MASS4 p.241 (referenced on the help  
> page).  This is done by the plot() method.
>
> Please do consult the references in the documentation: those who  
> very carefully documented these tools put them there because they  
> do give additional information.

This  reading and plotting of course was all done before contacting r- 
help.

The plots of the smooth curves will be used to communicate with most  
readers of our
publication.  However, reviewers for J. Geophysical Research and  
Atmospheric Environment
will require specifics of the fit performed.  We hope that the  
results will be used
to greatly improve our understanding of smog ozone formation, and for  
field use,
the pure spline fit is required, unencumbered by the ~2000 fitted  
points of our estimation
data set.

So, we have found that (my colleague Bob E's words)
> Both smooth.spline and ppr, with spline or gcvspline option,
> eventually come to a call to qsbart.f. Therefore, find out
> which parameters ppr is using in the qsbart call, and
> determine which parameters to smooth.spline cause the
> same arguments to qsbart.
>
>
> From ppr, qsbart is called from at
> the end of the chain (some wrappers left out):
> ppr -> smart1 -> {subfit, fulfit} -> onetrm -> oneone -> supsmu ->  
> spline -> qsbart
> Nearly all the code in the fortran routines deals with the
> supersmoother case.
>
> smooth.spline (stats/R/smspline.R) calls subroutine spline directly.

I have found that qsbart is  in this file  ..../library/stats/src/ 
sbart.c

and that the inputs and outputs seem to be clearly documented at the  
beginning
of the code:

> /* A Cubic B-spline Smoothing routine.
>
>    The algorithm minimises:
>
>         (1/n) * sum ws(i)^2 * (ys(i)-sz(i))^2 + lambda* int  
> ( s"(x) )^2 dx
>
>    lambda is a function of the spar which is assumed to be between  
> 0 and 1
>
> INPUT
> -----
>    penalt       A penalty > 1 to be used in the gcv criterion
>    dofoff       either `df.offset' for GCV or `df' (to be matched).
>    n            number of data points
>    ys(n)        vector of length n containing the observations
>    ws(n)        vector containing the weights given to each data point
>    xs(n)        vector containing the ordinates of the observations
>    ssw          `centered weighted sum of y^2'
>    nk           number of b-spline coefficients to be estimated
>                 nk <= n+2
>    knot(nk+4)   vector of knot points defining the cubic b-spline  
> basis.
>                 To obtain full cubic smoothing splines one might
>                 have (provided the xs-values are strictly increasing)
>    spar         penalised likelihood smoothing parameter
>    ispar        indicating if spar is supplied (ispar=1) or to be  
> estimated
>    lspar, uspar lower and upper values for spar search;  0.,1. are  
> good values
>    tol, eps     used in Golden Search routine
>    isetup       setup indicator [initially 0
>    icrit        indicator saying which cross validation score is to  
> be computed
>                 0: none ;  1: GCV ;  2: CV ;  3: 'df matching'
>    ld4          the leading dimension of abd (ie ld4=4)
>    ldnk         the leading dimension of p2ip (not referenced)
>
> OUTPUT
> ------
>    coef(nk)     vector of spline coefficients
>    sz(n)        vector of smoothed z-values
>    lev(n)       vector of leverages
>    crit         either ordinary or generalized CV score
>    spar         if ispar != 1
>    lspar         == lambda (a function of spar and the design)
>    iter         number of iterations needed for spar search (if  
> ispar != 1)
>    ier          error indicator
>                 ier = 0 ___  everything fine
>                 ier = 1 ___  spar too small or too big
>                         problem in cholesky decomposition

These are returned to the calling FORTRAN code (we should note the  
need to
respect FORTRAN<->C array and call conventions).

So it remains for us to find either the results of the last call to  
qsbart.f
made from  subroutine spline (n, x, y, w, smo, edf) in  the ppr.f
code "spline"
>       subroutine spline (n, x, y, w, smo, edf)
> c
> c------------------------------------------------------------------
> c
> c input:
> c    n : number of observations.
> c    x(n) : ordered abscissa values.
> c    y(n) : corresponding ordinate (response) values.
> c    w(n) : weight for each (x,y) observation.
> c output:
> c   smo(n) : smoothed ordinate (response) values.
> c   edf : equivalent degrees of freedom
> c
> c------------------------------------------------------------------

Note that coefficients are not passed back up from here.

The predict method in subroutine pppred(np,x,smod,y,sc) appears
to use an interpolation on fitted points, certainly a fast technique.

So it appears necessary to
(a) hack the code to return the last successfully used spline
parameters for each direction vector used.

or
(b) to use smooth.spline or some other technique to go back and
fit the fitted values all over again.  (Conceivably this would make
equivalent degrees of freedom more understandable to the non-ppr-user.)

smooth.spline (laudably!) returns the values:
> fit
> list for use by predict.smooth.spline, with components
> knot:
> the knot sequence (including the repeated boundary knots).
> nk:
> number of coefficients or number of “proper” knots plus 2.
> coef:
> coefficients for the spline basis used.
> min, range:
> numbers giving the corresponding quantities of x.

Now we appear to have the work of determining what was the
last, most successful use of sbart.c

Any help and correction would be gratefully accepted.

to wrap up the previous discussion ...
>
> On Tue, 27 Jun 2006, Robert Chatfield wrote:
>
>> The pursuit projection packages ppr is an excellent contribution  
>> to R. It is great for one-to-three ridge fits, often somewhat  
>> intuitive, and for multi-ridge fits, where it at least describes a  
>> lot of variance.
>>
>> Like many folk, I need to report the fits obtained from ppr to the  
>> greater, outside, non-R world.  It is fairly obvious how to use  
>> the terms alpha and beta to report on directionality and importance.

> In what sense do you claim it to be `undocumented'?  Methods should  
> do what the generic is documented to do, and that is the case here.
Fully documented inputs and outputs do indeed help.
Even better, a commented guide, deep in the algorithm, to the  
critical code that makes
and decides that we have best estimates at each direction.

> [...]
>
> -- 
> Brian D. Ripley,                  ripley@stats.ox.ac.uk
> Professor of Applied Statistics,  http://www.stats.ox.ac.uk/~ripley/
> University of Oxford,             Tel:  +44 1865 272861 (self)
> 1 South Parks Road,                     +44 1865 272866 (PA)
> Oxford OX1 3TG, UK                Fax:  +44 1865 272595

To conclude,  I am pleased that after all the statistics chasing  
after smog ozone
data (just check the various R datasets), that simplest-to-modest ppr  
studies have helped.

I have found the control over smoothness and sheer fitting ability
of ppr very helpful.  That I cannot push down degrees of freedom (df,  
gcvpen),  beyond a certain
point is a minor annoyance --  about which I have not enough  
information at this point.

Sometimes I need 2 or 3 directions, often I need only one,
I  even find that a spline fit on a conjectured, so-called "single  
index" variable,
(one direction, one term)  is all that I think I should credibly fit.
(Parsimony and physical insight are big-time review criteria in my  
field.)
Knowledge of the specific spline-fit parameters is crucial in these  
cases.
The plots will indicate that the listed equivalent degrees of freedom  
are appropriate.

It's really good idea to have all of this in one package.

Maybe we can work toward a confluence of techniques to address these  
needs.

Thanks again for ppr and the MASS book which led me to it.

Dr. Robert Chatfield
Earth Sciences, MS 245-5
NASA Ames Research Center
Moffett Field, CA 94035    USA

Ph: 650-604-5490  FAX 650-604-3625

http://geo.arc.nasa.gov/sgg/chatfield


	[[alternative HTML version deleted]]

