source: orange/orange/Orange/regression/lasso.py @ 9521:4a6752f76713

Revision 9521:4a6752f76713, 13.7 KB checked in by lan <lan@…>, 2 years ago (diff)

Updated printing model: function print_lasso_regression_model deleted; instead use
print model

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    >>> print c
18   
19      Variable  Coeff Est  Std Error          p
20     Intercept     22.533
21          CRIM     -0.044      0.030      0.510     
22            ZN      0.013      0.010      0.660     
23         INDUS     -0.003      0.023      0.980     
24          CHAS      2.318      1.304      0.200     
25           NOX     -7.530      2.803      0.370     
26            RM      4.231      0.819      0.000   ***
27           DIS     -0.710      0.130      0.070     .
28           RAD      0.074      0.029      0.510     
29           TAX     -0.004      0.002      0.560     
30       PTRATIO     -0.821      0.095      0.000   ***
31             B      0.007      0.002      0.170     
32         LSTAT     -0.503      0.085      0.000   ***
33    Signif. codes:  0 *** 0.001 ** 0.01 * 0.05 . 0.1 empty 1
34
35
36    For 1 variable the regression coefficient equals 0:
37    AGE
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       
336    @deprecated_keywords({"skipZero": "skip_zero"})
337    def to_string(self, skip_zero=True):
338        """Pretty-prints Lasso regression model,
339        i.e. estimated regression coefficients with standard errors
340        and significances. Standard errors are obtained using bootstrapping
341        method and significances by the permuation test
342
343        :param skip_zero: if True variables with estimated coefficient equal to 0
344            are omitted
345        :type skip_zero: boolean
346        """
347       
348        from string import join
349        labels = ('Variable', 'Coeff Est', 'Std Error', 'p')
350        lines = [join(['%10s' % l for l in labels], ' ')]
351
352        fmt = "%10s " + join(["%10.3f"]*3, " ") + " %5s"
353        fmt1 = "%10s %10.3f"
354
355        def get_star(p):
356            if p < 0.001: return  "*"*3
357            elif p < 0.01: return "*"*2
358            elif p < 0.05: return "*"
359            elif p < 0.1: return  "."
360            else: return " "
361
362        stars =  get_star(self.p_vals[0])
363        lines.append(fmt1 % ('Intercept', self.coef0))
364        skipped = []
365        for i in range(len(self.domain.attributes)):
366            if self.coefficients[i] == 0. and skip_zero:
367                skipped.append(self.domain.attributes[i].name)
368                continue           
369            stars = get_star(self.p_vals[i])
370            lines.append(fmt % (self.domain.attributes[i].name, 
371                         self.coefficients[i], self.std_errors_fixed_t[i], 
372                         self.p_vals[i], stars))
373        lines.append("Signif. codes:  0 *** 0.001 ** 0.01 * 0.05 . 0.1 empty 1")
374        lines.append("\n")
375        if skip_zero:
376            k = len(skipped)
377            if k == 0:
378                lines.append("All variables have non-zero regression coefficients. ")
379            else:
380                suff = "s" if k > 1 else ""
381                lines.append("For %d variable%s the regression coefficient equals 0: " \
382                      % (k, suff))
383                for var in skipped:
384                    lines.append(var)
385        return "\n".join(lines)
386
387    def __str__(self):
388        return self.to_string(skip_zero=True)
389
390deprecated_members({"muX": "mu_x",
391                    "stdErrorsFixedT": "std_errors_fixed_t",
392                    "pVals": "p_vals",
393                    "dictModel": "dict_model"},
394                   wrap_methods=["__init__"],
395                   in_place=True)(LassoRegression)
396
397if __name__ == "__main__":
398
399    import Orange
400   
401    table = Orange.data.Table("housing.tab")       
402
403    c = LassoRegressionLearner(table, t=len(table.domain))
404    print c
Note: See TracBrowser for help on using the repository browser.