# R/poslasso.R In gpDDE: General Profiling Method for Delay Differential Equation

```###Positive-lasso###
###Based on lars.pos provided by David Reiss, further modified by Kun Chen###
###Use positive=T, type="lasso"###
###Not fully modified for positive-lar and positive-stagewise.

#Example:
#source("poslasso.R")
#X <- matrix(ncol=20,nrow=100,rnorm(2000))
#beta <- as.vector(c(rep(1,5),rep(-1,5),rep(0,10)))
#E <- as.vector(rnorm(100,0,1))
#Y <- X%*%beta + E
#fit <- lars.pos(X,Y,positive=TRUE,type="lasso",intercept=FALSE,normalize=FALSE)
#coef(fit)
#plot(fit\$Cp)
#plot(fit)
###The last solution is the same as the LS-fit of the selected predictors
#lastsolution <- coef(fit)[nrow(coef(fit)),]
#as.vector(coef(lm(Y~X[,lastsolution!=0]-1)))

lars.pos <-
function (x, y, type = c("lasso", "lar", "forward.stagewise"),
trace = FALSE, Gram, eps = .Machine\$double.eps, max.steps,
use.Gram = TRUE, positive = FALSE, normalize=TRUE,intercept=FALSE)
{
call <- match.call()
type <- match.arg(type)
TYPE <- switch(type, lasso = "LASSO", lar = "LAR", forward.stagewise = "Forward Stagewise")
if (trace)
cat(paste(TYPE, "sequence\n"))
nm <- dim(x)
n <- nm
m <- nm
im <- inactive <- seq(m)
one <- rep(1, n)
vn <- dimnames(x)[]

if(intercept){
meanx <- drop(one %*% x)/n
x <- scale(x, meanx, FALSE)	# centers x
mu <- mean(y)
y <- drop(y - mu)
}else {
meanx <- rep(0,m)
mu <- 0
y <- drop(y)
}

#meanx <- (one %*% x)[1, ]/n
#x <- scale(x, meanx, FALSE)
#mu <- mean(y)
#y <- my.drop(y - mu)

if(normalize){
normx <- sqrt(drop(one %*% (x^2)))
nosignal<-normx/sqrt(n) < eps
if(any(nosignal))# ignore variables with too small a variance
{
ignores<-im[nosignal]
inactive<-im[-ignores]
normx[nosignal]<-eps*sqrt(n)
if(trace)
cat("LARS Step 0 :\t", sum(nosignal), "Variables with Variance < eps; dropped for good\n")
#
}else ignores <- NULL #singularities; augmented later as well
names(normx) <- NULL
x <- scale(x, FALSE, normx)	# scales x
}else {
normx <- rep(1,m)
ignores <- NULL
}

if (use.Gram & missing(Gram)) {
if (m > 500 && n < m)
cat("There are more than 500 variables and n<m;\nYou may wish to restart and set
use.Gram=FALSE\n")
if (trace)
cat("Computing X'X .....\n")
Gram <- t(x) %*% x
}

Cvec <- my.drop(t(y) %*% x)
ssy <- sum(y^2)
residuals <- y
if (missing(max.steps)) max.steps <- 8 * min(m, n - intercept)
beta <- matrix(0, max.steps + 1, m)
Gamrat <- NULL
arc.length <- NULL
R2 <- 1
first.in <- integer(m)
active <- NULL
actions <- as.list(seq(max.steps))
drops <- FALSE
Sign <- NULL
R <- NULL

k <- 0
pos_stop <- FALSE
while((k < max.steps) & (length(active) < min(m - length(ignores),n-intercept)) &
pos_stop==FALSE) {
action <- NULL
k <- k + 1

C <- Cvec[inactive]

if (!positive){
Cmax <- max(abs(C))
}else
{
Cmax <- max(C)
##if(Cmax < -eps) pos_stop <- TRUE
}

if (!any(drops)) {
if (!positive)
new <- abs(C) >= Cmax - eps
else new <- C >= Cmax - eps
C <- C[!new]
new <- inactive[new]
for (inew in new) {
if (use.Gram) {
R <- updateR.faster(Gram[inew, inew], R, my.drop(Gram[inew,
active]), Gram = TRUE, eps = eps)
}else {
R <- updateR.faster(x[, inew], R, x[, active],
Gram = FALSE, eps = eps)
}
if (attr(R, "rank") == length(active)) {
nR <- seq(length(active))
R <- R[nR, nR, drop = FALSE]
attr(R, "rank") <- length(active)
ignores <- c(ignores, inew)
action <- c(action, -inew)
if (trace)
cat("LARS Step", k, ":\t Variable", inew,
"\tcollinear; dropped for good\n")
}else {
if (first.in[inew] == 0)
first.in[inew] <- k
active <- c(active, inew)
if (!positive)
Sign <- c(Sign, sign(Cvec[inew]))
else Sign <- c(Sign, abs(sign(Cvec[inew])))
action <- c(action, inew)
if (trace)
cat("LARS Step", k, ":\t Variable", inew,
}
}
}
else action <- -dropid
Gi1 <- backsolve(R, lars::backsolvet(R, Sign))
dropouts <- NULL
if (type == "forward.stagewise") {
directions <- Gi1 * Sign
if (!all(directions > 0)) {
if (use.Gram) {
nnls.object <- nnls.lars.faster(active, Sign,
R, directions, Gram[active, active], trace = trace,
use.Gram = TRUE, eps = eps)
}
else {
nnls.object <- nnls.lars.faster(active, Sign,
R, directions, x[, active], trace = trace,
use.Gram = FALSE, eps = eps)
}
positive <- nnls.object\$positive
dropouts <- active[-positive]
action <- c(action, -dropouts)
active <- nnls.object\$active
Sign <- Sign[positive]
Gi1 <- nnls.object\$beta[positive] * Sign
R <- nnls.object\$R
C <- Cvec[-c(active, ignores)]
}
}
A <- 1/sqrt(sum(Gi1 * Sign))
w <- A * Gi1
if (!use.Gram)
u <- (x[, active, drop = FALSE] %*% w)[, 1]
if (length(active) >= min(n - intercept, m - length(ignores))) {
gamhat <- Cmax/A
}
else {
if (use.Gram) {
a <- my.drop(w %*% Gram[active, -c(active, ignores),
drop = FALSE])
}
else {
a <- (u %*% x[, -c(active, ignores), drop = FALSE])[1,
]
}
if (!positive)
gam <- c((Cmax - C)/(A - a), (Cmax + C)/(A +
a))
else gam <- (Cmax - C)/(A - a)
gamhat <- min(gam[gam > eps], Cmax/A)
}
if (type == "lasso") {
dropid <- NULL
b1 <- beta[k, active]
z1 <- -b1/w
zmin <- min(z1[z1 > eps], gamhat)
if (zmin < gamhat) {
gamhat <- zmin
drops <- z1 == zmin
}
else drops <- FALSE
}
beta[k + 1, ] <- beta[k, ]
beta[k + 1, active] <- beta[k + 1, active] + gamhat * w
if (use.Gram) {
Cvec <- Cvec - gamhat * Gram[, active, drop = FALSE] %*% w
}
else {
residuals <- residuals - gamhat * u
Cvec <- (t(residuals) %*% x)[1, ]
}
Gamrat <- c(Gamrat, gamhat/(Cmax/A))
arc.length <- c(arc.length, gamhat)
if (type == "lasso" && any(drops)) {
dropid <- seq(drops)[drops]
for (id in rev(dropid)) {
if (trace)
cat("Lasso Step", k + 1, ":\t Variable", active[id],
"\tdropped\n")
R <- lars::downdateR(R, id)
}
dropid <- active[drops]
beta[k + 1, dropid] <- 0
active <- active[!drops]
Sign <- Sign[!drops]
}
if (!is.null(vn))
names(action) <- vn[abs(action)]
actions[[k]] <- action
inactive <- im[-c(active, ignores)]

###Note that for positive-lasso, some predictors may never enter the model
###We shall stop when all the predictors left have negative correlations with y.
###The last step corresponds to LS fit of the selected predictors.
if (positive)
{
C <- Cvec[inactive]
Cmax <- max(C)
if(Cmax < -eps) pos_stop <- TRUE
}

}
beta <- beta[seq(k + 1), ]
dimnames(beta) <- list(paste(0:k), vn)
if (trace)
residuals <- y - x %*% t(beta)
beta <- scale(beta, FALSE, normx)

###Cp selection
actions=actions[seq(k)]
netdf=sapply(actions,function(x)sum(sign(x)))
df=cumsum(netdf)### This takes into account drops
if(intercept)df=c(Intercept=1,df+1)
else df=c(Null=0,df)
df.big=n-rev(df)
else
Cp <- RSS/sigma2 - n + 2 * df
attr(Cp,"sigma2")=sigma2
attr(Cp,"n")=n
##Cp <- ((n - k - 1) * RSS)/rev(RSS) - n + 2 * seq(k + 1)

object <- list(call = call, type = TYPE, R2 = R2, RSS = RSS,
Cp = Cp, actions = actions[seq(k)], entry = first.in,
Gamrat = Gamrat, arc.length = arc.length, Gram = if (use.Gram) Gram else NULL,
beta = beta, mu = mu, normx = normx, meanx = meanx)
class(object) <- c("lars", "lars.pos")
object
}

my.drop <-
function (x)
{
q = dim(x)
if (length(q) == 0)
x
else if (q == 1)
x[1, ]
else if (q == 1)
x[, 1]
else x
}

nnls.lars.faster <-
function (active, Sign, R, beta, Gram, eps = 1e-10, trace = FALSE,
use.Gram = TRUE)
{
if (!use.Gram)
x <- Gram
M <- m <- length(active)
im <- seq(m)
positive <- im
zero <- NULL
while (m > 1) {
zero.old <- c(m, zero)
R.old <- lars::downdateR(R, m)
beta0 <- backsolve(R.old, lars::backsolvet(R.old, Sign[-zero.old])) *
Sign[-zero.old]
beta.old <- c(beta0, rep(0, length(zero.old)))
if (all(beta0 > 0))
break
m <- m - 1
zero <- zero.old
positive <- im[-zero]
R <- R.old
beta <- beta.old
}
while (TRUE) {
while (!all(beta[positive] > 0)) {
alpha0 <- beta.old/(beta.old - beta)
alpha <- min(alpha0[positive][(beta <= 0)[positive]])
beta.old <- beta.old + alpha * (beta - beta.old)
dropouts <- match(alpha, alpha0[positive], 0)
for (i in rev(dropouts)) R <- lars::downdateR(R, i)
positive <- positive[-dropouts]
zero <- im[-positive]
beta0 <- backsolve(R, lars::backsolvet(R, Sign[positive])) *
Sign[positive]
beta <- beta.old * 0
beta[positive] <- beta0
}
if (use.Gram) {
w <- 1 - Sign * my.drop(Gram %*% (Sign * beta))
}
else {
jw <- x %*% (Sign * beta)
w <- 1 - Sign * my.drop(t(jw) %*% x)
}
if ((length(zero) == 0) || all(w[zero] <= 0))
break
if (use.Gram) {
positive]), Gram = TRUE, eps = eps)
}
else {
R <- updateR.faster(x[, add], R, x[, positive], Gram = FALSE,
eps = eps)
}
beta0 <- backsolve(R, lars::backsolvet(R, Sign[positive])) *
Sign[positive]
beta[positive] <- beta0
}
if (trace) {
dropouts <- active[-positive]
for (i in dropouts) {
cat("NNLS Step:\t Variable", i, "\tdropped\n")
}
}
list(active = active[positive], R = R, beta = beta, positive = positive)
}

updateR.faster <-
function (xnew, R = NULL, xold, eps = .Machine\$double.eps, Gram = FALSE)
{
xtx <- if (Gram)
xnew
else sum(xnew^2)
norm.xnew <- sqrt(xtx)
if (is.null(R)) {
R <- matrix(norm.xnew, 1, 1)
attr(R, "rank") <- 1
return(R)
}
Xtx <- if (Gram)
xold
else (t(xnew) %*% xold)[1, ]
r <- lars::backsolvet(R, Xtx)
rpp <- norm.xnew^2 - sum(r^2)
rank <- attr(R, "rank")
if (rpp <= eps)
rpp <- eps
else {
rpp <- sqrt(rpp)
rank <- rank + 1
}
R <- cbind(rbind(R, 0), c(r, rpp))
attr(R, "rank") <- rank
R
}

#
# lasso.adapt.pos <-function(x,y, pars, beta, pars.c = 1){
# # adaptive lasso from lars with BIC stopping rule
# # this one uses the "known variance" version of BIC with RSS/(full model mse)
# # must use a recent version of R so that normalize=FALSE can be used in lars
# #  standardize variables like lars does
#     w.pars <- abs(pars) * pars.c
#     w.beta <- abs(beta)                      # weights for adaptive lasso
#     beta.ind <- which(w.beta > 0)
#     w.beta <- w.beta[beta.ind]
#     x <- x[, c(1:length(pars), beta.ind + length(par))]
#
#     one <- rep(1, n)
#     meanx <- drop(one %*% x)/n
#     xc <- scale(x, meanx, FALSE)         # first subtracts mean
#     normx <- sqrt(drop(one %*% (xc^2)))
#     names(normx) <- NULL
#     xs <- scale(xc, FALSE, normx)        # now rescales with norm (not sd)
#
# xs=scale(xs,center=FALSE,scale=1/w)  # xs times the weights
# object=lars(xs,y,type="lasso",normalize=FALSE)
#
# # get min BIC
# sig2f=summary(out.ls)\$sigma^2        # full model mse
# step.bic2=which.min(bic2)            # step with min BIC
#
# fit=predict.lars(object,xs,s=step.bic2,type="fit",mode="step")\$fit
# coeff=predict.lars(object,xs,s=step.bic2,type="coef",mode="step")\$coefficients
# coeff=coeff*w/normx                  # get back in right scale
# st=sum(coeff !=0)                    # number nonzero
# mse=sum((y-fit)^2)/(n-st-1)          # 1 for the intercept
#
# # this next line just finds the variable id of coeff. not equal 0
# if(st>0) x.ind<-as.vector(which(coeff !=0)) else x.ind<-0
# intercept=as.numeric(mean(y)-meanx%*%coeff)
# return(list(fit=fit,st=st,mse=mse,x.ind=x.ind,coeff=coeff,intercept=intercept,object=object,
#             bic2=bic2,step.bic2=step.bic2))
# }
```

## Try the gpDDE package in your browser

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

gpDDE documentation built on May 2, 2019, 1:09 p.m.