setupLambda <- function(X, y, group, family, penalty, alpha, lambda.min, log.lambda, nlambda, group.multiplier) {
# Fit to unpenalized covariates
n <- length(y)
K <- table(group)
K1 <- if (min(group)==0) cumsum(K) else c(0, cumsum(K))
storage.mode(K1) <- "integer"
if (K1[1]!=0) {
fit <- glm(y~X[, group==0], family=family)
} else {
fit <- glm(y~1, family=family)
}
## Determine lambda.max
if (family=="gaussian") {
r <- fit$residuals
} else {
w <- fit$weights
if (max(w) < 1e-4) stop("Unpenalized portion of model is already saturated; exiting...", call.=FALSE)
r <- residuals(fit, "working")*w
}
if (strtrim(penalty, 2) == "gr") {
zmax <- .Call("maxgrad", X, r, K1, as.double(group.multiplier)) / n
} else {
zmax <- .Call("maxprod", X, r, K1, as.double(group.multiplier)) / n
}
lambda.max <- zmax/alpha
if (log.lambda) { # lambda sequence on log-scale
if (lambda.min==0) {
lambda <- c(exp(seq(log(lambda.max), log(.001*lambda.max), length=nlambda-1)), 0)
} else {
lambda <- exp(seq(log(lambda.max), log(lambda.min*lambda.max), length=nlambda))
}
} else { # lambda sequence on linear-scale
if (lambda.min==0) {
lambda <- c(seq(lambda.max, 0.001*lambda.max, length = nlambda-1), 0)
} else {
lambda <- seq(lambda.max, lambda.min*lambda.max, length = nlambda)
}
}
lambda
}
setupLambda.gBridge <- function(X, y, group, family, alpha, lambda.min, lambda.max, nlambda, gamma, group.multiplier) {
## Fit to unpenalized covariates
n <- length(y)
ind <- which(group!=0)
if (length(ind)!=length(group)) {
fit <- glm(y~X[, group==0], family=family)
} else {
fit <- glm(y~1, family=family)
}
## Guess lambda.max
if (missing(lambda.max)) {
if (family=="gaussian") {
z <- crossprod(X[, ind], fit$residuals) / n
a <- .35
} else {
z <- crossprod(X[, ind], fit$weights * residuals(fit, "working")) / n
a <- .2
}
lambda.max <- max(abs(z)/group.multiplier[group])*a^(1-gamma)/(gamma*alpha)
}
if (lambda.min==0) {
lambda <- c(exp(seq(log(lambda.max), log(.001*lambda.max), len=nlambda-1)), 0)
} else {
lambda <- exp(seq(log(lambda.max), log(lambda.min*lambda.max), len=nlambda))
}
return(rev(lambda))
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.