source: orange/orange/Orange/regression/linear.py @ 9243:cc2ddc8cc69d

Revision 9243:cc2ddc8cc69d, 19.2 KB checked in by ales_erjavec <ales.erjavec@…>, 2 years ago (diff)

Use statc for betai, fprob if scipy is not installed.
Fixed KeyboardInterupt handling in stepwise function.

Line 
1"""\
2==============================
3Linear regression (``linear``)
4==============================
5
6.. index:: regression, linear model
7
8.. _`Linear regression`: http://en.wikipedia.org/wiki/Linear_regression
9
10Example ::
11
12    >>> from Orange.regression import linear
13    >>> table = Orange.data.Table("housing")
14    >>> c = linear.LinearRegressionLearner(table)
15    >>> linear.print_linear_regression_model(c)
16   
17      Variable  Coeff Est  Std Error    t-value          p
18     Intercept     36.459      5.103      7.144      0.000   ***
19          CRIM     -0.108      0.033     -3.287      0.001    **
20            ZN      0.046      0.014      3.382      0.001   ***
21         INDUS      0.021      0.061      0.334      0.738     
22          CHAS      2.687      0.862      3.118      0.002    **
23           NOX    -17.767      3.820     -4.651      0.000   ***
24            RM      3.810      0.418      9.116      0.000   ***
25           AGE      0.001      0.013      0.052      0.958     
26           DIS     -1.476      0.199     -7.398      0.000   ***
27           RAD      0.306      0.066      4.613      0.000   ***
28           TAX     -0.012      0.004     -3.280      0.001    **
29       PTRATIO     -0.953      0.131     -7.283      0.000   ***
30             B      0.009      0.003      3.467      0.001   ***
31         LSTAT     -0.525      0.051    -10.347      0.000   ***
32    Signif. codes:  0 *** 0.001 ** 0.01 * 0.05 . 0.1 empty 1
33       
34    >>>
35
36
37.. autoclass:: LinearRegressionLearner
38    :members:
39
40.. autoclass:: LinearRegression
41    :members:
42
43Utility functions
44-----------------
45
46.. autofunction:: print_linear_regression_model
47
48.. autofunction:: stepwise
49
50
51"""
52
53
54import Orange
55from Orange.regression import base
56import numpy
57
58try:
59    from scipy import stats
60except ImportError:
61    import statc as stats
62
63from numpy import dot, sqrt
64from numpy.linalg import inv, pinv
65
66
67from Orange.misc import deprecated_members, deprecated_keywords
68
69class LinearRegressionLearner(base.BaseRegressionLearner):
70
71    """Fits the linear regression model, i.e. learns the regression parameters
72    The class is derived from
73    :class:`Orange.regression.base.BaseRegressionLearner`
74    which is used for preprocessing the data (continuization and imputation)
75    before fitting the regression parameters.
76
77    """   
78
79    def __init__(self, name='linear regression', intercept=True, \
80                 compute_stats=True, ridge_lambda=None,\
81                 imputer=None, continuizer=None, \
82                 use_vars=None, stepwise=False, add_sig=0.05,
83                 remove_sig=0.2, **kwds):
84        """
85        :param name: name of the linear model, default 'linear regression'
86        :type name: string
87        :param intercept: if True, the intercept beta0 is included
88            in the model
89        :type intercept: boolean
90        :param compute_stats: if True, statistical properties of
91            the estimators (standard error, t-scores, significances)
92            and statistical properties of the model
93            (sum of squares, R2, adjusted R2) are computed
94        :type compute_stats: boolean
95        :param ridge_lambda: if not None, the lambda parameter
96            in ridge regression
97        :type ridge_lambda: integer or None
98        :param use_vars: the list of independent varaiables included in
99            regression model. If None (default) all variables are used
100        :type use_vars: list of Orange.data.variable or None
101        :param stepwise: if True, _`stepwise regression`:
102            http://en.wikipedia.org/wiki/Stepwise_regression
103            based on F-test is performed. The significance parameters are
104            add_sig and remove_sig
105        :type stepwise: boolean
106        :param add_sig: lower bound of significance for which the variable
107            is included in regression model
108            default value = 0.05
109        :type add_sig: float
110        :param remove_sig: upper bound of significance for which
111            the variable is excluded from the regression model
112            default value = 0.2
113        :type remove_sig: float
114        """
115        self.name = name
116        self.intercept = intercept
117        self.compute_stats = compute_stats
118        self.ridge_lambda = ridge_lambda
119        self.set_imputer(imputer=imputer)
120        self.set_continuizer(continuizer=continuizer)
121        self.stepwise = stepwise
122        self.add_sig = add_sig
123        self.remove_sig = remove_sig
124        self.use_vars = use_vars
125        self.__dict__.update(kwds)
126       
127    def __call__(self, table, weight=None, verbose=0):
128        """
129        :param table: data instances.
130        :type table: :class:`Orange.data.Table`
131        :param weight: the weights for instances. Default: None, i.e.
132            all data instances are eqaully important in fitting
133            the regression parameters
134        :type weight: None or list of Orange.data.variable.Continuous
135            which stores weights for instances
136        """       
137        if not self.use_vars is None:
138            new_domain = Orange.data.Domain(self.use_vars,
139                                            table.domain.class_var)
140            new_domain.addmetas(table.domain.getmetas())
141            table = Orange.data.Table(new_domain, table)
142
143        # dicrete values are continuized       
144        table = self.continuize_table(table)
145         
146        # missing values are imputed
147        table = self.impute_table(table)
148
149        if self.stepwise:
150            use_vars = stepwise(table, weight, add_sig=self.add_sig,
151                                      remove_sig=self.remove_sig)
152            new_domain = Orange.data.Domain(use_vars, table.domain.class_var)
153            new_domain.addmetas(table.domain.getmetas())
154            table = Orange.data.Table(new_domain, table)
155
156        # convertion to numpy
157        A, y, w = table.to_numpy()
158        if A is None:
159            n, m = len(table), 0
160        else:
161            n, m = numpy.shape(A)
162     
163        if self.intercept:
164             if A is None:
165                 X = numpy.ones([n,1])
166             else:
167                 X = numpy.insert(A, 0, 1, axis=1) # adds a column of ones
168        else:
169             X = A
170             
171        domain = table.domain
172       
173        if numpy.std(y) < 10e-6: # almost constant variable
174            return Orange.regression.mean.MeanLearner(table)
175     
176        # set weights to the instances
177        W = numpy.identity(n)
178        if weight:
179            for i, ins in enumerate(table):
180                W[i, i] = float(ins[weight])
181
182        compute_stats = self.compute_stats
183        # adds some robustness by computing the pseudo inverse;
184        # normal inverse could fail due to singularity of the X.T * W * X
185        if self.ridge_lambda is None:
186            cov = pinv(dot(dot(X.T, W), X))
187        else:
188            cov = pinv(dot(dot(X.T, W), X) - self.ridge_lambda*numpy.eye(m+1))
189            compute_stats = False # TO DO: find inferential properties of the estimators
190        D = dot(dot(cov, X.T), W)
191        coefficients = dot(D, y)
192
193        mu_y, sigma_y = numpy.mean(y), numpy.std(y)
194        if A is not None:
195            cov_x = numpy.cov(X, rowvar=0)
196
197            # standardized coefficients
198            std_coefficients = (sqrt(cov_x.diagonal()) / sigma_y) \
199                                * coefficients
200        else:
201            std_coefficients = None
202
203        if compute_stats is False:
204            return LinearRegression(domain.class_var, domain, coefficients=coefficients,
205                                    std_coefficients=std_coefficients, intercept=self.intercept)
206           
207
208        fitted = dot(X, coefficients)
209        residuals = [ins.get_class() - fitted[i] \
210                     for i, ins in enumerate(table)]
211
212        # model summary       
213        # total sum of squares (total variance)
214        sst = numpy.sum((y - mu_y) ** 2)
215        # sum of squares due to regression (explained variance)
216        ssr = numpy.sum((fitted - mu_y)**2)
217        # error sum of squares (unexplaied variance)
218        sse = sst - ssr
219        # coefficient of determination
220        r2 = ssr / sst
221        r2adj = 1-(1-r2)*(n-1)/(n-m-1)
222        F = (ssr/m)/(sst-ssr/(n-m-1))
223        df = n-2 
224        sigma_square = sse/(n-m-1)
225        # standard error of the regression estimator, t-scores and p-values
226        std_error = sqrt(sigma_square*pinv(dot(X.T, X)).diagonal())
227        t_scores = coefficients/std_error
228        p_vals = [stats.betai(df*0.5,0.5,df/(df + t*t)) \
229                  for t in t_scores]
230
231        # dictionary of regression coefficients with standard errors
232        # and p-values
233        dict_model = {}
234        if self.intercept:
235            dict_model["Intercept"] = (coefficients[0],\
236                                      std_error[0], \
237                                      t_scores[0], \
238                                      p_vals[0])
239        for i, var in enumerate(domain.attributes):
240            j = i + 1 if self.intercept else i
241            dict_model[var.name] = (coefficients[j], \
242                                   std_error[j],\
243                                   t_scores[j],\
244                                   p_vals[j])
245       
246        return LinearRegression(domain.class_var, domain, coefficients, F,
247                 std_error=std_error, t_scores=t_scores, p_vals=p_vals, dict_model=dict_model,
248                 fitted=fitted, residuals=residuals, m=m, n=n, mu_y=mu_y,
249                 r2=r2, r2adj=r2adj, sst=sst, sse=sse, ssr=ssr,
250                 std_coefficients=std_coefficients, intercept=self.intercept)
251
252deprecated_members({"ridgeLambda": "ridge_lambda",
253                    "computeStats": "compute_stats",
254                    "useVars": "use_vars",
255                    "addSig": "add_sig",
256                    "removeSig": "remove_sig",
257                    }
258                   , ["__init__"],
259                   in_place=True)(LinearRegressionLearner)
260
261class LinearRegression(Orange.classification.Classifier):
262
263    """Linear regression predicts value of the response variable
264    based on the values of independent variables.
265
266    .. attribute:: F
267   
268        F-statistics of the model.
269
270    .. attribute:: coefficients
271
272        list of regression coefficients. If the intercept is included
273        the first item corresponds to the estimated intercept
274
275    .. attribute:: std_error
276
277        list of standard errors of the coefficient estimator.   
278
279    .. attribute:: t_scores
280
281        list of t-scores for the estimated regression coefficients   
282
283    .. attribute:: p_vals
284
285        list of p-values for the null hypothesis that the regression
286        coefficients equal 0 based on t-scores and two sided
287        alternative hypothesis   
288
289    .. attribute:: dict_model
290
291        dictionary of statistical properties of the model.
292        Keys - names of the independent variables (or "Intercept")
293        Values - tuples (coefficient, standard error,
294        t-value, p-value)
295
296    .. attribute:: fitted
297
298        estimated values of the dependent variable for all instances
299        from the training table
300
301    .. attribute:: residuals
302
303        differences between estimated and actual values of the
304        dependent variable for all instances from the training table
305
306    .. attribute:: m
307
308        number of independent variables   
309
310    .. attribute:: n
311
312        number of instances   
313
314    .. attribute:: mu_y
315
316        the sample mean of the dependent variable   
317
318    .. attribute:: r2
319
320        _`coefficient of determination`:
321        http://en.wikipedia.org/wiki/Coefficient_of_determination
322
323    .. attribute:: r2adj
324
325        adjusted coefficient of determination
326
327    .. attribute:: sst, sse, ssr
328
329        total sum of squares, explained sum of squares and
330        residual sum of squares respectively
331
332    .. attribute:: std_coefficients
333
334        standardized regression coefficients
335
336    """   
337
338
339   
340    def __init__(self, class_var=None, domain=None, coefficients=None, F=None,
341                 std_error=None, t_scores=None, p_vals=None, dict_model=None,
342                 fitted=None, residuals=None, m = None, n=None, mu_y=None,
343                 r2=None, r2adj=None, sst=None, sse=None, ssr=None,
344                 std_coefficients=None, intercept=None):
345        """
346        :param model: fitted linear regression model
347        :type model: :class:`LinearRegressionLearner`
348        """
349        self.class_var = class_var
350        self.domain = domain
351        self.coefficients = coefficients
352        self.F = F
353        self.std_error = std_error
354        self.t_scores = t_scores
355        self.p_vals = p_vals
356        self.dict_model = dict_model
357        self.fitted = fitted
358        self.residuals = residuals
359        self.m = m
360        self.n = n
361        self.mu_y = mu_y
362        self.r2 = r2
363        self.r2adj = r2adj
364        self.sst = sst
365        self.sse = sse
366        self.ssr = ssr
367        self.std_coefficients = std_coefficients
368        self.intercept = intercept
369
370    def __call__(self, instance, \
371                 resultType=Orange.classification.Classifier.GetValue):
372        """
373        :param instance: data instance for which the value of the response
374            variable will be predicted
375        :type instance:
376        """       
377        ins = Orange.data.Instance(self.domain, instance)
378        ins = numpy.array(ins.native())
379        if "?" in ins: # missing value -> corresponding coefficient omitted
380            def miss_2_0(x): return x if x != "?" else 0
381            ins = map(miss_2_0, ins)
382
383        if self.intercept:
384            if len(self.coefficients) > 1:
385                y_hat = self.coefficients[0] + \
386                       dot(self.coefficients[1:], ins[:-1])
387            else:
388                if len(ins) == 1:
389                    print ins
390                    y_hat = self.mu_y
391                else:
392                    y_hat = dot(self.coefficients, ins[:-1])
393        else:
394            y_hat = dot(self.coefficients, ins[:-1])
395#        y_hat = Orange.data.Value(y_hat)
396        y_hat = self.class_var(y_hat)
397        dist = Orange.statistics.distribution.Continuous(self.class_var)
398        dist[y_hat] = 1.0
399        if resultType == Orange.classification.Classifier.GetValue:
400            return y_hat
401        if resultType == Orange.classification.Classifier.GetProbabilities:
402            return dist
403        return (y_hat, dist)
404
405
406def print_linear_regression_model(lr):
407    """Pretty-prints linear regression model,
408    i.e. estimated regression coefficients with standard errors, t-scores
409    and significances.
410
411    :param lr: a linear regression model object.
412    :type lr: :class:`LinearRegression`   
413
414    """
415    from string import join
416    m = lr   
417    labels = ('Variable', 'Coeff Est', 'Std Error', 't-value', 'p')
418    print join(['%10s' % l for l in labels], ' ')
419
420    fmt = "%10s " + join(["%10.3f"]*4, " ") + " %5s"
421    if not lr.p_vals:
422        raise ValueError("Model does not contain model statistics.")
423    def get_star(p):
424        if p < 0.001: return  "*"*3
425        elif p < 0.01: return "*"*2
426        elif p < 0.05: return "*"
427        elif p < 0.1: return  "."
428        else: return " "
429   
430    if m.intercept == True:
431        stars =  get_star(m.p_vals[0])
432        print fmt % ('Intercept', m.coefficients[0], \
433                     m.std_error[0], m.t_scores[0], m.p_vals[0], stars)
434        for i in range(len(m.domain.attributes)):
435            stars = get_star(m.p_vals[i+1])
436            print fmt % (m.domain.attributes[i].name,\
437                         m.coefficients[i+1], m.std_error[i+1],\
438                         m.t_scores[i+1], m.p_vals[i+1], stars)
439    else:
440        for i in range(len(m.domain.attributes)):
441            stars = get_star(m.p_vals[i])
442            print fmt % (m.domain.attributes[i].name,\
443                         m.coefficients[i], m.std_error[i],\
444                         m.t_scores[i], m.p_vals[i], stars)
445    print "Signif. codes:  0 *** 0.001 ** 0.01 * 0.05 . 0.1 empty 1"
446
447
448
449def compare_models(c1, c2):
450    """ Compares if classifiaction model c1 is significantly better
451    than model c2. The comparison is based on F-test, the p-value
452    is returned.
453
454    :param c1, c2: linear regression model objects.
455    :type lr: :class:`LinearRegression`     
456
457    """
458    if c1 == None or c2 == None:
459        return 1.0
460    p1, p2, n = c1.m, c2.m, c1.n
461    RSS1, RSS2 = c1.sse, c2.sse
462    if RSS1 <= RSS2 or p2 <= p1 or n <= p2 or RSS2 <= 0:
463        return 1.0
464    F = ((RSS1-RSS2)/(p2-p1))/(RSS2/(n-p2))
465    return stats.fprob(int(p2-p1), int(n-p2), F)
466
467
468@deprecated_keywords({"addSig": "add_sig", "removeSig": "remove_sig"})
469def stepwise(table, weight, add_sig=0.05, remove_sig=0.2):
470    """ Performs _`stepwise linear regression`:
471    http://en.wikipedia.org/wiki/Stepwise_regression
472    on table and returns the list of remaing independent variables
473    which fit a significant linear regression model.coefficients
474
475    :param table: data instances.
476    :type table: :class:`Orange.data.Table`
477    :param weight: the weights for instances. Default: None, i.e. all data
478        instances are eqaully important in fitting the regression parameters
479    :type weight: None or list of Orange.data.variable.Continuous
480        which stores the weights
481    :param add_sig: lower bound of significance for which the variable
482        is included in regression model
483        default value = 0.05
484    :type add_sig: float
485    :param remove_sig: upper bound of significance for which the variable
486        is excluded from the regression model
487        default value = 0.2
488    :type remove_sig: float
489    """
490
491   
492    inc_vars = []
493    not_inc_vars = table.domain.attributes
494
495    changed_model = True
496    while changed_model:
497        changed_model = False
498        # remove all unsignificant conditions (learn several models,
499        # where each time one variable is removed and check significance)
500        c0 = LinearRegressionLearner(table, use_vars=inc_vars)
501        reduced_model = [] # reduced model
502        for ati in range(len(inc_vars)):
503            try:
504                reduced_model.append(LinearRegressionLearner(table, weight,
505                        use_vars=inc_vars[:ati] + inc_vars[(ati + 1):]))
506            except Exception:
507                reduced_model.append(None)
508       
509        sigs = [compare_models(r, c0) for r in reduced_model]
510        if sigs and max(sigs) > remove_sig:
511            # remove that variable, start again
512            crit_var = inc_vars[sigs.index(max(sigs))]
513            not_inc_vars.append(crit_var)
514            inc_vars.remove(crit_var)
515            changed_model = True
516            continue
517
518        # add all significant conditions (go through all attributes in
519        # not_inc_vars, is there one that significantly improves the model?
520        extended_model = []
521        for ati in range(len(not_inc_vars)):
522            try:
523                extended_model.append(LinearRegressionLearner(table,
524                        weight, use_vars=inc_vars + [not_inc_vars[ati]]))
525            except Exception:
526                extended_model.append(None)
527             
528        sigs = [compare_models(c0, r) for r in extended_model]
529        if sigs and min(sigs) < add_sig:
530            best_var = not_inc_vars[sigs.index(min(sigs))]
531            inc_vars.append(best_var)
532            not_inc_vars.remove(best_var)
533            changed_model = True
534    return inc_vars
535
536
537if __name__ == "__main__":
538
539    import Orange
540    from Orange.regression import linear
541
542    table = Orange.data.Table("housing.tab")
543    c = LinearRegressionLearner(table)
544    print_linear_regression_model(c)
Note: See TracBrowser for help on using the repository browser.