source: orange/Orange/regression/linear.py @ 10557:5af235500ac2

Revision 10557:5af235500ac2, 21.6 KB checked in by Lan Zagar <lan.zagar@…>, 2 years ago (diff)

Better solution for intercept (see #1115).

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