Nothing
## (c) Natalya Pya (2024). Provided under GPL 2.
## routines for univariate local SCOP-spline construction, LSOP-spline,
## in collaboration with Jens Lichter and Thomas Kneib, University of Gottingen
#####################################################
### Local monotone increasing LSCOP-spline construction, increasing up until the change point, xc,
######################################################
## building B-spline bases functions over the whole range of x
## and making m knots at xc change point ...
smooth.construct.lmpi.smooth.spec<- function(object, data, knots)
## construction of the monotone increasing smooth
## xc specifies a change point
## scop-spline up until xc and unconstrained p-spline from xc,
## using ceiling(q/2) basis functions for both parts
{
m <- object$p.order[1]
if (is.na(m)) m <- 2 ## default for cubic spline
if (m<1) stop("silly m supplied")
if (object$bs.dim<0) object$bs.dim <- 10 ## default
q <- object$bs.dim
## q1 <- q2 <- ceiling(q/2) ## basis dimention for each part
xc <- object$xt$xc ## a change point
rg <- range(x)
share <- (xc-rg[1])/(rg[2]-rg[1]) ## share of the x values up untill the change point
q1 <- max(ceiling(q*share),5) ## basis dimention for the constrained part
q2 <- max((q-q1),5)
n <- length(x)
nk <- q1+q2+1 ## number of knots
if (nk<=0) stop("k too small for m")
x <- data[[object$term]] ## the data
xk <- knots[[object$term]] ## will be NULL if none supplied
if (!is.null(xk)) ## if not NULL
stop(paste("'lmpi' smooth currenly does not work with user-supplied knots"))
if (is.null(xc))
stop(paste("a change point 'xc' is not supplied"))
## getting equally spaced knots from both sides of the change point...
xk <- rep(NA,q1+q2+1+1)
xk[(m+2):(q1+1)] <- seq(rg[1],xc,length=q1-m) ## inner knots for the constrained part (to the left from the change point)
xk[(q1+1):(q1+q2-m+1)] <- seq(xc,rg[2],length=q2-m+1) ## inner knots for the unconstrained part (to the right from the change point), add one knot more than for the constrained part, as otherwise there 5 basis functions for the constrained part and only 3 for the unconstrained
for (i in 1:(m+1)) ## outer knots on the left-hand-side of the x range
xk[i] <- xk[m+2]-(m+2-i)*(xk[m+3]-xk[m+2])
for (i in (q1+q2-m+2):(q1+q2+2)) ## outer knots on the right-hand-side of the x range
xk[i] <- xk[q1+q2-m]+(i-q1-q2+m)*(xk[q1+2]-xk[q1+1])
xk <- c(xk[which(xk < xc)],rep(xc,m),xk[which(xk > xc)]) ## including xc change point m times
if (!is.null(object$point.con[[1]])) ## a point constraint is supplied?
stop(paste("'lmpi' smooth does not work with a point constraint; use 'miso' for a start-at-zero constraint, or 'mifo' for a finish-at-zero constraint"))
## get model matrix...
X <- splineDesign(xk,x,ord=m+2)
Sig1 <- matrix(1,q1,q1) ## coef summation matrix
Sig1[upper.tri(Sig1)] <- 0
q.t <- ncol(X) ## number of coefficients of the final lscop-spline, which is (q-2)
Sig <- matrix(0,q.t,q.t)
Sig[1:q1,1:q1] <- Sig1
Sig[(q1+1):q.t,(q1+1):q.t] <- diag(1,q.t-q1)
X <- X%*%Sig
## object$X <- X
## applying sum-to-zero (centering) constraint...
cmx <- colMeans(X)
X <- sweep(X,2,cmx) ## subtract cmx from columns
object$X <- X # the final model matrix
object$cmX <- cmx
object$P <- list()
object$S <- list()
object$Sigma <- Sig
if (!object$fixed)
{
# P <- diff(diag(q.t-1),difference=1)
# P <- rbind(rep(0,q.t-1),P) ## adding 1st row of zeros
# P <- cbind(rep(0,q.t-1),P) ## adding first column of zeros
## making 1st order difference penalty for the constrained part and 2nd order diff-s for the unconstrained
P <- matrix(0,q.t-3,q.t)
d1 <- diff(diag(q1),difference=1)
P[2:(q1),2:(q1+1)] <- d1
d1 <- diff(diag(q.t-q1-1),difference=2)
P[(q1+1):(q.t-3),(q1+2):q.t] <- d1
object$P[[1]] <- P
object$S[[1]] <- crossprod(P)
## create a block diagonal penalty matrix, penalizing separately constrained and unconstrained parts
## block diagonal penalty showed slightly worse mse when compared on two simul. examples ...
## penalty for the constrained first part
# P <- matrix(0,q.t-2,q.t) # matrix(0,q.t-1,q.t)
# d1 <- diff(diag(q1-1),difference=1)
# P[2:(q1-1),2:q1] <- d1
# # P[q1:(q.t-1),(q1+1):q.t] <-diag(1,(q.t-q1))
# P[q1:(q.t-2),(q1+1):(q.t-1)] <-diag(1,(q.t-q1-1))
# object$P[[1]] <- P
# object$S[[1]] <- crossprod(P)
# ## block for the unconstrained second part...
# P <- matrix(0,q.t-2,q.t) # matrix(0,q.t-1,q.t)
# d1 <- diff(diag(q.t-q1),difference=2)
# P[(q1+1):(q.t-2),(q1+1):q.t] <- d1
# object$P[[2]] <- P
# object$S[[2]] <- crossprod(P)
}
object$p.ident <- c(FALSE,rep(TRUE,q1-1),rep(FALSE,q.t-q1)) ## p.ident is an indicator of which coefficients must be positive (exponentiated)
object$rank <- ncol(object$X) # penalty rank
object$null.space.dim <- 2 ## ##m+1 # dim. of unpenalized space, 2 as the basis of a straight line is two-dimensional
object$C <- matrix(0, 0, ncol(X)) # to have no other constraints
object$knots <- xk;
object$m <- m;
object$df<-ncol(object$X) # maximum DoF
object$q1 <- q1
## get model matrix for 1st and 2nd derivatives of the smooth...
h <- xk[m+3]-xk[m+2] ## distance between two adjacent knots
object$Xdf1 <- splineDesign(xk,x,ord=m+1)[,1:(q.t-1)]/h ## ord is by one less for the 1st derivative
object$Xdf2 <- splineDesign(xk,x,ord=m)[,1:(q.t-2)]/h^2 ## ord is by two less for the 2nd derivative
class(object)<-"lmpi.smooth" # Give object a class
object
}
## Prediction matrix for the `lmpi` smooth class...
Predict.matrix.lmpi.smooth<-function(object,data)
## prediction method function for the `mpi' smooth class
{ m <- object$m # spline order, m+1=3 default for cubic spline
q <- object$df
q1 <- object$q1 ## basis dimention for the constrained part
x <- data[[object$term]]
Sig1 <- matrix(1,q1,q1) ## coef summation matrix
Sig1[upper.tri(Sig1)] <-0
Sig <- matrix(0,q,q)
Sig[1:q1,1:q1] <- Sig1
Sig[(q1+1):q,(q1+1):q] <- diag(1,q-q1)
## find spline basis inner knot range...
ll <- object$knots[m+2];ul <- object$knots[length(object$knots)-m-1]
n <- length(x)
ind <- x<=ul & x>=ll ## data in range
if (sum(ind)==n) ## all in range
X <- spline.des(object$knots,x,m+2)$design
else { ## some extrapolation needed
## matrix mapping coefs to value and slope at end points...
D <- spline.des(object$knots,c(ll,ll,ul,ul),m+2,c(0,1,0,1))$design
X <- matrix(0,n,ncol(D)) ## full predict matrix
if (sum(ind)> 0) X[ind,] <- spline.des(object$knots,x[ind],m+2)$design ## interior rows
## Now add rows for linear extrapolation...
ind <- x < ll
if (sum(ind)>0) X[ind,] <- cbind(1,x[ind]-ll)%*%D[1:2,]
ind <- x > ul
if (sum(ind)>0) X[ind,] <- cbind(1,x[ind]-ul)%*%D[3:4,]
}
X <- X%*%Sig
X <- sweep(X,2,object$cmX)
X
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.