From de1fa6c60881841896bed5a5fd0238db564288a2 Mon Sep 17 00:00:00 2001 From: Patrick Breheny Date: Thu, 18 May 2017 08:30:43 -0500 Subject: [PATCH] Changed convergence criterion --- NEWS | 4 +- inst/tests/basic-functionality.R | 148 ++++++++++++------------------- src/gdfit_binomial.c | 15 ++-- src/gdfit_cox.c | 19 ++-- src/gdfit_poisson.c | 20 ++--- src/grpreg_init.c | 21 ----- src/lcdfit_binomial.c | 24 ++--- src/lcdfit_cox.c | 24 ++--- src/lcdfit_gaussian.c | 34 +++---- src/lcdfit_poisson.c | 22 ++--- src/maxprod.c | 4 +- 11 files changed, 144 insertions(+), 191 deletions(-) diff --git a/NEWS b/NEWS index 39d96d4..5a2fa3b 100644 --- a/NEWS +++ b/NEWS @@ -1,10 +1,8 @@ -Pending: - * Change: Convergence criterion now based on RMSD of linear predictors - 3.1-0 () * New: Additional tests and support for coersion of various types with respect to both X and y * Internal: new SSR-BEDPP feature screening rule for group lasso + * Change: Convergence criterion now based on RMSD of linear predictors * Change: 'Lung' and 'Birthwt' data sets now use factor representation of group, as character vectors are inherently ambiguous with respect to order * Change: max.iter now based on total number of iterations for entire path diff --git a/inst/tests/basic-functionality.R b/inst/tests/basic-functionality.R index 7252370..59367c5 100644 --- a/inst/tests/basic-functionality.R +++ b/inst/tests/basic-functionality.R @@ -1,5 +1,4 @@ #source("~/dev/.grpreg.setup.R") -set.seed(1) .test = "grpreg() reproduces simple linear regression" n <- 5 @@ -12,22 +11,22 @@ nlam=100 par(mfcol=c(3,2)) gel <- coef(fit <- grpreg(X, y, group, penalty="gel", nlambda=nlam, lambda.min=0))[,nlam] plot(fit, main=fit$penalty) -check(gel, reg, tolerance=.01) +check(gel, reg) cMCP <- coef(fit <- grpreg(X, y, group, penalty="cMCP", nlambda=nlam, lambda.min=0))[,nlam] plot(fit, main=fit$penalty) -check(cMCP, reg, tolerance=.01) +check(cMCP, reg) bridge <- coef(fit <- gBridge(X, y, group, nlambda=nlam, lambda.min=0))[,1] plot(fit, main=fit$penalty) -check(bridge, reg, tolerance=.01) +check(bridge, reg) grLasso <- coef(fit <- grpreg(X, y, group, penalty="grLasso", nlambda=nlam, lambda.min=0))[,nlam] plot(fit, main=fit$penalty) -check(grLasso, reg, tolerance=.01) +check(grLasso, reg) grMCP <- coef(fit <- grpreg(X, y, group, penalty="grMCP", nlambda=nlam, lambda.min=0))[,nlam] plot(fit, main=fit$penalty) -check(grMCP, reg, tolerance=.01) +check(grMCP, reg) grSCAD <- coef(fit <- grpreg(X, y, group, penalty="grSCAD", nlambda=nlam, lambda.min=0))[,nlam] plot(fit, main=fit$penalty) -check(grSCAD, reg, tolerance=.01) +check(grSCAD, reg) .test = "grpreg reproduces linear regression" n <- 50 @@ -39,25 +38,25 @@ fit.mle <- lm(y~X) reg <- coef(fit.mle) nlam=100 par(mfcol=c(3,2)) -gel <- coef(fit <- grpreg(X, y, group, penalty="gel", nlambda=nlam, lambda.min=0))[,nlam] +gel <- coef(fit <- grpreg(X, y, group, penalty="gel", nlambda=nlam, lambda.min=0, eps=1e-10))[,nlam] plot(fit, main=fit$penalty) -check(gel, reg, tolerance=.01) -cMCP <- coef(fit <- grpreg(X, y, group, penalty="cMCP", lambda.min=0))[,100] +check(gel, reg) +cMCP <- coef(fit <- grpreg(X, y, group, penalty="cMCP", lambda.min=0, eps=1e-10))[,100] plot(fit, main=fit$penalty) -check(cMCP, reg, tolerance=.01) -bridge <- coef(fit <- gBridge(X, y, group, lambda.min=0))[,1] +check(cMCP, reg) +bridge <- coef(fit <- gBridge(X, y, group, lambda.min=0, eps=1e-10))[,1] plot(fit, main=fit$penalty) -check(bridge, reg, tolerance=.01) -grLasso <- coef(fit <- grpreg(X, y, group, penalty="grLasso", lambda.min=0))[,100] +check(bridge, reg) +grLasso <- coef(fit <- grpreg(X, y, group, penalty="grLasso", lambda.min=0, eps=1e-10))[,100] plot(fit, main=fit$penalty) -check(grLasso, reg, tolerance=.01) -grMCP <- coef(fit <- grpreg(X, y, group, penalty="grMCP", lambda.min=0))[,100] +check(grLasso, reg) +grMCP <- coef(fit <- grpreg(X, y, group, penalty="grMCP", lambda.min=0, eps=1e-10))[,100] plot(fit, main=fit$penalty) -check(grMCP, reg, tolerance=.01) -grSCAD <- coef(fit <- grpreg(X, y, group, penalty="grSCAD", lambda.min=0))[,100] +check(grMCP, reg) +grSCAD <- coef(fit <- grpreg(X, y, group, penalty="grSCAD", lambda.min=0, eps=1e-10))[,100] plot(fit, main=fit$penalty) -check(grSCAD, reg, tolerance=.01) -check(predict(fit, X)[,100], predict(fit.mle), tolerance=.001) +check(grSCAD, reg) +check(predict(fit, X)[,100], predict(fit.mle)) .test = "grpreg reproduces simple logistic regression" n <- 20 @@ -68,24 +67,24 @@ group <- 1 reg <- glm(y~X, family="binomial")$coef nlam <- 100 par(mfcol=c(3,2)) -gel <- coef(fit <- grpreg(X, y, group, penalty="gel", nlambda=nlam, lambda.min=0, family="binomial"))[,nlam] +gel <- coef(fit <- grpreg(X, y, group, penalty="gel", nlambda=nlam, lambda.min=0, family="binomial", eps=1e-10))[,nlam] plot(fit, main=fit$penalty) -check(gel, reg, tolerance=.01) -cMCP <- coef(fit <- grpreg(X, y, group, penalty="cMCP", nlambda=nlam, lambda.min=0, family="binomial", gamma=12))[,nlam] +check(gel, reg) +cMCP <- coef(fit <- grpreg(X, y, group, penalty="cMCP", nlambda=nlam, lambda.min=0, family="binomial", gamma=9, eps=1e-10))[,nlam] plot(fit, main=fit$penalty) -check(cMCP, reg, tolerance=.01) -bridge <- coef(fit <- gBridge(X, y, group, lambda.min=0, nlambda=nlam, family="binomial"))[,1] +check(cMCP, reg) +bridge <- coef(fit <- gBridge(X, y, group, lambda.min=0, nlambda=nlam, family="binomial", eps=1e-10))[,1] plot(fit, main=fit$penalty) -check(bridge, reg, tolerance=.01) -grLasso <- coef(fit <- grpreg(X, y, group, penalty="grLasso", nlambda=nlam, lambda.min=0, family="binomial"))[,nlam] +check(bridge, reg) +grLasso <- coef(fit <- grpreg(X, y, group, penalty="grLasso", nlambda=nlam, lambda.min=0, family="binomial", eps=1e-10))[,nlam] plot(fit, main=fit$penalty) -check(grLasso, reg, tolerance=.01) -grMCP <- coef(fit <- grpreg(X, y, group, penalty="grMCP", nlambda=nlam, lambda.min=0, family="binomial"))[,nlam] +check(grLasso, reg) +grMCP <- coef(fit <- grpreg(X, y, group, penalty="grMCP", nlambda=nlam, lambda.min=0, family="binomial", eps=1e-10))[,nlam] plot(fit, main=fit$penalty) -check(grMCP, reg, tolerance=.01) -grSCAD <- coef(fit <- grpreg(X, y, group, penalty="grSCAD", nlambda=nlam, lambda.min=0, family="binomial"))[,nlam] +check(grMCP, reg) +grSCAD <- coef(fit <- grpreg(X, y, group, penalty="grSCAD", nlambda=nlam, lambda.min=0, family="binomial", eps=1e-10))[,nlam] plot(fit, main=fit$penalty) -check(grSCAD, reg, tolerance=.01) +check(grSCAD, reg) .test = "grpreg() reproduces logistic regression" n <- 50 @@ -96,26 +95,26 @@ y <- runif(n) > .5 fit.mle <- glm(y~X, family="binomial") reg <- coef(fit.mle) par(mfcol=c(3,2)) -gel <- coef(fit <- grpreg(X, y, group, penalty="gel", family="binomial"))[,100] +gel <- coef(fit <- grpreg(X, y, group, penalty="gel", family="binomial", eps=1e-10, lambda.min=0))[,100] plot(fit, main=fit$penalty) -check(gel, reg, tolerance=.01) -cMCP <- coef(fit <- grpreg(X, y, group, penalty="cMCP", family="binomial"))[,100] +check(gel, reg) +cMCP <- coef(fit <- grpreg(X, y, group, penalty="cMCP", family="binomial", eps=1e-10, lambda.min=0))[,100] plot(fit, main=fit$penalty) -check(cMCP, reg, tolerance=.01) -bridge <- coef(fit <- gBridge(X, y, group, family="binomial"))[,1] +check(cMCP, reg) +bridge <- coef(fit <- gBridge(X, y, group, family="binomial", eps=1e-10, lambda.min=0))[,1] plot(fit, main=fit$penalty) -check(bridge, reg, tolerance=.01) -grLasso <- coef(fit <- grpreg(X, y, group, penalty="grLasso", family="binomial"))[,100] +check(bridge, reg) +grLasso <- coef(fit <- grpreg(X, y, group, penalty="grLasso", family="binomial", eps=1e-10, lambda.min=0))[,100] plot(fit, main=fit$penalty) -check(grLasso, reg, tolerance=.01) -grMCP <- coef(fit <- grpreg(X, y, group, penalty="grMCP", family="binomial", gamma=2))[,100] +check(grLasso, reg) +grMCP <- coef(fit <- grpreg(X, y, group, penalty="grMCP", family="binomial", gamma=2, eps=1e-10, lambda.min=0))[,100] plot(fit, main=fit$penalty) -check(grMCP, reg, tolerance=.01) -grSCAD <- coef(fit <- grpreg(X, y, group, penalty="grSCAD", family="binomial", gamma=2.1))[,100] +check(grMCP, reg) +grSCAD <- coef(fit <- grpreg(X, y, group, penalty="grSCAD", family="binomial", gamma=2.1, eps=1e-10, lambda.min=0))[,100] plot(fit, main=fit$penalty) -check(grSCAD, reg, tolerance=.01) -check(predict(fit, X)[,100], predict(fit.mle), tolerance=.001) -check(predict(fit, X, type="response")[,100], predict(fit.mle, type="response"), tolerance=.001) +check(grSCAD, reg) +check(predict(fit, X)[,100], predict(fit.mle)) +check(predict(fit, X, type="response")[,100], predict(fit.mle, type="response")) .test = "grpreg() reproduces poisson regression" n <- 50 @@ -126,52 +125,23 @@ y <- sample(1:5, n, replace=TRUE) fit.mle <- glm(y~X, family="poisson") reg <- coef(fit.mle) par(mfcol=c(3,2)) -gel <- coef(fit <- grpreg(X, y, group, penalty="gel", family="poisson"))[,100] +gel <- coef(fit <- grpreg(X, y, group, penalty="gel", family="poisson", eps=1e-10, lambda.min=0))[,100] plot(fit, main=fit$penalty) -check(gel, reg, tolerance=.01) -cMCP <- coef(fit <- grpreg(X, y, group, penalty="cMCP", family="poisson"))[,100] +check(gel, reg) +cMCP <- coef(fit <- grpreg(X, y, group, penalty="cMCP", family="poisson", eps=1e-10, lambda.min=0))[,100] plot(fit, main=fit$penalty) -check(cMCP, reg, tolerance=.01) -bridge <- coef(fit <- gBridge(X, y, group, family="poisson"))[,1] +check(cMCP, reg) +bridge <- coef(fit <- gBridge(X, y, group, family="poisson", eps=1e-10, lambda.min=0))[,1] plot(fit, main=fit$penalty) -check(bridge, reg, tolerance=.01) -grLasso <- coef(fit <- grpreg(X, y, group, penalty="grLasso", family="poisson"))[,100] +check(bridge, reg) +grLasso <- coef(fit <- grpreg(X, y, group, penalty="grLasso", family="poisson", eps=1e-10, lambda.min=0))[,100] plot(fit, main=fit$penalty) -check(grLasso, reg, tolerance=.01) -grMCP <- coef(fit <- grpreg(X, y, group, penalty="grMCP", family="poisson", gamma=2))[,100] +check(grLasso, reg) +grMCP <- coef(fit <- grpreg(X, y, group, penalty="grMCP", family="poisson", gamma=2, eps=1e-10, lambda.min=0))[,100] plot(fit, main=fit$penalty) -check(grMCP, reg, tolerance=.01) -grSCAD <- coef(fit <- grpreg(X, y, group, penalty="grSCAD", family="poisson", gamma=2.1))[,100] +check(grMCP, reg) +grSCAD <- coef(fit <- grpreg(X, y, group, penalty="grSCAD", family="poisson", gamma=2.1, eps=1e-10, lambda.min=0))[,100] plot(fit, main=fit$penalty) -check(grSCAD, reg, tolerance=.01) -check(predict(fit, X)[,100], predict(fit.mle), tolerance=.001) -check(predict(fit, X, type="response")[,100], predict(fit.mle, type="response"), tolerance=.001) - - - - -.test = "grpreg exactly reproduces linear regression" -n <- 50 -p <- 10 -X <- matrix(rnorm(n*p),ncol=p) -y <- rnorm(n) -group <- rep(0:3,4:1) -group <- 1:10 -reg <- coef(fit.mle <- lm(y~X)) -gel <- coef(fit <- grpreg(X, y, group, penalty="gel", lambda.min=0, eps=1e-10), lambda=0) -grl <- coef(fit <- grpreg(X, y, group, penalty="grLasso", lambda.min=0, eps=1e-10), lambda=0) -check(reg, grl) -check(reg, gel) - -.test = "grpreg exactly reproduces linear regression" -n <- 50 -p <- 10 -X <- ncvreg::std(matrix(rnorm(n*p),ncol=p)) -y <- rnorm(n) -group <- rep(0:3,4:1) -group <- 1:10 -reg <- coef(fit.mle <- lm(y~X)) -gel <- coef(fit <- grpreg(X, y, group, penalty="gel", lambda.min=0, eps=1e-10), lambda=0) -grl <- coef(fit <- grpreg(X, y, group, penalty="grLasso", lambda.min=0, eps=1e-10), lambda=0) -check(reg, grl) -check(reg, gel) +check(grSCAD, reg) +check(predict(fit, X)[,100], predict(fit.mle)) +check(predict(fit, X, type="response")[,100], predict(fit.mle, type="response")) diff --git a/src/gdfit_binomial.c b/src/gdfit_binomial.c index 4705f5a..0fbcb95 100644 --- a/src/gdfit_binomial.c +++ b/src/gdfit_binomial.c @@ -4,7 +4,6 @@ #include #include #include -int checkConvergence(double *beta, double *beta_old, double eps, int l, int J); double crossprod(double *x, double *y, int n, int j); double sum(double *x, int n); double norm(double *x, int p); @@ -15,7 +14,7 @@ double MCP(double theta, double l, double a); double dMCP(double theta, double l, double a); // Group descent update -- binomial -double gd_binomial(double *b, double *x, double *r, double *eta, int g, int *K1, int n, int l, int p, const char *penalty, double lam1, double lam2, double gamma, SEXP df, double *a, double maxChange) { +void gd_binomial(double *b, double *x, double *r, double *eta, int g, int *K1, int n, int l, int p, const char *penalty, double lam1, double lam2, double gamma, SEXP df, double *a, double *maxChange) { // Calculate z int K = K1[g+1] - K1[g]; @@ -34,7 +33,7 @@ double gd_binomial(double *b, double *x, double *r, double *eta, int g, int *K1, for (int j=K1[g]; j maxChange) maxChange = fabs(shift); + if (fabs(shift) > maxChange[0]) maxChange[0] = fabs(shift); for (int i=0; i 0) REAL(df)[l] += K * len / z_norm; Free(z); - return(maxChange); } SEXP gdfit_binomial(SEXP X_, SEXP y_, SEXP penalty_, SEXP K1_, SEXP K0_, SEXP lambda, SEXP alpha_, SEXP eps_, SEXP max_iter_, SEXP gamma_, SEXP group_multiplier, SEXP dfmax_, SEXP gmax_, SEXP warn_, SEXP user_) { @@ -101,7 +99,7 @@ SEXP gdfit_binomial(SEXP X_, SEXP y_, SEXP penalty_, SEXP K1_, SEXP K0_, SEXP la int *e = Calloc(J, int); for (int g=0; g maxChange) maxChange = fabs(shift); b[l*p+j] = shift + a[j]; for (int i=0; i #include #include -int checkConvergence(double *beta, double *beta_old, double eps, int l, int J); double crossprod(double *x, double *y, int n, int j); double norm(double *x, int p); double S(double z, double l); @@ -14,7 +13,7 @@ double Fs(double z, double l1, double l2, double gamma); // Group descent update -- cox void gd_cox(double *b, double *x, double *r, double *eta, double v, int g, int *K1, int n, int l, int p, const char *penalty, double lam1, - double lam2, double gamma, SEXP df, double *a) { + double lam2, double gamma, SEXP df, double *a, double *maxChange) { // Calculate z int K = K1[g+1] - K1[g]; @@ -32,6 +31,7 @@ void gd_cox(double *b, double *x, double *r, double *eta, double v, int g, for (int j=K1[g]; j maxChange[0]) maxChange[0] = fabs(shift); for (int i=0; i maxChange) maxChange = fabs(shift); b[l*p+j] = shift + a[j]; for (int i=0; i #include #include -int checkConvergence(double *beta, double *beta_old, double eps, int l, int J); double crossprod(double *x, double *y, int n, int j); double sum(double *x, int n); double max(double *x, int n); @@ -16,7 +15,7 @@ double MCP(double theta, double l, double a); double dMCP(double theta, double l, double a); // Group descent update -- poisson -void gd_poisson(double *b, double *x, double *r, double v, double *eta, int g, int *K1, int n, int l, int p, const char *penalty, double lam1, double lam2, double gamma, SEXP df, double *a) { +void gd_poisson(double *b, double *x, double *r, double v, double *eta, int g, int *K1, int n, int l, int p, const char *penalty, double lam1, double lam2, double gamma, SEXP df, double *a, double *maxChange) { // Calculate z int K = K1[g+1] - K1[g]; @@ -34,6 +33,7 @@ void gd_poisson(double *b, double *x, double *r, double v, double *eta, int g, i for (int j=K1[g]; j maxChange[0]) maxChange[0] = fabs(shift); for (int i=0; i maxChange) maxChange = fabs(shift); b[l*p+j] = shift + a[j]; for (int i=0; i eps) { - converged = 0; - break; - } - } else if (beta[l*J+j]==0 & beta_old[j]!=0) { - converged = 0; - break; - } else if (beta[l*J+j]!=0 & beta_old[j]==0) { - converged = 0; - break; - } - } - return(converged); -} - // Cross product of the jth column of x with y double crossprod(double *x, double *y, int n, int j) { double val = 0; diff --git a/src/lcdfit_binomial.c b/src/lcdfit_binomial.c index b6eb952..487aa99 100644 --- a/src/lcdfit_binomial.c +++ b/src/lcdfit_binomial.c @@ -4,7 +4,6 @@ #include #include #include -int checkConvergence(double *beta, double *beta_old, double eps, int l, int J); double crossprod(double *x, double *y, int n, int j); double sum(double *x, int n); double norm(double *x, int p); @@ -15,11 +14,12 @@ double MCP(double theta, double l, double a); double dMCP(double theta, double l, double a); // Groupwise local coordinate descent updates -- binomial -void gLCD_binomial(double *b, const char *penalty, double *x, double *r, double *eta, int g, int *K1, int n, int l, int p, double lam1, double lam2, double gamma, double tau, SEXP df, double *a, double delta, int *e) { +void gLCD_binomial(double *b, const char *penalty, double *x, double *r, double *eta, int g, int *K1, int n, int l, int p, double lam1, double lam2, double gamma, double tau, SEXP df, double *a, double delta, int *e, double *maxChange) { // Calculate v int K = K1[g+1] - K1[g]; double v = 0.25; + double shift; // Make initial local approximation double sG = 0; // Sum of inner penalties for group @@ -34,7 +34,9 @@ void gLCD_binomial(double *b, const char *penalty, double *x, double *r, double if (sG < delta) { for (int j=K1[g]; j maxChange[0]) maxChange[0] = fabs(shift); + for (int i=0; i maxChange[0]) maxChange[0] = fabs(shift); for (int i=0; i maxChange) maxChange = fabs(shift); b[l*p+j] = shift + a[j]; for (int i=0; i #include #include -int checkConvergence(double *beta, double *beta_old, double eps, int l, int J); double crossprod(double *x, double *y, int n, int j); double wcrossprod(double *X, double *y, double *w, int n, int j); double wsqsum(double *X, double *w, int n, int j); @@ -16,12 +15,13 @@ double dMCP(double theta, double l, double a); // Groupwise local coordinate descent updates -- Cox void gLCD_cox(double *b, const char *penalty, double *X, double *r, double *eta, double *h, int g, int *K1, int n, int l, int p, double lam1, double lam2, - double gamma, double tau, SEXP df, double *a, double delta, int *e) { + double gamma, double tau, SEXP df, double *a, double delta, int *e, double *maxChange) { // Pre-calculcate v int K = K1[g+1] - K1[g]; double xwr, u; double *v = Calloc(K, double); + double shift; for (int j=K1[g]; j maxChange[0]) maxChange[0] = fabs(shift); + for (int i=0; i maxChange[0]) maxChange[0] = fabs(shift); for (int i=0; i maxChange) maxChange = fabs(shift); b[l*p+j] = shift + a[j]; for (int i=0; i #include #include -int checkConvergence(double *beta, double *beta_old, double eps, int l, int J); double crossprod(double *x, double *y, int n, int j); double norm(double *x, int p); double S(double z, double l); @@ -15,11 +14,12 @@ double dMCP(double theta, double l, double a); double gLoss(double *r, int n); // Groupwise local coordinate descent updates -void gLCD_gaussian(double *b, const char *penalty, double *x, double *r, int g, int *K1, int n, int l, int p, double lam1, double lam2, double gamma, double tau, SEXP df, double *a, double delta, int *e) -{ +void gLCD_gaussian(double *b, const char *penalty, double *x, double *r, int g, int *K1, int n, int l, int p, double lam1, double lam2, double gamma, double tau, SEXP df, double *a, double delta, int *e, double *maxChange) { + // Make initial local approximation int K = K1[g+1] - K1[g]; double sG = 0; // Sum of inner penalties for group + double shift; if (strcmp(penalty, "gel")==0) for (int j=K1[g]; j maxChange[0]) maxChange[0] = fabs(shift); + for (int i=0; i maxChange[0]) maxChange[0] = fabs(shift); for (int i=0; i maxChange) maxChange = fabs(shift); b[l*p+j] = shift + a[j]; for (int i=0; i maxChange[0]) maxChange[0] = fabs(shift); for (int i=0; i maxChange[0]) maxChange[0] = fabs(shift); if (shift != 0) { for (int i=0; i maxChange) maxChange = fabs(shift); b[l*p+j] = shift + a[j]; for (int i=0; i