IDEMSInternational / R-Instat

A statistics software package powered by R
http://r-instat.org/
GNU General Public License v3.0
38 stars 103 forks source link

Periodic splines for our Markov chain fitting? #5911

Open rdstern opened 4 years ago

rdstern commented 4 years ago

@dannyparsons and @jkmusyoka I keep re-discovering periodic splines that I suggest we should consider as an alternative to the current circular functions for our Markov chain modelling.

a) There is a package called splines, with a function periodicSpline, which sounds what we need. Unfortunately it has been removed from CRAN. I assume we could still extract the function from there if needed. (There is a package called splines2 that is available. It has lots, but not -it seems - periodic splines. It says is is actually "splines too" and not intended as an update for the splines package. b) There is a one-function package, called pbs (periodic b spline). It has not been updated since 2013, but perhaps doesn't need to be! c) there is the splinesfun function in the stats package. It also refers to the now sort of defunct splines package. It has a "periodic" option, which might be what we need?

wenjie2wang commented 3 years ago

@rdstern You may be interested in splines2::mSpline() from the latest version of splines2, which returns periodic M-splines if periodic=TRUE.

rdstern commented 3 years ago

I had missed acting on this helpful message from the maintainer of the splines2 package.

I propose we now include splines2 in R-Instat instead of the splines package. Can we detect where (if) we currently use splines. It has been withdrawn from CRAN.

We also include an example of the use of this package in geom_smooth

dannyparsons commented 2 years ago

We're currently using pbs::pbs for periodic splines and splines::ns for natural splines. So we definitiely need to replace splines::ns and happy to consider splines2 as well.

@wenjie2wang thanks for the info on splines2. Are you aware of pbs::pbs() and the difference between this and your package?

wenjie2wang commented 2 years ago

@wenjie2wang thanks for the info on splines2. Are you aware of pbs::pbs() and the difference between this and your package?

The pbs::pbs() constructs periodic splines based on B-splines, while mSpline() with periodic = TRUE provides periodic splines based on M-splines, which are essentially different in scales. See the following example for illustration.

library(pbs)
library(splines2)

x <- seq.int(0, 1, 0.01)
bound_knots <- c(0, 1)

mat1 <- pbs(x, knots = c(0.25, 0.5, 0.75), degree = 3, intercept = TRUE,
            Boundary.knots = bound_knots)
mat2 <- mSpline(x, knots = c(0.25, 0.5, 0.75), degree = 3, intercept = TRUE,
                Boundary.knots = bound_knots, periodic = TRUE)

op <- par(mfrow = c(1, 2))
matplot(x, mat1, type = "l")
matplot(x, mat2, type = "l")

In addition, I found that pbs::pbs() could be too restricted when checking the number of knots.

mat1 <- pbs(x = x, knots = c(0.25, 0.75), degree = 3, intercept = TRUE,
            Boundary.knots = bound_knots)
#> Error in getPeriodBsplineKnots(sort(c(Boundary.knots, knots)), degree = degree): number of internal knots(no boudary) must be greater than or equal to 2+ degree
mat2 <- mSpline(x, knots = c(0.25, 0.75), degree = 3, intercept = TRUE,
                Boundary.knots = bound_knots, periodic = TRUE)
par(op)
matplot(x, mat2, type = "l")

Created on 2021-11-13 by the reprex package (v2.0.1)