# R/unireg.R In uniReg: Unimodal Penalized Spline Regression using B-Splines

```negloglikFREQ <- function(lambda,tBSB,tySB,sigmaest,tbetaO,tbetaObeta,Dtilde,beta0,Om,rangO,Cm,constr){
################################################################################
## Calculates the value of the negative (restricted) log-likelihood of lambda. #
################################################################################
Einv <- tBSB + exp(lambda)*Om
E <- solve(Einv)
e <- as.numeric((tySB + exp(lambda)*tbetaO)%*%E)
if(constr!="none"){
prob2 <- pmvnorm(lower=0, mean=as.numeric(Cm%*%e), sigma=Cm%*%E%*%t(Cm))[[1]]
term2 <- -log(prob2)
Omtilde <- Dtilde + exp(lambda)*Om
prob3 <- pmvnorm(lower=0, mean=as.numeric(Cm%*%beta0), sigma=Cm%*%solve(Omtilde)%*%t(Cm))[[1]]
term3 <- log(prob3)
}else{term2 <- 0
term3 <- 0
}
term1 <- 0.5*log(det(mean(sigmaest)*Einv)) - (rangO/2)*lambda - 0.5*t(e)%*%Einv%*%e + 0.5*exp(lambda)*tbetaObeta #

negloglik <- term1+term2+term3

if(term3==-Inf){negloglik <- .Machine\$double.xmax}
if(negloglik==Inf){negloglik <- .Machine\$double.xmax}
if(negloglik==-Inf){negloglik <- -.Machine\$double.xmax}
return(negloglik)
}

unisplinem <- function(m,tBSB,tySB,sigmaest,tbetaO,tbetaObeta,Dtilde,rangO,B,beta0,Om,constr,inverse,penalty,tuning){
######################################################################################
# Optimizes the tuning factor lambda and calculates the constrained (with mode m) or #
# unconstrained estimate of the coefficient vector.                                  #
######################################################################################
d <- dim(B)[2]
Cm <- inverse*unimat(d,m)          # Cm is the (transposed) matrix of constraints on the beta coefficients if the mode is m

if(penalty=="none"){lambdaopt <- 0
}else{
# if tuning=TRUE, lambda ist optimized with constr=constr, otherwise with constr="none"
if(!tuning){
opt <- optimize(negloglikFREQ,interval=c(3,10),tBSB=tBSB,tySB=tySB,sigmaest=sigmaest,tbetaO=tbetaO,tbetaObeta=tbetaObeta,Dtilde=Dtilde,beta0=beta0,Om=Om,rangO=rangO,Cm=Cm,constr="none")
lambdaopt <- exp(opt\$minimum)
}else{
opt <- optimize(negloglikFREQ,interval=c(3,10),tBSB=tBSB,tySB=tySB,sigmaest=sigmaest,tbetaO=tbetaO,tbetaObeta=tbetaObeta,Dtilde=Dtilde,beta0=beta0,Om=Om,rangO=rangO,Cm=Cm,constr=constr)
lambdaopt <- exp(opt\$minimum)
}
}

# constructing E_{lambda}^{-1}, E_{lambda} and e_{lambda} with the optimal lambda
Einv <- tBSB+ lambdaopt*Om
tySBlBetaO <- as.numeric(tySB + lambdaopt*tbetaO)
e <- solve(Einv,tySBlBetaO,system=Einv)

# matrices for solve.QP.compact
Aneu <- unimatind(d,m)
Amat <- inverse*Aneu\$Amat
Aind <- Aneu\$Aind

sc <- norm(Einv,"2")
# constrained coefficients are estimated with solve.QP.compact; otherwise the vector e is the unconstrained solution
if(constr!="none"){coeff <- try(solve.QP.compact(Dmat=Einv/sc,dvec=tySBlBetaO/sc,Amat=Amat,Aind=Aind)\$solution, silent=TRUE)
if(class(coeff)=="try-error" || any(is.nan(coeff))){
coeff <- solve.QP.compact(Dmat=t(solve(chol(Einv/sc))),dvec=tySBlBetaO/sc,Amat=Amat,Aind=Aind, factorized=TRUE)\$solution}
}else{coeff <- e}

fit <- as.numeric(B%*%coeff)
return(list(coef=coeff,fitted.values=fit,lambdaopt=lambdaopt))
}

unireg <- function(x, y, w=NULL, sigmasq=NULL, a=min(x), b=max(x), g=10, k=3, constr=c("unimodal","none","invuni","isotonic","antitonic"),
penalty=c("diff", "none", "sigEmax", "self", "diag"), Om=NULL, beta0=NULL, coinc=NULL, tuning=TRUE, abstol=0.01,vari=5,ordpen=2,m=1:(g+k+1),
allfits=FALSE, nCores=1){
#######################################################
# Applies the specified spline regression to x and y. #
#######################################################
n <- length(x)
if(!is.vector(x,mode="numeric")){stop("x should be a numeric vector.")}
if(!is.vector(y,mode="numeric")){stop("y should be a numeric vector.")}
if(!is.null(w) && !is.vector(w,mode="numeric")){stop("w should be NULL or a numeric vector.")}
if(!is.null(sigmasq) && !is.vector(sigmasq,mode="numeric")){stop("sigmasq should be NULL or a numeric vector.")}
if(!is.null(sigmasq) && length(sigmasq)!=1 && length(sigmasq)!=n){stop("sigmasq a numeric vector of length 1 or same length as x.")}
if(!(is.vector(a,mode="numeric") && length(a)==1 && is.finite(a))){stop("a should be a finite numeric vector of length 1.")}
if(!(is.vector(b,mode="numeric") && length(b)==1 && is.finite(b))){stop("b should be a finite numeric vector of length 1.")}
if(!(g%%1==0 && g>=0 && is.finite(g))){stop("g should be a finite whole number >=0.")}
if(!(k%%1==0 && k>=0 && is.finite(k))){stop("k should be a finite whole number >=0.")}
constr <- match.arg(constr)
penalty <- match.arg(penalty)
if(!is.null(Om) && !(is.matrix(Om) && isTRUE(all.equal(dim(Om),c(g+k+1,g+k+1))))){stop("Om should be NULL or a (g+k+1)x(g+k+1) matrix.")}
if(!is.null(beta0) && !is.vector(beta0,mode="numeric")){stop("beta0 should be NULL or a numeric vector.")}
if(!is.null(coinc) && !is.logical(coinc)){stop("coinc should be NULL, TRUE or FALSE.")}
if(!is.logical(tuning)){stop("tuning should be TRUE or FALSE.")}
if(!(is.vector(abstol,mode="numeric") && length(abstol)==1 && is.finite(abstol)) && !is.null(abstol)){stop("abstol should be NULL or a finite numeric vector of length 1.")}

if(!is.null(sigmasq) && !(length(sigmasq)==1 || all(sigmasq==sigmasq[1])) && !is.null(abstol)){stop("Cannot estimate sigmasq in case of heteroscedasticity. Set abstol to NULL.")}
if(is.null(sigmasq) && is.null(abstol)){stop("Either sigmasq or abstol have to be non-NULL.")}
if(is.null(sigmasq) && !is.null(abstol) && !all(table(x)>=2)){stop("Provide a startvalue for sigmasq or at least 2 repeated measurements for each x-value.")}
if(b<=a){stop("[a,b] is not a proper interval.")}
if(penalty=="none" && (g+k+1)>length(unique(x))){warning("Parameters not estimable. Reduce g+k or increase number ob observation points.")}

if(penalty=="self" && (is.null(Om) || is.null(beta0))){warning("Om and beta0 have to be specified, if penalty='self'.")}
if(penalty!="self" && (!is.null(Om) || !is.null(beta0))){warning("Om and beta0 will be ignored, if penalty!='self'.")}
if(penalty!="self" && is.logical(coinc)){warning("coinc has no influence, if penalty!='self'.")}
if(penalty=="self" && is.null(coinc)){warning("coinc has to be TRUE or FALSE, if penalty='self'.")}

if(constr=="none" && !tuning){warning("tuning has no influence, if constr='none'.")}
if(penalty=="none" && !tuning){warning("tuning has no influence, if penalty='none'.")}

if(!is.null(w)){
if(length(w)!=n){stop("w should be a vector of length n.")
}else{w <- n*w/sum(w)}
}else{w <- rep(1,n)}

orderx <- order(x)
y <- y[orderx]
w <- w[orderx]
x <- sort(x)
dose <- unique(x)
nod <- length(dose)
nd <- as.numeric(table(x))
d <- g+k+1

constraint <- constr
inverse <- 1
if(constr=="invuni"){
constr <- "unimodal"
inverse <- -1
}
if(constr=="antitonic"){
constr <- "isotonic"
inverse <- -1
}

if(is.null(sigmasq)){
wssd <- numeric()
for(u in seq_along(dose)){
yd <- y[x==dose[u]]
wd <- w[x==dose[u]]
wmd <- sum(wd*yd)/sum(wd)
wssd[u] <- sum(wd*(yd-wmd)^2)
}
varest <- sum(wssd)/sum(w)
}else{varest <- sigmasq}

yold <- y
scalfactor <- 0.5*(max(yold) - min(yold))
shiftfactor <- min(yold) + scalfactor
y <- (yold-shiftfactor)/scalfactor

if(length(varest)!=1){
varest <- varest[orderx]
if(!all(varest==varest[1])){
w <- 1/varest
sigmaest <- 1/(scalfactor)^2
}else{sigmaest <- varest[1]/(scalfactor)^2}
}else{sigmaest <- varest/(scalfactor)^2}

if(penalty!="self"){
if(penalty=="none"){coinc <- TRUE; ordpen=2; Dtilde <- diag(0,d); beta0 <- rep(0,d)
}else if(penalty=="diff"){coinc <- FALSE; ordpen=ordpen; Dtilde <- diag(1/vari,d); beta0 <- rep(0,d)
}else if(penalty=="sigEmax"){coinc <- TRUE; ordpen=1; Dtilde <- diag(1/vari,d)
}else if(penalty=="diag"){coinc <- TRUE; ordpen=0; Dtilde <- diag(0,d); beta0 <- rep(0,d)
}
# penalty matrix
Dq <- diag(d)
if(ordpen>=1){Dq <- diff(Dq,difference=ordpen)
}else{Dtilde <- diag(0,d)} # wie bei "diag"
Om <- t(Dq)%*%Dq

# determination of knots depending on the interval [a,b], number g of inner knots and degree k of the spline
# if coinc=T there are k coincident knots at the boundaries a and b
knotseq <- equiknots(a,b,g,k,coinc)

# calculating beta0 for sigEmax
# we fit a sigmoidEmax model and evaluate it at the knot averages
if(penalty=="sigEmax"){
bounds <- defBnds(mD = 8)\$sigEmax
bounds[2,] <- c(1,10)
sigE <- fitMod(x, y, model="sigEmax", bnds=bounds)
knotloc <- knotave(knotseq,d=k)
beta0 <- predict(sigE, newdata=data.frame(x = knotloc), predType="full-model")
}
}else{Dtilde <- diag(0,d)
beta0 <- (beta0-shiftfactor)/scalfactor
knotseq <- equiknots(a,b,g,k,coinc)}

# B-spline design matrix
B <- splineDesign(knots=knotseq, x=x, ord = k+1, outer.ok = T)

tbetaObeta <- t(beta0)%*%Om%*%beta0
tbetaO <- t(beta0)%*%Om
rangO <- qr(Om)\$rank
tBB <- t(B)%*%diag(w)%*%B
tyB <- t(y)%*%diag(w)%*%B

variter <- 0
repeat{
tBSB <- tBB/sigmaest
tySB <- tyB/sigmaest

if(constr=="unimodal"){
vals <- numeric(length(m))
if(nCores>1){
cl <- makeCluster(rep("localhost", nCores), type = "SOCK")
regs <- parLapply(cl=cl, X=m, fun=unisplinem, tBSB=tBSB,tySB=tySB,sigmaest=sigmaest,tbetaO=tbetaO,tbetaObeta=tbetaObeta,Dtilde=Dtilde,rangO=rangO,B=B,beta0=beta0,Om=Om,constr=constr,inverse=inverse,penalty=penalty,tuning=tuning)
stopCluster(cl)
}else{
regs <- lapply(X=m,FUN=unisplinem,tBSB=tBSB,tySB=tySB,sigmaest=sigmaest,tbetaO=tbetaO,tbetaObeta=tbetaObeta,Dtilde=Dtilde,rangO=rangO,B=B,beta0=beta0,Om=Om,constr=constr,inverse=inverse,penalty=penalty,tuning=tuning)
}
for (i in seq_along(m)){vals[i] <- sum(w*(y - regs[[i]]\$fitted.values)^2)}
mini <- which.min(vals)
coef.unimod <- regs[[mini]]\$coef
fit.unimod <- regs[[mini]]\$fitted.values
lambdaopt <- regs[[mini]]\$lambdaopt
}else{ # constr = isotonic or none
erg <- unisplinem(m=d,tBSB=tBSB,tySB=tySB,sigmaest=sigmaest,tbetaO=tbetaO,tbetaObeta=tbetaObeta,Dtilde=Dtilde,rangO=rangO,B=B,beta0=beta0,Om=Om,constr=constr,inverse=inverse,penalty=penalty,tuning=tuning)
coef.unimod <- erg\$coef
fit.unimod <- erg\$fitted.values
lambdaopt <- erg\$lambdaopt
}
Hperm <- tBSB%*%solve(tBSB+ lambdaopt*Om)
ed <- sum(diag(Hperm))

if(!is.null(sigmasq) && is.null(abstol)){break}

resd <- y-fit.unimod
varest <- sum(w*(resd-sum(resd*w)/sum(w))^2)/sum(w)
variter <- variter+1
if(abs(varest - sigmaest)<abstol || variter>=10){
sigmaest <- varest
break}
sigmaest <- varest
}
if(allfits){
allcoefs <- matrix(0,nrow=0,ncol=d)
for(i in seq_along(m)){
allcoefs <- rbind(allcoefs, scalfactor*regs[[i]]\$coef + shiftfactor)
}
}else{
allcoefs <- NULL
}

res <- list(x=x,y=yold,w=w,a=a,b=b,g=g,degree=k,knotsequence=knotseq,constr=constraint,penalty=penalty,Om=Om,beta0=beta0,coinc=coinc,tuning=tuning,abstol=abstol,vari=vari,ordpen=ordpen,
coef=scalfactor*coef.unimod+shiftfactor,fitted.values=scalfactor*fit.unimod+shiftfactor,lambdaopt=lambdaopt,sigmasq=sigmaest*scalfactor^2,variter=variter,ed=ed,modes=m,allcoefs=allcoefs)
class(res) <- "unireg"
res
}
```

## Try the uniReg package in your browser

Any scripts or data that you put into this service are public.

uniReg documentation built on May 2, 2019, 6:50 a.m.