harmonization.elastic_net
1import numpy as np 2import warnings 3 4from sklearn.model_selection import KFold 5from harmonization._glmnet import elnet 6 7 8def lasso_path(X, y, nlam=100, fit_intercept=False, pos=False, standardize=False, weights=None, penalty=None, criterion=None): 9 """return full path for Lasso""" 10 return elastic_net_path(X, y, rho=1.0, nlam=nlam, weights=weights, penalty=penalty, criterion=criterion, 11 fit_intercept=fit_intercept, pos=pos, standardize=standardize) 12 13 14def lasso_crossval(X, y, nlam=100, fit_intercept=False, pos=False, standardize=False, 15 weights=None, penalty=None, n_splits=3, use_min_se=True): 16 """return k-fold cross validation for Lasso""" 17 18 # In this call we find a list of lambda values 19 # we then crossval using those same values on different splits down later. 20 lmu, a0, ca, ia, nin, rsq, alm, nlp, jerr = elastic_net(X, 21 y, 22 1., 23 nlam=nlam, 24 fit_intercept=fit_intercept, 25 pos=pos, 26 weights=weights, 27 vp=penalty, 28 standardize=standardize) 29 alm = alm[:lmu] 30 cv_error = np.zeros((n_splits, lmu)) 31 32 for idx, (train_idx, test_idx) in enumerate(KFold(n_splits=n_splits, shuffle=True).split(X)): 33 beta, a0, _, _ = elastic_net_path(X[train_idx], 34 y[train_idx], 35 rho=1.0, 36 ulam=alm, 37 nlam=lmu, 38 penalty=penalty, 39 fit_intercept=fit_intercept, 40 pos=pos, 41 standardize=standardize) 42 43 y_pred = np.dot(X[test_idx], beta) + a0 44 error = np.mean((y[test_idx, None] - y_pred)**2, axis=0) 45 # seems like the lambda sequence can sometimes be truncated somehow? 46 cv_error[idx, :len(error)] = error 47 48 mean_cv = np.mean(cv_error, axis=0) 49 best_idx = np.argmin(mean_cv) 50 51 use_min_se = False 52 if use_min_se: 53 std_cv = np.std(cv_error, axis=0, ddof=1) 54 semin = mean_cv[best_idx] + std_cv[best_idx] 55 56 # we want the smallest cv + std_cv as the new model instead of just the smallest one 57 # that's the length of the vector where everything is below the upper bound on cv 58 best_idx = np.sum(mean_cv <= semin) 59 60 best_lambda = alm[best_idx] 61 62 # redo a full fit with the best lambda value we found 63 beta_best, a0_best, Xhat_best, alm_best = elastic_net_path(X, 64 y, 65 rho=1.0, 66 ulam=best_lambda, 67 nlam=1, 68 penalty=penalty, 69 fit_intercept=fit_intercept, 70 pos=pos, 71 standardize=standardize) 72 73 return Xhat_best.squeeze(), beta_best.squeeze(), a0_best, alm_best 74 75 76def elastic_net_path(X, y, rho, nlam=100, ulam=None, criterion=None, variance=None, 77 fit_intercept=False, pos=False, standardize=False, weights=None, penalty=None): 78 """return full path for ElasticNet""" 79 80 lmu, a0, ca, ia, nin, rsq, alm, nlp, jerr = elastic_net(X, 81 y, 82 rho, 83 nlam=nlam, 84 fit_intercept=fit_intercept, 85 pos=pos, 86 ulam=ulam, 87 weights=weights, 88 vp=penalty, 89 standardize=standardize) 90 nobs, nx = X.shape 91 92 if jerr > 0: # jerr < 0 is a non fatal error 93 warnings.warn('Non zero error code {}'.format(jerr)) 94 95 a0 = a0[:lmu] 96 ca = ca[:nx, :lmu] 97 ia = ia[:nx] 98 nin = nin[:lmu] 99 rsq = rsq[:lmu] 100 alm = alm[:lmu] 101 102 if len(nin) == 0: 103 ninmax = None 104 else: 105 ninmax = max(nin) 106 107 if ninmax is None: 108 return (np.zeros([nx, nlam], dtype=np.float32), 109 np.zeros((nlam), dtype=np.float32), 110 np.zeros([nobs, nlam], dtype=np.float32), 111 np.zeros((nlam), dtype=np.float32)) 112 113 ca = ca[:ninmax] 114 ja = ia[:ninmax] - 1 # ia is 1-indexed in fortran 115 oja = np.argsort(ja) 116 ja1 = ja[oja] 117 beta = np.zeros([nx, lmu], dtype=np.float32) 118 beta[ja1] = ca[oja] 119 120 yhat = np.dot(X, beta) + a0 121 122 # no criterion - just return the whole path 123 # else we choose the best value of lambda and just return that 124 if criterion is None: 125 return beta, a0, yhat, alm 126 127 n = y.shape[0] 128 p = X.shape[1] 129 130 if criterion == 'aic' or criterion == 'aicc': 131 w = 2 132 elif criterion == 'bic': 133 w = np.log(n) 134 elif criterion == 'ric': 135 w = 2 * np.log(p) 136 else: 137 raise ValueError('Criterion {} is not supported!'.format(criterion)) 138 139 mse = np.mean((y[:, None] - yhat)**2, axis=0, dtype=np.float32) 140 df_mu = np.sum(beta != 0, axis=0, dtype=np.float32) 141 142 # criterion from burnham and anderson 143 criterion_value = np.array(n * np.log(mse) + w * df_mu, dtype=np.float64) 144 145 # criterion according to 2.15 and 2.16 of https://projecteuclid.org/download/pdfview_1/euclid.aos/1194461726 146 # if variance is None: 147 # variance = np.var(y) 148 # criterion_value = mse / variance + w * df_mu / n 149 150 if criterion == 'aicc': 151 with np.errstate(divide='ignore'): 152 aicc = (2 * df_mu * (df_mu + 1)) / (n - df_mu - 1) 153 criterion_value += aicc 154 155 criterion_value[df_mu == 0] = 1e300 156 criterion_value[np.isnan(criterion_value)] = 1e300 157 best_idx = np.argmin(criterion_value, axis=0) 158 159 # Seems like the memory is not correctly freed internally when we call from fortran, 160 # so we copy everything once to not keep references to it during the parallel processing 161 # or whatever, computer sciency stuff, it works :/ 162 beta = beta[:, best_idx].copy() 163 yhat = yhat[:, best_idx].copy() 164 a0 = a0[best_idx].copy() 165 alm = alm[best_idx].copy() 166 167 return beta, a0, yhat, alm 168 169 170def elastic_net(X, y, rho, pos=False, thr=1e-7, weights=None, vp=None, copy=True, ulam=None, jd=np.zeros(1), 171 standardize=False, nlam=100, maxit=100000, flmin=1e-4, fit_intercept=False, custom_path=False): 172 """ 173 Raw-output wrapper for elastic net linear regression. 174 X is D 175 y is X 176 rho for lasso/elastic net tradeoff 177 """ 178 179 # X/y is overwritten in the fortran function at every loop, so we must copy it each time 180 if copy: 181 X = np.array(X, copy=True, dtype=np.float64, order='F') 182 y = np.array(y, copy=True, dtype=np.float64, order='F').ravel() 183 184 box_constraints = np.zeros((2, X.shape[1]), order='F') 185 box_constraints[1] = 9.9e35 186 187 if not pos: 188 box_constraints[0] = -9.9e35 189 190 # Uniform weighting if no weights are specified. 191 if weights is None: 192 weights = np.ones(X.shape[0], order='F') 193 else: 194 weights = weights.copy(order='F') 195 196 # Uniform penalties if none were specified. 197 if vp is None: 198 vp = np.ones(X.shape[1], order='F') 199 else: 200 vp = vp.copy(order='F') 201 202 if (X.shape[1] > 500) or (X.shape[1] >= 2*X.shape[0]): 203 # use naive algorithm 204 ka = 2 205 else: 206 # use covariance algorithm 207 ka = 1 208 209 if custom_path: 210 ulam = np.logspace(2, -5, num=nlam) 211 212 # User activated values are used only if flmin >= 1. 213 if ulam is not None: 214 flmin = 2. 215 216 nvars = X.shape[1] 217 ne = nvars + 1 218 nx = int(min(ne * 2 + 20, X.shape[1])) 219 220 if y.ndim == 1: 221 lmu, a0, ca, ia, nin, rsq, alm, nlp, jerr = elnet(rho, 222 X, 223 y, 224 weights, 225 jd, 226 vp, 227 box_constraints, 228 nx, 229 flmin, 230 ulam, 231 thr, 232 nlam=nlam, 233 ka=ka, 234 # ne=ne, 235 isd=standardize, 236 maxit=maxit, 237 intr=fit_intercept) 238 else: 239 raise ValueError('y needs to be 1D but is {}'.format(y.ndim)) 240 241 return lmu, a0, ca, ia, nin, rsq, alm, nlp, jerr
def
lasso_path( X, y, nlam=100, fit_intercept=False, pos=False, standardize=False, weights=None, penalty=None, criterion=None):
9def lasso_path(X, y, nlam=100, fit_intercept=False, pos=False, standardize=False, weights=None, penalty=None, criterion=None): 10 """return full path for Lasso""" 11 return elastic_net_path(X, y, rho=1.0, nlam=nlam, weights=weights, penalty=penalty, criterion=criterion, 12 fit_intercept=fit_intercept, pos=pos, standardize=standardize)
return full path for Lasso
def
lasso_crossval( X, y, nlam=100, fit_intercept=False, pos=False, standardize=False, weights=None, penalty=None, n_splits=3, use_min_se=True):
15def lasso_crossval(X, y, nlam=100, fit_intercept=False, pos=False, standardize=False, 16 weights=None, penalty=None, n_splits=3, use_min_se=True): 17 """return k-fold cross validation for Lasso""" 18 19 # In this call we find a list of lambda values 20 # we then crossval using those same values on different splits down later. 21 lmu, a0, ca, ia, nin, rsq, alm, nlp, jerr = elastic_net(X, 22 y, 23 1., 24 nlam=nlam, 25 fit_intercept=fit_intercept, 26 pos=pos, 27 weights=weights, 28 vp=penalty, 29 standardize=standardize) 30 alm = alm[:lmu] 31 cv_error = np.zeros((n_splits, lmu)) 32 33 for idx, (train_idx, test_idx) in enumerate(KFold(n_splits=n_splits, shuffle=True).split(X)): 34 beta, a0, _, _ = elastic_net_path(X[train_idx], 35 y[train_idx], 36 rho=1.0, 37 ulam=alm, 38 nlam=lmu, 39 penalty=penalty, 40 fit_intercept=fit_intercept, 41 pos=pos, 42 standardize=standardize) 43 44 y_pred = np.dot(X[test_idx], beta) + a0 45 error = np.mean((y[test_idx, None] - y_pred)**2, axis=0) 46 # seems like the lambda sequence can sometimes be truncated somehow? 47 cv_error[idx, :len(error)] = error 48 49 mean_cv = np.mean(cv_error, axis=0) 50 best_idx = np.argmin(mean_cv) 51 52 use_min_se = False 53 if use_min_se: 54 std_cv = np.std(cv_error, axis=0, ddof=1) 55 semin = mean_cv[best_idx] + std_cv[best_idx] 56 57 # we want the smallest cv + std_cv as the new model instead of just the smallest one 58 # that's the length of the vector where everything is below the upper bound on cv 59 best_idx = np.sum(mean_cv <= semin) 60 61 best_lambda = alm[best_idx] 62 63 # redo a full fit with the best lambda value we found 64 beta_best, a0_best, Xhat_best, alm_best = elastic_net_path(X, 65 y, 66 rho=1.0, 67 ulam=best_lambda, 68 nlam=1, 69 penalty=penalty, 70 fit_intercept=fit_intercept, 71 pos=pos, 72 standardize=standardize) 73 74 return Xhat_best.squeeze(), beta_best.squeeze(), a0_best, alm_best
return k-fold cross validation for Lasso
def
elastic_net_path( X, y, rho, nlam=100, ulam=None, criterion=None, variance=None, fit_intercept=False, pos=False, standardize=False, weights=None, penalty=None):
77def elastic_net_path(X, y, rho, nlam=100, ulam=None, criterion=None, variance=None, 78 fit_intercept=False, pos=False, standardize=False, weights=None, penalty=None): 79 """return full path for ElasticNet""" 80 81 lmu, a0, ca, ia, nin, rsq, alm, nlp, jerr = elastic_net(X, 82 y, 83 rho, 84 nlam=nlam, 85 fit_intercept=fit_intercept, 86 pos=pos, 87 ulam=ulam, 88 weights=weights, 89 vp=penalty, 90 standardize=standardize) 91 nobs, nx = X.shape 92 93 if jerr > 0: # jerr < 0 is a non fatal error 94 warnings.warn('Non zero error code {}'.format(jerr)) 95 96 a0 = a0[:lmu] 97 ca = ca[:nx, :lmu] 98 ia = ia[:nx] 99 nin = nin[:lmu] 100 rsq = rsq[:lmu] 101 alm = alm[:lmu] 102 103 if len(nin) == 0: 104 ninmax = None 105 else: 106 ninmax = max(nin) 107 108 if ninmax is None: 109 return (np.zeros([nx, nlam], dtype=np.float32), 110 np.zeros((nlam), dtype=np.float32), 111 np.zeros([nobs, nlam], dtype=np.float32), 112 np.zeros((nlam), dtype=np.float32)) 113 114 ca = ca[:ninmax] 115 ja = ia[:ninmax] - 1 # ia is 1-indexed in fortran 116 oja = np.argsort(ja) 117 ja1 = ja[oja] 118 beta = np.zeros([nx, lmu], dtype=np.float32) 119 beta[ja1] = ca[oja] 120 121 yhat = np.dot(X, beta) + a0 122 123 # no criterion - just return the whole path 124 # else we choose the best value of lambda and just return that 125 if criterion is None: 126 return beta, a0, yhat, alm 127 128 n = y.shape[0] 129 p = X.shape[1] 130 131 if criterion == 'aic' or criterion == 'aicc': 132 w = 2 133 elif criterion == 'bic': 134 w = np.log(n) 135 elif criterion == 'ric': 136 w = 2 * np.log(p) 137 else: 138 raise ValueError('Criterion {} is not supported!'.format(criterion)) 139 140 mse = np.mean((y[:, None] - yhat)**2, axis=0, dtype=np.float32) 141 df_mu = np.sum(beta != 0, axis=0, dtype=np.float32) 142 143 # criterion from burnham and anderson 144 criterion_value = np.array(n * np.log(mse) + w * df_mu, dtype=np.float64) 145 146 # criterion according to 2.15 and 2.16 of https://projecteuclid.org/download/pdfview_1/euclid.aos/1194461726 147 # if variance is None: 148 # variance = np.var(y) 149 # criterion_value = mse / variance + w * df_mu / n 150 151 if criterion == 'aicc': 152 with np.errstate(divide='ignore'): 153 aicc = (2 * df_mu * (df_mu + 1)) / (n - df_mu - 1) 154 criterion_value += aicc 155 156 criterion_value[df_mu == 0] = 1e300 157 criterion_value[np.isnan(criterion_value)] = 1e300 158 best_idx = np.argmin(criterion_value, axis=0) 159 160 # Seems like the memory is not correctly freed internally when we call from fortran, 161 # so we copy everything once to not keep references to it during the parallel processing 162 # or whatever, computer sciency stuff, it works :/ 163 beta = beta[:, best_idx].copy() 164 yhat = yhat[:, best_idx].copy() 165 a0 = a0[best_idx].copy() 166 alm = alm[best_idx].copy() 167 168 return beta, a0, yhat, alm
return full path for ElasticNet
def
elastic_net( X, y, rho, pos=False, thr=1e-07, weights=None, vp=None, copy=True, ulam=None, jd=array([0.]), standardize=False, nlam=100, maxit=100000, flmin=0.0001, fit_intercept=False, custom_path=False):
171def elastic_net(X, y, rho, pos=False, thr=1e-7, weights=None, vp=None, copy=True, ulam=None, jd=np.zeros(1), 172 standardize=False, nlam=100, maxit=100000, flmin=1e-4, fit_intercept=False, custom_path=False): 173 """ 174 Raw-output wrapper for elastic net linear regression. 175 X is D 176 y is X 177 rho for lasso/elastic net tradeoff 178 """ 179 180 # X/y is overwritten in the fortran function at every loop, so we must copy it each time 181 if copy: 182 X = np.array(X, copy=True, dtype=np.float64, order='F') 183 y = np.array(y, copy=True, dtype=np.float64, order='F').ravel() 184 185 box_constraints = np.zeros((2, X.shape[1]), order='F') 186 box_constraints[1] = 9.9e35 187 188 if not pos: 189 box_constraints[0] = -9.9e35 190 191 # Uniform weighting if no weights are specified. 192 if weights is None: 193 weights = np.ones(X.shape[0], order='F') 194 else: 195 weights = weights.copy(order='F') 196 197 # Uniform penalties if none were specified. 198 if vp is None: 199 vp = np.ones(X.shape[1], order='F') 200 else: 201 vp = vp.copy(order='F') 202 203 if (X.shape[1] > 500) or (X.shape[1] >= 2*X.shape[0]): 204 # use naive algorithm 205 ka = 2 206 else: 207 # use covariance algorithm 208 ka = 1 209 210 if custom_path: 211 ulam = np.logspace(2, -5, num=nlam) 212 213 # User activated values are used only if flmin >= 1. 214 if ulam is not None: 215 flmin = 2. 216 217 nvars = X.shape[1] 218 ne = nvars + 1 219 nx = int(min(ne * 2 + 20, X.shape[1])) 220 221 if y.ndim == 1: 222 lmu, a0, ca, ia, nin, rsq, alm, nlp, jerr = elnet(rho, 223 X, 224 y, 225 weights, 226 jd, 227 vp, 228 box_constraints, 229 nx, 230 flmin, 231 ulam, 232 thr, 233 nlam=nlam, 234 ka=ka, 235 # ne=ne, 236 isd=standardize, 237 maxit=maxit, 238 intr=fit_intercept) 239 else: 240 raise ValueError('y needs to be 1D but is {}'.format(y.ndim)) 241 242 return lmu, a0, ca, ia, nin, rsq, alm, nlp, jerr
Raw-output wrapper for elastic net linear regression. X is D y is X rho for lasso/elastic net tradeoff