source: orange/orange/Orange/regression/lasso.py @ 9290:c96f47743746

Revision 9290:c96f47743746, 13.6 KB checked in by ales_erjavec <ales.erjavec@…>, 2 years ago (diff)

Normalize X by standard deviation before computing the lasso solution.
Added 's' parameter as an alternative way to specify 't'.

Line 
1"""\
2============================
3Lasso regression (``lasso``)
4============================
5
6.. index:: regression
7
8.. _`Lasso regression. Regression shrinkage and selection via the lasso`:
9    http://www-stat.stanford.edu/~tibs/lasso/lasso.pdf
10
11
12Example ::
13
14    >>> from Orange.regression import lasso
15    >>> table = Orange.data.Table("housing")
16    >>> c = lasso.LassoRegressionLearner(table)
17    >>> linear.print_lasso_regression_model(c)
18   
19      Variable  Coeff Est  Std Error          p
20     Intercept     22.533
21          CRIM     -0.049      0.282      0.770     
22            ZN      0.106      0.055      0.030     *
23         INDUS     -0.111      0.442      0.920     
24          CHAS      1.757      0.669      0.180     
25           NOX      0.318      0.483      0.680     
26            RM      1.643      0.461      0.480     
27           AGE      0.062      0.051      0.230     
28           DIS      0.627      0.538      0.930     
29           RAD      1.260      0.472      0.070     .
30           TAX     -0.074      0.027      0.120     
31       PTRATIO      1.331      0.464      0.050     .
32             B      0.017      0.007      0.080     .
33         LSTAT     -0.209      0.323      0.650     
34    Signif. codes:  0 *** 0.001 ** 0.01 * 0.05 . 0.1 empty 1
35
36
37    All variables have non-zero regression coefficients.
38       
39    >>>
40
41
42.. autoclass:: LassoRegressionLearner
43    :members:
44
45.. autoclass:: LassoRegression
46    :members:
47
48Utility functions
49-----------------
50
51.. autofunction:: center
52
53.. autofunction:: get_bootstrap_sample
54
55.. autofunction:: permute_responses
56
57.. autofunction:: print_lasso_regression_model
58
59
60"""
61
62import Orange
63import numpy
64
65from Orange.regression import base
66
67from Orange.misc import deprecated_members, deprecated_keywords
68
69def center(X):
70    """Centers the data, i.e. subtracts the column means.
71    Returns the centered data and the mean.
72
73    :param X: the data arry
74    :type table: :class:`numpy.array`
75    """
76    mu = X.mean(axis=0)
77    return X - mu, mu
78
79def standardize(X):
80    """Standardizes the data, i.e. subtracts the column means and divide by
81    standard deviation.
82    Returns the centered data, the mean, and deviations.
83
84    :param X: the data arry
85    :type table: :class:`numpy.array`
86    """
87    mu = numpy.mean(X, axis=0)
88    std = numpy.std(X, axis=0)
89    return (X - mu) / std, mu, std
90
91def get_bootstrap_sample(table):
92    """Generates boostrap sample from an Orange Example Table
93    and stores it in a new :class:`Orange.data.Table` object
94
95    :param table: the original data sample
96    :type table: :class:`Orange.data.Table`
97    """
98    n = len(table)
99    bootTable = Orange.data.Table(table.domain)
100    for i in range(n):
101        id = numpy.random.randint(0,n)
102        bootTable.append(table[id])
103    return bootTable
104
105def permute_responses(table):
106    """ Permutes values of the class (response) variable.
107    The independence between independent variables and the response
108    is obtained but the distribution of the response variable is kept.
109
110    :param table: the original data sample
111    :type table: :class:`Orange.data.Table`
112    """
113    n = len(table)
114    perm = numpy.random.permutation(n)
115    permTable = Orange.data.Table(table.domain, table)
116    for i, ins in enumerate(table):
117        permTable[i].set_class(table[perm[i]].get_class())
118    return permTable
119
120class LassoRegressionLearner(base.BaseRegressionLearner):
121    """Fits the lasso regression model, i.e. learns the regression parameters
122    The class is derived from
123    :class:`Orange.regression.base.BaseRegressionLearner`
124    which is used for preprocessing the data (continuization and imputation)
125    before fitting the regression parameters
126
127    """
128
129    def __init__(self, name='lasso regression', t=1, s=None, tol=0.001, \
130                 n_boot=100, n_perm=100, imputer=None, continuizer=None):
131        """
132        :param name: name of the linear model, default 'lasso regression'
133        :type name: string
134       
135        :param t: tuning parameter, upper bound for the L1-norm of the
136            regression coefficients
137        :type t: float
138       
139        :param s: An alternative way to specify the tuning parameter ``t``.
140            Here ``t`` is taken to be t = s * sum(abs(B)) where B are the
141            coefficients of an ordinary least square linear fit. ``t`` parameter is ignored if ``s`` is specified (by default it
142            is None).
143        :type s: float
144       
145        :param tol: tolerance parameter, regression coefficients
146            (absoulute value) under tol are set to 0,
147            default=0.001
148        :type tol: float
149       
150        :param n_boot: number of bootstrap samples used for non-parametric
151            estimation of standard errors
152        :type n_boot: int
153       
154        :param n_perm: number of permuations used for non-parametric
155            estimation of p-values
156        :type n_perm: int
157       
158        """
159
160        self.name = name
161        self.t = t
162        self.s = s
163        self.tol = tol
164        self.n_boot = n_boot
165        self.n_perm = n_perm
166        self.set_imputer(imputer=imputer)
167        self.set_continuizer(continuizer=continuizer)
168       
169       
170    def __call__(self, table, weight=None):
171        """
172        :param table: data instances.
173        :type table: :class:`Orange.data.Table`
174        :param weight: the weights for instances. Default: None, i.e.
175            all data instances are eqaully important in fitting
176            the regression parameters
177        :type weight: None or list of Orange.data.variable.Continuous
178            which stores weights for instances
179       
180        """ 
181        # dicrete values are continuized       
182        table = self.continuize_table(table)
183        # missing values are imputed
184        table = self.impute_table(table)
185
186        domain = table.domain
187        X, y, w = table.to_numpy()
188        n, m = numpy.shape(X)
189       
190        X, mu_x, sigma_x = standardize(X)
191        y, coef0 = center(y)
192       
193        t = self.t
194       
195        if self.s is not None:
196            beta_full, rss, _, _ = numpy.linalg.lstsq(X, y)
197            t = self.s * numpy.sum(numpy.abs(beta_full))
198            print "t =", t
199           
200        import scipy.optimize
201           
202        # objective function to be minimized
203        objective = lambda beta: numpy.linalg.norm(y - numpy.dot(X, beta))
204        # initial guess for the regression parameters
205        beta_init = numpy.random.random(m)
206        # constraints for the regression coefficients
207        cnstr = lambda beta: t - numpy.sum(numpy.abs(beta))
208        # optimal solution
209        coefficients = scipy.optimize.fmin_cobyla(objective, beta_init,\
210                                                       cnstr, disp=0)
211
212        # set small coefficients to 0
213        def set_2_0(c): return c if abs(c) > self.tol else 0
214        coefficients = numpy.array(map(set_2_0, coefficients))
215        coefficients /= sigma_x
216       
217        # bootstrap estimator of standard error of the coefficient estimators
218        # assumption: fixed t
219        if self.n_boot > 0:
220            coeff_b = [] # bootstrapped coefficients
221            for i in range(self.n_boot):
222                tmp_table = get_bootstrap_sample(table)
223                l = LassoRegressionLearner(t=t, n_boot=0, n_perm=0)
224                c = l(tmp_table)
225                coeff_b.append(c.coefficients)
226            std_errors_fixed_t = numpy.std(coeff_b, axis=0)
227        else:
228            std_errors_fixed_t = [float("nan")] * m
229
230        # permutation test to obtain the significance of the regression
231        #coefficients
232        if self.n_perm > 0:
233            coeff_p = []
234            for i in range(self.n_perm):
235                tmp_table = permute_responses(table)
236                l = LassoRegressionLearner(t=t, n_boot=0, n_perm=0)
237                c = l(tmp_table)
238                coeff_p.append(c.coefficients)
239            p_vals = \
240                   numpy.sum(abs(numpy.array(coeff_p))>\
241                             abs(numpy.array(coefficients)), \
242                             axis=0)/float(self.n_perm)
243        else:
244            p_vals = [float("nan")] * m
245
246        # dictionary of regression coefficients with standard errors
247        # and p-values
248        dict_model = {}
249        for i, var in enumerate(domain.attributes):
250            dict_model[var.name] = (coefficients[i], std_errors_fixed_t[i], p_vals[i])           
251       
252        return LassoRegression(domain=domain, class_var=domain.class_var,
253                               coef0=coef0, coefficients=coefficients,
254                               std_errors_fixed_t=std_errors_fixed_t,
255                               p_vals=p_vals,
256                               dict_model= dict_model,
257                               mu_x=mu_x)
258
259deprecated_members({"nBoot": "n_boot",
260                    "nPerm": "n_perm"}, 
261                   wrap_methods=["__init__"],
262                   in_place=True)(LassoRegressionLearner)
263
264class LassoRegression(Orange.classification.Classifier):
265    """Lasso regression predicts value of the response variable
266    based on the values of independent variables.
267
268    .. attribute:: coef0
269
270        intercept (sample mean of the response variable)   
271
272    .. attribute:: coefficients
273
274        list of regression coefficients.
275
276    .. attribute:: std_errors_fixed_t
277
278        list of standard errors of the coefficient estimator for the fixed
279        tuning parameter t. The standard errors are estimated using
280        bootstrapping method.
281
282    .. attribute:: p_vals
283
284        list of p-values for the null hypothesis that the regression
285        coefficients equal 0 based on non-parametric permutation test
286
287    .. attribute:: dict_model
288
289        dictionary of statistical properties of the model.
290        Keys - names of the independent variables
291        Values - tuples (coefficient, standard error, p-value)
292
293    .. attribute:: mu_x
294
295        the sample mean of the all independent variables   
296
297    """ 
298    def __init__(self, domain=None, class_var=None, coef0=None,
299                 coefficients=None, std_errors_fixed_t=None, p_vals=None,
300                 dict_model=None, mu_x=None):
301        self.domain = domain
302        self.class_var = class_var
303        self.coef0 = coef0
304        self.coefficients = coefficients
305        self.std_errors_fixed_t = std_errors_fixed_t
306        self.p_vals = p_vals
307        self.dict_model = dict_model
308        self.mu_x = mu_x
309
310    @deprecated_keywords({"resultType": "result_type"})
311    def __call__(self, instance, result_type=Orange.core.GetValue):
312        """
313        :param instance: data instance for which the value of the response
314            variable will be predicted
315        :type instance:
316        """ 
317        ins = Orange.data.Instance(self.domain, instance)
318        if "?" in ins: # missing value -> corresponding coefficient omitted
319            def miss_2_0(x): return x if x != "?" else 0
320            ins = map(miss_2_0, ins)
321            ins = numpy.array(ins)[:-1] - self.mu_x
322        else:
323            ins = numpy.array(ins.native())[:-1] - self.mu_x
324
325        y_hat = numpy.dot(self.coefficients, ins) + self.coef0
326        y_hat = self.class_var(y_hat)
327        dist = Orange.statistics.distribution.Continuous(self.class_var)
328        dist[y_hat] = 1.0
329        if result_type == Orange.core.GetValue:
330            return y_hat
331        if result_type == Orange.core.GetProbabilities:
332            return dist
333        else:
334            return (y_hat, dist)   
335
336deprecated_members({"muX": "mu_x",
337                    "stdErrorsFixedT": "std_errors_fixed_t",
338                    "pVals": "p_vals",
339                    "dictModel": "dict_model"},
340                   wrap_methods=["__init__"],
341                   in_place=True)(LassoRegression)
342
343
344@deprecated_keywords({"skipZero": "skip_zero"})
345def print_lasso_regression_model(lr, skip_zero=True):
346    """Pretty-prints Lasso regression model,
347    i.e. estimated regression coefficients with standard errors
348    and significances. Standard errors are obtained using bootstrapping
349    method and significances by the permuation test
350
351    :param lr: a Lasso regression model object.
352    :type lr: :class:`LassoRegression`
353    :param skip_zero: if True variables with estimated coefficient equal to 0
354        are omitted
355    :type skip_zero: boolean
356    """
357   
358    from string import join
359    m = lr
360    labels = ('Variable', 'Coeff Est', 'Std Error', 'p')
361    print join(['%10s' % l for l in labels], ' ')
362
363    fmt = "%10s " + join(["%10.3f"]*3, " ") + " %5s"
364    fmt1 = "%10s %10.3f"
365
366    def get_star(p):
367        if p < 0.001: return  "*"*3
368        elif p < 0.01: return "*"*2
369        elif p < 0.05: return "*"
370        elif p < 0.1: return  "."
371        else: return " "
372
373    stars =  get_star(m.p_vals[0])
374    print fmt1 % ('Intercept', m.coef0)
375    skipped = []
376    for i in range(len(m.domain.attributes)):
377        if m.coefficients[i] == 0. and skip_zero:
378            skipped.append(m.domain.attributes[i].name)
379            continue           
380        stars = get_star(m.p_vals[i])
381        print fmt % (m.domain.attributes[i].name, \
382                     m.coefficients[i], m.std_errors_fixed_t[i], \
383                     m.p_vals[i], stars)
384    print "Signif. codes:  0 *** 0.001 ** 0.01 * 0.05 . 0.1 empty 1"
385    print "\n"
386    if skip_zero:
387        k = len(skipped)
388        if k == 0:
389            print "All variables have non-zero regression coefficients. "
390        else:
391            suff = "s" if k > 1 else ""
392            print "For %d variable%s the regression coefficient equals 0: " \
393                  % (k, suff)
394            for var in skipped:
395                print var
396
397
398if __name__ == "__main__":
399
400    import Orange
401   
402    table = Orange.data.Table("housing.tab")       
403
404    c = LassoRegressionLearner(table, t=len(table.domain))
405    print_lasso_regression_model(c)
Note: See TracBrowser for help on using the repository browser.