chjackson / survextrap

Survival Extrapolation With a Flexible Parametric Model and External Data
https://chjackson.github.io/survextrap/
GNU General Public License v3.0
7 stars 3 forks source link

computing uniform weights analytically #2

Closed irtimmins closed 1 year ago

irtimmins commented 1 year ago

Hi Chris,

with the mspline_uniform_weights() function, just an idea to suggest: given that B-splines have uniform weights by default, we can simply say that the scaling factor that takes M-splines to B-splines will be equivalent to the uniform weights. So uniform weights can be calculated analytically straight from the set of knots - current variance minimising approach may not be needed?

Here's an example attached


# example set of knots for M-splines.

degree <- 3
bknots <- c(0, 15)
iknots <- c(2,4,5,8,9,12)
# number of M-splines
K <- length(iknots)+degree+1

# define knot sequence with repetition of boundary knots
# (see Ramsay, 1988, www.jstor.org/stable/2245395)

knot.seq <- c(rep(bknots[1], degree+1), iknots, rep(bknots[2], degree+1))

# scaling factor from M-splines to B-splines is (knot.seq[i+degree+1]-knot.seq[i])

p_const <- rep(0,K)
for(i in 1:K){
  p_const[i] <- (knot.seq[i+degree+1]-knot.seq[i]) 
}
p_const <- p_const/sum(p_const)
p_const
[1] 0.03333 0.06667 0.08333 0.13333 0.11667 0.13333 0.16667 0.11667 0.10000 0.05000

# this should match weights estimated from variance minimising:

mspline_uniform_weights(iknots=iknots, bknots=bknots, degree=degree)
[1] 0.03329 0.06693 0.08301 0.13372 0.11640 0.13356 0.16633 0.11698 0.09975 0.05002
chjackson commented 1 year ago

Brilliant - well spotted, thank you! I shall put this in. Just wondering if the 4 in the code above should be degree+1 to generalise this to different degrees?

irtimmins commented 1 year ago

Yes that's right, thanks.

One small extra thing is that the sum of all the weights (knot.seq[i+degree+1]-knot.seq[i]) is (degree+1)*(bknots[2]-bknots[1]) since it telescopes with the only the repeated boundary knots remaining. So can write the uniform weights in single step.

Can make a small edit:

p_const <- rep(0,K)
rescale <- (degree+1)*(bknots[2]-bknots[1])
for(i in 1:K){
  p_const[i] <- (knot.seq[i+degree+1]-knot.seq[i])/(rescale) 
}
p_const
[1] 0.03333 0.06667 0.08333 0.13333 0.11667 0.13333 0.16667 0.11667 0.10000 0.05000