R/picasso.poisson.R In picasso: Pathwise Calibrated Sparse Shooting Algorithm

Documented in coef.poissonplot.poissonpredict.poissonprint.poisson

```picasso.poisson <- function(X,
Y,
lambda = NULL,
nlambda = NULL,
lambda.min.ratio = NULL,
method = "l1",
gamma = 3,
standardize = TRUE,
intercept = FALSE,
prec = 1e-4,
max.ite = 1e4,
verbose = FALSE)
{
n = nrow(X)
d = ncol(X)
if (!isTRUE(all(Y == floor(Y))) || !isTRUE(all(Y >= 0)))
stop("The response must only contain non-negative integer values for poisson regression")

if (sum(Y) <= 0)
stop("The response vector is an all-zero vector. The problem is ill-conditioned.")

begt = Sys.time()
if (verbose)
cat("Sparse poisson regression. \n")
if (n == 0 || d == 0) {
cat("No data input.\n")
return(NULL)
}
if (method != "l1" && method != "mcp" && method != "scad"){
cat(" Wrong \"method\" input. \n \"method\"
should be one of \"l1\", \"mcp\", \"scad\".\n",
method,"is not supported in this version. \n")
return(NULL)
}

if (standardize){
xx = rep(0.0, n*d)
xm = rep(0.0, d)
xinvc.vec = rep(0.0, d)
str = .C("standardize_design", as.double(X), as.double(xx),
as.double(xm), as.double(xinvc.vec),
as.integer(n), as.integer(d), PACKAGE="picasso")
xx = matrix(unlist(str[2]), nrow=n, ncol=d, byrow=FALSE)
xm = matrix(unlist(str[3]), nrow=1)
xinvc.vec = unlist(str[4])
} else {
xinvc.vec = rep(1,d)
xx = X
}
yy = Y
avr_y = mean(yy)

if (!is.null(lambda))
nlambda = length(lambda)
if (is.null(lambda)){
if (is.null(nlambda))
nlambda = 100

lambda.max = max(abs(crossprod(xx,(yy-avr_y)/n)))

if (is.null(lambda.min.ratio)){
lambda.min = 0.05*lambda.max
} else {
lambda.min = min(lambda.min.ratio*lambda.max, lambda.max)
}

if (lambda.min>=lambda.max)
cat("lambda.min is too small. \n")
lambda = exp(seq(log(lambda.max), log(lambda.min), length = nlambda))
}
if (method == "l1" || method == "mcp" || method == "scad") {
if (method=="l1") {
method.flag = 1
}
method.flag = 3
if (gamma<=2) {
cat("gamma > 2 is required for SCAD. Set to default value 3. \n")
gamma = 3
}
}
if (method=="mcp") {
method.flag = 2
if (gamma<=1) {
cat("gamma > 1 is required for MCP. Set to default value 3. \n")
gamma = 3
}
}

out = poisson_solver(yy, xx, lambda, nlambda, gamma,
n, d, max.ite, prec, intercept = FALSE, verbose,
method.flag)
}

df = rep(0,nlambda)
for (i in 1:nlambda)
df[i] = sum(out\$beta[[i]]!=0)

est = list()
intcpt = matrix(0, nrow=1, ncol=nlambda)
beta1 = matrix(0, nrow=d, ncol=nlambda)

if (standardize==TRUE){
for (k in 1:nlambda){
tmp.beta = out\$beta[[k]]
beta1[,k] = xinvc.vec*tmp.beta
intcpt[k] = -as.numeric(xm[1,]%*%beta1[,k])+out\$intcpt[k]
}
} else {
for (k in 1:nlambda){
beta1[,k]=out\$beta[[k]]
intcpt[k] = out\$intcpt[k]
}
}
runt = Sys.time()-begt
est\$runt = out\$runt
est\$beta = Matrix(beta1)
est\$intercept = intcpt
est\$lambda = lambda
est\$nlambda = nlambda
est\$df = df
est\$method = method

est\$ite =out\$ite
est\$verbose = verbose
est\$runtime = runt
class(est) = "poisson"
return(est)
}

print.poisson <- function(x, ...)
{
cat("\n Poisson regression options summary: \n")
cat(x\$nlambda, " lambdas used:\n")
print(signif(x\$lambda,digits=3))
cat("Regularization =", x\$method, "\n")
cat("Degree of freedom:",min(x\$df),"----->",max(x\$df),"\n")
if(units.difftime(x\$runtime)=="secs") unit="secs"
if(units.difftime(x\$runtime)=="mins") unit="mins"
if(units.difftime(x\$runtime)=="hours") unit="hours"
cat("Runtime:",x\$runtime," ",unit,"\n")
}

plot.poisson <- function(x, ...)
{
matplot(x\$lambda, t(x\$beta), type="l", main="Regularization Path",
xlab="Regularization Parameter", ylab="Coefficient")
}

coef.poisson <- function(object, lambda.idx = c(1:3), beta.idx = c(1:3), ...)
{
lambda.n = length(lambda.idx)
beta.n = length(beta.idx)
cat("\n Values of estimated coefficients: \n")
cat(" index     ")
for (i in 1:lambda.n){
cat("",formatC(lambda.idx[i], digits=5, width=10),"")
}
cat("\n")
cat(" lambda    ")
for (i in 1:lambda.n){
cat("",formatC(object\$lambda[lambda.idx[i]], digits=4, width=10),"")
}
cat("\n")
cat(" intercept ")
for (i in 1:lambda.n){
cat("",formatC(object\$intercept[i], digits=4, width=10),"")
}
cat("\n")
for (i in 1:beta.n){
cat(" beta",formatC(beta.idx[i], digits=5, width=-5))
for (j in 1:lambda.n){
cat("",formatC(object\$beta[beta.idx[i],lambda.idx[j]],digits=4,width=10),"")
}
cat("\n")
}
}

predict.poisson <- function(object, newdata, lambda.idx = c(1:3), p.pred.idx = c(1:5), ...)
{
pred.n = nrow(newdata)
lambda.n = length(lambda.idx)
p.pred.n = length(p.pred.idx)
intcpt = matrix(rep(object\$intercept[lambda.idx],pred.n), nrow=pred.n,
ncol=lambda.n, byrow=T)
res = newdata%*%object\$beta[,lambda.idx] + intcpt

p.pred = res
cat("\n Values of predicted Poisson parameter: \n")
cat("   index   ")
for (i in 1:lambda.n){
cat("",formatC(lambda.idx[i], digits=5, width=10),"")
}
cat("\n")
cat("   lambda  ")
for (i in 1:lambda.n){
cat("",formatC(object\$lambda[lambda.idx[i]], digits=4, width=10),"")
}
cat("\n")
for (i in 1:p.pred.n){
cat("    Y",formatC(p.pred.idx[i], digits=5, width=-5))
for (j in 1:lambda.n){
cat("",formatC(p.pred[p.pred.idx[i],j], digits=4, width=10),"")
}
cat("\n")
}
return(p.pred)
}
```

Try the picasso package in your browser

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

picasso documentation built on May 1, 2019, 6:28 p.m.