Hello:
I hope that someone more knowledgeable will confirm or correct
what I'm about to say. I don't have time now to check this by studying
the "pspline" code.
From reading the pspline{survival} help page, I believe it uses
standard B-splines, "penalising the integrated second derivative", with
"theta" being the penalty or weight used to multiply the integral of the
square of the second derivative before adding it to the AIC to be
minimized. The help page also suggests that the basis functions (i.e.,
the knots) are evenly spaced, though I'm not completely sure about any
of this.
A basic rule for B-splines is as follows:
nbasis = order + n.knots
where n.knots = the number of interior knots, nbasis = the number of
basis functions, and order = the order of the spline = one more than the
degree of the polynomial segments. With the default cubic splines,
degree = 3, so order = 4.
In a fit, nbasis = the number of coefficients estimated, I think.
Consider the following modification of one of the examples in
help('pspline'):
coef(coxph(Surv(time, status) ~ ph.ecog + pspline(age,3), cancer))
ph.ecog ps(age)2 ps(age)3 ps(age)4 ps(age)5 ps(age)6 ps(age)7
ps(age)8
0.4480154 0.3221152 0.6470688 0.9342460 1.0289846 0.9639187 0.9004191
1.0119748
ps(age)9 ps(age)10 ps(age)11
1.2422085 1.6403940 2.0785682
This fit includes one coefficient for ph.ecog and 10 for
pspline(age,3). I infer from this that there are 10 basis functions.
Since degree = 3 (the default), order = 4, so the number of interior
knots = 6 ( = 10-4). To find where they are located, I would guess that
they are evenly spaced in range(cancer$age) = (39, 82):
seq(39, 82, length=8)
[1] 39.00000 45.14286 51.28571 57.42857 63.57143 69.71429 75.85714 82.00000
HOWEVER, these are only guesses; I hope someone else will either
confirm what I've said or correct any errors.
It would help to have a function "as.fd.coxph.penal" in the "fda"
package, as "fda" includes functions to help clarify this question:
"knots.fd" extracts the knots from an object of class "fd", answering
Federrico's question directly. Unfortunately, I don't have time to
write such. It also includes a function "TaylorSpline" for converting a
standard B-spline into the coefficients of a Taylor series approximation
for each segment of the spline. Chapter 3 in Ramsay, Hooker, and
Graves (2009) Functional Data Analysis with R and Matlab (Springer)
includes code for plotting a B-spline basis set, and the "fda" package
includes script files for working all but one of the 76 examples in that
book. You can find that script file without obtaining the book as follows:
fdaDir <- system.file(package='fda')
fdaDirs <- dir(fdaDir, full=TRUE)
fdaScripts <- grep('scripts$', fdaDirs, value=TRUE)
fdaCh <- dir(fdaScripts, full=TRUE)
fdaCh3 <- grep('fdarm-ch03.R$', fdaCh, value=TRUE)
You can then make a local copy of "fdaCh3" and work through it
line by line to see plots of the B-spline basis. This would help
Federico more if "fda" included a "as.fd.coxph.penal" function, because
then he could plot the spline basis used in his example.
Hope this helps.
Spencer
On 11/17/2010 11:39 AM, David Winsemius wrote:
On Nov 17, 2010, at 1:05 PM, Federico Calboli wrote:
Hi All,
I am trying to figure out how to get the position of the knots in a
pspline used in a cox model.
As far as I understand it, the term "knot" should be used with natural
splines, but not for penalized smoothing splines. See p 124 of
Therneau and Gramsch for better description and illustrations. If you
want to see the weighting of a particular pspline then this example
may help:
nn=rnorm(20)
matrix(pspline(nn), nrow=20)[order(nn), ]
# sorts the spline basis weights so the overlap is "visible"
my.model = coxph(Surv(agein, ageout, status) ~ pspline(x), mydata) #
x being continuous
How do I find out where the knot of the spline are? I would like to
know to figure out how many cases are there between each knot.
Best,
Federico
--
Federico C. F. Calboli
Department of Epidemiology and Biostatistics
Imperial College, St. Mary's Campus
Norfolk Place, London W2 1PG
Tel +44 (0)20 75941602 Fax +44 (0)20 75943193
f.calboli [.a.t] imperial.ac.uk
f.calboli [.a.t] gmail.com
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide
http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
