source: orange/Orange/regression/linear.py @ 10777:546aa4c30fdc

Revision 10777:546aa4c30fdc, 21.3 KB checked in by Ales Erjavec <ales.erjavec@…>, 2 years ago (diff)

BUGFIX: Fixed 'to_string' method to work even if model statistics are not available.

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