Nothing
R/TSMCP.R
In MTAFT: Data-Driven Estimation for Multi-Threshold Accelerate Failure Time Model
Defines functions
TSMCP
cumcp
lmres
plus.knots
plus.recip
plus.penalty
plus.not.in.loop
plus.new.eta
plus.check.eta
plus.single.step
plot.plus
coef.plus
predict.plus
print.plus
plus
Documented in
TSMCP
####This R script mainly implements the two-stage method of Li and Jin (2018)
########################################################
# plus(), the main plus function; documented
###
plus <- function(x,y, method = c("lasso", "mc+", "scad", "general"), m=2, gamma,v,t,
monitor=FALSE, normalize = TRUE, intercept = TRUE,
Gram, use.Gram = FALSE, eps=1e-15, max.steps=500, lam)
{
## fill missing values
if (missing(gamma)) gamma <- 0
if (missing(v)) v <- 0
if (missing(t)) t <- 0
if (missing(Gram)) Gram <- 0
if (missing(lam)) lam <- -1
lam.min <- min(lam)
## check data
data.dim.error <- FALSE
if (! is.matrix(x)) data.dim.error <- TRUE
else {
n <- dim(x)[1]
if ( (n<=1) | (n!=length(y)) ) data.dim.error <- TRUE
p <- dim(x)[2]
if (p<= 1) data.dim.error <- TRUE
}
if (data.dim.error == TRUE){
warning("Data dimention error.", "\n")
return(list(data.dim.error = TRUE))
}
## standardize x
orig.x <- x
if (intercept) x <- t(t(x)-apply(x,2,mean)) # center to mean zero
if (normalize) {
normalize.factor <- sqrt(apply(x^2,2,mean))
x <- t(t(x)/normalize.factor)
}
## compute Gram matrix if necessary
## if ((!use.Gram) & (p<=200)) use.Gram <- TRUE
if ((use.Gram) & (p > 1000)) use.Gram <- FALSE
if (use.Gram)
if (! is.matrix(Gram)) Gram <- t(x)%*%x
else if ((dim(Gram)[1]!=p) || (dim(Gram)[2]!=p)) Gram <- t(x)%*%x
## determine method and compute penalty
if (length(method)==1)
if (method == "lasso") m <- 1
else if (method == "mc+") m <- 2
else if (method == "scad") m <- 3
if (m==1) method <- "LASSO"
else if (m==2) method <- "MC+"
else if (m==3) method <- "SCAD"
else method <- "PLUS"
if (monitor == TRUE) {
cat("Sequence of", method, "moves:", "\n","\n")
}
if (m<4) { # provide default penalty for lasso, MCP and SCAD
if ((m>1) & (use.Gram) & (gamma == 0)) {
max.corr <- max(abs(Gram[ row(Gram) != col(Gram) ]))/n
if ((1 - max.corr) >= (1e-5) ) gamma <- 2/(1- max.corr)
}
if ((m==3) & (gamma <=1)) gamma <- 0
tmp <- plus.penalty(m,a=gamma)
m <- tmp$m
gamma <- tmp$gamma
v <- tmp$v
t <- tmp$t
}
## declare variables
names <- dimnames(x)
b <- matrix(0, max.steps+1,p)
etatil <- matrix(0,2,p+1)
tau <- rep(0, max.steps+1)
my.env <-environment()
my.env$x <-x
my.env$y <- y
my.env$tx <- t(x)/n
my.env$eta <- b
my.env$g.prime <- rep(0,p)
my.env$grad <- rep(0,p)
my.env$etatil2 <- etatil
my.env$etatil2[1,] <- 1:(p+1)
if (use.Gram) my.env$Sigma <- Gram/n
else {
my.env$Sigma <- matrix(0,p,100) # e$a.set.var will match a.set
my.env$var.list <- 0
}
my.env$use.Gram <- use.Gram
my.env$ties <- rep(FALSE,p)
my.env$t.fun <- plus.knots(t,m)
z <- my.env$tx %*%y
my.env$tau1 <- 1/max(abs(z))
## initialization: for the initial step k=0, tau^{(0)} is needed, b[1,] <- 0 done already
tau[1] <- my.env$tau1
etatil[1,] <- 1:(p+1) # the first etatil
etatil[2,1:p] <- sign(z)
etatil <- etatil[ , c(abs(z) >= max(abs(z)) - eps,TRUE) ]
k <-1
exit.while <- FALSE
last.valid.s <- rep(0,p)
## the main loop
while (exit.while==FALSE) { # segment k: model eta[k+1,] begins with b[k,] and ends with b[k+1,]
tmp <- plus.single.step(k, etatil, z, b[k,], tau[k], m, v, t, eps=eps,e=my.env)
## tmp holds return(new.eta, etatil, s, tau, b, singular.Q, forced.stop, full.path)
etatil <- tmp$etatil # etatil = 0 if forced.stop, length(etatil)=2 if Ise.termination
exit.while <- tmp$forced.stop # cannot find a valid new.eta
full.path <- tmp$full.path
if (exit.while == FALSE) { # save data from tmp if not forced to stop
eta[k+1,] <- tmp$new.eta
if (monitor == TRUE) { # print info to monitor the iterations
var.change <- sign(abs(sign(eta[k+1,]))-abs(sign(eta[k,])))
if (sum(abs(var.change))==1) {
if(sum(var.change)==1)
action <-"added "
else
action <-"removed "
var_name <- (order(-abs(var.change))[1])
if (!is.null(names))
var_name <- names[[2]][var_name]
cat("Step",k,": ",action,var_name,"\n")
}
}
last.valid.s <- tmp$s # save the s for the last valid eta
tau[k+1]<- tmp$tau
b[k+1,] <- tmp$b # b[k+1] is at the end of segment eta[k+1,]
b[k+1,eta[k+1,]==0] <- 0 # more accurate b=0, b[k,eta[k+1,]==0] <- 0 is also ok
k <- k+1
if ((k > max.steps)|(full.path)|(tau[k]*lam.min > 1)) { # stop with a valid eta
exit.while <- TRUE
if ((full.path) & (m>1)) k <- k-1
}
} # if (exit.while == FALSE)
} # end the main loop
## tmp holds return(new.eta, etatil, s, tau, b, singular.Q, forced.stop, full.path)
singular.Q <- tmp$singular.Q
forced.stop <- tmp$forced.stop
lam.path <- 1/tau[1:k]
beta.path <- lam.path * b[1:k,]
if (normalize) {
beta.path <- t(t(beta.path)/normalize.factor)
}
lam.path[k] <- max(abs(my.env$tx %*%(y-x%*%beta.path[k,])))
if (full.path) lam.path[k] <- 0
eta <- eta[1:k,]
total.steps <- k-1
## calculate the estimator beta at the given lam or the set of values of lam.path
x <- orig.x
if (sum(lam) == -1) lam <- sort(lam.path, decreasing=TRUE)
tmp <- plus.hit.points(lam,lam.path)
total.hits <- tmp$total.hits
if (total.hits > 0) {
beta <- tmp$w1*beta.path[tmp$k1,] + (1-tmp$w1)*beta.path[tmp$k2,]
lam <- lam[1:total.hits]
dim <- apply(abs(beta)>0, 1, sum)
if (total.hits == 1)
r.square <- 1- sum((y - x%*%beta)^2)/sum(y^2)
else
r.square <- 1- apply((x%*%t(beta) - y)^2, 2, sum)/sum(y^2)
obj<-list(x=x,y=y,eta=eta, beta.path=beta.path, lam.path=lam.path,
beta=beta, lam=lam, dim=dim, r.square = r.square, total.hits=total.hits,
method = method, gamma = gamma, total.steps=total.steps, max.steps=max.steps,
full.path=full.path, forced.stop=forced.stop, singular.Q=singular.Q)
}
else
obj<-list(x=x,y=y,eta=eta, beta.path=beta.path, lam.path=lam.path, total.hits=total.hits,
method = method, gamma = gamma, total.steps=total.steps, max.steps=max.steps,
full.path=full.path, forced.stop=forced.stop, singular.Q=singular.Q)
class(obj) <- "plus"
return(obj)
}
########################################################
# print.plus(), to print plus() moves; documented
###
print.plus <- function(x, print.moves = 20, ...) {
#cat("\n")
#cat("Sequence of", x$method, "moves:", "\n","\n")
print.moves <- round(print.moves)
if (print.moves < 1) print.moves <- 20
names <- dimnames(x$x)
m <- 1
k <- 1
exit.while <- FALSE
while (exit.while == FALSE) {
var.change <- sign(abs(sign(x$eta[k+1,]))-abs(sign(x$eta[k,])))
if (sum(abs(var.change))==1) {
if(sum(var.change)==1)
action <-"added "
else
action <-"removed "
var_name <- (order(-abs(var.change))[1])
if (!is.null(names))
var_name <- names[[2]][var_name]
#cat("Step",k,": ",action,var_name,"\n")
m <- m+1
if (m>print.moves) exit.while <- TRUE
}
k <- k+1
if (k > x$total.steps) exit.while <- TRUE
} # end while
}
########################################################
# predict.plus(), to extract coefficients and preditions; documented
###
predict.plus <- function(object,lam,newx, ...) {
if (missing(newx)) {
flag <- FALSE
warning("Warning message: ")
warning("no newx argument; newy not produced", "\n")
}
else flag <- TRUE
if (missing(lam)) {
lam <- sort(object$lam.path, decreasing=TRUE)
warning("Warning message: ")
warning("no lam argument; object$lam.path used", "\n")
}
else if (max(lam) < min(object$lam.path)) {
lam <- sort(object$lam.path, decreasing=TRUE)
warning("Warning message: ")
warning("lam not reached by the plus path; object$lam.path used", "\n")
}
#cat("\n")
tmp <- plus.hit.points(lam,object$lam.path)
total.hits <- tmp$total.hits
beta <- tmp$w1* object$beta.path[tmp$k1,] + (1-tmp$w1)* object$beta.path[tmp$k2,]
lam <- lam[1:total.hits]
dim <- apply(abs(beta)>0, 1, sum)
if (total.hits == 1) {
r.square <- 1- sum((object$y - object$x%*%beta)^2)/sum(object$y^2)
if (flag) {
if (is.matrix(newx)) newy <- newx%*% beta
else newy <- sum(beta * newx)
}
}
else {
r.square <- 1- apply((object$x%*%t(beta) - object$y)^2, 2, sum)/sum(object$y^2)
if (flag) {
if (is.matrix(newx)) newy <- beta %*% t(newx)
else newy <- beta %*% newx
}
}
if (flag)
return(list(lambda=lam, coefficients=beta, dimension=dim, r.square = r.square, step = tmp$k2-1,
method = object$method, newy = newy))
else
return(list(lambda=lam, coefficients=beta, dimension=dim, r.square = r.square, step = tmp$k2-1,
method = object$method))
}
########################################################
# coef.plus(), to extract coefficients; documented
###
coef.plus <- function(object,lam, ...){
if (missing(lam)) {
lam <- sort(object$lam.path, decreasing=TRUE)
warning("Warning message: ")
warning("no lam argument; object$lam.path used", "\n")
}
else if (max(lam) < min(object$lam.path)) {
lam <- sort(object$lam.path, decreasing=TRUE)
warning("Warning message: ")
warning("lam not reached by the plus path; x$lam.path used", "\n")
}
#cat("\n")
tmp <- plus.hit.points(lam,object$lam.path)
beta <- tmp$w1* object$beta.path[tmp$k1,] + (1-tmp$w1)* object$beta.path[tmp$k2,]
return(beta)
}
########################################################
# plot.plus(), to plot coefficients, preditions, penalty level, dimension or R-square
# againt plus step or penalty level; documented
###
plot.plus <- function(x, xvar=c("lam","step"), yvar=c("coef","newy","lam","dim","R-sq"),
newx, step.interval, lam.interval, predictors, ...) {
# determin xvar in the plot
if (length(unique(c(xvar,"lam")))==length(unique(c(xvar)))) xvar <- "lambda"
else if (length(unique(c(xvar,"step")))==length(unique(c(xvar)))) xvar <- "step"
else xvar <- "lambda"
# determin yvar in the plot
if (length(unique(c(yvar,"coef")))==length(unique(c(yvar)))) yvar <- "coefficients"
else if (length(unique(c(yvar,"newy")))==length(unique(c(yvar)))) yvar <- "newy"
else if (length(unique(c(yvar,"lam")))==length(unique(c(yvar)))) yvar <- "lambda"
else if (length(unique(c(yvar,"dim")))==length(unique(c(yvar)))) yvar <- "dimension"
else if (length(unique(c(yvar,"R-sq")))==length(unique(c(yvar)))) yvar <- "R-square"
else yvar <- "coefficients"
if (yvar != "newy") newx <- matrix(0,2,dim(x$x)[2])
if ((yvar == "newy") & missing(newx)) {
warning("Warning message: ")
warning("no newx argument", "\n")
return()
}
if (missing(predictors)) predictors <- 1:(dim(x$x)[2])
## plot
if (xvar == "step") {
xtmp <- 0:(x$total.steps)
plot.set <- rep(TRUE, length(xtmp))
if (! missing(step.interval)) {
plot.set[xtmp > max(step.interval)] <- FALSE
plot.set[xtmp < min(step.interval)-1] <- FALSE
xtmp <- xtmp[plot.set]
}
if (yvar == "coefficients") ytmp <- x$beta.path[plot.set, ]
if ((yvar == "newy") & is.matrix(newx)) ytmp <- x$beta.path[plot.set, ]%*%t(newx)
if ((yvar == "newy") & (! is.matrix(newx))) ytmp <- x$beta.path[plot.set, ]%*%newx
if (yvar == "lambda") ytmp <- x$lam.path[plot.set]
if (yvar == "dimension") ytmp <- apply(abs(sign(x$eta[plot.set,])),1, sum)
if (yvar == "R-square")
ytmp <- 1- apply((x$x%*%t(x$beta.path[plot.set, ]) - x$y)^2, 2, sum)/sum(x$y^2)
}
if (xvar == "lambda") {
xtmp <- sort(x$lam.path,decreasing=TRUE)
tmp <- predict(x,xtmp,newx)
plot.set <- rep(TRUE, length(xtmp))
if (! missing(lam.interval)) {
plot.set[xtmp > max(lam.interval)] <- FALSE
plot.set[xtmp < min(lam.interval)] <- FALSE
xtmp <- - xtmp[plot.set]
}
else xtmp <- -xtmp
if (yvar == "coefficients") ytmp <- tmp$coefficients[plot.set, ]
if ((yvar == "newy") & is.matrix(newx)) ytmp <- tmp$newy[plot.set,]
if ((yvar == "newy") & (! is.matrix(newx))) ytmp <- tmp$newy[plot.set]
if (yvar == "lambda") ytmp <- tmp$lambda[plot.set]
if (yvar == "dimension") ytmp <- tmp$dimension[plot.set]
if (yvar == "R-square") ytmp <- tmp$r.square[plot.set]
}
xmax <- max(xtmp)
xmin <- min(xtmp)
ymax <- max(ytmp)
ymin <- min(ytmp)
if (xvar == "lambda") {
plot(c(xmin,xmax),c(ymin,ymax), xaxt = "n", xlab = xvar, ylab=yvar, type="l",lty=0, main=x$method)
axis(1, at = axTicks(1), labels = abs(axTicks(1)) )
}
else plot(c(xmin,xmax),c(ymin,ymax), xlab = xvar, ylab=yvar, type="l",lty=0, main=x$method)
if ( ((yvar=="coefficients")|(yvar=="newy")) & is.matrix(ytmp)) {
index <- 1: (dim(ytmp)[2])
if (yvar=="coefficients") {
f.set <- rep(FALSE, max(length(index),max(predictors)))
f.set[predictors] <- TRUE
f.set <- f.set[1:length(index)]
}
else f.set <- rep(TRUE, length(index))
axis(4, at = ytmp[dim(ytmp)[1],f.set], labels = index[f.set])
for (j in 1: (dim(ytmp)[2]))
if (f.set[j]) lines(xtmp, ytmp[,j])
}
else lines(xtmp, ytmp)
}
########################################################
# internal functions used by plus()
###
########################################################
# one single step in the iterative algorithm
###
plus.single.step <- function(k, etatil, z, b, tau, m, v, t, eps=1e-15, e) {
## input
# k: number of parallelepiped indicators already found, including 0, k-1 steps completed
# etatil, 2 X (|C|+1): label and new indicator for max crossing with dummy variable p+1
# NOTE: |C| > 0 required here
# z: t(x)y/n; b, tau: old turning point b at old tau; m, v, t: penalty; eps: numerical zero,
## out put: return(new.eta, etatil, s, tau,b, singular.Q, forced.stop, full.path)
# including info at the end of segment new.eta, e.g. etatil for the eta beyond new.eta
n <- dim(e$x)[1] # no dummy
p <- dim(e$x)[2]
# this is in step k, eta^{(0)},...,eta^{(k-1)} have already been found
sides <- length(etatil)/2-1 # number of boundaries for possible crossing
singular.Q <- FALSE # this and next values are needed in case plus.check.eta not called
forced.stop <- TRUE
if (sides == 1) { # one-at-a-time scenario
new.eta <- plus.new.eta(e$eta[k,],etatil,rep(TRUE,2),p) # T for crossing the boundary and dummy
## always a new eta in theory, only check when k is a multiplier of 50 to prevent numerical loop
if (k != round(k/50)*50) check.new.eta <- TRUE
else check.new.eta <- plus.not.in.loop(k, new.eta, e=e)
if (check.new.eta==TRUE) { # do not run plus.check.eta if eta is not new, try xi = 1 first
tmp <- plus.check.eta(e$eta[k,], new.eta, etatil, 1, z, v, eps=eps,e=e)
## tmp holds return(xi, singular.Q, s, invalid.cross)
singular.Q <- tmp$singular.Q # save if invalid cross is due to singular Q, an almost nonevent
if ( (m>1) & ((tmp$invalid.cross) == TRUE) ) { # try xi = -1
tmp <- plus.check.eta(e$eta[k,], new.eta, etatil, -1, z, v, eps=eps,e=e)
if (singular.Q == FALSE) singular.Q <- tmp$singular.Q
}
forced.stop <- tmp$invalid.cross # due to singular Q or not
} # if (check.new.eta==T)
} # end of the one-at-a-time case
else { # sides > 1, first check ties
etatil <- etatil[, c((! e$ties[ etatil[1,1:sides] ]),TRUE) ] # remove known ties
sides <- length(etatil)/2-1
if (sides > 1) { # remove new ties
if (e$use.Gram) Sgm.B <- e$Sigma[etatil[1,1:sides], etatil[1,1:sides]]
else Sgm.B <- e$tx[etatil[1,1:sides],] %*% e$x[,etatil[1,1:sides]]
for (i in 1:(sides-1)) for (j in (i+1):sides)
if (Sgm.B[i,j] > 1- eps) {
e$ties[j] <- TRUE
}
etatil <- etatil[, c(! e$ties[ etatil[1,1:sides] ],TRUE) ] # remove known ties
sides <- length(etatil)/2-1
}
if (sides == 1) {# go back to the one-at-time case
new.eta <- plus.new.eta(e$eta[k,],etatil,rep(TRUE,2),p) # T for crossing the boundary and dummy
## always a new eta in theory, only check when k is a multiplier of 50 to prevent numerical loop
if (k != round(k/50)*50) check.new.eta <- TRUE
else check.new.eta <- plus.not.in.loop(k, new.eta, e=e)
if (check.new.eta==TRUE) { # do not run plus.check.eta if eta is not new, try xi = 1 first
tmp <- plus.check.eta(e$eta[k,], new.eta, etatil, 1, z, v, eps=eps,e=e)
## tmp holds return(xi, singular.Q, s, invalid.cross)
singular.Q <- tmp$singular.Q # save if invalid cross is due to singular Q, an almost nonevent
if ( (m>1) & ((tmp$invalid.cross) == TRUE) ) { # try xi = -1
tmp <- plus.check.eta(e$eta[k,], new.eta, etatil, -1, z, v, eps=eps,e=e)
if (singular.Q == FALSE) singular.Q <- tmp$singular.Q
}
forced.stop <- tmp$invalid.cross # due to singular Q or not
} # if (check.new.eta==T)
} # end of the one-at-a-time case
else if (sides > 1) { # general case, (more than) two-at-a-time
M <- matrix(FALSE, 2^sides, sides + 1) # generate all possible candidates for new eta
M[, sides+1] <- TRUE # last T for dummy
M[ (2^(sides-1)+1):(2^sides), 1 ] <- TRUE # T for cross the indicated side
for (i in (2:sides) ) M[,i] <- M[(2^(sides-i)+1):(3*2^(sides-i)),i-1]
M <- M[order(apply(M,1,sum)),] # closer eta (smaller number of crossings) has higher priority
# initial while loop look for new eta
continue.while <- TRUE
j <- 1
while (continue.while == TRUE) {
new.eta <- plus.new.eta(e$eta[k,],etatil,M[j+1,],p) # skip the first M[1,]=F...FT
new.is.new <- plus.not.in.loop(k, new.eta, e=e)
if (new.is.new == TRUE) { # try xi =1
tmp <- plus.check.eta(e$eta[k,], new.eta, etatil, 1, z, v, eps=eps, e=e)
## tmp holds return(xi, singular.Q, s, invalid.cross)
## invalid.cross == T iff either (singular.Q == T) or (xi and s fail to cross to new.eta)
if (singular.Q == FALSE) singular.Q <- tmp$singular.Q
if ( (m>1) & ((tmp$invalid.cross) == TRUE)) { # try xi = -1
tmp <- plus.check.eta(e$eta[k,], new.eta, etatil, -1, z, v, eps=eps, e=e)
if (singular.Q == FALSE) singular.Q <- tmp$singular.Q # singular.Q should be T if happens once
}
forced.stop <- tmp$invalid.cross # forced.stop is true if all "not new" or invalid
continue.while <- tmp$invalid.cross # exit while loop for valid cross to new eta
} # end if (new.is.new == T)
j <- j+1 # continue.while does not change when new.is.new == F
if (j==2^sides) continue.while <- FALSE
} # end while, if forced.stop == F, then tmp holds info for a valid new eta
} # end of (more than) two-at-a-time
else if (sides < 1) new.eta <- rep(m+2,p)
}
full.path <- FALSE
if (forced.stop == TRUE) { # terminate the program, need value for the 9 out of 10 var
## dummy output for return(new.eta, etatil, s, tau,b, singular.Q, forced.stop, full.path)
## values have already been assigned to new.eta, singular.Q, forced.stop, full.path
etatil <- c(p+1,m+2) # etatil <- as.matrix(c(p+1,m+2))
s <- rep(0,p)
tau <- 0
b <- rep(0,p)
}
if (forced.stop==FALSE) { # valid new eta has been found as new.eta
## tmp holds return(xi, singular.Q, s, invalid.cross), singular.Q = invalid.cross = F
xi <- tmp$xi
s <- tmp$s
## need values for return(etatil, tau,b, full.path); new.eta, singular.Q, forced.stop done
a.set <- (new.eta!=0)
if ( k - 50*round(k/50) == 1) # directly calculate the gradient to prevent accumulation of error
# e$grad <- tau*z - as.matrix(e$Sigma[,e$a.set.var])%*%b[a.set]
if (sum(a.set)<2)
e$grad <- tau*z - (e$Sigma[,e$a.set.var])*b[a.set] # gradient
else
e$grad <- tau*z - (e$Sigma[,e$a.set.var])%*%b[a.set]
## compute DeltaJ and new etatil
DeltaJ <- rep(-1,p) # -1 represents infinity
new.etatil <- rep(m+2,p) # m+2 is the dummy value for etatil calculation
ind <- (xi*s>=eps) & (new.eta != 0) & (new.eta < m)
if (sum(ind)>0) { # move towards higher eta value
DeltaJ[ind] <- xi*(e$t.fun[new.eta[ind]+1+m]-b[ind])*plus.recip(s[ind],eps=eps)
new.etatil[ind] <- new.eta[ind] + 1
}
ind <- (xi*s<= -eps) & (new.eta != 0) & (new.eta > -m)
if (sum(ind)>0) { # move towards lower eta value
DeltaJ[ind] <- xi*(e$t.fun[new.eta[ind]+m]-b[ind])*plus.recip(s[ind],eps=eps)
new.etatil[ind] <- new.eta[ind] - 1
}
ind <- (xi*e$g.prime>= eps) & (new.eta == 0)
if (sum(ind)>0) { # new variable to with positive b
DeltaJ[ind] <- xi*(1-e$grad[ind])*plus.recip(e$g.prime[ind],eps=eps)
new.etatil[ind] <- 1
}
ind <- (xi*e$g.prime<= -eps) & (new.eta == 0)
if (sum(ind)>0) { # new variable to with negative b
DeltaJ[ind] <- xi*(-1-e$grad[ind])*plus.recip(e$g.prime[ind],eps=eps)
new.etatil[ind] <- -1
}
# compute Delta, new.etatil, full.path
old.etatil <- rep(m+2,p+1) # again m+2 is the dummy value
old.etatil[etatil[1,]] <- etatil[2,]
old.etatil <- old.etatil[1:p] # old etatil in long format
if ( all(DeltaJ == (-1)) ) { # all infinity
Delta <- -1 # lse has attained for m > 1 and will be the next for lasso
full.path <- TRUE
index.new.etatil <- rep(FALSE,p) # no new sides, all done
}
else { # some stopping time is finite
Delta <- min(DeltaJ[DeltaJ != (-1)])
# index.new.etatil <- (DeltaJ < (Delta + eps)) & (DeltaJ !=(-1)) # does not work well
index.new.etatil <- DeltaJ==Delta # sides hit, numerically delicate
ind <- (abs(s) < eps) & (new.eta != 0) & (old.etatil !=m+2)
if (sum(ind)>0) { # sides not moving
new.etatil[ind] <- old.etatil[ind]
index.new.etatil[ind] <- TRUE
}
ind <- (abs(e$g.prime) < eps) & (new.eta == 0) & (old.etatil !=m+2)
if (sum(ind)>0) { # sides not moving
new.etatil[ind] <- old.etatil[ind]
index.new.etatil[ind] <- TRUE
}
} # end of else
# calculate the new etatil
index.new.etatil <- c(index.new.etatil,TRUE) # add dummy
e$etatil2[2,1:p] <- new.etatil # other entries of e$etatil2 never change value
etatil <- e$etatil2[,index.new.etatil]
# etatil is m+2 at dummy column p+1 when full.path == T
# still need to compute return(tau,b)
if (Delta != -1) {
tau <- tau + xi*Delta
b <- b + xi*Delta*s # update b, b[k+1,] <- b[k,] + ...
e$grad <- e$grad + xi*Delta*e$g.prime
if (tau * eps > 1/5) { # end at perfect fit with zero gradient numerically
full.path <- TRUE
index.new.etatil <- rep(FALSE,p) #
}
if (tau < e$tau1 - eps) { # shall not return to the beginning
forced.stop <- TRUE
}
}
else {
tau <- xi
b <- s
}
} # end if forced.stop == F
return(list(new.eta=new.eta, etatil=etatil, s=s, tau=tau, b=b, singular.Q=singular.Q,
forced.stop = forced.stop, full.path = full.path))
}
########################################################
# check the validity of a specific parallelepiped
###
plus.check.eta <- function(old.eta, new.eta, etatil, xi, z, v, eps=1e-15, e) {
### test a possible new eta and xi combination
## input
# old.eta, 1 X p: old parallelepiped indicator; new.eta: the new one to be tried
# etatil, 2 X (|C|+1): label and new indicator for max crossing, with a dummy variable p+1
# xi: the sign of new.tau - old.tau
# z: t(x)y/n; v: 2nd derivative of penalty; eps: numerical zero; e: environment including x, y
## out put: return(xi, s, singular.Q, invalid.cross)
# singular.Q: error for solve(Q, z[a.set])
# invalid.cross: either singular.Q or fail to cross to new.eta
p <- dim(e$x)[2]
a.set <- rep(TRUE,p)
a.set[new.eta==0] <- FALSE
## update e$Sigma, e$var.list and e$a.set.var
## e$Sigma[,e$a.set.var] will equal to t(x)%*%x[,a.set]/n when done
if (e$use.Gram) e$a.set.var <- a.set # e$Sigma is complete
else if ( sum(e$var.list) == 0 ) { # the initial case of e$Sigma = 0
e$var.list <- order(!a.set)[1:sum(a.set)]
e$a.set.var <- 1:sum(a.set)
e$Sigma[, e$a.set.var] <- e$tx %*% e$x[,a.set]
}
else { # incrementally add new columns to e$Sigma if necessary
new.var <- a.set
new.var[e$var.list] <- FALSE # new.var are those in a.set but nor in e$var.list
if (sum(new.var)>0) { # need to add
e$var.list <- c(e$var.list, order(!new.var)[1:sum(new.var)])
if (length(e$var.list) > (dim(e$Sigma)[2])) { # need to make e$Sigma larger
tmp.Sigma <- e$Sigma
e$Sigma <- matrix(0,p,dim(tmp.Sigma)[2]+100)
e$Sigma[,1:(dim(tmp.Sigma)[2])] <- tmp.Sigma
}
e$Sigma[,(length(e$var.list) - sum(new.var)+1):length(e$var.list)] <- e$tx %*% e$x[,new.var]
}
e$a.set.var <- (1:length(e$var.list))[a.set[e$var.list]] # the equivalent labels of a.set
e$a.set.var <- e$a.set.var[order(e$var.list[ e$a.set.var])]
}
## calculate for return(xi, singular.Q, s, invalid.cross), xi is given
if (sum(a.set)==1) Q <- e$Sigma[a.set,e$a.set.var] - v[ abs(new.eta[a.set]) ]
else Q <- e$Sigma[a.set,e$a.set.var] - diag( v[ abs(new.eta[a.set]) ] )
singular.Q <- FALSE
# LINPACK=T generated an incorrect calculation, due to the use of single precision
# s1 <- try(solve(Q, z[a.set], LINPACK = (length(Q)>100)), silent=T)
# s1 <- try(solve(Q, z[a.set], symmetry=T), silent=T) # option has no effect
s1 <- try(solve(Q, z[a.set]), silent=TRUE)
if (inherits(s1,"try-error")==TRUE) {
s1 <- rep(0, sum(a.set))
singular.Q <- TRUE
}
s <- rep(0,p)
if (singular.Q==FALSE) {
s[ a.set ] <- s1
## check the new s
etatil2 <- c(old.eta,0) # etatil in full length
etatil2[ etatil[1,] ] <- etatil[2,]
etatil2 <- etatil2[1:p] # remove dummy
# e$g.prime <- z - as.matrix((e$Sigma[,e$a.set.var]))%*%s1 # Eq. (2.13)
if (length(s1)==1) e$g.prime <- z - (e$Sigma[,e$a.set.var])*s1 # Eq. (2.13)
else e$g.prime <- z - (e$Sigma[,e$a.set.var])%*%s1
flag <- rep(TRUE,p) # T means violation
flag[old.eta==etatil2] <- FALSE # noncritical j are fine
# check according to (2.15) with eps = numerical zero
flag[ (old.eta != new.eta) & (new.eta != 0) & (xi*(new.eta-old.eta)*s > -eps) ] <- FALSE
flag[ (old.eta == new.eta) & (new.eta != 0) & (etatil2 != old.eta) & (xi*(etatil2-old.eta)*s < eps) ] <- FALSE
flag[ (old.eta != 0) & (new.eta == 0) & (xi*old.eta*e$g.prime < eps)] <- FALSE
flag[ (old.eta == 0) & (new.eta == 0) & (etatil2 != 0) & (xi*etatil2*e$g.prime < eps)] <- FALSE
invalid.cross <- sum(flag) > 0 # T here iff the new eta and xi pair is invalid and singular.Q==FALSE
}
else invalid.cross <- TRUE # invalid.cross is T if singular.Q == T
return(list(xi=xi, singular.Q=singular.Q, s=s, invalid.cross=invalid.cross))
}
########################################################
# find the segment and weight in the path for the hitting point
###
plus.hit.points <- function (lam, lam.path) {
if (is.null(lam) || !is.numeric(lam) || is.nan(sum(lam)))
{
lam <- sort(lam.path, decreasing = TRUE)
warning("Warning message: ")
warning("invalid lam; lam.path used\n")
}
if (max(lam)< min(lam.path))
{
lam <- sort(lam.path, decreasing = TRUE)
warning("Warning message: ")
warning("lam not reached by the plus path; lam.path used\n")
}
lam <- mapply(function(x) min(x,max(lam.path)),lam)
if (min(lam) < min(lam.path))
{
lam <- lam[lam>=min(lam.path)]
warning("Warning message: ")
warning("some lam not reached by the plus path and dropped\n")
}
lam.lesser.path <- mapply(function(x) x<=lam.path,lam)
lam.greater.path <- mapply(function(x) x>=lam.path,lam)
hit <- (lam.lesser.path[1:(length(lam.path)-1),] &
lam.greater.path[2:length(lam.path),])
k1 <- apply(as.matrix(hit),2,function(bool.vec) which(bool.vec)[1])
k2 <- k1+1
w1 <- (lam - lam.path[k2])/(lam.path[k1]-lam.path[k2])
return(list("k1"=k1,"k2"=k2,"w1"=w1,"total.hits"=length(lam)))
}
########################################################
# create a new eta indicating a new parallelepiped
###
plus.new.eta <- function(eta, etatil, cross,p) {
eta <- c(eta,0) # put the dummy in the zero interval
eta[ etatil[1,cross] ] <- etatil[2, cross]
eta <- eta[1:p] # remove dummy
return(eta)
}
########################################################
# check if a newly created eta is indeed new
###
plus.not.in.loop <- function( k, new.eta,e) {
## input: eta, (k+1)Xp; new.eta, 1Xp
## output is T iff new.eta is not any of eta[i,], i\le k
new.is.new <- TRUE
for (i in 1:k)
if( all(e$eta[i,]==new.eta)) return(FALSE)
return(new.is.new)
}
########################################################
# compute the knots and second derivative of a penalty
###
plus.penalty <- function(m,a=0) {
### lasso for m=1, MCP for m=2, SCAD for m=3; default gamma = a = 3.7
### vectors v and t with dummy are used for lasso to maintain vector attribute
m <- round(m)
if (m < 1) m <- 2
if (m > 3) m <- 2
if (m == 2) {
if (a==0) a <- 3.7
v <- c(1/a,0)
t <- c(0,a)
}
if (m == 1) {
a <- 0
v <- c(0,0)
t <- c(0,0)
}
if (m==3) {
if (a==0) a <- 3.7
v <- c(0,1/(a-1),0)
t <- c(0,1,a)
}
return(list(m=m, gamma=a, v=v,t=t))
}
########################################################
# stablized reciprocals
###
plus.recip <- function(Q, eps=1e-15){
### compute reciprocal of reals
Q[ (Q>=0) & (Q <eps) ] <- eps
Q[ (Q<0) & (Q> - eps) ] <- - eps
return(1/Q)
}
########################################################
# generate knots for both sides of zero
###
plus.knots <- function(t=0,m=0) {
### compute the function t() as in (2.6) without infinity and - infinity
### note that t_1=t(1)=t(0) in (2.6)
if (m==0) { # default MC penalty
m <- 2
t <- c(0,3)
}
t2 <- rep(0,2*m)
t2[1:m] <- - t[m:1]
t2[(m+1):(2*m)] <- t[1:m]
return(c(t2,0))
}
##### The above are the plus function
lmres <- function(Y, X, delta) {
X <- as.matrix(X)
o <- order(Y)
Y1 <- Y[o]
X1 <- X[o, ]
delta1 <- delta[o]
Wn1 <- delta1
n1 <- length(Y)
Wn1[1] <- delta1[1]/n1
tt <- 1
for (i in 2:n1) {
j <- i - 1
tt <- tt * ((n1 - j)/(n1 - j + 1))^delta1[j]
Wn1[i] <- delta1[i]/(n1 - i + 1) * tt
}
W <- sqrt(Wn1)
Y2 <- W * Y1
X2 <- W * X1
X2 <- as.matrix(X2)
id <- which(Wn1 > 0)
if (length(id) >= dim(X)[2]) {
Y3 <- Y2[id]
X3 <- X2[id, ]
W3 <- W[id]
X3 <- as.matrix(X3)
lmr <- lm(Y3 ~ W3 + X3 - 1)
return(list(length(Y) * sum((lmr$residuals)^2), lmr$coefficients,
sqrt(length(Y)) * Y3, sqrt(length(Y)) * X3, sqrt(length(Y)) *
W3))
}
if (length(id) < dim(X)[2])
return(NULL)
}
cumcp <- function(Y, X, delta) {
X <- as.matrix(X)
p <- dim(X)[2]
n <- dim(X)[1]
cp <- rep(0, n)
delta1 <- which(delta > 0)
n1 <- length(which(delta > 0))
if (n1 > 2 * p) {
for (i in delta1[p + 1]:delta1[n1 - p]){
cpi <- try(lmres(Y[1:i],
X[1:i, ],
delta[1:i])[[1]] + lmres(Y[(i + 1):n], X[(i + 1):n,],
delta[(i + 1):n])[[1]],silent = TRUE)
if(!is(cpi,"try-error")){
cp[i] <- cpi
}
}
cpm <- min(cp[which(cp > 0)])
return(which(cp == cpm))
} else return(0)
}
#' TSMCP: Two stage multiple change points detection for AFT model.
#'
#' This function first formulates the threshold problem as a group model selection problem
#' so that a concave 2-norm group selection method can be applied using the `grpreg`
#' package in R, and then finalizes it via a refining method.
#'
#' @param Y the censored logarithm of the failure time.
#' @param X the design matrix without the intercept.
#' @param delta the censoring indicator.
#' @param c the length of each segment in the splitting stage, defined as
#' \code{ceiling(c * sqrt(length(Y)))}.
#' @param penalty Penalty type (default is "scad").
#'
#' @return An object with the following components:
#' \describe{
#' \item{cp}{the change points.}
#' \item{coef}{the estimated coefficients.}
#' \item{sigma}{the variance of the error.}
#' \item{residuals}{the residuals.}
#' \item{Yn}{weighted Y by Kaplan-Meier weight.}
#' \item{Xn}{weighted Xn by Kaplan-Meier weight.}
#' }
#'
#' @export
#'
#' @references
#' Li, Jialiang, and Baisuo Jin. 2018. “Multi-Threshold Accelerated Failure Time Model.”
#' The Annals of Statistics 46 (6A): 2657–82.
#'
#' @seealso
#' \code{\link{grpreg}}
#'
#' @examples
#' \donttest{
#' library(grpreg)
#' # Generate simulated data with 500 samples and normal error distribution
#' dataset <- MTAFT_simdata(n = 500, err = "normal")
#' Y <- dataset[, 1]
#' delta <- dataset[, 2]
#' Tq <- dataset[, 3]
#' X <- dataset[, -c(1:3)]
#' n1 = sum(delta)
#' c=seq(0.5,1.5,0.1)
#' m=ceiling(c*sqrt(n1))
#' bicy= rep(NA,length(c))
#' tsmc=NULL
#' p = ncol(X)
#' for(i in 1:length(c)){
#' tsm=try(TSMCP(Y,X,delta,c[i],penalty = "scad"),silent=TRUE)
#' if(is(tsm,"try-error")) next()
#' bicy[i]=log(n)*((length(tsm[[1]])+1)*(p+1))+n*log(tsm[[3]])
#' tsmc[[i]]=tsm
#' }
#'
#' if((any(!is.na(bicy)))){
#' tsmcp=tsmc[[which(bicy==min(bicy))[1]]]
#' thre.LJ = Tq[tsmcp[[1]]]
#' thre.num.Lj = length(thre.LJ)
#' thre.LJ
#' thre.num.Lj
#' }
#' }
TSMCP <- function(Y, X, delta, c,penalty = "scad") {
n <- length(Y)
X <- as.matrix(X)
p <- dim(X)[2]
id1 <- which(delta == 1)
n1 <- length(id1)
m <- ceiling(c * sqrt(n1))
q <- floor(n1/m)
x <- NULL
y <- NULL
de <- NULL
tt <- lmres(Y[1:id1[n1 - (q - 1) * m]], X[1:id1[n1 - (q - 1) * m],], delta[1:id1[n1 - (q - 1) * m]])
y[[1]] <- tt[[3]]
x[[1]] <- cbind(tt[[5]], tt[[4]])
for (i in 2:(q - 1)) {
tt <- lmres(Y[id1[n1 - (q - i + 1) * m + 1]:id1[n1 - (q - i) *
m]], X[id1[n1 - (q - i + 1) * m + 1]:id1[n1 - (q - i) * m],
], delta[id1[n1 - (q - i + 1) * m + 1]:id1[n1 - (q - i) * m]])
y[[i]] <- tt[[3]]
x[[i]] <- cbind(tt[[5]], tt[[4]])
}
i <- q
tt <- lmres(Y[id1[n1 - (q - i + 1) * m + 1]:n], X[id1[n1 - (q - i +
1) * m + 1]:n, ], delta[id1[n1 - (q - i + 1) * m + 1]:n])
y[[i]] <- tt[[3]]
x[[i]] <- cbind(tt[[5]], tt[[4]])
K_temp <- matrix(0, nrow = q, ncol = q, byrow = TRUE)
for (i in 1:q) K_temp[i, 1:i] <- rep(1, i)
X_temp1 <- lapply(1:length(x), function(j, mat, list) kronecker(K_temp[j,
, drop = FALSE], x[[j]]),
mat = K_temp, list = x)
Xn <- do.call("rbind", X_temp1)
Yn <- NULL
for (i in 1:q) {
Yn <- c(Yn, y[[i]])
}
group <- kronecker(c(1:q), rep(1, p + 1))
fit <- grpreg(Xn, Yn, group, penalty = "grSCAD")
mcp.coef <- grpreg::select(fit, "BIC")$beta[-1]
mcp.coef.s <- sum(abs(mcp.coef))
mcp.coef.v.m <- abs(matrix(c(mcp.coef), q, (p + 1), byrow = T))
mcp.coef.m <- c(apply(mcp.coef.v.m, 1, max))
mcp.cp <- which(mcp.coef.m != 0)
if (q > 10) {
if (length(mcp.cp) > 1) {
for (i in 2:length(mcp.cp)) {
if (mcp.cp[i] - mcp.cp[i - 1] == 1)
mcp.cp[i] <- 0
}
}
mcp.cp1 <- mcp.cp[mcp.cp > 1 & mcp.cp < q]
d1 <- length(mcp.cp1)
if (d1 == 0) {
mcpcss.cp <- NULL
adcpcss.cp <- mcpcss.cp
}
if (d1 >= 1) {
mcpcss.cp <- NULL
for (i in 1:d1) {
t <- mcp.cp1[i]
id <- c(id1[n1 - (q - t + 2) * m + 1]:id1[n1 - (q - t -
1) * m])
cp.t <- cumcp(Y[id], X[id, ], delta[id])[[1]]
if (cp.t > 0)
mcpcss.cp <- c(mcpcss.cp, cp.t + id1[n1 - (q - t + 2) *
m + 1] - 1)
}
id2 <- c(mcpcss.cp, n)
tt1 <- NULL
tt1 <- c(tt1, table(delta[1:mcpcss.cp[1]])[2])
for (i in 2:length(id2)) {
tt1 <- c(tt1, table(delta[(id2[i - 1] + 1):id2[i]])[2])
}
tt2 <- which(tt1 < p + 2)
tt2[which(tt2 == length(id2))] <- length(id2) - 1
if (length(tt2) > 0) {
mcpcss.cp <- mcpcss.cp[-tt2]
}
tt <- which(abs(diff(mcpcss.cp)) < 2 * m)
if (length(tt) > 0)
mcpcss.cp <- mcpcss.cp[-tt]
}
}
if (q <= 10) {
mcp.cp1 <- mcp.cp[mcp.cp > 1 & mcp.cp < q]
d1 <- length(mcp.cp1)
if (d1 == 0) {
mcpcss.cp <- NULL
adcpcss.cp <- mcpcss.cp
}
if (d1 >= 1) {
mcpcss.cp <- NULL
for (i in 1:d1) {
t <- mcp.cp1[i]
if (t == 2) {
id <- c(1:id1[n1 - (q - t - 1) * m])
cp.t <- cumcp(Y[id], X[id, ], delta[id])[[1]]
if (cp.t > 0)
mcpcss.cp <- c(mcpcss.cp, cp.t)
}
if (t == q) {
id <- c(id1[n1 - (q - t + 2) * m + 1]:n)
cp.t <- cumcp(Y[id], X[id, ], delta[id])[[1]]
if (cp.t > 0)
mcpcss.cp <- c(mcpcss.cp, cp.t + id1[n1 - (q - t +
2) * m + 1] - 1)
}
if (t < q & t > 2) {
id <- c(id1[n1 - (q - t + 2) * m + 1]:id1[n1 - (q - t -
1) * m])
cp.t <- cumcp(Y[id], X[id, ], delta[id])[[1]]
if (cp.t > 0)
mcpcss.cp <- c(mcpcss.cp, cp.t + id1[n1 - (q - t +
2) * m + 1] - 1)
}
}
id2 <- c(mcpcss.cp, n)
tt1 <- NULL
tt1 <- c(tt1, table(delta[1:mcpcss.cp[1]])[2])
for (i in 2:length(id2)) {
tt1 <- c(tt1, table(delta[(id2[i - 1] + 1):id2[i]])[2])
}
tt2 <- which(tt1 < p + 2)
tt2[which(tt2 == length(id2))] <- length(id2) - 1
if (length(tt2) > 0) {
mcpcss.cp <- mcpcss.cp[-tt2]
}
tt <- which(abs(diff(mcpcss.cp)) < 2 * m)
if (length(tt) > 0)
mcpcss.cp <- mcpcss.cp[-tt]
}
}
if (length(mcpcss.cp) >= 1 && mcpcss.cp[1] != 0) {
id2 <- c(mcpcss.cp, n)
x <- NULL
y <- NULL
de <- NULL
tt <- lmres(Y[1:mcpcss.cp[1]], X[1:mcpcss.cp[1], ], delta[1:mcpcss.cp[1]])
y[[1]] <- tt[[3]]
x[[1]] <- cbind(tt[[5]], tt[[4]])
for (i in 2:length(id2)) {
tt <- lmres(Y[(id2[i - 1] + 1):id2[i]], X[(id2[i - 1] + 1):id2[i],
], delta[(id2[i - 1] + 1):id2[i]])
y[[i]] <- tt[[3]]
x[[i]] <- cbind(tt[[5]], tt[[4]])
}
K_temp <- matrix(0, nrow = length(id2), ncol = length(id2), byrow = TRUE)
for (i in 1:length(id2)) K_temp[i, 1:i] <- rep(1, i)
X_temp1 <- lapply(1:length(x), function(j, mat, list) kronecker(K_temp[j,
, drop = FALSE], x[[j]]), mat = K_temp, list = x)
Xn <- do.call("rbind", X_temp1)
Yn <- NULL
for (i in 1:length(id2)) {
Yn <- c(Yn, y[[i]])
}
#//...................................................................................................
object <- plus(Xn, Yn, method = penalty, gamma = 2.4, intercept = F,normalize = F, eps = 1e-30)
# step 1 estimate coef using BIC.
bic <- log(dim(Xn)[1]) * object$dim + dim(Xn)[1] * log(as.vector((1 -object$r.square) * sum(Yn^2))/length(Yn))
step.bic <- which.min(bic)
mcp.coef <- coef(object, lam = object$lam[step.bic])
#//...................................................................................................
#reg.fit <- ncvreg(X = Xn, y = drop(Yn), family = "gaussian",penalty = penalty)
#coef.beta <- reg.fit$beta # extract coefficients at all values of lambda, including the intercept
#reg.df <- apply(abs(coef.beta[-1, , drop = FALSE]) > 1e-10,2, sum)
#BIC <- reg.fit$loss + log(length(Yn)) * reg.df
#lambda.ind <- which.min(BIC)
#mcp.coef <- coef.beta[-1, lambda.ind]
#..................
ep <- (Yn - Xn %*% mcp.coef)
sigmaep <- sum(ep^2)/length(Y)
mcp.coef.v.m <- matrix(c(mcp.coef), length(id2), p + 1, byrow = T)
mcp.coef.m <- c(apply(abs(mcp.coef.v.m), 1, max))
mcp.m.id <- which(mcp.coef.m == 0)
if (length(mcp.m.id) > 0) {
mcp.coef <- c(t(mcp.coef.v.m[-mcp.m.id, ]))
mcpcss.cp <- mcpcss.cp[-c(mcp.m.id - 1)]
}
}
if (mcpcss.cp == 0 || length(mcpcss.cp) == 0) {
id2 <- c(mcpcss.cp, n)
x <- NULL
y <- NULL
de <- NULL
tt <- lmres(Y, X, delta)
Yn <- tt[[3]]
Xn <- cbind(tt[[5]], tt[[4]])
object <- plus(Xn, Yn, method = "scad", gamma = 2.4, intercept = F,
normalize = F, eps = 1e-30)
# step 1 estimate coef using BIC.
bic <- log(dim(Xn)[1]) * object$dim + dim(Xn)[1] * log(as.vector((1 -
object$r.square) * sum(Yn^2))/length(Yn))
step.bic <- which.min(bic)
mcp.coef <- coef(object, lam = object$lam[step.bic])
ep <- (Yn - Xn %*% mcp.coef)
sigmaep <- sum(ep^2)/length(Y)
mcp.coef.v.m <- matrix(c(mcp.coef), length(id2), p + 1, byrow = T)
mcp.coef.m <- c(apply(abs(mcp.coef.v.m), 1, max))
mcp.m.id <- which(mcp.coef.m == 0)
if (length(mcp.m.id) > 0) {
mcp.coef <- c(t(mcp.coef.v.m[-mcp.m.id, ]))
mcpcss.cp <- mcpcss.cp[-c(mcp.m.id - 1)]
}
}
return(list(cp = mcpcss.cp, coef = mcp.coef, sigma = sigmaep, residuals = ep,
Yn = Yn, Xn = Xn))
}
Try the MTAFT package in your browser
Any scripts or data that you put into this service are public.
MTAFT documentation built on Nov. 14, 2023, 1:08 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions TSMCP cumcp lmres plus.knots plus.recip plus.penalty plus.not.in.loop plus.new.eta plus.check.eta plus.single.step plot.plus coef.plus predict.plus print.plus plus
Documented in TSMCP
####This R script mainly implements the two-stage method of Li and Jin (2018)
########################################################
# plus(), the main plus function; documented
###
plus <- function(x,y, method = c("lasso", "mc+", "scad", "general"), m=2, gamma,v,t,
monitor=FALSE, normalize = TRUE, intercept = TRUE,
Gram, use.Gram = FALSE, eps=1e-15, max.steps=500, lam)
{
## fill missing values
if (missing(gamma)) gamma <- 0
if (missing(v)) v <- 0
if (missing(t)) t <- 0
if (missing(Gram)) Gram <- 0
if (missing(lam)) lam <- -1
lam.min <- min(lam)
## check data
data.dim.error <- FALSE
if (! is.matrix(x)) data.dim.error <- TRUE
else {
n <- dim(x)[1]
if ( (n<=1) | (n!=length(y)) ) data.dim.error <- TRUE
p <- dim(x)[2]
if (p<= 1) data.dim.error <- TRUE
}
if (data.dim.error == TRUE){
warning("Data dimention error.", "\n")
return(list(data.dim.error = TRUE))
}
## standardize x
orig.x <- x
if (intercept) x <- t(t(x)-apply(x,2,mean)) # center to mean zero
if (normalize) {
normalize.factor <- sqrt(apply(x^2,2,mean))
x <- t(t(x)/normalize.factor)
}
## compute Gram matrix if necessary
## if ((!use.Gram) & (p<=200)) use.Gram <- TRUE
if ((use.Gram) & (p > 1000)) use.Gram <- FALSE
if (use.Gram)
if (! is.matrix(Gram)) Gram <- t(x)%*%x
else if ((dim(Gram)[1]!=p) || (dim(Gram)[2]!=p)) Gram <- t(x)%*%x
## determine method and compute penalty
if (length(method)==1)
if (method == "lasso") m <- 1
else if (method == "mc+") m <- 2
else if (method == "scad") m <- 3
if (m==1) method <- "LASSO"
else if (m==2) method <- "MC+"
else if (m==3) method <- "SCAD"
else method <- "PLUS"
if (monitor == TRUE) {
cat("Sequence of", method, "moves:", "\n","\n")
}
if (m<4) { # provide default penalty for lasso, MCP and SCAD
if ((m>1) & (use.Gram) & (gamma == 0)) {
max.corr <- max(abs(Gram[ row(Gram) != col(Gram) ]))/n
if ((1 - max.corr) >= (1e-5) ) gamma <- 2/(1- max.corr)
}
if ((m==3) & (gamma <=1)) gamma <- 0
tmp <- plus.penalty(m,a=gamma)
m <- tmp$m
gamma <- tmp$gamma
v <- tmp$v
t <- tmp$t
}
## declare variables
names <- dimnames(x)
b <- matrix(0, max.steps+1,p)
etatil <- matrix(0,2,p+1)
tau <- rep(0, max.steps+1)
my.env <-environment()
my.env$x <-x
my.env$y <- y
my.env$tx <- t(x)/n
my.env$eta <- b
my.env$g.prime <- rep(0,p)
my.env$grad <- rep(0,p)
my.env$etatil2 <- etatil
my.env$etatil2[1,] <- 1:(p+1)
if (use.Gram) my.env$Sigma <- Gram/n
else {
my.env$Sigma <- matrix(0,p,100) # e$a.set.var will match a.set
my.env$var.list <- 0
}
my.env$use.Gram <- use.Gram
my.env$ties <- rep(FALSE,p)
my.env$t.fun <- plus.knots(t,m)
z <- my.env$tx %*%y
my.env$tau1 <- 1/max(abs(z))
## initialization: for the initial step k=0, tau^{(0)} is needed, b[1,] <- 0 done already
tau[1] <- my.env$tau1
etatil[1,] <- 1:(p+1) # the first etatil
etatil[2,1:p] <- sign(z)
etatil <- etatil[ , c(abs(z) >= max(abs(z)) - eps,TRUE) ]
k <-1
exit.while <- FALSE
last.valid.s <- rep(0,p)
## the main loop
while (exit.while==FALSE) { # segment k: model eta[k+1,] begins with b[k,] and ends with b[k+1,]
tmp <- plus.single.step(k, etatil, z, b[k,], tau[k], m, v, t, eps=eps,e=my.env)
## tmp holds return(new.eta, etatil, s, tau, b, singular.Q, forced.stop, full.path)
etatil <- tmp$etatil # etatil = 0 if forced.stop, length(etatil)=2 if Ise.termination
exit.while <- tmp$forced.stop # cannot find a valid new.eta
full.path <- tmp$full.path
if (exit.while == FALSE) { # save data from tmp if not forced to stop
eta[k+1,] <- tmp$new.eta
if (monitor == TRUE) { # print info to monitor the iterations
var.change <- sign(abs(sign(eta[k+1,]))-abs(sign(eta[k,])))
if (sum(abs(var.change))==1) {
if(sum(var.change)==1)
action <-"added "
else
action <-"removed "
var_name <- (order(-abs(var.change))[1])
if (!is.null(names))
var_name <- names[[2]][var_name]
cat("Step",k,": ",action,var_name,"\n")
}
}
last.valid.s <- tmp$s # save the s for the last valid eta
tau[k+1]<- tmp$tau
b[k+1,] <- tmp$b # b[k+1] is at the end of segment eta[k+1,]
b[k+1,eta[k+1,]==0] <- 0 # more accurate b=0, b[k,eta[k+1,]==0] <- 0 is also ok
k <- k+1
if ((k > max.steps)|(full.path)|(tau[k]*lam.min > 1)) { # stop with a valid eta
exit.while <- TRUE
if ((full.path) & (m>1)) k <- k-1
}
} # if (exit.while == FALSE)
} # end the main loop
## tmp holds return(new.eta, etatil, s, tau, b, singular.Q, forced.stop, full.path)
singular.Q <- tmp$singular.Q
forced.stop <- tmp$forced.stop
lam.path <- 1/tau[1:k]
beta.path <- lam.path * b[1:k,]
if (normalize) {
beta.path <- t(t(beta.path)/normalize.factor)
}
lam.path[k] <- max(abs(my.env$tx %*%(y-x%*%beta.path[k,])))
if (full.path) lam.path[k] <- 0
eta <- eta[1:k,]
total.steps <- k-1
## calculate the estimator beta at the given lam or the set of values of lam.path
x <- orig.x
if (sum(lam) == -1) lam <- sort(lam.path, decreasing=TRUE)
tmp <- plus.hit.points(lam,lam.path)
total.hits <- tmp$total.hits
if (total.hits > 0) {
beta <- tmp$w1*beta.path[tmp$k1,] + (1-tmp$w1)*beta.path[tmp$k2,]
lam <- lam[1:total.hits]
dim <- apply(abs(beta)>0, 1, sum)
if (total.hits == 1)
r.square <- 1- sum((y - x%*%beta)^2)/sum(y^2)
else
r.square <- 1- apply((x%*%t(beta) - y)^2, 2, sum)/sum(y^2)
obj<-list(x=x,y=y,eta=eta, beta.path=beta.path, lam.path=lam.path,
beta=beta, lam=lam, dim=dim, r.square = r.square, total.hits=total.hits,
method = method, gamma = gamma, total.steps=total.steps, max.steps=max.steps,
full.path=full.path, forced.stop=forced.stop, singular.Q=singular.Q)
}
else
obj<-list(x=x,y=y,eta=eta, beta.path=beta.path, lam.path=lam.path, total.hits=total.hits,
method = method, gamma = gamma, total.steps=total.steps, max.steps=max.steps,
full.path=full.path, forced.stop=forced.stop, singular.Q=singular.Q)
class(obj) <- "plus"
return(obj)
}
########################################################
# print.plus(), to print plus() moves; documented
###
print.plus <- function(x, print.moves = 20, ...) {
#cat("\n")
#cat("Sequence of", x$method, "moves:", "\n","\n")
print.moves <- round(print.moves)
if (print.moves < 1) print.moves <- 20
names <- dimnames(x$x)
m <- 1
k <- 1
exit.while <- FALSE
while (exit.while == FALSE) {
var.change <- sign(abs(sign(x$eta[k+1,]))-abs(sign(x$eta[k,])))
if (sum(abs(var.change))==1) {
if(sum(var.change)==1)
action <-"added "
else
action <-"removed "
var_name <- (order(-abs(var.change))[1])
if (!is.null(names))
var_name <- names[[2]][var_name]
#cat("Step",k,": ",action,var_name,"\n")
m <- m+1
if (m>print.moves) exit.while <- TRUE
}
k <- k+1
if (k > x$total.steps) exit.while <- TRUE
} # end while
}
########################################################
# predict.plus(), to extract coefficients and preditions; documented
###
predict.plus <- function(object,lam,newx, ...) {
if (missing(newx)) {
flag <- FALSE
warning("Warning message: ")
warning("no newx argument; newy not produced", "\n")
}
else flag <- TRUE
if (missing(lam)) {
lam <- sort(object$lam.path, decreasing=TRUE)
warning("Warning message: ")
warning("no lam argument; object$lam.path used", "\n")
}
else if (max(lam) < min(object$lam.path)) {
lam <- sort(object$lam.path, decreasing=TRUE)
warning("Warning message: ")
warning("lam not reached by the plus path; object$lam.path used", "\n")
}
#cat("\n")
tmp <- plus.hit.points(lam,object$lam.path)
total.hits <- tmp$total.hits
beta <- tmp$w1* object$beta.path[tmp$k1,] + (1-tmp$w1)* object$beta.path[tmp$k2,]
lam <- lam[1:total.hits]
dim <- apply(abs(beta)>0, 1, sum)
if (total.hits == 1) {
r.square <- 1- sum((object$y - object$x%*%beta)^2)/sum(object$y^2)
if (flag) {
if (is.matrix(newx)) newy <- newx%*% beta
else newy <- sum(beta * newx)
}
}
else {
r.square <- 1- apply((object$x%*%t(beta) - object$y)^2, 2, sum)/sum(object$y^2)
if (flag) {
if (is.matrix(newx)) newy <- beta %*% t(newx)
else newy <- beta %*% newx
}
}
if (flag)
return(list(lambda=lam, coefficients=beta, dimension=dim, r.square = r.square, step = tmp$k2-1,
method = object$method, newy = newy))
else
return(list(lambda=lam, coefficients=beta, dimension=dim, r.square = r.square, step = tmp$k2-1,
method = object$method))
}
########################################################
# coef.plus(), to extract coefficients; documented
###
coef.plus <- function(object,lam, ...){
if (missing(lam)) {
lam <- sort(object$lam.path, decreasing=TRUE)
warning("Warning message: ")
warning("no lam argument; object$lam.path used", "\n")
}
else if (max(lam) < min(object$lam.path)) {
lam <- sort(object$lam.path, decreasing=TRUE)
warning("Warning message: ")
warning("lam not reached by the plus path; x$lam.path used", "\n")
}
#cat("\n")
tmp <- plus.hit.points(lam,object$lam.path)
beta <- tmp$w1* object$beta.path[tmp$k1,] + (1-tmp$w1)* object$beta.path[tmp$k2,]
return(beta)
}
########################################################
# plot.plus(), to plot coefficients, preditions, penalty level, dimension or R-square
# againt plus step or penalty level; documented
###
plot.plus <- function(x, xvar=c("lam","step"), yvar=c("coef","newy","lam","dim","R-sq"),
newx, step.interval, lam.interval, predictors, ...) {
# determin xvar in the plot
if (length(unique(c(xvar,"lam")))==length(unique(c(xvar)))) xvar <- "lambda"
else if (length(unique(c(xvar,"step")))==length(unique(c(xvar)))) xvar <- "step"
else xvar <- "lambda"
# determin yvar in the plot
if (length(unique(c(yvar,"coef")))==length(unique(c(yvar)))) yvar <- "coefficients"
else if (length(unique(c(yvar,"newy")))==length(unique(c(yvar)))) yvar <- "newy"
else if (length(unique(c(yvar,"lam")))==length(unique(c(yvar)))) yvar <- "lambda"
else if (length(unique(c(yvar,"dim")))==length(unique(c(yvar)))) yvar <- "dimension"
else if (length(unique(c(yvar,"R-sq")))==length(unique(c(yvar)))) yvar <- "R-square"
else yvar <- "coefficients"
if (yvar != "newy") newx <- matrix(0,2,dim(x$x)[2])
if ((yvar == "newy") & missing(newx)) {
warning("Warning message: ")
warning("no newx argument", "\n")
return()
}
if (missing(predictors)) predictors <- 1:(dim(x$x)[2])
## plot
if (xvar == "step") {
xtmp <- 0:(x$total.steps)
plot.set <- rep(TRUE, length(xtmp))
if (! missing(step.interval)) {
plot.set[xtmp > max(step.interval)] <- FALSE
plot.set[xtmp < min(step.interval)-1] <- FALSE
xtmp <- xtmp[plot.set]
}
if (yvar == "coefficients") ytmp <- x$beta.path[plot.set, ]
if ((yvar == "newy") & is.matrix(newx)) ytmp <- x$beta.path[plot.set, ]%*%t(newx)
if ((yvar == "newy") & (! is.matrix(newx))) ytmp <- x$beta.path[plot.set, ]%*%newx
if (yvar == "lambda") ytmp <- x$lam.path[plot.set]
if (yvar == "dimension") ytmp <- apply(abs(sign(x$eta[plot.set,])),1, sum)
if (yvar == "R-square")
ytmp <- 1- apply((x$x%*%t(x$beta.path[plot.set, ]) - x$y)^2, 2, sum)/sum(x$y^2)
}
if (xvar == "lambda") {
xtmp <- sort(x$lam.path,decreasing=TRUE)
tmp <- predict(x,xtmp,newx)
plot.set <- rep(TRUE, length(xtmp))
if (! missing(lam.interval)) {
plot.set[xtmp > max(lam.interval)] <- FALSE
plot.set[xtmp < min(lam.interval)] <- FALSE
xtmp <- - xtmp[plot.set]
}
else xtmp <- -xtmp
if (yvar == "coefficients") ytmp <- tmp$coefficients[plot.set, ]
if ((yvar == "newy") & is.matrix(newx)) ytmp <- tmp$newy[plot.set,]
if ((yvar == "newy") & (! is.matrix(newx))) ytmp <- tmp$newy[plot.set]
if (yvar == "lambda") ytmp <- tmp$lambda[plot.set]
if (yvar == "dimension") ytmp <- tmp$dimension[plot.set]
if (yvar == "R-square") ytmp <- tmp$r.square[plot.set]
}
xmax <- max(xtmp)
xmin <- min(xtmp)
ymax <- max(ytmp)
ymin <- min(ytmp)
if (xvar == "lambda") {
plot(c(xmin,xmax),c(ymin,ymax), xaxt = "n", xlab = xvar, ylab=yvar, type="l",lty=0, main=x$method)
axis(1, at = axTicks(1), labels = abs(axTicks(1)) )
}
else plot(c(xmin,xmax),c(ymin,ymax), xlab = xvar, ylab=yvar, type="l",lty=0, main=x$method)
if ( ((yvar=="coefficients")|(yvar=="newy")) & is.matrix(ytmp)) {
index <- 1: (dim(ytmp)[2])
if (yvar=="coefficients") {
f.set <- rep(FALSE, max(length(index),max(predictors)))
f.set[predictors] <- TRUE
f.set <- f.set[1:length(index)]
}
else f.set <- rep(TRUE, length(index))
axis(4, at = ytmp[dim(ytmp)[1],f.set], labels = index[f.set])
for (j in 1: (dim(ytmp)[2]))
if (f.set[j]) lines(xtmp, ytmp[,j])
}
else lines(xtmp, ytmp)
}
########################################################
# internal functions used by plus()
###
########################################################
# one single step in the iterative algorithm
###
plus.single.step <- function(k, etatil, z, b, tau, m, v, t, eps=1e-15, e) {
## input
# k: number of parallelepiped indicators already found, including 0, k-1 steps completed
# etatil, 2 X (|C|+1): label and new indicator for max crossing with dummy variable p+1
# NOTE: |C| > 0 required here
# z: t(x)y/n; b, tau: old turning point b at old tau; m, v, t: penalty; eps: numerical zero,
## out put: return(new.eta, etatil, s, tau,b, singular.Q, forced.stop, full.path)
# including info at the end of segment new.eta, e.g. etatil for the eta beyond new.eta
n <- dim(e$x)[1] # no dummy
p <- dim(e$x)[2]
# this is in step k, eta^{(0)},...,eta^{(k-1)} have already been found
sides <- length(etatil)/2-1 # number of boundaries for possible crossing
singular.Q <- FALSE # this and next values are needed in case plus.check.eta not called
forced.stop <- TRUE
if (sides == 1) { # one-at-a-time scenario
new.eta <- plus.new.eta(e$eta[k,],etatil,rep(TRUE,2),p) # T for crossing the boundary and dummy
## always a new eta in theory, only check when k is a multiplier of 50 to prevent numerical loop
if (k != round(k/50)*50) check.new.eta <- TRUE
else check.new.eta <- plus.not.in.loop(k, new.eta, e=e)
if (check.new.eta==TRUE) { # do not run plus.check.eta if eta is not new, try xi = 1 first
tmp <- plus.check.eta(e$eta[k,], new.eta, etatil, 1, z, v, eps=eps,e=e)
## tmp holds return(xi, singular.Q, s, invalid.cross)
singular.Q <- tmp$singular.Q # save if invalid cross is due to singular Q, an almost nonevent
if ( (m>1) & ((tmp$invalid.cross) == TRUE) ) { # try xi = -1
tmp <- plus.check.eta(e$eta[k,], new.eta, etatil, -1, z, v, eps=eps,e=e)
if (singular.Q == FALSE) singular.Q <- tmp$singular.Q
}
forced.stop <- tmp$invalid.cross # due to singular Q or not
} # if (check.new.eta==T)
} # end of the one-at-a-time case
else { # sides > 1, first check ties
etatil <- etatil[, c((! e$ties[ etatil[1,1:sides] ]),TRUE) ] # remove known ties
sides <- length(etatil)/2-1
if (sides > 1) { # remove new ties
if (e$use.Gram) Sgm.B <- e$Sigma[etatil[1,1:sides], etatil[1,1:sides]]
else Sgm.B <- e$tx[etatil[1,1:sides],] %*% e$x[,etatil[1,1:sides]]
for (i in 1:(sides-1)) for (j in (i+1):sides)
if (Sgm.B[i,j] > 1- eps) {
e$ties[j] <- TRUE
}
etatil <- etatil[, c(! e$ties[ etatil[1,1:sides] ],TRUE) ] # remove known ties
sides <- length(etatil)/2-1
}
if (sides == 1) {# go back to the one-at-time case
new.eta <- plus.new.eta(e$eta[k,],etatil,rep(TRUE,2),p) # T for crossing the boundary and dummy
## always a new eta in theory, only check when k is a multiplier of 50 to prevent numerical loop
if (k != round(k/50)*50) check.new.eta <- TRUE
else check.new.eta <- plus.not.in.loop(k, new.eta, e=e)
if (check.new.eta==TRUE) { # do not run plus.check.eta if eta is not new, try xi = 1 first
tmp <- plus.check.eta(e$eta[k,], new.eta, etatil, 1, z, v, eps=eps,e=e)
## tmp holds return(xi, singular.Q, s, invalid.cross)
singular.Q <- tmp$singular.Q # save if invalid cross is due to singular Q, an almost nonevent
if ( (m>1) & ((tmp$invalid.cross) == TRUE) ) { # try xi = -1
tmp <- plus.check.eta(e$eta[k,], new.eta, etatil, -1, z, v, eps=eps,e=e)
if (singular.Q == FALSE) singular.Q <- tmp$singular.Q
}
forced.stop <- tmp$invalid.cross # due to singular Q or not
} # if (check.new.eta==T)
} # end of the one-at-a-time case
else if (sides > 1) { # general case, (more than) two-at-a-time
M <- matrix(FALSE, 2^sides, sides + 1) # generate all possible candidates for new eta
M[, sides+1] <- TRUE # last T for dummy
M[ (2^(sides-1)+1):(2^sides), 1 ] <- TRUE # T for cross the indicated side
for (i in (2:sides) ) M[,i] <- M[(2^(sides-i)+1):(3*2^(sides-i)),i-1]
M <- M[order(apply(M,1,sum)),] # closer eta (smaller number of crossings) has higher priority
# initial while loop look for new eta
continue.while <- TRUE
j <- 1
while (continue.while == TRUE) {
new.eta <- plus.new.eta(e$eta[k,],etatil,M[j+1,],p) # skip the first M[1,]=F...FT
new.is.new <- plus.not.in.loop(k, new.eta, e=e)
if (new.is.new == TRUE) { # try xi =1
tmp <- plus.check.eta(e$eta[k,], new.eta, etatil, 1, z, v, eps=eps, e=e)
## tmp holds return(xi, singular.Q, s, invalid.cross)
## invalid.cross == T iff either (singular.Q == T) or (xi and s fail to cross to new.eta)
if (singular.Q == FALSE) singular.Q <- tmp$singular.Q
if ( (m>1) & ((tmp$invalid.cross) == TRUE)) { # try xi = -1
tmp <- plus.check.eta(e$eta[k,], new.eta, etatil, -1, z, v, eps=eps, e=e)
if (singular.Q == FALSE) singular.Q <- tmp$singular.Q # singular.Q should be T if happens once
}
forced.stop <- tmp$invalid.cross # forced.stop is true if all "not new" or invalid
continue.while <- tmp$invalid.cross # exit while loop for valid cross to new eta
} # end if (new.is.new == T)
j <- j+1 # continue.while does not change when new.is.new == F
if (j==2^sides) continue.while <- FALSE
} # end while, if forced.stop == F, then tmp holds info for a valid new eta
} # end of (more than) two-at-a-time
else if (sides < 1) new.eta <- rep(m+2,p)
}
full.path <- FALSE
if (forced.stop == TRUE) { # terminate the program, need value for the 9 out of 10 var
## dummy output for return(new.eta, etatil, s, tau,b, singular.Q, forced.stop, full.path)
## values have already been assigned to new.eta, singular.Q, forced.stop, full.path
etatil <- c(p+1,m+2) # etatil <- as.matrix(c(p+1,m+2))
s <- rep(0,p)
tau <- 0
b <- rep(0,p)
}
if (forced.stop==FALSE) { # valid new eta has been found as new.eta
## tmp holds return(xi, singular.Q, s, invalid.cross), singular.Q = invalid.cross = F
xi <- tmp$xi
s <- tmp$s
## need values for return(etatil, tau,b, full.path); new.eta, singular.Q, forced.stop done
a.set <- (new.eta!=0)
if ( k - 50*round(k/50) == 1) # directly calculate the gradient to prevent accumulation of error
# e$grad <- tau*z - as.matrix(e$Sigma[,e$a.set.var])%*%b[a.set]
if (sum(a.set)<2)
e$grad <- tau*z - (e$Sigma[,e$a.set.var])*b[a.set] # gradient
else
e$grad <- tau*z - (e$Sigma[,e$a.set.var])%*%b[a.set]
## compute DeltaJ and new etatil
DeltaJ <- rep(-1,p) # -1 represents infinity
new.etatil <- rep(m+2,p) # m+2 is the dummy value for etatil calculation
ind <- (xi*s>=eps) & (new.eta != 0) & (new.eta < m)
if (sum(ind)>0) { # move towards higher eta value
DeltaJ[ind] <- xi*(e$t.fun[new.eta[ind]+1+m]-b[ind])*plus.recip(s[ind],eps=eps)
new.etatil[ind] <- new.eta[ind] + 1
}
ind <- (xi*s<= -eps) & (new.eta != 0) & (new.eta > -m)
if (sum(ind)>0) { # move towards lower eta value
DeltaJ[ind] <- xi*(e$t.fun[new.eta[ind]+m]-b[ind])*plus.recip(s[ind],eps=eps)
new.etatil[ind] <- new.eta[ind] - 1
}
ind <- (xi*e$g.prime>= eps) & (new.eta == 0)
if (sum(ind)>0) { # new variable to with positive b
DeltaJ[ind] <- xi*(1-e$grad[ind])*plus.recip(e$g.prime[ind],eps=eps)
new.etatil[ind] <- 1
}
ind <- (xi*e$g.prime<= -eps) & (new.eta == 0)
if (sum(ind)>0) { # new variable to with negative b
DeltaJ[ind] <- xi*(-1-e$grad[ind])*plus.recip(e$g.prime[ind],eps=eps)
new.etatil[ind] <- -1
}
# compute Delta, new.etatil, full.path
old.etatil <- rep(m+2,p+1) # again m+2 is the dummy value
old.etatil[etatil[1,]] <- etatil[2,]
old.etatil <- old.etatil[1:p] # old etatil in long format
if ( all(DeltaJ == (-1)) ) { # all infinity
Delta <- -1 # lse has attained for m > 1 and will be the next for lasso
full.path <- TRUE
index.new.etatil <- rep(FALSE,p) # no new sides, all done
}
else { # some stopping time is finite
Delta <- min(DeltaJ[DeltaJ != (-1)])
# index.new.etatil <- (DeltaJ < (Delta + eps)) & (DeltaJ !=(-1)) # does not work well
index.new.etatil <- DeltaJ==Delta # sides hit, numerically delicate
ind <- (abs(s) < eps) & (new.eta != 0) & (old.etatil !=m+2)
if (sum(ind)>0) { # sides not moving
new.etatil[ind] <- old.etatil[ind]
index.new.etatil[ind] <- TRUE
}
ind <- (abs(e$g.prime) < eps) & (new.eta == 0) & (old.etatil !=m+2)
if (sum(ind)>0) { # sides not moving
new.etatil[ind] <- old.etatil[ind]
index.new.etatil[ind] <- TRUE
}
} # end of else
# calculate the new etatil
index.new.etatil <- c(index.new.etatil,TRUE) # add dummy
e$etatil2[2,1:p] <- new.etatil # other entries of e$etatil2 never change value
etatil <- e$etatil2[,index.new.etatil]
# etatil is m+2 at dummy column p+1 when full.path == T
# still need to compute return(tau,b)
if (Delta != -1) {
tau <- tau + xi*Delta
b <- b + xi*Delta*s # update b, b[k+1,] <- b[k,] + ...
e$grad <- e$grad + xi*Delta*e$g.prime
if (tau * eps > 1/5) { # end at perfect fit with zero gradient numerically
full.path <- TRUE
index.new.etatil <- rep(FALSE,p) #
}
if (tau < e$tau1 - eps) { # shall not return to the beginning
forced.stop <- TRUE
}
}
else {
tau <- xi
b <- s
}
} # end if forced.stop == F
return(list(new.eta=new.eta, etatil=etatil, s=s, tau=tau, b=b, singular.Q=singular.Q,
forced.stop = forced.stop, full.path = full.path))
}
########################################################
# check the validity of a specific parallelepiped
###
plus.check.eta <- function(old.eta, new.eta, etatil, xi, z, v, eps=1e-15, e) {
### test a possible new eta and xi combination
## input
# old.eta, 1 X p: old parallelepiped indicator; new.eta: the new one to be tried
# etatil, 2 X (|C|+1): label and new indicator for max crossing, with a dummy variable p+1
# xi: the sign of new.tau - old.tau
# z: t(x)y/n; v: 2nd derivative of penalty; eps: numerical zero; e: environment including x, y
## out put: return(xi, s, singular.Q, invalid.cross)
# singular.Q: error for solve(Q, z[a.set])
# invalid.cross: either singular.Q or fail to cross to new.eta
p <- dim(e$x)[2]
a.set <- rep(TRUE,p)
a.set[new.eta==0] <- FALSE
## update e$Sigma, e$var.list and e$a.set.var
## e$Sigma[,e$a.set.var] will equal to t(x)%*%x[,a.set]/n when done
if (e$use.Gram) e$a.set.var <- a.set # e$Sigma is complete
else if ( sum(e$var.list) == 0 ) { # the initial case of e$Sigma = 0
e$var.list <- order(!a.set)[1:sum(a.set)]
e$a.set.var <- 1:sum(a.set)
e$Sigma[, e$a.set.var] <- e$tx %*% e$x[,a.set]
}
else { # incrementally add new columns to e$Sigma if necessary
new.var <- a.set
new.var[e$var.list] <- FALSE # new.var are those in a.set but nor in e$var.list
if (sum(new.var)>0) { # need to add
e$var.list <- c(e$var.list, order(!new.var)[1:sum(new.var)])
if (length(e$var.list) > (dim(e$Sigma)[2])) { # need to make e$Sigma larger
tmp.Sigma <- e$Sigma
e$Sigma <- matrix(0,p,dim(tmp.Sigma)[2]+100)
e$Sigma[,1:(dim(tmp.Sigma)[2])] <- tmp.Sigma
}
e$Sigma[,(length(e$var.list) - sum(new.var)+1):length(e$var.list)] <- e$tx %*% e$x[,new.var]
}
e$a.set.var <- (1:length(e$var.list))[a.set[e$var.list]] # the equivalent labels of a.set
e$a.set.var <- e$a.set.var[order(e$var.list[ e$a.set.var])]
}
## calculate for return(xi, singular.Q, s, invalid.cross), xi is given
if (sum(a.set)==1) Q <- e$Sigma[a.set,e$a.set.var] - v[ abs(new.eta[a.set]) ]
else Q <- e$Sigma[a.set,e$a.set.var] - diag( v[ abs(new.eta[a.set]) ] )
singular.Q <- FALSE
# LINPACK=T generated an incorrect calculation, due to the use of single precision
# s1 <- try(solve(Q, z[a.set], LINPACK = (length(Q)>100)), silent=T)
# s1 <- try(solve(Q, z[a.set], symmetry=T), silent=T) # option has no effect
s1 <- try(solve(Q, z[a.set]), silent=TRUE)
if (inherits(s1,"try-error")==TRUE) {
s1 <- rep(0, sum(a.set))
singular.Q <- TRUE
}
s <- rep(0,p)
if (singular.Q==FALSE) {
s[ a.set ] <- s1
## check the new s
etatil2 <- c(old.eta,0) # etatil in full length
etatil2[ etatil[1,] ] <- etatil[2,]
etatil2 <- etatil2[1:p] # remove dummy
# e$g.prime <- z - as.matrix((e$Sigma[,e$a.set.var]))%*%s1 # Eq. (2.13)
if (length(s1)==1) e$g.prime <- z - (e$Sigma[,e$a.set.var])*s1 # Eq. (2.13)
else e$g.prime <- z - (e$Sigma[,e$a.set.var])%*%s1
flag <- rep(TRUE,p) # T means violation
flag[old.eta==etatil2] <- FALSE # noncritical j are fine
# check according to (2.15) with eps = numerical zero
flag[ (old.eta != new.eta) & (new.eta != 0) & (xi*(new.eta-old.eta)*s > -eps) ] <- FALSE
flag[ (old.eta == new.eta) & (new.eta != 0) & (etatil2 != old.eta) & (xi*(etatil2-old.eta)*s < eps) ] <- FALSE
flag[ (old.eta != 0) & (new.eta == 0) & (xi*old.eta*e$g.prime < eps)] <- FALSE
flag[ (old.eta == 0) & (new.eta == 0) & (etatil2 != 0) & (xi*etatil2*e$g.prime < eps)] <- FALSE
invalid.cross <- sum(flag) > 0 # T here iff the new eta and xi pair is invalid and singular.Q==FALSE
}
else invalid.cross <- TRUE # invalid.cross is T if singular.Q == T
return(list(xi=xi, singular.Q=singular.Q, s=s, invalid.cross=invalid.cross))
}
########################################################
# find the segment and weight in the path for the hitting point
###
plus.hit.points <- function (lam, lam.path) {
if (is.null(lam) || !is.numeric(lam) || is.nan(sum(lam)))
{
lam <- sort(lam.path, decreasing = TRUE)
warning("Warning message: ")
warning("invalid lam; lam.path used\n")
}
if (max(lam)< min(lam.path))
{
lam <- sort(lam.path, decreasing = TRUE)
warning("Warning message: ")
warning("lam not reached by the plus path; lam.path used\n")
}
lam <- mapply(function(x) min(x,max(lam.path)),lam)
if (min(lam) < min(lam.path))
{
lam <- lam[lam>=min(lam.path)]
warning("Warning message: ")
warning("some lam not reached by the plus path and dropped\n")
}
lam.lesser.path <- mapply(function(x) x<=lam.path,lam)
lam.greater.path <- mapply(function(x) x>=lam.path,lam)
hit <- (lam.lesser.path[1:(length(lam.path)-1),] &
lam.greater.path[2:length(lam.path),])
k1 <- apply(as.matrix(hit),2,function(bool.vec) which(bool.vec)[1])
k2 <- k1+1
w1 <- (lam - lam.path[k2])/(lam.path[k1]-lam.path[k2])
return(list("k1"=k1,"k2"=k2,"w1"=w1,"total.hits"=length(lam)))
}
########################################################
# create a new eta indicating a new parallelepiped
###
plus.new.eta <- function(eta, etatil, cross,p) {
eta <- c(eta,0) # put the dummy in the zero interval
eta[ etatil[1,cross] ] <- etatil[2, cross]
eta <- eta[1:p] # remove dummy
return(eta)
}
########################################################
# check if a newly created eta is indeed new
###
plus.not.in.loop <- function( k, new.eta,e) {
## input: eta, (k+1)Xp; new.eta, 1Xp
## output is T iff new.eta is not any of eta[i,], i\le k
new.is.new <- TRUE
for (i in 1:k)
if( all(e$eta[i,]==new.eta)) return(FALSE)
return(new.is.new)
}
########################################################
# compute the knots and second derivative of a penalty
###
plus.penalty <- function(m,a=0) {
### lasso for m=1, MCP for m=2, SCAD for m=3; default gamma = a = 3.7
### vectors v and t with dummy are used for lasso to maintain vector attribute
m <- round(m)
if (m < 1) m <- 2
if (m > 3) m <- 2
if (m == 2) {
if (a==0) a <- 3.7
v <- c(1/a,0)
t <- c(0,a)
}
if (m == 1) {
a <- 0
v <- c(0,0)
t <- c(0,0)
}
if (m==3) {
if (a==0) a <- 3.7
v <- c(0,1/(a-1),0)
t <- c(0,1,a)
}
return(list(m=m, gamma=a, v=v,t=t))
}
########################################################
# stablized reciprocals
###
plus.recip <- function(Q, eps=1e-15){
### compute reciprocal of reals
Q[ (Q>=0) & (Q <eps) ] <- eps
Q[ (Q<0) & (Q> - eps) ] <- - eps
return(1/Q)
}
########################################################
# generate knots for both sides of zero
###
plus.knots <- function(t=0,m=0) {
### compute the function t() as in (2.6) without infinity and - infinity
### note that t_1=t(1)=t(0) in (2.6)
if (m==0) { # default MC penalty
m <- 2
t <- c(0,3)
}
t2 <- rep(0,2*m)
t2[1:m] <- - t[m:1]
t2[(m+1):(2*m)] <- t[1:m]
return(c(t2,0))
}
##### The above are the plus function
lmres <- function(Y, X, delta) {
X <- as.matrix(X)
o <- order(Y)
Y1 <- Y[o]
X1 <- X[o, ]
delta1 <- delta[o]
Wn1 <- delta1
n1 <- length(Y)
Wn1[1] <- delta1[1]/n1
tt <- 1
for (i in 2:n1) {
j <- i - 1
tt <- tt * ((n1 - j)/(n1 - j + 1))^delta1[j]
Wn1[i] <- delta1[i]/(n1 - i + 1) * tt
}
W <- sqrt(Wn1)
Y2 <- W * Y1
X2 <- W * X1
X2 <- as.matrix(X2)
id <- which(Wn1 > 0)
if (length(id) >= dim(X)[2]) {
Y3 <- Y2[id]
X3 <- X2[id, ]
W3 <- W[id]
X3 <- as.matrix(X3)
lmr <- lm(Y3 ~ W3 + X3 - 1)
return(list(length(Y) * sum((lmr$residuals)^2), lmr$coefficients,
sqrt(length(Y)) * Y3, sqrt(length(Y)) * X3, sqrt(length(Y)) *
W3))
}
if (length(id) < dim(X)[2])
return(NULL)
}
cumcp <- function(Y, X, delta) {
X <- as.matrix(X)
p <- dim(X)[2]
n <- dim(X)[1]
cp <- rep(0, n)
delta1 <- which(delta > 0)
n1 <- length(which(delta > 0))
if (n1 > 2 * p) {
for (i in delta1[p + 1]:delta1[n1 - p]){
cpi <- try(lmres(Y[1:i],
X[1:i, ],
delta[1:i])[[1]] + lmres(Y[(i + 1):n], X[(i + 1):n,],
delta[(i + 1):n])[[1]],silent = TRUE)
if(!is(cpi,"try-error")){
cp[i] <- cpi
}
}
cpm <- min(cp[which(cp > 0)])
return(which(cp == cpm))
} else return(0)
}
#' TSMCP: Two stage multiple change points detection for AFT model.
#'
#' This function first formulates the threshold problem as a group model selection problem
#' so that a concave 2-norm group selection method can be applied using the `grpreg`
#' package in R, and then finalizes it via a refining method.
#'
#' @param Y the censored logarithm of the failure time.
#' @param X the design matrix without the intercept.
#' @param delta the censoring indicator.
#' @param c the length of each segment in the splitting stage, defined as
#' \code{ceiling(c * sqrt(length(Y)))}.
#' @param penalty Penalty type (default is "scad").
#'
#' @return An object with the following components:
#' \describe{
#' \item{cp}{the change points.}
#' \item{coef}{the estimated coefficients.}
#' \item{sigma}{the variance of the error.}
#' \item{residuals}{the residuals.}
#' \item{Yn}{weighted Y by Kaplan-Meier weight.}
#' \item{Xn}{weighted Xn by Kaplan-Meier weight.}
#' }
#'
#' @export
#'
#' @references
#' Li, Jialiang, and Baisuo Jin. 2018. “Multi-Threshold Accelerated Failure Time Model.”
#' The Annals of Statistics 46 (6A): 2657–82.
#'
#' @seealso
#' \code{\link{grpreg}}
#'
#' @examples
#' \donttest{
#' library(grpreg)
#' # Generate simulated data with 500 samples and normal error distribution
#' dataset <- MTAFT_simdata(n = 500, err = "normal")
#' Y <- dataset[, 1]
#' delta <- dataset[, 2]
#' Tq <- dataset[, 3]
#' X <- dataset[, -c(1:3)]
#' n1 = sum(delta)
#' c=seq(0.5,1.5,0.1)
#' m=ceiling(c*sqrt(n1))
#' bicy= rep(NA,length(c))
#' tsmc=NULL
#' p = ncol(X)
#' for(i in 1:length(c)){
#' tsm=try(TSMCP(Y,X,delta,c[i],penalty = "scad"),silent=TRUE)
#' if(is(tsm,"try-error")) next()
#' bicy[i]=log(n)*((length(tsm[[1]])+1)*(p+1))+n*log(tsm[[3]])
#' tsmc[[i]]=tsm
#' }
#'
#' if((any(!is.na(bicy)))){
#' tsmcp=tsmc[[which(bicy==min(bicy))[1]]]
#' thre.LJ = Tq[tsmcp[[1]]]
#' thre.num.Lj = length(thre.LJ)
#' thre.LJ
#' thre.num.Lj
#' }
#' }
TSMCP <- function(Y, X, delta, c,penalty = "scad") {
n <- length(Y)
X <- as.matrix(X)
p <- dim(X)[2]
id1 <- which(delta == 1)
n1 <- length(id1)
m <- ceiling(c * sqrt(n1))
q <- floor(n1/m)
x <- NULL
y <- NULL
de <- NULL
tt <- lmres(Y[1:id1[n1 - (q - 1) * m]], X[1:id1[n1 - (q - 1) * m],], delta[1:id1[n1 - (q - 1) * m]])
y[[1]] <- tt[[3]]
x[[1]] <- cbind(tt[[5]], tt[[4]])
for (i in 2:(q - 1)) {
tt <- lmres(Y[id1[n1 - (q - i + 1) * m + 1]:id1[n1 - (q - i) *
m]], X[id1[n1 - (q - i + 1) * m + 1]:id1[n1 - (q - i) * m],
], delta[id1[n1 - (q - i + 1) * m + 1]:id1[n1 - (q - i) * m]])
y[[i]] <- tt[[3]]
x[[i]] <- cbind(tt[[5]], tt[[4]])
}
i <- q
tt <- lmres(Y[id1[n1 - (q - i + 1) * m + 1]:n], X[id1[n1 - (q - i +
1) * m + 1]:n, ], delta[id1[n1 - (q - i + 1) * m + 1]:n])
y[[i]] <- tt[[3]]
x[[i]] <- cbind(tt[[5]], tt[[4]])
K_temp <- matrix(0, nrow = q, ncol = q, byrow = TRUE)
for (i in 1:q) K_temp[i, 1:i] <- rep(1, i)
X_temp1 <- lapply(1:length(x), function(j, mat, list) kronecker(K_temp[j,
, drop = FALSE], x[[j]]),
mat = K_temp, list = x)
Xn <- do.call("rbind", X_temp1)
Yn <- NULL
for (i in 1:q) {
Yn <- c(Yn, y[[i]])
}
group <- kronecker(c(1:q), rep(1, p + 1))
fit <- grpreg(Xn, Yn, group, penalty = "grSCAD")
mcp.coef <- grpreg::select(fit, "BIC")$beta[-1]
mcp.coef.s <- sum(abs(mcp.coef))
mcp.coef.v.m <- abs(matrix(c(mcp.coef), q, (p + 1), byrow = T))
mcp.coef.m <- c(apply(mcp.coef.v.m, 1, max))
mcp.cp <- which(mcp.coef.m != 0)
if (q > 10) {
if (length(mcp.cp) > 1) {
for (i in 2:length(mcp.cp)) {
if (mcp.cp[i] - mcp.cp[i - 1] == 1)
mcp.cp[i] <- 0
}
}
mcp.cp1 <- mcp.cp[mcp.cp > 1 & mcp.cp < q]
d1 <- length(mcp.cp1)
if (d1 == 0) {
mcpcss.cp <- NULL
adcpcss.cp <- mcpcss.cp
}
if (d1 >= 1) {
mcpcss.cp <- NULL
for (i in 1:d1) {
t <- mcp.cp1[i]
id <- c(id1[n1 - (q - t + 2) * m + 1]:id1[n1 - (q - t -
1) * m])
cp.t <- cumcp(Y[id], X[id, ], delta[id])[[1]]
if (cp.t > 0)
mcpcss.cp <- c(mcpcss.cp, cp.t + id1[n1 - (q - t + 2) *
m + 1] - 1)
}
id2 <- c(mcpcss.cp, n)
tt1 <- NULL
tt1 <- c(tt1, table(delta[1:mcpcss.cp[1]])[2])
for (i in 2:length(id2)) {
tt1 <- c(tt1, table(delta[(id2[i - 1] + 1):id2[i]])[2])
}
tt2 <- which(tt1 < p + 2)
tt2[which(tt2 == length(id2))] <- length(id2) - 1
if (length(tt2) > 0) {
mcpcss.cp <- mcpcss.cp[-tt2]
}
tt <- which(abs(diff(mcpcss.cp)) < 2 * m)
if (length(tt) > 0)
mcpcss.cp <- mcpcss.cp[-tt]
}
}
if (q <= 10) {
mcp.cp1 <- mcp.cp[mcp.cp > 1 & mcp.cp < q]
d1 <- length(mcp.cp1)
if (d1 == 0) {
mcpcss.cp <- NULL
adcpcss.cp <- mcpcss.cp
}
if (d1 >= 1) {
mcpcss.cp <- NULL
for (i in 1:d1) {
t <- mcp.cp1[i]
if (t == 2) {
id <- c(1:id1[n1 - (q - t - 1) * m])
cp.t <- cumcp(Y[id], X[id, ], delta[id])[[1]]
if (cp.t > 0)
mcpcss.cp <- c(mcpcss.cp, cp.t)
}
if (t == q) {
id <- c(id1[n1 - (q - t + 2) * m + 1]:n)
cp.t <- cumcp(Y[id], X[id, ], delta[id])[[1]]
if (cp.t > 0)
mcpcss.cp <- c(mcpcss.cp, cp.t + id1[n1 - (q - t +
2) * m + 1] - 1)
}
if (t < q & t > 2) {
id <- c(id1[n1 - (q - t + 2) * m + 1]:id1[n1 - (q - t -
1) * m])
cp.t <- cumcp(Y[id], X[id, ], delta[id])[[1]]
if (cp.t > 0)
mcpcss.cp <- c(mcpcss.cp, cp.t + id1[n1 - (q - t +
2) * m + 1] - 1)
}
}
id2 <- c(mcpcss.cp, n)
tt1 <- NULL
tt1 <- c(tt1, table(delta[1:mcpcss.cp[1]])[2])
for (i in 2:length(id2)) {
tt1 <- c(tt1, table(delta[(id2[i - 1] + 1):id2[i]])[2])
}
tt2 <- which(tt1 < p + 2)
tt2[which(tt2 == length(id2))] <- length(id2) - 1
if (length(tt2) > 0) {
mcpcss.cp <- mcpcss.cp[-tt2]
}
tt <- which(abs(diff(mcpcss.cp)) < 2 * m)
if (length(tt) > 0)
mcpcss.cp <- mcpcss.cp[-tt]
}
}
if (length(mcpcss.cp) >= 1 && mcpcss.cp[1] != 0) {
id2 <- c(mcpcss.cp, n)
x <- NULL
y <- NULL
de <- NULL
tt <- lmres(Y[1:mcpcss.cp[1]], X[1:mcpcss.cp[1], ], delta[1:mcpcss.cp[1]])
y[[1]] <- tt[[3]]
x[[1]] <- cbind(tt[[5]], tt[[4]])
for (i in 2:length(id2)) {
tt <- lmres(Y[(id2[i - 1] + 1):id2[i]], X[(id2[i - 1] + 1):id2[i],
], delta[(id2[i - 1] + 1):id2[i]])
y[[i]] <- tt[[3]]
x[[i]] <- cbind(tt[[5]], tt[[4]])
}
K_temp <- matrix(0, nrow = length(id2), ncol = length(id2), byrow = TRUE)
for (i in 1:length(id2)) K_temp[i, 1:i] <- rep(1, i)
X_temp1 <- lapply(1:length(x), function(j, mat, list) kronecker(K_temp[j,
, drop = FALSE], x[[j]]), mat = K_temp, list = x)
Xn <- do.call("rbind", X_temp1)
Yn <- NULL
for (i in 1:length(id2)) {
Yn <- c(Yn, y[[i]])
}
#//...................................................................................................
object <- plus(Xn, Yn, method = penalty, gamma = 2.4, intercept = F,normalize = F, eps = 1e-30)
# step 1 estimate coef using BIC.
bic <- log(dim(Xn)[1]) * object$dim + dim(Xn)[1] * log(as.vector((1 -object$r.square) * sum(Yn^2))/length(Yn))
step.bic <- which.min(bic)
mcp.coef <- coef(object, lam = object$lam[step.bic])
#//...................................................................................................
#reg.fit <- ncvreg(X = Xn, y = drop(Yn), family = "gaussian",penalty = penalty)
#coef.beta <- reg.fit$beta # extract coefficients at all values of lambda, including the intercept
#reg.df <- apply(abs(coef.beta[-1, , drop = FALSE]) > 1e-10,2, sum)
#BIC <- reg.fit$loss + log(length(Yn)) * reg.df
#lambda.ind <- which.min(BIC)
#mcp.coef <- coef.beta[-1, lambda.ind]
#..................
ep <- (Yn - Xn %*% mcp.coef)
sigmaep <- sum(ep^2)/length(Y)
mcp.coef.v.m <- matrix(c(mcp.coef), length(id2), p + 1, byrow = T)
mcp.coef.m <- c(apply(abs(mcp.coef.v.m), 1, max))
mcp.m.id <- which(mcp.coef.m == 0)
if (length(mcp.m.id) > 0) {
mcp.coef <- c(t(mcp.coef.v.m[-mcp.m.id, ]))
mcpcss.cp <- mcpcss.cp[-c(mcp.m.id - 1)]
}
}
if (mcpcss.cp == 0 || length(mcpcss.cp) == 0) {
id2 <- c(mcpcss.cp, n)
x <- NULL
y <- NULL
de <- NULL
tt <- lmres(Y, X, delta)
Yn <- tt[[3]]
Xn <- cbind(tt[[5]], tt[[4]])
object <- plus(Xn, Yn, method = "scad", gamma = 2.4, intercept = F,
normalize = F, eps = 1e-30)
# step 1 estimate coef using BIC.
bic <- log(dim(Xn)[1]) * object$dim + dim(Xn)[1] * log(as.vector((1 -
object$r.square) * sum(Yn^2))/length(Yn))
step.bic <- which.min(bic)
mcp.coef <- coef(object, lam = object$lam[step.bic])
ep <- (Yn - Xn %*% mcp.coef)
sigmaep <- sum(ep^2)/length(Y)
mcp.coef.v.m <- matrix(c(mcp.coef), length(id2), p + 1, byrow = T)
mcp.coef.m <- c(apply(abs(mcp.coef.v.m), 1, max))
mcp.m.id <- which(mcp.coef.m == 0)
if (length(mcp.m.id) > 0) {
mcp.coef <- c(t(mcp.coef.v.m[-mcp.m.id, ]))
mcpcss.cp <- mcpcss.cp[-c(mcp.m.id - 1)]
}
}
return(list(cp = mcpcss.cp, coef = mcp.coef, sigma = sigmaep, residuals = ep,
Yn = Yn, Xn = Xn))
}
Try the MTAFT package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.