#' @title Generate design matrix for periodic restricted cubic spline
#'
#' @description
#' Generate design matrix for periodic restricted cubic spline.
#'
#' @param x numerical x values to transform to new basis
#' @param knots vector with locations of the knots of the spline
#' @param nk number of knots, used only if the knots are not specified, overridden otherwise
#' @param xmin value of the (theoretical) minimum of x
#' @param xmax value of the (theoretical) maximum of x
#'
#' @examples
#' # load example data; see help("viral_east_mediteranean")
#' data("viral_east_mediteranean")
#'
#' # calculate location of knots to use
#' Knots <-
#' Hmisc::rcspline.eval(x = viral_east_mediteranean$EpiWeek,
#' nk = 5, knots.only = TRUE)
#'
#' # model viral infections vs weeks
#' model <- glm(RSV ~ rcs_per(EpiWeek, knots = Knots), data = viral_east_mediteranean)
#'
#' # plot model (with many points, to make it smooth)
#' plot_per_mod(Model = model, XvarName = "EpiWeek")
#'
#' @export
rcs_per <- function(x,
knots = NULL,
nk = 5,
xmin = min(x, na.rm=TRUE),
xmax = max(x, na.rm=TRUE)){
# derive the knot locations as in the rms package, if they are not provided
if(is.null(knots)) {
knots <- rcspline.eval(x, nk = nk, knots.only = TRUE)
}
# check if the number of knots if at least 4
nk <- length(knots)
if (nk < 4) stop("To use the periodic RCS you must specify at least 4 knots")
b.x.all <- b_rcs(x, knots) # matrix with the expansion of the splines, n*(k-2)
b.xmax.all <- b_rcs(xmax, knots) # value of the spline for x=xmax, vector 1*(k-2)
b.prime.x.all <- b_rcs_prime(x, knots) # matrix with the first derivative of the expansion of the splines, n*(k-2)
b.prime.xmax.all <- b_rcs_prime(xmax, knots) # vector with the value of the first derivative for x=xmax, vector 1*(k-2)
# terms to add in gamma1
gamma1.c1 <- (xmin-xmax) / b.xmax.all[nk-2] # constant
gamma.f1 <- gamma1.c1 * b.x.all[, nk-2] # vector n*1
gamma.f2 <- gamma1.c1 * b.prime.xmax.all[nk-2]
denom.beta.kMinus3 <- b.prime.xmax.all[nk-3] - b.xmax.all[nk-3] / b.xmax.all[nk-2] * b.prime.xmax.all[nk-2]
beta.j <- matrix(NA, ncol=nk-3, nrow = length(x))
for(j in 1:(nk-3)) {beta.j[,j] <- b.x.all[,j] - b.xmax.all[j] / b.xmax.all[nk-2] * b.x.all[,nk-2]}
if(nk > 4) {
beta.c2.j <- matrix(NA, ncol = nk-4, nrow = length(x))
for(j in 1:(nk-4)) {beta.c2.j[,j] <- b.prime.xmax.all[,j] - b.xmax.all[j] / b.xmax.all[nk-2] * b.prime.xmax.all[nk-2]}
}
num.beta.kMinus3 <- beta.j[, nk-3]
# construct dataframe with the result
if (nk > 4) {
result <-
cbind(x+gamma.f1-gamma.f2/denom.beta.kMinus3*num.beta.kMinus3,
beta.j[,1:(nk-4)]-beta.c2.j[,1:(nk-4)]*num.beta.kMinus3/denom.beta.kMinus3
)
} else {
result <- x+gamma.f1-gamma.f2/denom.beta.kMinus3*num.beta.kMinus3
}
# saves the knot locations as attributes
attr(result, "knots") <- knots
# return result
return(result)
} # end rcs_per
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.