| Title: | Pretraining the Lasso |
|---|---|
| Description: | Implements lasso pretraining for generalized linear models for data with grouped observations; data without grouped observations, but with multinomial responses; data with multiple Gaussian responses; and data with repeated measurements over time. The method fits a lasso model to an initial dataset and then uses its support and coefficient values to guide model fitting for a related task, which can improve performance by leveraging larger datasets and preserve interpretability by identifying features which are generally shared, and which are task specific. See Craig et al. (2026) <doi:10.1093/jrsssb/qkaf050> for details. Model fitting in this package is done with cv.glmnet() from the 'glmnet' package, and syntax closely follows that of 'glmnet' (Friedman, Hastie and Tibshirani (2010) <doi:10.18637/jss.v033.i01>). |
| Authors: | Erin Craig [aut, cre, cph] (ORCID: <https://orcid.org/0000-0002-4893-6671>), Rob Tibshirani [aut, cph] |
| Maintainer: | Erin Craig <[email protected]> |
| License: | GPL-2 |
| Version: | 1.0.2 |
| Built: | 2026-10-28 11:33:52 UTC |
| Source: | https://github.com/erincr/ptlasso |
No required arguments; used primarily for documentation. Simply calls makedata with a reasonable set of features.
binomial.example.data( k = 3, class.sizes = rep(100, k), n = sum(class.sizes), scommon = 5, sindiv = rep(5, k), p = 2 * (sum(sindiv) + scommon), beta.common = list(c(-0.5, 0.5, 0.3, -0.9, 0.1), c(-0.3, 0.9, 0.1, -0.1, 0.2), c(0.1, 0.2, -0.1, 0.2, 0.3)), beta.indiv = lapply(1:k, function(i) 0.9 * beta.common[[i]]), intercepts = rep(0, k), sigma = NULL )binomial.example.data( k = 3, class.sizes = rep(100, k), n = sum(class.sizes), scommon = 5, sindiv = rep(5, k), p = 2 * (sum(sindiv) + scommon), beta.common = list(c(-0.5, 0.5, 0.3, -0.9, 0.1), c(-0.3, 0.9, 0.1, -0.1, 0.2), c(0.1, 0.2, -0.1, 0.2, 0.3)), beta.indiv = lapply(1:k, function(i) 0.9 * beta.common[[i]]), intercepts = rep(0, k), sigma = NULL )
k |
Default: 3. |
class.sizes |
Default: rep(100, k). |
n |
Default: sum(class.sizes). |
scommon |
Default: 5. |
sindiv |
Default: rep(5, k). |
p |
Default: 2*(sum(sindiv) + scommon). |
beta.common |
Default: list(c(-.5, .5, .3, -.9, .1), c(-.3, .9, .1, -.1, .2), c(0.1, .2, -.1, .2, .3)). |
beta.indiv |
Default: lapply(1:k, function(i) 0.9 * beta.common[[i]]). |
intercepts |
Default: rep(0,k). |
sigma |
Default: NULL. |
A list for data with 5 groups and a binomial outcome, n=300 and p=40:
x |
Simulated features, size n x p. |
y |
Outcomes y, length n. |
groups |
Vector of length n, indicating which observations belong to which group. |
snr |
Gaussian outcome only: signal to noise ratio. |
mu |
Gaussian outcome only: the value of y before noise is added. |
Erin Craig and Rob Tibshirani
Maintainer: Erin Craig <[email protected]>
cv.ptLasso, ptLasso.
out = binomial.example.data() x = out$x; y=out$y; groups = out$groupout = binomial.example.data() x = out$x; y=out$y; groups = out$group
Get the coefficients from a fitted cv.ptLasso model.
## S3 method for class 'cv.ptLasso' coef( object, model = c("all", "individual", "overall", "pretrain"), alpha = NULL, ... )## S3 method for class 'cv.ptLasso' coef( object, model = c("all", "individual", "overall", "pretrain"), alpha = NULL, ... )
object |
fitted |
model |
string indicating which coefficients to retrieve. Must be one of "all", "individual", "overall" or "pretrain". |
alpha |
value between 0 and 1, indicating which alpha to use. If |
... |
other arguments to be passed to the |
Model coefficients. If model = "overall", this returns the output of coef for the overall model. If model is "individual", this returns a list containing the results of coef for each group-specific (or response-specific) model. If model = "pretrain", this returns the pretrained model coefficients: a list with one entry per value of alpha if alpha = NULL, otherwise the coefficients for the chosen alpha. If model = "all", this returns a list containing all (overall, individual and pretrain) coefficients. (For time series data, there is no overall model, so only individual and pretrain coefficients are returned.) Note that the pretrained coefficients are from the second stage of pretraining only: they do not include the offset from the overall model. The full pretrained linear predictor for a group (or response) is (1 - alpha) times the overall model's linear predictor (at fitoverall.lambda), plus the linear predictor from these coefficients. For time series data, the offset instead comes from the pretrained model for the previous time point.
Erin Craig and Rob Tibshirani
Maintainer: Erin Craig <[email protected]>
cv.ptLasso, ptLasso.
set.seed(1234) out = gaussian.example.data(k=2, class.sizes = c(50, 50)) x = out$x; y=out$y; groups = out$group; cvfit = cv.ptLasso(x, y, groups = groups, family = "gaussian", type.measure = "mse") # Get all model coefficients. names(coef(cvfit)) coef(cvfit, model = "overall") # Overall model only length(coef(cvfit, model = "individual")) # List of coefficients for each group model length(coef(cvfit, model = "pretrain", alpha = .5)) # List of coefficients for each group modelset.seed(1234) out = gaussian.example.data(k=2, class.sizes = c(50, 50)) x = out$x; y=out$y; groups = out$group; cvfit = cv.ptLasso(x, y, groups = groups, family = "gaussian", type.measure = "mse") # Get all model coefficients. names(coef(cvfit)) coef(cvfit, model = "overall") # Overall model only length(coef(cvfit, model = "individual")) # List of coefficients for each group model length(coef(cvfit, model = "pretrain", alpha = .5)) # List of coefficients for each group model
Get the coefficients from a fitted ptLasso model.
## S3 method for class 'ptLasso' coef(object, model = c("all", "individual", "overall", "pretrain"), ...)## S3 method for class 'ptLasso' coef(object, model = c("all", "individual", "overall", "pretrain"), ...)
object |
fitted |
model |
string indicating which coefficients to retrieve. Must be one of "all", "individual", "overall" or "pretrain". |
... |
other arguments to be passed to the |
Model coefficients. If model = "overall", this function returns the output of coef. If model is "individual" or "pretrain", this function returns a list containing the results of coef for each group-specific model. If model = "all", this returns a list containing all (overall, individual and pretrain) coefficients. Note that the pretrained coefficients are from the second stage of pretraining only: they do not include the offset from the overall model. The full pretrained linear predictor for a group (or response) is (1 - alpha) times the overall model's linear predictor (at fitoverall.lambda), plus the linear predictor from these coefficients. For time series data, the offset instead comes from the pretrained model for the previous time point.
Erin Craig and Rob Tibshirani
Maintainer: Erin Craig <[email protected]>
ptLasso.
# Train data out = gaussian.example.data() x = out$x; y=out$y; groups = out$group; fit = ptLasso(x, y, groups = groups, family = "gaussian", type.measure = "mse") # Get all model coefficients. names(coef(fit)) coef(fit, model = "overall") # Overall model only length(coef(fit, model = "individual")) # List of coefficients for each group model length(coef(fit, model = "pretrain")) # List of coefficients for each group model# Train data out = gaussian.example.data() x = out$x; y=out$y; groups = out$group; fit = ptLasso(x, y, groups = groups, family = "gaussian", type.measure = "mse") # Get all model coefficients. names(coef(fit)) coef(fit, model = "overall") # Overall model only length(coef(fit, model = "individual")) # List of coefficients for each group model length(coef(fit, model = "pretrain")) # List of coefficients for each group model
Cross-validation for ptLasso.
cv.ptLasso( x, y, groups = NULL, alphalist = seq(0, 1, length = 11), family = c("default", "gaussian", "multinomial", "binomial", "cox"), use.case = c("inputGroups", "targetGroups", "multiresponse", "timeSeries"), type.measure = c("default", "mse", "mae", "auc", "deviance", "class", "C"), nfolds = 10, foldid = NULL, verbose = FALSE, fitoverall = NULL, fitind = NULL, s = "lambda.min", gamma = "gamma.min", alphahat.choice = "overall", group.intercepts = TRUE, ... )cv.ptLasso( x, y, groups = NULL, alphalist = seq(0, 1, length = 11), family = c("default", "gaussian", "multinomial", "binomial", "cox"), use.case = c("inputGroups", "targetGroups", "multiresponse", "timeSeries"), type.measure = c("default", "mse", "mae", "auc", "deviance", "class", "C"), nfolds = 10, foldid = NULL, verbose = FALSE, fitoverall = NULL, fitind = NULL, s = "lambda.min", gamma = "gamma.min", alphahat.choice = "overall", group.intercepts = TRUE, ... )
x |
|
y |
|
groups |
A vector of length nobs indicating to which group each observation belongs. For data with k groups, groups should be coded as integers 1 through k. Only for 'use.case = "inputGroups"'. |
alphalist |
A vector of values of the pretraining hyperparameter alpha. Defaults to |
family |
Response type as in |
use.case |
The type of grouping observed in the data. Can be one of "inputGroups", "targetGroups", "multiresponse" or "timeSeries". |
type.measure |
Measure computed in |
nfolds |
Number of folds for CV (default is 10). Although |
foldid |
An optional vector of values between 1 and |
verbose |
If |
fitoverall |
An optional cv.glmnet object specifying the overall model. This should have been trained on the full training data, with the argument keep = TRUE. |
fitind |
An optional list of cv.glmnet objects specifying the individual models. These should have been trained on the training data, with the argumnet keep = TRUE. |
s |
The choice of lambda to be used by all models when estimating the CV performance for each choice of alpha. Defaults to "lambda.min". May be "lambda.1se", or a numeric value. (Use caution when supplying a numeric value: the same lambda will be used for all models.) |
gamma |
For use only when |
alphahat.choice |
When choosing alphahat, we may prefer the best performance using all data ( |
group.intercepts |
For 'use.case = "inputGroups"' only. If 'TRUE', fit the overall model with a separate intercept for each group. If 'FALSE', ignore the grouping and fit one overall intercept. Default is 'TRUE'. |
... |
Additional arguments to be passed to the 'cv.glmnet' function. Notable choices include |
This function runs ptLasso once for each requested choice of alpha, and returns the cross validated performance.
An object of class "cv.ptLasso", which is a list with the ingredients of the cross-validation fit.
call |
The call that produced this object. |
alphahat |
Value of |
varying.alphahat |
Vector of values of |
alphalist |
Vector of all alphas that were compared. |
errall |
CV performance for the overall model. |
errpre |
CV performance for the pretrained models (one for each |
errind |
CV performance for the individual model. |
fit |
List of |
fitoverall |
The fitted overall model used for the first stage of pretraining. |
fitoverall.lambda |
The value of |
fitind |
A list containing one individual model for each group. |
use.case |
The use case: "inputGroups" or "targetGroups". |
family |
The family used. |
type.measure |
The type.measure used. |
ptLasso and plot.cv.ptLasso.
# Getting started. First, we simulate data: we need covariates x, response y and group IDs. set.seed(1234) n = 80 p = 20 x = matrix(rnorm(n*p), n, p) y = rnorm(n) groups = sort(rep(1:5, n/5)) xtest = matrix(rnorm(n*p), n, p) ytest = rnorm(n) groupstest = sort(rep(1:5, n/5)) # Model fitting cvfit = cv.ptLasso(x, y, groups = groups, family = "gaussian", nfolds=3, type.measure = "mse") cvfit plot(cvfit) # to see CV performance as a function of alpha predict(cvfit, xtest, groupstest, s="lambda.min") # to predict with held out data predict(cvfit, xtest, groupstest, s="lambda.min", ytest=ytest) # to also measure performance # By default, we used s = "lambda.min" to compute CV performance. # We could instead use s = "lambda.1se": cvfit = cv.ptLasso(x, y, groups = groups, family = "gaussian", nfolds=3, type.measure = "mse", s = "lambda.1se") # We could have used the glmnet option relax = TRUE: cvfit = cv.ptLasso(x, y, groups = groups, family = "gaussian", nfolds=3, type.measure = "mse", relax = TRUE) # And, as we did with lambda, we may want to specify the choice of gamma to compute CV performance: cvfit = cv.ptLasso(x, y, groups = groups, family = "gaussian", nfolds=3, type.measure = "mse", relax = TRUE, gamma = "gamma.1se") # Note that the first stage of pretraining uses "lambda.1se" and "gamma.1se" by default. # This behavior can be modified by specifying overall.lambda and overall.gamma; # see the documentation for ptLasso for more information. # Now, we are ready to simulate slightly more realistic data. # This continuous outcome example has k = 5 groups, where each group has 200 observations. # There are scommon = 10 features shared across all groups, and # sindiv = 10 features unique to each group. # n = 1000 and p = 120 (60 informative features and 60 noise features). # The coefficients of the common features differ across groups (beta.common). # In group 1, these coefficients are rep(1, 10); in group 2 they are rep(2, 10), etc. # Each group has 10 unique features, the coefficients of which are all 3 (beta.indiv). # The intercept in all groups is 0. # The variable sigma = 20 indicates that we add noise to y according to 20 * rnorm(n). set.seed(1234) k=5 class.sizes=rep(200, k) scommon=10; sindiv=rep(10, k) n=sum(class.sizes); p=2*(sum(sindiv) + scommon) beta.common=3*(1:k); beta.indiv=rep(3, k) intercepts=rep(0, k) sigma=20 out = gaussian.example.data(k=k, class.sizes=class.sizes, scommon=scommon, sindiv=sindiv, n=n, p=p, beta.common=beta.common, beta.indiv=beta.indiv, intercepts=intercepts, sigma=20) x = out$x; y=out$y; groups = out$group outtest = gaussian.example.data(k=k, class.sizes=class.sizes, scommon=scommon, sindiv=sindiv, n=n, p=p, beta.common=beta.common, beta.indiv=beta.indiv, intercepts=intercepts, sigma=20) xtest=outtest$x; ytest=outtest$y; groupstest=outtest$groups cvfit = cv.ptLasso(x, y, groups = groups, family = "gaussian", type.measure = "mse") cvfit # plot(cvfit) # to see CV performance as a function of alpha predict(cvfit, xtest, groupstest, ytest=ytest, s="lambda.min") # Now, we repeat with a binomial outcome. # This example has k = 3 groups, where each group has 100 observations. # There are scommon = 5 features shared across all groups, and # sindiv = 5 features unique to each group. # n = 300 and p = 40 (20 informative features and 20 noise features). # The coefficients of the common features differ across groups (beta.common), # as do the coefficients specific to each group (beta.indiv). set.seed(1234) k=3 class.sizes=rep(100, k) scommon=5; sindiv=rep(5, k) n=sum(class.sizes); p=2*(sum(sindiv) + scommon) beta.common=list(c(-.5, .5, .3, -.9, .1), c(-.3, .9, .1, -.1, .2), c(0.1, .2, -.1, .2, .3)) beta.indiv = lapply(1:k, function(i) 0.9 * beta.common[[i]]) out = binomial.example.data(k=k, class.sizes=class.sizes, scommon=scommon, sindiv=sindiv, n=n, p=p, beta.common=beta.common, beta.indiv=beta.indiv) x = out$x; y=out$y; groups = out$group outtest = binomial.example.data(k=k, class.sizes=class.sizes, scommon=scommon, sindiv=sindiv, n=n, p=p, beta.common=beta.common, beta.indiv=beta.indiv) xtest=outtest$x; ytest=outtest$y; groupstest=outtest$groups cvfit = cv.ptLasso(x, y, groups = groups, family = "binomial", type.measure = "auc", nfolds=3, verbose=TRUE, alpha = c(0, .5, 1), alphahat.choice="mean") cvfit # plot(cvfit) # to see CV performance as a function of alpha predict(cvfit, xtest, groupstest, ytest=ytest, s="lambda.1se") ## Not run: ### Model fitting with parallel = TRUE require(doMC) registerDoMC(cores = 4) cvfit = cv.ptLasso(x, y, groups = groups, family = "binomial", type.measure = "auc", parallel=TRUE) ## End(Not run) # Multiresponse pretraining # Now let's consider the case of a multiresponse outcome. We'll start by simulating data: set.seed(1234) n = 1000; ntrain = 500; p = 500 sigma = 2 x = matrix(rnorm(n*p), n, p) beta1 = c(rep(1, 5), rep(0.5, 5), rep(0, p - 10)) beta2 = c(rep(1, 5), rep(0, 5), rep(0.5, 5), rep(0, p - 15)) mu = cbind(x %*% beta1, x %*% beta2) y = cbind(mu[, 1] + sigma * rnorm(n), mu[, 2] + sigma * rnorm(n)) cat("SNR for the two tasks:", round(diag(var(mu)/var(y-mu)), 2), fill=TRUE) cat("Correlation between two tasks:", cor(y[, 1], y[, 2]), fill=TRUE) xtest = x[-(1:ntrain), ] ytest = y[-(1:ntrain), ] x = x[1:ntrain, ] y = y[1:ntrain, ] # Now, we can fit a ptLasso model: fit = cv.ptLasso(x, y, type.measure = "mse", use.case = "multiresponse") plot(fit) # to see the cv curve. predict(fit, xtest) # to predict with new data predict(fit, xtest, ytest=ytest) # if ytest is included, we also measure performance # By default, we used s = "lambda.min" to compute CV performance. # We could instead use s = "lambda.1se": cvfit = cv.ptLasso(x, y, type.measure = "mse", s = "lambda.1se", use.case = "multiresponse") # We could also use the glmnet option relax = TRUE: cvfit = cv.ptLasso(x, y, type.measure = "mse", relax = TRUE, use.case = "multiresponse") # And, as we did with lambda, we may want to specify the choice of gamma to compute CV performance: cvfit = cv.ptLasso(x, y, type.measure = "mse", relax = TRUE, gamma = "gamma.1se", use.case = "multiresponse") # Time series pretraining # Now suppose we have time series data with a binomial outcome measured at 3 different time points. set.seed(1234) n = 600; ntrain = 300; p = 50 x = matrix(rnorm(n*p), n, p) beta1 = c(rep(0.25, 10), rep(0, p-10)) beta2 = beta1 + c(rep(0.1, 10), runif(5, min = -0.25, max = 0), rep(0, p-15)) beta3 = beta1 + c(rep(0.2, 10), runif(5, min = -0.25, max = 0), runif(5, min = 0, max = 0.1), rep(0, p-20)) y1 = rbinom(n, 1, prob = 1/(1 + exp(-x %*% beta1))) y2 = rbinom(n, 1, prob = 1/(1 + exp(-x %*% beta2))) y3 = rbinom(n, 1, prob = 1/(1 + exp(-x %*% beta3))) y = cbind(y1, y2, y3) xtest = x[-(1:ntrain), ] ytest = y[-(1:ntrain), ] x = x[1:ntrain, ] y = y[1:ntrain, ] cvfit = cv.ptLasso(x, y, use.case="timeSeries", family="binomial", type.measure="auc") plot(cvfit, plot.alphahat = TRUE) predict(cvfit, xtest, ytest=ytest) # The glmnet option relax = TRUE: cvfit = cv.ptLasso(x, y, type.measure = "auc", family = "binomial", relax = TRUE, use.case = "timeSeries")# Getting started. First, we simulate data: we need covariates x, response y and group IDs. set.seed(1234) n = 80 p = 20 x = matrix(rnorm(n*p), n, p) y = rnorm(n) groups = sort(rep(1:5, n/5)) xtest = matrix(rnorm(n*p), n, p) ytest = rnorm(n) groupstest = sort(rep(1:5, n/5)) # Model fitting cvfit = cv.ptLasso(x, y, groups = groups, family = "gaussian", nfolds=3, type.measure = "mse") cvfit plot(cvfit) # to see CV performance as a function of alpha predict(cvfit, xtest, groupstest, s="lambda.min") # to predict with held out data predict(cvfit, xtest, groupstest, s="lambda.min", ytest=ytest) # to also measure performance # By default, we used s = "lambda.min" to compute CV performance. # We could instead use s = "lambda.1se": cvfit = cv.ptLasso(x, y, groups = groups, family = "gaussian", nfolds=3, type.measure = "mse", s = "lambda.1se") # We could have used the glmnet option relax = TRUE: cvfit = cv.ptLasso(x, y, groups = groups, family = "gaussian", nfolds=3, type.measure = "mse", relax = TRUE) # And, as we did with lambda, we may want to specify the choice of gamma to compute CV performance: cvfit = cv.ptLasso(x, y, groups = groups, family = "gaussian", nfolds=3, type.measure = "mse", relax = TRUE, gamma = "gamma.1se") # Note that the first stage of pretraining uses "lambda.1se" and "gamma.1se" by default. # This behavior can be modified by specifying overall.lambda and overall.gamma; # see the documentation for ptLasso for more information. # Now, we are ready to simulate slightly more realistic data. # This continuous outcome example has k = 5 groups, where each group has 200 observations. # There are scommon = 10 features shared across all groups, and # sindiv = 10 features unique to each group. # n = 1000 and p = 120 (60 informative features and 60 noise features). # The coefficients of the common features differ across groups (beta.common). # In group 1, these coefficients are rep(1, 10); in group 2 they are rep(2, 10), etc. # Each group has 10 unique features, the coefficients of which are all 3 (beta.indiv). # The intercept in all groups is 0. # The variable sigma = 20 indicates that we add noise to y according to 20 * rnorm(n). set.seed(1234) k=5 class.sizes=rep(200, k) scommon=10; sindiv=rep(10, k) n=sum(class.sizes); p=2*(sum(sindiv) + scommon) beta.common=3*(1:k); beta.indiv=rep(3, k) intercepts=rep(0, k) sigma=20 out = gaussian.example.data(k=k, class.sizes=class.sizes, scommon=scommon, sindiv=sindiv, n=n, p=p, beta.common=beta.common, beta.indiv=beta.indiv, intercepts=intercepts, sigma=20) x = out$x; y=out$y; groups = out$group outtest = gaussian.example.data(k=k, class.sizes=class.sizes, scommon=scommon, sindiv=sindiv, n=n, p=p, beta.common=beta.common, beta.indiv=beta.indiv, intercepts=intercepts, sigma=20) xtest=outtest$x; ytest=outtest$y; groupstest=outtest$groups cvfit = cv.ptLasso(x, y, groups = groups, family = "gaussian", type.measure = "mse") cvfit # plot(cvfit) # to see CV performance as a function of alpha predict(cvfit, xtest, groupstest, ytest=ytest, s="lambda.min") # Now, we repeat with a binomial outcome. # This example has k = 3 groups, where each group has 100 observations. # There are scommon = 5 features shared across all groups, and # sindiv = 5 features unique to each group. # n = 300 and p = 40 (20 informative features and 20 noise features). # The coefficients of the common features differ across groups (beta.common), # as do the coefficients specific to each group (beta.indiv). set.seed(1234) k=3 class.sizes=rep(100, k) scommon=5; sindiv=rep(5, k) n=sum(class.sizes); p=2*(sum(sindiv) + scommon) beta.common=list(c(-.5, .5, .3, -.9, .1), c(-.3, .9, .1, -.1, .2), c(0.1, .2, -.1, .2, .3)) beta.indiv = lapply(1:k, function(i) 0.9 * beta.common[[i]]) out = binomial.example.data(k=k, class.sizes=class.sizes, scommon=scommon, sindiv=sindiv, n=n, p=p, beta.common=beta.common, beta.indiv=beta.indiv) x = out$x; y=out$y; groups = out$group outtest = binomial.example.data(k=k, class.sizes=class.sizes, scommon=scommon, sindiv=sindiv, n=n, p=p, beta.common=beta.common, beta.indiv=beta.indiv) xtest=outtest$x; ytest=outtest$y; groupstest=outtest$groups cvfit = cv.ptLasso(x, y, groups = groups, family = "binomial", type.measure = "auc", nfolds=3, verbose=TRUE, alpha = c(0, .5, 1), alphahat.choice="mean") cvfit # plot(cvfit) # to see CV performance as a function of alpha predict(cvfit, xtest, groupstest, ytest=ytest, s="lambda.1se") ## Not run: ### Model fitting with parallel = TRUE require(doMC) registerDoMC(cores = 4) cvfit = cv.ptLasso(x, y, groups = groups, family = "binomial", type.measure = "auc", parallel=TRUE) ## End(Not run) # Multiresponse pretraining # Now let's consider the case of a multiresponse outcome. We'll start by simulating data: set.seed(1234) n = 1000; ntrain = 500; p = 500 sigma = 2 x = matrix(rnorm(n*p), n, p) beta1 = c(rep(1, 5), rep(0.5, 5), rep(0, p - 10)) beta2 = c(rep(1, 5), rep(0, 5), rep(0.5, 5), rep(0, p - 15)) mu = cbind(x %*% beta1, x %*% beta2) y = cbind(mu[, 1] + sigma * rnorm(n), mu[, 2] + sigma * rnorm(n)) cat("SNR for the two tasks:", round(diag(var(mu)/var(y-mu)), 2), fill=TRUE) cat("Correlation between two tasks:", cor(y[, 1], y[, 2]), fill=TRUE) xtest = x[-(1:ntrain), ] ytest = y[-(1:ntrain), ] x = x[1:ntrain, ] y = y[1:ntrain, ] # Now, we can fit a ptLasso model: fit = cv.ptLasso(x, y, type.measure = "mse", use.case = "multiresponse") plot(fit) # to see the cv curve. predict(fit, xtest) # to predict with new data predict(fit, xtest, ytest=ytest) # if ytest is included, we also measure performance # By default, we used s = "lambda.min" to compute CV performance. # We could instead use s = "lambda.1se": cvfit = cv.ptLasso(x, y, type.measure = "mse", s = "lambda.1se", use.case = "multiresponse") # We could also use the glmnet option relax = TRUE: cvfit = cv.ptLasso(x, y, type.measure = "mse", relax = TRUE, use.case = "multiresponse") # And, as we did with lambda, we may want to specify the choice of gamma to compute CV performance: cvfit = cv.ptLasso(x, y, type.measure = "mse", relax = TRUE, gamma = "gamma.1se", use.case = "multiresponse") # Time series pretraining # Now suppose we have time series data with a binomial outcome measured at 3 different time points. set.seed(1234) n = 600; ntrain = 300; p = 50 x = matrix(rnorm(n*p), n, p) beta1 = c(rep(0.25, 10), rep(0, p-10)) beta2 = beta1 + c(rep(0.1, 10), runif(5, min = -0.25, max = 0), rep(0, p-15)) beta3 = beta1 + c(rep(0.2, 10), runif(5, min = -0.25, max = 0), runif(5, min = 0, max = 0.1), rep(0, p-20)) y1 = rbinom(n, 1, prob = 1/(1 + exp(-x %*% beta1))) y2 = rbinom(n, 1, prob = 1/(1 + exp(-x %*% beta2))) y3 = rbinom(n, 1, prob = 1/(1 + exp(-x %*% beta3))) y = cbind(y1, y2, y3) xtest = x[-(1:ntrain), ] ytest = y[-(1:ntrain), ] x = x[1:ntrain, ] y = y[1:ntrain, ] cvfit = cv.ptLasso(x, y, use.case="timeSeries", family="binomial", type.measure="auc") plot(cvfit, plot.alphahat = TRUE) predict(cvfit, xtest, ytest=ytest) # The glmnet option relax = TRUE: cvfit = cv.ptLasso(x, y, type.measure = "auc", family = "binomial", relax = TRUE, use.case = "timeSeries")
No required arguments; used primarily for documentation. Simply calls makedata with a reasonable set of features.
gaussian.example.data( k = 5, class.sizes = rep(200, k), n = sum(class.sizes), scommon = 10, sindiv = rep(10, k), p = 2 * (sum(sindiv) + scommon), beta.common = 3 * (1:k), beta.indiv = rep(3, k), intercepts = rep(0, k), sigma = 20 )gaussian.example.data( k = 5, class.sizes = rep(200, k), n = sum(class.sizes), scommon = 10, sindiv = rep(10, k), p = 2 * (sum(sindiv) + scommon), beta.common = 3 * (1:k), beta.indiv = rep(3, k), intercepts = rep(0, k), sigma = 20 )
k |
Default: 5. |
class.sizes |
Default: rep(200, k). |
n |
Default: sum(class.sizes). |
scommon |
Default: 10. |
sindiv |
Default: rep(10, k). |
p |
Default: 2*(sum(sindiv) + scommon). |
beta.common |
Default: 3*(1:k). |
beta.indiv |
Default: rep(3, k). |
intercepts |
Default: rep(0, k). |
sigma |
Default: 20. |
A list for data with 5 groups and a gaussian outcome, n=1000 and p=120:
x |
Simulated features, size n x p. |
y |
Outcomes y, length n. |
groups |
Vector of length n, indicating which observations belong to which group. |
snr |
Gaussian outcome only: signal to noise ratio. |
mu |
Gaussian outcome only: the value of y before noise is added. |
Erin Craig and Rob Tibshirani
Maintainer: Erin Craig <[email protected]>
cv.ptLasso, ptLasso.
out = gaussian.example.data() x = out$x; y=out$y; groups = out$groupout = gaussian.example.data() x = out$x; y=out$y; groups = out$group
Get the indices of nonzero coefficients from the individual models in a fitted ptLasso or cv.ptLasso object, excluding the intercept.
get.individual.support( fit, s = "lambda.min", gamma = "gamma.min", commonOnly = FALSE, groups = 1:length(fit$fitind) )get.individual.support( fit, s = "lambda.min", gamma = "gamma.min", commonOnly = FALSE, groups = 1:length(fit$fitind) )
fit |
fitted |
s |
the choice of lambda to use. May be "lambda.min", "lambda.1se" or a numeric value. Default is "lambda.min". |
gamma |
for use only when 'relax = TRUE' was specified during training. The choice of 'gamma' to use. May be "gamma.min" or "gamma.1se". Default is "gamma.min". |
commonOnly |
whether to return the features that are chosen by more than half of the group- or response-specific models (TRUE) or the features that are chosen by any of the group-specific models (FALSE). Default is FALSE. |
groups |
which groups or responses to include when computing the support. Default is to include all groups/responses. |
This returns a vector containing the indices of nonzero coefficients (excluding the intercept).
Erin Craig and Rob Tibshirani
Maintainer: Erin Craig <[email protected]>
ptLasso, cv.ptLasso.
# Train data set.seed(1234) out = gaussian.example.data() x = out$x; y=out$y; groups = out$group; fit = ptLasso(x, y, groups = groups, family = "gaussian", type.measure = "mse") get.individual.support(fit) # only return features common to all groups get.individual.support(fit, commonOnly = TRUE) # group 1 only get.individual.support(fit, groups = 1) cvfit = cv.ptLasso(x, y, groups = groups, alphalist = c(0, .5, 1), family = "gaussian", type.measure = "mse") get.individual.support(cvfit) # group 1 only get.individual.support(cvfit, groups = 1)# Train data set.seed(1234) out = gaussian.example.data() x = out$x; y=out$y; groups = out$group; fit = ptLasso(x, y, groups = groups, family = "gaussian", type.measure = "mse") get.individual.support(fit) # only return features common to all groups get.individual.support(fit, commonOnly = TRUE) # group 1 only get.individual.support(fit, groups = 1) cvfit = cv.ptLasso(x, y, groups = groups, alphalist = c(0, .5, 1), family = "gaussian", type.measure = "mse") get.individual.support(cvfit) # group 1 only get.individual.support(cvfit, groups = 1)
Get the indices of nonzero coefficients from the overall model in a fitted ptLasso or cv.ptLasso object, excluding the intercept.
get.overall.support(fit, s = "lambda.min", gamma = "gamma.min")get.overall.support(fit, s = "lambda.min", gamma = "gamma.min")
fit |
fitted |
s |
the choice of lambda to use. May be "lambda.min", "lambda.1se" or a numeric value. Default is "lambda.min". |
gamma |
for use only when 'relax = TRUE' was specified during training. The choice of 'gamma' to use. May be "gamma.min" or "gamma.1se". Default is "gamma.min". |
This returns a vector containing the indices of nonzero coefficients (excluding the intercept).
Erin Craig and Rob Tibshirani
Maintainer: Erin Craig <[email protected]>
ptLasso, cv.ptLasso.
# Train data set.seed(1234) out = gaussian.example.data(k=2, class.sizes = c(50, 50)) x = out$x; y=out$y; groups = out$group; fit = ptLasso(x, y, groups = groups, family = "gaussian", type.measure = "mse") get.overall.support(fit, s="lambda.min") get.overall.support(fit, s="lambda.1se") cvfit = cv.ptLasso(x, y, groups = groups, family = "gaussian", type.measure = "mse") get.overall.support(cvfit, s="lambda.min") get.overall.support(cvfit, s="lambda.1se")# Train data set.seed(1234) out = gaussian.example.data(k=2, class.sizes = c(50, 50)) x = out$x; y=out$y; groups = out$group; fit = ptLasso(x, y, groups = groups, family = "gaussian", type.measure = "mse") get.overall.support(fit, s="lambda.min") get.overall.support(fit, s="lambda.1se") cvfit = cv.ptLasso(x, y, groups = groups, family = "gaussian", type.measure = "mse") get.overall.support(cvfit, s="lambda.min") get.overall.support(cvfit, s="lambda.1se")
Get the indices of nonzero coefficients from the pretrained models in a fitted ptLasso or cv.ptLasso object, excluding the intercept.
get.pretrain.support( fit, s = "lambda.min", gamma = "gamma.min", commonOnly = FALSE, includeOverall = TRUE, groups = 1:length(fit$fitind) )get.pretrain.support( fit, s = "lambda.min", gamma = "gamma.min", commonOnly = FALSE, includeOverall = TRUE, groups = 1:length(fit$fitind) )
fit |
fitted |
s |
the choice of lambda to use. May be "lambda.min", "lambda.1se" or a numeric value. Default is "lambda.min". |
gamma |
for use only when 'relax = TRUE' was specified during training. The choice of 'gamma' to use. May be "gamma.min" or "gamma.1se". Default is "gamma.min". |
commonOnly |
whether to return the features that are chosen by more than half of the group- or response-specific models (TRUE) or the features that are chosen by any of the group-specific models (FALSE). Default is FALSE. |
includeOverall |
whether to return the features that are chosen by the overall model and not the group-specific models (TRUE) or the features that are chosen by the overall model or the group-specific models (FALSE). Default is TRUE. Not used when 'use.case = "timeSeries"'. |
groups |
which groups or responses to include when computing the support. Default is to include all groups/responses. |
If a ptLasso object is supplied, this returns a vector containing the indices of nonzero coefficients (excluding the intercept). If a cv.ptLasso object is supplied, this returns a list of results - one for each value of alpha.
ptLasso, cv.ptLasso.
# Train data set.seed(1234) out = gaussian.example.data(k=2, class.sizes = c(50, 50)) x = out$x; y=out$y; groups = out$group; fit = ptLasso(x, y, groups = groups, family = "gaussian", type.measure = "mse") get.pretrain.support(fit) # only return features common to all groups get.pretrain.support(fit, commonOnly = TRUE) # group 1 only, don't include the overall model support get.pretrain.support(fit, groups = 1, includeOverall = FALSE) # group 1 only, include the overall model support get.pretrain.support(fit, groups = 1, includeOverall = TRUE) cvfit = cv.ptLasso(x, y, groups = groups, family = "gaussian", type.measure = "mse") get.pretrain.support(cvfit) get.pretrain.support(cvfit, groups = 1)# Train data set.seed(1234) out = gaussian.example.data(k=2, class.sizes = c(50, 50)) x = out$x; y=out$y; groups = out$group; fit = ptLasso(x, y, groups = groups, family = "gaussian", type.measure = "mse") get.pretrain.support(fit) # only return features common to all groups get.pretrain.support(fit, commonOnly = TRUE) # group 1 only, don't include the overall model support get.pretrain.support(fit, groups = 1, includeOverall = FALSE) # group 1 only, include the overall model support get.pretrain.support(fit, groups = 1, includeOverall = TRUE) cvfit = cv.ptLasso(x, y, groups = groups, family = "gaussian", type.measure = "mse") get.pretrain.support(cvfit) get.pretrain.support(cvfit, groups = 1)
Simulate input grouped data for testing with ptLasso.
makedata( n, p, k, scommon, sindiv, class.sizes, beta.common, beta.indiv, intercepts = rep(0, k), sigma = 0, outcome = c("gaussian", "binomial", "multinomial"), mult.classes = 3 )makedata( n, p, k, scommon, sindiv, class.sizes, beta.common, beta.indiv, intercepts = rep(0, k), sigma = 0, outcome = c("gaussian", "binomial", "multinomial"), mult.classes = 3 )
n |
Total number of observations to simulate. |
p |
Total number of features to simulate. |
k |
Number of groups. |
scommon |
Number of features shared by all groups. |
sindiv |
Vector of length k. The i^th entry indicates the number of features specific to group i. |
class.sizes |
Vector of length k. The i^th entry indicates the number of observations in group i. |
beta.common |
The coefficients for the common features. This can be a vector of length k, in which case, the i^th entry is the coefficient for all scommon features for group i. This can alternatively be a list of length k (one for each group). Each entry of this list should be a vector of length scommon, containing the coefficients for the scommon features. |
beta.indiv |
The coefficients for the individual features, in the same form as beta.common. |
intercepts |
A vector of length k, indicating the intercept for each group. Default is 0. |
sigma |
Only used for the Gaussian outcome. Should be a number greater than or equal to 0, used to modify the amount of noise added. Default is 0. |
outcome |
May be '"gaussian"', '"binomial"' or '"multinomial"'. |
mult.classes |
Number of classes to simulate for the multinomial setting. |
A list:
x |
Simulated features, size n x p. |
y |
Outcomes y, length n. |
groups |
Vector of length n, indicating which observations belong to which group. |
snr |
Gaussian outcome only: signal to noise ratio. |
mu |
Gaussian outcome only: the value of y before noise is added. |
Erin Craig and Rob Tibshirani
Maintainer: Erin Craig <[email protected]>
cv.ptLasso, ptLasso.
# Data with a binary outcome: k = 3 class.sizes = rep(100, k) n = sum(class.sizes) scommon = 5 sindiv = rep(5, k) p = 2*(sum(sindiv) + scommon) beta.common = lapply(1:k, function(i) c(-.5, .5, .3, -.9, .1)) beta.indiv = lapply(1:k, function(i) 0.9 * beta.common[[i]]) out = makedata(n=n, p=p, k=k, scommon=scommon, sindiv=sindiv, beta.common=beta.common, beta.indiv=beta.indiv, class.sizes=class.sizes, outcome="binomial") x = out$x; y=out$y; groups = out$group# Data with a binary outcome: k = 3 class.sizes = rep(100, k) n = sum(class.sizes) scommon = 5 sindiv = rep(5, k) p = 2*(sum(sindiv) + scommon) beta.common = lapply(1:k, function(i) c(-.5, .5, .3, -.9, .1)) beta.indiv = lapply(1:k, function(i) 0.9 * beta.common[[i]]) out = makedata(n=n, p=p, k=k, scommon=scommon, sindiv=sindiv, beta.common=beta.common, beta.indiv=beta.indiv, class.sizes=class.sizes, outcome="binomial") x = out$x; y=out$y; groups = out$group
Simulate target grouped data for testing with ptLasso.
makedata.targetgroups( n, p, scommon, sindiv, class.sizes, shift.common, shift.indiv )makedata.targetgroups( n, p, scommon, sindiv, class.sizes, shift.common, shift.indiv )
n |
Total number of observations to simulate. |
p |
Total number of features to simulate. |
scommon |
Number of features shared by all groups. |
sindiv |
Vector of length k. The i^th entry indicates the number of features specific to group i. |
class.sizes |
Vector of length k. The i^th entry indicates the number of observations in group i. |
shift.common |
A list of length k (one for each group). Each entry of this list should be a vector of length scommon, containing the shifts for the scommon features. The i^th entry of this list will be added to the first scommon columns of x for observations in group i. |
shift.indiv |
The shifts for the individual features, in the same form as shift.common. |
A list:
x |
Simulated features, size n x p. |
y |
Outcomes y, length n. |
Erin Craig and Rob Tibshirani
Maintainer: Erin Craig <[email protected]>
cv.ptLasso, ptLasso.
k = 5 class.sizes = rep(50, k) n = sum(class.sizes) scommon = 3 sindiv = rep(3, k) p = 3*(sum(sindiv) + scommon) shift.common = lapply(seq(-.1, .1, length.out = k), function(i) rep(i, scommon)) shift.indiv = lapply(1:k, function(i) -shift.common[[i]]) out = makedata.targetgroups(n=n, p=p, scommon=scommon, sindiv=sindiv, class.sizes=class.sizes, shift.common=shift.common, shift.indiv=shift.indiv) x = out$x; y=out$yk = 5 class.sizes = rep(50, k) n = sum(class.sizes) scommon = 3 sindiv = rep(3, k) p = 3*(sum(sindiv) + scommon) shift.common = lapply(seq(-.1, .1, length.out = k), function(i) rep(i, scommon)) shift.indiv = lapply(1:k, function(i) -shift.common[[i]]) out = makedata.targetgroups(n=n, p=p, scommon=scommon, sindiv=sindiv, class.sizes=class.sizes, shift.common=shift.common, shift.indiv=shift.indiv) x = out$x; y=out$y
alpha values used.A plot is produced, and the plot object is returned invisibly.
## S3 method for class 'cv.ptLasso' plot(x, plot.alphahat = FALSE, ...)## S3 method for class 'cv.ptLasso' plot(x, plot.alphahat = FALSE, ...)
x |
Fitted |
plot.alphahat |
If |
... |
Other graphical parameters to plot. |
Invisibly returns a gtable object (from gridExtra::grid.arrange) containing the combined plot. Called mainly for its side effect of drawing the plot.
ptLasso, cv.ptLasso and predict.cv.ptLasso.
set.seed(1234) out = gaussian.example.data(k=2, class.sizes = c(50, 50)) x = out$x; y=out$y; groups = out$group cvfit = cv.ptLasso(x, y, groups = groups, family = "gaussian", type.measure = "mse") plot(cvfit)set.seed(1234) out = gaussian.example.data(k=2, class.sizes = c(50, 50)) x = out$x; y=out$y; groups = out$group cvfit = cv.ptLasso(x, y, groups = groups, family = "gaussian", type.measure = "mse") plot(cvfit)
A plot is produced, and nothing is returned.
## S3 method for class 'ptLasso' plot(x, ...)## S3 method for class 'ptLasso' plot(x, ...)
x |
Fitted |
... |
Other parameters to pass to the plot function. |
No return value (returns NULL invisibly); called for its side effect of drawing the plot.
Erin Craig and Rob Tibshirani
Maintainer: Erin Craig <[email protected]>
ptLasso, cv.ptLasso and predict.cv.ptLasso.
set.seed(1234) out = gaussian.example.data() x = out$x; y=out$y; groups = out$group fit = ptLasso(x, y, groups = groups, alpha = 0.5, family = "gaussian", type.measure = "mse") plot(fit)set.seed(1234) out = gaussian.example.data() x = out$x; y=out$y; groups = out$group fit = ptLasso(x, y, groups = groups, alpha = 0.5, family = "gaussian", type.measure = "mse") plot(fit)
Return predictions and performance measures for a test set.
## S3 method for class 'cv.ptLasso' predict( object, xtest, groupstest = NULL, ytest = NULL, alpha = NULL, alphatype = c("fixed", "varying"), type = c("link", "response", "class"), s = "lambda.min", gamma = "gamma.min", return.link = FALSE, ... )## S3 method for class 'cv.ptLasso' predict( object, xtest, groupstest = NULL, ytest = NULL, alpha = NULL, alphatype = c("fixed", "varying"), type = c("link", "response", "class"), s = "lambda.min", gamma = "gamma.min", return.link = FALSE, ... )
object |
Fitted |
xtest |
Input matrix, matching the form used by |
groupstest |
A vector indicating to which group each observation belongs. Coding should match that used for model training. Will be NULL for target grouped data. |
ytest |
Response variable. Optional. If included, |
alpha |
The chosen alpha to use for prediction. May be a vector containing one value of alpha for each group. If NULL, this will rely on the choice of "alphatype". |
alphatype |
Choice of '"fixed"' or '"varying"'. If '"fixed"', use the alpha that achieved best cross-validated performance. If '"varying"', each group uses the alpha that optimized the group-specific cross-validated performance. |
type |
Type of prediction required. Type '"link"' gives the linear predictors for '"binomial", '"multinomial"' or '"cox"' models; for '"gaussian"' models it gives the fitted values. Type '"response"' gives the fitted probabilities for '"binomial"' or '"multinomial"', and the fitted relative-risk for '"cox"'; for '"gaussian"' type '"response"' is equivalent to type '"link"'. Note that for '"binomial"' models, results are returned only for the class corresponding to the second level of the factor response. Type '"class"' applies only to '"binomial"' or '"multinomial"' models, and produces the class label corresponding to the maximum probability. |
s |
Value of the penalty parameter 'lambda' at which predictions are required. Will use the same lambda for all models; can be a numeric value, '"lambda.min"' or '"lambda.1se"'. Default is '"lambda.min"'. |
gamma |
For use only when 'relax = TRUE' was specified during training. Value of the penalty parameter 'gamma' at which predictions are required. Will use the same gamma for all models; can be a numeric value, '"gamma.min"' or '"gamma.1se"'. Default is '"gamma.min"'. |
return.link |
If |
... |
other arguments to be passed to the |
A list containing the requested predictions. If ytest is included, will also return error measures.
call |
The call that produced this object. |
alpha |
The value(s) of alpha used to generate predictions. |
yhatoverall |
Predictions from the overall model. |
yhatind |
Predictions from the individual models. |
yhatpre |
Predictions from the pretrained models. |
supoverall |
Indices of the features selected by the overall model. |
supind |
Union of the indices of the features selected by the individual models. |
suppre.common |
Features selected in the first stage of pretraining. |
suppre.individual |
Union of the indices of the features selected by the pretrained models, without the features selected in the first stage. |
type.measure |
If |
erroverall |
If |
errind |
If |
errpre |
If |
linkoverall |
If |
linkind |
If |
linkpre |
If |
Erin Craig and Rob Tibshirani
Maintainer:
Erin Craig <[email protected]>
ptLasso, cv.ptLasso and predict.cv.ptLasso.
#### Gaussian example set.seed(1234) out = gaussian.example.data(k=2, class.sizes = c(50, 50)) x = out$x; y=out$y; groups = out$group; outtest = gaussian.example.data(k=2, class.sizes = c(50, 50)) xtest=outtest$x; ytest=outtest$y; groupstest=outtest$groups; # Model fitting # By default, use the single value of alpha that had the best CV performance on the entire dataset: cvfit = cv.ptLasso(x, y, groups = groups, family = "gaussian", type.measure = "mse") pred = predict(cvfit, xtest, groupstest, ytest=ytest, s="lambda.min") pred # For each group, use the value of alpha that had the best CV performance for that group: pred = predict(cvfit, xtest, groupstest, ytest=ytest, s="lambda.min", alphatype = "varying") pred # Specify a single value of alpha and use lambda.1se. pred = predict(cvfit, xtest, groupstest, ytest=ytest, s="lambda.1se", alphatype = "varying", alpha = .3) pred # Specify a vector of choices for alpha: pred = predict(cvfit, xtest, groupstest, ytest=ytest, s="lambda.min", alphatype = "varying", alpha = c(.1, .5)) pred#### Gaussian example set.seed(1234) out = gaussian.example.data(k=2, class.sizes = c(50, 50)) x = out$x; y=out$y; groups = out$group; outtest = gaussian.example.data(k=2, class.sizes = c(50, 50)) xtest=outtest$x; ytest=outtest$y; groupstest=outtest$groups; # Model fitting # By default, use the single value of alpha that had the best CV performance on the entire dataset: cvfit = cv.ptLasso(x, y, groups = groups, family = "gaussian", type.measure = "mse") pred = predict(cvfit, xtest, groupstest, ytest=ytest, s="lambda.min") pred # For each group, use the value of alpha that had the best CV performance for that group: pred = predict(cvfit, xtest, groupstest, ytest=ytest, s="lambda.min", alphatype = "varying") pred # Specify a single value of alpha and use lambda.1se. pred = predict(cvfit, xtest, groupstest, ytest=ytest, s="lambda.1se", alphatype = "varying", alpha = .3) pred # Specify a vector of choices for alpha: pred = predict(cvfit, xtest, groupstest, ytest=ytest, s="lambda.min", alphatype = "varying", alpha = c(.1, .5)) pred
Return predictions and performance measures for a test set.
## S3 method for class 'ptLasso' predict( object, xtest, groupstest = NULL, ytest = NULL, type = c("link", "response", "class"), s = "lambda.min", gamma = "gamma.min", return.link = FALSE, ... )## S3 method for class 'ptLasso' predict( object, xtest, groupstest = NULL, ytest = NULL, type = c("link", "response", "class"), s = "lambda.min", gamma = "gamma.min", return.link = FALSE, ... )
object |
Fitted |
xtest |
Input matrix, matching the form used by |
groupstest |
A vector indicating to which group each observation belongs. Coding should match that used for model training. Will be NULL for target grouped data. |
ytest |
Response variable. Optional. If included, |
type |
Type of prediction required. Type '"link"' gives the linear predictors for '"binomial", '"multinomial"' or '"cox"' models; for '"gaussian"' models it gives the fitted values. Type '"response"' gives the fitted probabilities for '"binomial"' or '"multinomial"', and the fitted relative-risk for '"cox"'; for '"gaussian"' type '"response"' is equivalent to type '"link"'. Note that for '"binomial"' models, results are returned only for the class corresponding to the second level of the factor response. Type '"class"' applies only to '"binomial"' or '"multinomial"' models, and produces the class label corresponding to the maximum probability. |
s |
Value of the penalty parameter 'lambda' at which predictions are required. Will use the same lambda for all models; can be a numeric value, '"lambda.min"' or '"lambda.1se"'. Default is '"lambda.min"'. |
gamma |
For use only when 'relax = TRUE' was specified during training. Value of the penalty parameter 'gamma' at which predictions are required. Will use the same gamma for all models; can be a numeric value, '"gamma.min"' or '"gamma.1se"'. Default is '"gamma.min"'. |
return.link |
If |
... |
other arguments to be passed to the |
A list containing the requested predictions. If ytest is included, will also return error measures.
call |
The call that produced this object. |
alpha |
The value(s) of alpha used to generate predictions. Will be the same alpha used to in model training. |
yhatoverall |
Predictions from the overall model. |
yhatind |
Predictions from the individual models. |
yhatpre |
Predictions from the pretrained models. |
supoverall |
Indices of the features selected by the overall model. |
supind |
Union of the indices of the features selected by the individual models. |
suppre.common |
Features selected in the first stage of pretraining. |
suppre.individual |
Union of the indices of the features selected by the pretrained models, without the features selected in the first stage. |
type.measure |
If |
erroverall |
If |
errind |
If |
errpre |
If |
linkoverall |
If |
linkind |
If |
linkpre |
If |
# Gaussian example set.seed(1234) out = gaussian.example.data() x = out$x; y=out$y; groups = out$group outtest = gaussian.example.data() xtest=outtest$x; ytest=outtest$y; groupstest=outtest$groups fit = ptLasso(x, y, groups = groups, alpha = 0.5, family = "gaussian", type.measure = "mse") pred = predict(fit, xtest, groupstest, ytest=ytest) pred# Gaussian example set.seed(1234) out = gaussian.example.data() x = out$x; y=out$y; groups = out$group outtest = gaussian.example.data() xtest=outtest$x; ytest=outtest$y; groupstest=outtest$groups fit = ptLasso(x, y, groups = groups, alpha = 0.5, family = "gaussian", type.measure = "mse") pred = predict(fit, xtest, groupstest, ytest=ytest) pred
Print the cv.ptLasso object.
## S3 method for class 'cv.ptLasso' print(x, ...)## S3 method for class 'cv.ptLasso' print(x, ...)
x |
fitted |
... |
other arguments to pass to the print function. |
Invisibly returns x. Called for its side effect of printing a summary to the console.
Erin Craig and Rob Tibshirani
Maintainer:
Erin Craig <[email protected]>
ptLasso, cv.ptLasso and predict.cv.ptLasso.
out = gaussian.example.data(k=2, class.sizes = c(50, 50)) x = out$x; y=out$y; groups = out$group; cvfit = cv.ptLasso(x, y, groups = groups, family = "gaussian", type.measure = "mse") print(cvfit)out = gaussian.example.data(k=2, class.sizes = c(50, 50)) x = out$x; y=out$y; groups = out$group; cvfit = cv.ptLasso(x, y, groups = groups, family = "gaussian", type.measure = "mse") print(cvfit)
Print the predict.cv.ptLasso object.
## S3 method for class 'predict.cv.ptLasso' print(x, ...)## S3 method for class 'predict.cv.ptLasso' print(x, ...)
x |
output of predict called with a ptLasso object. |
... |
other arguments to pass to the print function. |
Invisibly returns x. Called for its side effect of printing a summary to the console.
Erin Craig and Rob Tibshirani
Maintainer:
Erin Craig <[email protected]>
cv.ptLasso and predict.cv.ptLasso.
# Train data out = gaussian.example.data(k=2, class.sizes = c(50, 50)) x = out$x; y=out$y; groups = out$group; # Test data outtest = gaussian.example.data(k=2, class.sizes = c(50, 50)) xtest=outtest$x; ytest=outtest$y; groupstest=outtest$groups cvfit = cv.ptLasso(x, y, groups = groups, nfolds = 3, family = "gaussian", type.measure = "mse") pred = predict(cvfit, xtest, groupstest, ytest=ytest, s="lambda.min") print(pred) # If ytest is not supplied, just prints the pretrained predictions. pred = predict(cvfit, xtest, groupstest, s="lambda.min") print(pred)# Train data out = gaussian.example.data(k=2, class.sizes = c(50, 50)) x = out$x; y=out$y; groups = out$group; # Test data outtest = gaussian.example.data(k=2, class.sizes = c(50, 50)) xtest=outtest$x; ytest=outtest$y; groupstest=outtest$groups cvfit = cv.ptLasso(x, y, groups = groups, nfolds = 3, family = "gaussian", type.measure = "mse") pred = predict(cvfit, xtest, groupstest, ytest=ytest, s="lambda.min") print(pred) # If ytest is not supplied, just prints the pretrained predictions. pred = predict(cvfit, xtest, groupstest, s="lambda.min") print(pred)
Print the predict.ptLasso object.
## S3 method for class 'predict.ptLasso' print(x, ...)## S3 method for class 'predict.ptLasso' print(x, ...)
x |
output of predict called with a ptLasso object. |
... |
other arguments to pass to the print function. |
Invisibly returns x. Called for its side effect of printing a summary to the console.
Erin Craig and Rob Tibshirani
Maintainer:
Erin Craig <[email protected]>
ptLasso and predict.ptLasso.
# Train data out = gaussian.example.data(k=2, class.sizes = c(50, 50)) x = out$x; y=out$y; groups = out$group; # Test data outtest = gaussian.example.data(k=2, class.sizes = c(50, 50)) xtest=outtest$x; ytest=outtest$y; groupstest=outtest$groups fit = ptLasso(x, y, groups = groups, nfolds = 3, family = "gaussian", type.measure = "mse") pred = predict(fit, xtest, groupstest, ytest=ytest, s="lambda.min") print(pred) # If ytest is not supplied, just prints the pretrained predictions. pred = predict(fit, xtest, groupstest, s="lambda.min") print(pred)# Train data out = gaussian.example.data(k=2, class.sizes = c(50, 50)) x = out$x; y=out$y; groups = out$group; # Test data outtest = gaussian.example.data(k=2, class.sizes = c(50, 50)) xtest=outtest$x; ytest=outtest$y; groupstest=outtest$groups fit = ptLasso(x, y, groups = groups, nfolds = 3, family = "gaussian", type.measure = "mse") pred = predict(fit, xtest, groupstest, ytest=ytest, s="lambda.min") print(pred) # If ytest is not supplied, just prints the pretrained predictions. pred = predict(fit, xtest, groupstest, s="lambda.min") print(pred)
Print the ptLasso object.
## S3 method for class 'ptLasso' print(x, ...)## S3 method for class 'ptLasso' print(x, ...)
x |
fitted |
... |
other arguments to pass to the print function. |
Invisibly returns x. Called for its side effect of printing a summary to the console.
Erin Craig and Rob Tibshirani
Maintainer:
Erin Craig <[email protected]>
ptLasso, cv.ptLasso and predict.cv.ptLasso.
out = gaussian.example.data() x = out$x; y=out$y; groups = out$group; fit = ptLasso(x, y, groups = groups, family = "gaussian", type.measure = "mse") print(fit)out = gaussian.example.data() x = out$x; y=out$y; groups = out$group; fit = ptLasso(x, y, groups = groups, family = "gaussian", type.measure = "mse") print(fit)
Fits a pretrained lasso model using the glmnet package, for a fixed choice of the pretraining hyperparameter alpha. Additionally fits an "overall" model (using all data) and "individual" models (use each individual group). Can fit input-grouped data with Gaussian, multinomial, binomial or Cox outcomes, and target-grouped data, which necessarily has a multinomial outcome. Many ptLasso arguments are passed directly to glmnet, and therefore the glmnet documentation is another good reference for ptLasso.
ptLasso( x, y, groups = NULL, alpha = 0.5, family = c("default", "gaussian", "multinomial", "binomial", "cox"), type.measure = c("default", "mse", "mae", "auc", "deviance", "class", "C"), use.case = c("inputGroups", "targetGroups", "multiresponse", "timeSeries"), overall.lambda = c("lambda.1se", "lambda.min"), overall.gamma = "gamma.1se", foldid = NULL, nfolds = 10, standardize = TRUE, verbose = FALSE, weights = NULL, penalty.factor = rep(1, p), fitoverall = NULL, fitind = NULL, en.alpha = 1, group.intercepts = TRUE, ... )ptLasso( x, y, groups = NULL, alpha = 0.5, family = c("default", "gaussian", "multinomial", "binomial", "cox"), type.measure = c("default", "mse", "mae", "auc", "deviance", "class", "C"), use.case = c("inputGroups", "targetGroups", "multiresponse", "timeSeries"), overall.lambda = c("lambda.1se", "lambda.min"), overall.gamma = "gamma.1se", foldid = NULL, nfolds = 10, standardize = TRUE, verbose = FALSE, weights = NULL, penalty.factor = rep(1, p), fitoverall = NULL, fitind = NULL, en.alpha = 1, group.intercepts = TRUE, ... )
x |
input matrix, of dimension nobs x nvars; each row is an observation vector. Can be in sparse matrix format (inherit from class '"sparseMatrix"' as in package 'Matrix'). Requirement: 'nvars >1'; in other words, 'x' should have 2 or more columns. If 'use.case = "timeSeries"', then x may be a list of matrices with identical dimensions, one for each point in time. |
y |
response variable. Quantitative for 'family="gaussian"'. For 'family="binomial"' should be either a factor with two levels, or a two-column matrix of counts or proportions (the second column is treated as the target class; for a factor, the last level in alphabetical order is the target class). For 'family="multinomial"', can be a 'nc>=2' level factor, or a matrix with 'nc' columns of counts or proportions. For either '"binomial"' or '"multinomial"', if 'y' is presented as a vector, it will be coerced into a factor. For 'family="cox"', preferably a 'Surv' object from the survival package: see Detail section for more information. For 'use.case = "multiresponse"' or 'use.case = "timeSeries"', 'y' is a matrix of responses. |
groups |
A vector of length nobs indicating to which group each observation belongs. Only for 'use.case = "inputGroups"'. |
alpha |
The pretrained lasso hyperparameter, with |
family |
Either a character string representing one of the built-in families, or else a 'glm()' family object. For more information, see Details section below or the documentation for response type (see above). |
type.measure |
loss to use for cross-validation within each individual, overall, or pretrained lasso model. Currently five options, not all available for all models. The default is 'type.measure="deviance"', which uses squared-error for gaussian models (a.k.a 'type.measure="mse"' there), deviance for logistic and poisson regression, and partial-likelihood for the Cox model. 'type.measure="class"' applies to binomial and multinomial logistic regression only, and gives misclassification error. 'type.measure="auc"' is for two-class logistic regression only, and gives area under the ROC curve. 'type.measure="mse"' or 'type.measure="mae"' (mean absolute error) can be used by all models except the '"cox"'; they measure the deviation from the fitted mean to the response. 'type.measure="C"' is Harrel's concordance measure, only available for 'cox' models. |
use.case |
The type of grouping observed in the data. Can be one of "inputGroups", "targetGroups", "multiresponse" or "timeSeries". |
overall.lambda |
The choice of lambda to be used by the overall model to define the offset and penalty factor for pretrained lasso. Defaults to "lambda.1se", could alternatively be "lambda.min". This choice of lambda will be used to compute the offset and penalty factor (1) during model training and (2) during prediction. In the predict function, another lambda must be specified for the individual models, the second stage of pretraining and the overall model. |
overall.gamma |
For use only when the option |
foldid |
An optional vector of values between 1 and |
nfolds |
Number of folds for CV (default is 10). Although |
standardize |
Should the predictors be standardized before fitting (default is TRUE). |
verbose |
If |
weights |
observation weights. Default is 1 for each observation. |
penalty.factor |
Separate penalty factors can be applied to each coefficient. This is a number that multiplies 'lambda' to allow differential shrinkage. Can be 0 for some variables, which implies no shrinkage, and that variable is always included in the model. Default is 1 for all variables (and implicitly infinity for variables listed in 'exclude'). For more information, see |
fitoverall |
An optional cv.glmnet object specifying the overall model. This should have been trained on the full training data, with the argument 'keep = TRUE'. |
fitind |
An optional list of cv.glmnet objects specifying the individual models. These should have been trained on the original training data, with the argument 'keep = TRUE'. |
en.alpha |
The elasticnet mixing parameter, with 0 <= en.alpha <= 1. The penalty is defined as (1-alpha)/2||beta||_2^2+alpha||beta||_1. 'alpha=1' is the lasso penalty, and 'alpha=0' the ridge penalty. Default is 'en.alpha = 1' (lasso). |
group.intercepts |
For 'use.case = "inputGroups"' only. If 'TRUE', fit the overall model with a separate intercept for each group. If 'FALSE', ignore the grouping and fit one overall intercept. Default is 'TRUE'. |
... |
Additional arguments to be passed to the cv.glmnet functions. Notable choices include |
An object of class "ptLasso", which is a list with the ingredients of the fitted models.
call |
The call that produced this object. |
k |
The number of groups (for 'use.case = "inputGroups"'). |
nresps |
The number of responses (for 'use.case = "multiresponse"' and 'use.case = "timeseries"'). |
alpha |
The value of alpha used for pretraining. |
group.levels |
IDs for all of the groups used in training (for 'use.case = "inputGroups"'). |
group.legend |
Mapping from user-supplied group ids to numeric group ids. For internal use (e.g. with predict). (For 'use.case = "inputGroups"') |
fitoverall |
A fitted |
fitpre |
A list of fitted (pretrained) |
fitind |
A list of fitted |
fitoverall.lambda |
Lambda used with fitoverall, to compute the offset for pretraining. |
fitoverall.gamma |
Gamma used with fitoverall when 'relax = TRUE', to compute the offset for pretraining. |
Friedman, J., Hastie, T., & Tibshirani, R. (2010). Regularization paths for generalized linear models via coordinate descent. Journal of Statistical Software, 33(1), 1-22.
# Getting started. First, we simulate data: we need covariates x, response y and group IDs. set.seed(1234) n=100 p=10 n.groups=2 x = matrix(rnorm(n*p), n, p) y = rnorm(n) groups = sort(rep(1:n.groups, n/n.groups)) xtest = matrix(rnorm(n*p), n, p) ytest = rnorm(n) groupstest = sort(rep(1:n.groups, n/n.groups)) # Now, we can fit a ptLasso model: fit = ptLasso(x, y, groups = groups, alpha = 0.5, family = "gaussian", nfolds = 3, type.measure = "mse") plot(fit) # to see all of the cv.glmnet models trained predict(fit, xtest, groupstest) # to predict on new data predict(fit, xtest, groupstest, ytest=ytest) # if ytest is included, we also measure performance # When we trained our model, we used "lambda.1se" in the first stage of pretraining by default. # This is a necessary choice to make during model training; we need to select the model # we want to use to define the offset and penalty factor for the second stage of pretraining. # We could instead have used "lambda.min": fit = ptLasso(x, y, groups = groups, alpha = 0.5, family = "gaussian", type.measure = "mse", nfolds = 3, overall.lambda = "lambda.min") # We can use the 'relax' option to fit relaxed lasso models: fit = ptLasso(x, y, groups = groups, alpha = 0.5, family = "gaussian", type.measure = "mse", relax = TRUE) # As we did for lambda, we may want to specify the choice of gamma for stage one # of pretraining. (The default is "gamma.1se".) fit = ptLasso(x, y, groups = groups, alpha = 0.5, family = "gaussian", type.measure = "mse", relax = TRUE, overall.gamma = "gamma.min") # In practice, we may want to try many values of alpha. # alpha may range from 0 (the overall model with fine tuning) to 1 (the individual models). # To choose alpha, you may either (1) run ptLasso with different values of alpha # and measure performance with a validation set, or (2) use cv.ptLasso. # Now, we are ready to simulate slightly more realistic data. # This continuous outcome example has k = 5 groups, where each group has 200 observations. # There are scommon = 10 features shared across all groups, and # sindiv = 10 features unique to each group. # n = 1000 and p = 120 (60 informative features and 60 noise features). # The coefficients of the common features differ across groups (beta.common). # In group 1, these coefficients are rep(1, 10); in group 2 they are rep(2, 10), etc. # Each group has 10 unique features, the coefficients of which are all 3 (beta.indiv). # The intercept in all groups is 0. # The variable sigma = 20 indicates that we add noise to y according to 20 * rnorm(n). set.seed(1234) k=5 class.sizes=rep(200, k) scommon=10; sindiv=rep(10, k) n=sum(class.sizes); p=2*(sum(sindiv) + scommon) beta.common=3*(1:k); beta.indiv=rep(3, k) intercepts=rep(0, k) sigma=20 out = gaussian.example.data(k=k, class.sizes=class.sizes, scommon=scommon, sindiv=sindiv, n=n, p=p, beta.common=beta.common, beta.indiv=beta.indiv, intercepts=intercepts, sigma=20) x = out$x; y=out$y; groups = out$group outtest = gaussian.example.data(k=k, class.sizes=class.sizes, scommon=scommon, sindiv=sindiv, n=n, p=p, beta.common=beta.common, beta.indiv=beta.indiv, intercepts=intercepts, sigma=20) xtest=outtest$x; ytest=outtest$y; groupstest=outtest$groups fit = ptLasso(x, y, groups = groups, alpha = 0.5, family = "gaussian", type.measure = "mse") plot(fit) # to see all of the cv.glmnet models trained predict(fit, xtest, groupstest, ytest=ytest) # Now, we repeat with a binomial outcome. # This example has k = 3 groups, where each group has 100 observations. # There are scommon = 5 features shared across all groups, and # sindiv = 5 features unique to each group. # n = 300 and p = 40 (20 informative features and 20 noise features). # The coefficients of the common features differ across groups (beta.common), # as do the coefficients specific to each group (beta.indiv). set.seed(1234) k=3 class.sizes=rep(100, k) scommon=5; sindiv=rep(5, k) n=sum(class.sizes); p=2*(sum(sindiv) + scommon) beta.common=list(c(-.5, .5, .3, -.9, .1), c(-.3, .9, .1, -.1, .2), c(0.1, .2, -.1, .2, .3)) beta.indiv = lapply(1:k, function(i) 0.9 * beta.common[[i]]) out = binomial.example.data(k=k, class.sizes=class.sizes, scommon=scommon, sindiv=sindiv, n=n, p=p, beta.common=beta.common, beta.indiv=beta.indiv) x = out$x; y=out$y; groups = out$group outtest = binomial.example.data(k=k, class.sizes=class.sizes, scommon=scommon, sindiv=sindiv, n=n, p=p, beta.common=beta.common, beta.indiv=beta.indiv) xtest=outtest$x; ytest=outtest$y; groupstest=outtest$groups fit = ptLasso(x, y, groups = groups, alpha = 0.5, family = "binomial", type.measure = "auc") plot(fit) # to see all of the cv.glmnet models trained predict(fit, xtest, groupstest, ytest=ytest) ## Not run: ### Model fitting with parallel = TRUE require(doMC) registerDoMC(cores = 4) fit = ptLasso(x, y, groups = groups, family = "gaussian", type.measure = "mse", parallel=TRUE) ## End(Not run) # Multiresponse pretraining: # Now let's consider the case of a multiresponse outcome. We'll start by simulating data: set.seed(1234) n = 1000; ntrain = 500; p = 500 sigma = 2 x = matrix(rnorm(n*p), n, p) beta1 = c(rep(1, 5), rep(0.5, 5), rep(0, p - 10)) beta2 = c(rep(1, 5), rep(0, 5), rep(0.5, 5), rep(0, p - 15)) mu = cbind(x %*% beta1, x %*% beta2) y = cbind(mu[, 1] + sigma * rnorm(n), mu[, 2] + sigma * rnorm(n)) cat("SNR for the two tasks:", round(diag(var(mu)/var(y-mu)), 2), fill=TRUE) cat("Correlation between two tasks:", cor(y[, 1], y[, 2]), fill=TRUE) xtest = x[-(1:ntrain), ] ytest = y[-(1:ntrain), ] x = x[1:ntrain, ] y = y[1:ntrain, ] fit = ptLasso(x, y, type.measure = "mse", use.case = "multiresponse") plot(fit) # to see all of the cv.glmnet models trained predict(fit, xtest) # to predict with new data predict(fit, xtest, ytest=ytest) # if ytest is included, we also measure performance # By default, we used lambda = "lambda.min" to measure performance. # We could instead use lambda = "lambda.1se": predict(fit, xtest, ytest=ytest, s="lambda.1se") # We could also use the glmnet option relax = TRUE: fit = ptLasso(x, y, type.measure = "mse", relax = TRUE, use.case = "multiresponse") # Time series pretraining # Now suppose we have time series data with a binomial outcome measured at 3 different time points. set.seed(1234) n = 600; ntrain = 300; p = 50 x = matrix(rnorm(n*p), n, p) beta1 = c(rep(0.5, 10), rep(0, p-10)) beta2 = beta1 + c(rep(0, 10), runif(5, min = 0, max = 0.5), rep(0, p-15)) beta3 = beta1 + c(rep(0, 10), runif(5, min = 0, max = 0.5), rep(.5, 5), rep(0, p-20)) y1 = rbinom(n, 1, prob = 1/(1 + exp(-x %*% beta1))) y2 = rbinom(n, 1, prob = 1/(1 + exp(-x %*% beta2))) y3 = rbinom(n, 1, prob = 1/(1 + exp(-x %*% beta3))) y = cbind(y1, y2, y3) xtest = x[-(1:ntrain), ] ytest = y[-(1:ntrain), ] x = x[1:ntrain, ] y = y[1:ntrain, ] fit = ptLasso(x, y, use.case="timeSeries", family="binomial", type.measure = "auc") plot(fit) predict(fit, xtest, ytest=ytest) # The glmnet option relax = TRUE: fit = ptLasso(x, y, type.measure = "auc", family = "binomial", relax = TRUE, use.case = "timeSeries") plot(fit) predict(fit, xtest, ytest=ytest)# Getting started. First, we simulate data: we need covariates x, response y and group IDs. set.seed(1234) n=100 p=10 n.groups=2 x = matrix(rnorm(n*p), n, p) y = rnorm(n) groups = sort(rep(1:n.groups, n/n.groups)) xtest = matrix(rnorm(n*p), n, p) ytest = rnorm(n) groupstest = sort(rep(1:n.groups, n/n.groups)) # Now, we can fit a ptLasso model: fit = ptLasso(x, y, groups = groups, alpha = 0.5, family = "gaussian", nfolds = 3, type.measure = "mse") plot(fit) # to see all of the cv.glmnet models trained predict(fit, xtest, groupstest) # to predict on new data predict(fit, xtest, groupstest, ytest=ytest) # if ytest is included, we also measure performance # When we trained our model, we used "lambda.1se" in the first stage of pretraining by default. # This is a necessary choice to make during model training; we need to select the model # we want to use to define the offset and penalty factor for the second stage of pretraining. # We could instead have used "lambda.min": fit = ptLasso(x, y, groups = groups, alpha = 0.5, family = "gaussian", type.measure = "mse", nfolds = 3, overall.lambda = "lambda.min") # We can use the 'relax' option to fit relaxed lasso models: fit = ptLasso(x, y, groups = groups, alpha = 0.5, family = "gaussian", type.measure = "mse", relax = TRUE) # As we did for lambda, we may want to specify the choice of gamma for stage one # of pretraining. (The default is "gamma.1se".) fit = ptLasso(x, y, groups = groups, alpha = 0.5, family = "gaussian", type.measure = "mse", relax = TRUE, overall.gamma = "gamma.min") # In practice, we may want to try many values of alpha. # alpha may range from 0 (the overall model with fine tuning) to 1 (the individual models). # To choose alpha, you may either (1) run ptLasso with different values of alpha # and measure performance with a validation set, or (2) use cv.ptLasso. # Now, we are ready to simulate slightly more realistic data. # This continuous outcome example has k = 5 groups, where each group has 200 observations. # There are scommon = 10 features shared across all groups, and # sindiv = 10 features unique to each group. # n = 1000 and p = 120 (60 informative features and 60 noise features). # The coefficients of the common features differ across groups (beta.common). # In group 1, these coefficients are rep(1, 10); in group 2 they are rep(2, 10), etc. # Each group has 10 unique features, the coefficients of which are all 3 (beta.indiv). # The intercept in all groups is 0. # The variable sigma = 20 indicates that we add noise to y according to 20 * rnorm(n). set.seed(1234) k=5 class.sizes=rep(200, k) scommon=10; sindiv=rep(10, k) n=sum(class.sizes); p=2*(sum(sindiv) + scommon) beta.common=3*(1:k); beta.indiv=rep(3, k) intercepts=rep(0, k) sigma=20 out = gaussian.example.data(k=k, class.sizes=class.sizes, scommon=scommon, sindiv=sindiv, n=n, p=p, beta.common=beta.common, beta.indiv=beta.indiv, intercepts=intercepts, sigma=20) x = out$x; y=out$y; groups = out$group outtest = gaussian.example.data(k=k, class.sizes=class.sizes, scommon=scommon, sindiv=sindiv, n=n, p=p, beta.common=beta.common, beta.indiv=beta.indiv, intercepts=intercepts, sigma=20) xtest=outtest$x; ytest=outtest$y; groupstest=outtest$groups fit = ptLasso(x, y, groups = groups, alpha = 0.5, family = "gaussian", type.measure = "mse") plot(fit) # to see all of the cv.glmnet models trained predict(fit, xtest, groupstest, ytest=ytest) # Now, we repeat with a binomial outcome. # This example has k = 3 groups, where each group has 100 observations. # There are scommon = 5 features shared across all groups, and # sindiv = 5 features unique to each group. # n = 300 and p = 40 (20 informative features and 20 noise features). # The coefficients of the common features differ across groups (beta.common), # as do the coefficients specific to each group (beta.indiv). set.seed(1234) k=3 class.sizes=rep(100, k) scommon=5; sindiv=rep(5, k) n=sum(class.sizes); p=2*(sum(sindiv) + scommon) beta.common=list(c(-.5, .5, .3, -.9, .1), c(-.3, .9, .1, -.1, .2), c(0.1, .2, -.1, .2, .3)) beta.indiv = lapply(1:k, function(i) 0.9 * beta.common[[i]]) out = binomial.example.data(k=k, class.sizes=class.sizes, scommon=scommon, sindiv=sindiv, n=n, p=p, beta.common=beta.common, beta.indiv=beta.indiv) x = out$x; y=out$y; groups = out$group outtest = binomial.example.data(k=k, class.sizes=class.sizes, scommon=scommon, sindiv=sindiv, n=n, p=p, beta.common=beta.common, beta.indiv=beta.indiv) xtest=outtest$x; ytest=outtest$y; groupstest=outtest$groups fit = ptLasso(x, y, groups = groups, alpha = 0.5, family = "binomial", type.measure = "auc") plot(fit) # to see all of the cv.glmnet models trained predict(fit, xtest, groupstest, ytest=ytest) ## Not run: ### Model fitting with parallel = TRUE require(doMC) registerDoMC(cores = 4) fit = ptLasso(x, y, groups = groups, family = "gaussian", type.measure = "mse", parallel=TRUE) ## End(Not run) # Multiresponse pretraining: # Now let's consider the case of a multiresponse outcome. We'll start by simulating data: set.seed(1234) n = 1000; ntrain = 500; p = 500 sigma = 2 x = matrix(rnorm(n*p), n, p) beta1 = c(rep(1, 5), rep(0.5, 5), rep(0, p - 10)) beta2 = c(rep(1, 5), rep(0, 5), rep(0.5, 5), rep(0, p - 15)) mu = cbind(x %*% beta1, x %*% beta2) y = cbind(mu[, 1] + sigma * rnorm(n), mu[, 2] + sigma * rnorm(n)) cat("SNR for the two tasks:", round(diag(var(mu)/var(y-mu)), 2), fill=TRUE) cat("Correlation between two tasks:", cor(y[, 1], y[, 2]), fill=TRUE) xtest = x[-(1:ntrain), ] ytest = y[-(1:ntrain), ] x = x[1:ntrain, ] y = y[1:ntrain, ] fit = ptLasso(x, y, type.measure = "mse", use.case = "multiresponse") plot(fit) # to see all of the cv.glmnet models trained predict(fit, xtest) # to predict with new data predict(fit, xtest, ytest=ytest) # if ytest is included, we also measure performance # By default, we used lambda = "lambda.min" to measure performance. # We could instead use lambda = "lambda.1se": predict(fit, xtest, ytest=ytest, s="lambda.1se") # We could also use the glmnet option relax = TRUE: fit = ptLasso(x, y, type.measure = "mse", relax = TRUE, use.case = "multiresponse") # Time series pretraining # Now suppose we have time series data with a binomial outcome measured at 3 different time points. set.seed(1234) n = 600; ntrain = 300; p = 50 x = matrix(rnorm(n*p), n, p) beta1 = c(rep(0.5, 10), rep(0, p-10)) beta2 = beta1 + c(rep(0, 10), runif(5, min = 0, max = 0.5), rep(0, p-15)) beta3 = beta1 + c(rep(0, 10), runif(5, min = 0, max = 0.5), rep(.5, 5), rep(0, p-20)) y1 = rbinom(n, 1, prob = 1/(1 + exp(-x %*% beta1))) y2 = rbinom(n, 1, prob = 1/(1 + exp(-x %*% beta2))) y3 = rbinom(n, 1, prob = 1/(1 + exp(-x %*% beta3))) y = cbind(y1, y2, y3) xtest = x[-(1:ntrain), ] ytest = y[-(1:ntrain), ] x = x[1:ntrain, ] y = y[1:ntrain, ] fit = ptLasso(x, y, use.case="timeSeries", family="binomial", type.measure = "auc") plot(fit) predict(fit, xtest, ytest=ytest) # The glmnet option relax = TRUE: fit = ptLasso(x, y, type.measure = "auc", family = "binomial", relax = TRUE, use.case = "timeSeries") plot(fit) predict(fit, xtest, ytest=ytest)