source: orange/Orange/regression/linear.py @ 10643:79cc41d29eb7

Revision 10643:79cc41d29eb7, 21.1 KB checked in by Lan Zagar <lan.zagar@…>, 2 years ago (diff)

Corrections to computations of statistics in linear regression.

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, std
137from numpy.linalg import inv, pinv
138
139
140from Orange.utils 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 equally 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 self.use_vars is not 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        # discrete 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        domain = table.domain
230
231        # convert to numpy
232        X, y, w = table.to_numpy()
233        n, m = numpy.shape(X)
234
235        if self.intercept:
236            X = numpy.insert(X, 0, 1, axis=1) # adds a column of ones
237
238        if weight:
239            weights = numpy.sqrt([float(ins[weight]) for ins in table])
240            X = weights.reshape(n, 1) * X
241            y = weights * y
242
243        cov = dot(X.T, X)
244
245        if self.ridge_lambda:
246            stride = cov.shape[0] + 1
247            cov.flat[self.intercept * stride::stride] += self.ridge_lambda
248
249        # adds some robustness by computing the pseudo inverse;
250        # normal inverse could fail due to the singularity of X.T * X
251        invcov = pinv(cov)
252        D = dot(invcov, X.T)
253        coefficients = dot(D, y)
254
255        mu_y, sigma_y = numpy.mean(y), numpy.std(y)
256        if m > 0:
257            # standardized coefficients
258            std_coefficients = std(X, axis=0, ddof=1) / sigma_y * coefficients
259        else:
260            std_coefficients = None
261
262        # TODO: find inferential properties of the estimators for ridge
263        if self.compute_stats is False or self.ridge_lambda:
264            return LinearRegression(domain.class_var, domain,
265                coefficients=coefficients, std_coefficients=std_coefficients,
266                intercept=self.intercept)
267
268        fitted = dot(X, coefficients)
269        residuals = [ins.get_class() - fitted[i]
270                     for i, ins in enumerate(table)]
271
272        # model summary
273        df_reg = n - m - self.intercept
274        # total sum of squares (total variance)
275        sst = numpy.sum((y - mu_y) ** 2)
276        # regression sum of squares (explained variance)
277        ssr = numpy.sum((fitted - mu_y) ** 2)
278        # residual sum of squares
279        sse = numpy.sum((y - fitted) ** 2)
280        # coefficient of determination
281        r2 = ssr / sst
282        r2 = 1 - sse / sst
283        r2adj = 1 - (1 - r2) * (n - 1) / df_reg
284        F = (ssr / m) / ((sst - ssr) / df_reg) if m else 0
285        sigma_square = sse / df_reg
286        # standard error of the regression estimator, t-scores and p-values
287        std_error = sqrt(sigma_square * invcov.diagonal())
288        t_scores = coefficients / std_error
289        df_res = n - 2
290        p_vals = [stats.betai(df_res * 0.5, 0.5, df_res / (df_res + t * t))
291                  for t in t_scores]
292
293        # dictionary of regression coefficients with standard errors
294        # and p-values
295        dict_model = {}
296        if self.intercept:
297            dict_model["Intercept"] = (coefficients[0], std_error[0],
298                                       t_scores[0], p_vals[0])
299        for i, var in enumerate(domain.features):
300            j = i + 1 if self.intercept else i
301            dict_model[var.name] = (coefficients[j], std_error[j],
302                                    t_scores[j], p_vals[j])
303
304        return LinearRegression(domain.class_var, domain, coefficients, F,
305                 std_error=std_error, t_scores=t_scores, p_vals=p_vals,
306                 dict_model=dict_model, fitted=fitted, residuals=residuals,
307                 m=m, n=n, mu_y=mu_y, r2=r2, r2adj=r2adj, sst=sst, sse=sse,
308                 ssr=ssr, std_coefficients=std_coefficients,
309                 intercept=self.intercept)
310
311deprecated_members({"ridgeLambda": "ridge_lambda",
312                    "computeStats": "compute_stats",
313                    "useVars": "use_vars",
314                    "addSig": "add_sig",
315                    "removeSig": "remove_sig",
316                    }
317                   , ["__init__"],
318                   in_place=True)(LinearRegressionLearner)
319
320class LinearRegression(Orange.classification.Classifier):
321
322    """Linear regression predicts value of the response variable
323    based on the values of independent variables.
324
325    .. attribute:: F
326
327        F-statistics of the model.
328
329    .. attribute:: coefficients
330
331        Regression coefficients stored in list. If the intercept is included
332        the first item corresponds to the estimated intercept.
333
334    .. attribute:: std_error
335
336        Standard errors of the coefficient estimator, stored in list.
337
338    .. attribute:: t_scores
339
340        List of t-scores for the estimated regression coefficients.
341
342    .. attribute:: p_vals
343
344        List of p-values for the null hypothesis that the regression
345        coefficients equal 0 based on t-scores and two sided
346        alternative hypothesis.
347
348    .. attribute:: dict_model
349
350        Statistical properties of the model in a dictionary:
351        Keys - names of the independent variables (or "Intercept")
352        Values - tuples (coefficient, standard error,
353        t-value, p-value)
354
355    .. attribute:: fitted
356
357        Estimated values of the dependent variable for all instances
358        from the training table.
359
360    .. attribute:: residuals
361
362        Differences between estimated and actual values of the
363        dependent variable for all instances from the training table.
364
365    .. attribute:: m
366
367        Number of independent (predictor) variables.
368
369    .. attribute:: n
370
371        Number of instances.
372
373    .. attribute:: mu_y
374
375        Sample mean of the dependent variable.
376
377    .. attribute:: r2
378
379        `Coefficient of determination
380        <http://en.wikipedia.org/wiki/Coefficient_of_determination>`_.
381
382
383    .. attribute:: r2adj
384
385        Adjusted coefficient of determination.
386
387    .. attribute:: sst, sse, ssr
388
389        Total sum of squares, explained sum of squares and
390        residual sum of squares respectively.
391
392    .. attribute:: std_coefficients
393
394        Standardized regression coefficients.
395
396    """
397
398    def __init__(self, class_var=None, domain=None, coefficients=None, F=None,
399                 std_error=None, t_scores=None, p_vals=None, dict_model=None,
400                 fitted=None, residuals=None, m=None, n=None, mu_y=None,
401                 r2=None, r2adj=None, sst=None, sse=None, ssr=None,
402                 std_coefficients=None, intercept=None):
403        """
404        :param model: fitted linear regression model
405        :type model: :class:`LinearRegressionLearner`
406        """
407        self.class_var = class_var
408        self.domain = domain
409        self.coefficients = coefficients
410        self.F = F
411        self.std_error = std_error
412        self.t_scores = t_scores
413        self.p_vals = p_vals
414        self.dict_model = dict_model
415        self.fitted = fitted
416        self.residuals = residuals
417        self.m = m
418        self.n = n
419        self.mu_y = mu_y
420        self.r2 = r2
421        self.r2adj = r2adj
422        self.sst = sst
423        self.sse = sse
424        self.ssr = ssr
425        self.std_coefficients = std_coefficients
426        self.intercept = intercept
427
428    def __call__(self, instance,
429                 result_type=Orange.classification.Classifier.GetValue):
430        """
431        :param instance: data instance for which the value of the response
432                         variable will be predicted
433        :type instance: :obj:`~Orange.data.Instance`
434        """
435        ins = Orange.data.Instance(self.domain, instance)
436        ins = numpy.array(ins.native())
437        if "?" in ins: # missing value -> corresponding coefficient omitted
438            def miss_2_0(x): return x if x != "?" else 0
439            ins = map(miss_2_0, ins)
440
441        if self.intercept:
442            if len(self.coefficients) > 1:
443                y_hat = self.coefficients[0] + dot(self.coefficients[1:],
444                                                   ins[:-1])
445            else:
446                if len(ins) == 1:
447                    print ins
448                    y_hat = self.mu_y
449                else:
450                    y_hat = dot(self.coefficients, ins[:-1])
451        else:
452            y_hat = dot(self.coefficients, ins[:-1])
453#        y_hat = Orange.data.Value(y_hat)
454        y_hat = self.class_var(y_hat)
455        dist = Orange.statistics.distribution.Continuous(self.class_var)
456        dist[y_hat] = 1.
457        if result_type == Orange.classification.Classifier.GetValue:
458            return y_hat
459        if result_type == Orange.classification.Classifier.GetProbabilities:
460            return dist
461        return (y_hat, dist)
462
463    def to_string(self):
464        """Pretty-prints linear regression model,
465        i.e. estimated regression coefficients with standard errors, t-scores
466        and significances.
467
468        """
469        from string import join
470        labels = ('Variable', 'Coeff Est', 'Std Error', 't-value', 'p')
471        lines = [join(['%10s' % l for l in labels], ' ')]
472
473        fmt = "%10s " + join(["%10.3f"] * 4, " ") + " %5s"
474        if not self.p_vals:
475            raise ValueError("Model does not contain model statistics.")
476        def get_star(p):
477            if p < 0.001: return  "*" * 3
478            elif p < 0.01: return "*" * 2
479            elif p < 0.05: return "*"
480            elif p < 0.1: return  "."
481            else: return " "
482
483        if self.intercept == True:
484            stars =  get_star(self.p_vals[0])
485            lines.append(fmt % ('Intercept', self.coefficients[0],
486                                self.std_error[0], self.t_scores[0],
487                                self.p_vals[0], stars))
488            for i in range(len(self.domain.attributes)):
489                stars = get_star(self.p_vals[i + 1])
490                lines.append(fmt % (self.domain.attributes[i].name,
491                             self.coefficients[i + 1], self.std_error[i + 1],
492                             self.t_scores[i + 1], self.p_vals[i + 1], stars))
493        else:
494            for i in range(len(self.domain.attributes)):
495                stars = get_star(self.p_vals[i])
496                lines.append(fmt % (self.domain.attributes[i].name,
497                             self.coefficients[i], self.std_error[i],
498                             self.t_scores[i], self.p_vals[i], stars))
499        lines.append("Signif. codes:  0 *** 0.001 ** 0.01 * 0.05 . 0.1 empty 1")
500        return "\n".join(lines)
501
502    def __str__(self):
503        return self.to_string()
504
505
506def compare_models(c1, c2):
507    """ Compares if classifiaction model c1 is significantly better
508    than model c2. The comparison is based on F-test, the p-value
509    is returned.
510
511    :param c1, c2: linear regression model objects.
512    :type lr: :class:`LinearRegression`
513
514    """
515    if c1 == None or c2 == None:
516        return 1.0
517    p1, p2, n = c1.m, c2.m, c1.n
518    RSS1, RSS2 = c1.sse, c2.sse
519    if RSS1 <= RSS2 or p2 <= p1 or n <= p2 or RSS2 <= 0:
520        return 1.0
521    F = ((RSS1 - RSS2) / (p2 - p1)) / (RSS2 / (n - p2))
522    return stats.fprob(int(p2 - p1), int(n - p2), F)
523
524
525@deprecated_keywords({"addSig": "add_sig", "removeSig": "remove_sig"})
526def stepwise(table, weight, add_sig=0.05, remove_sig=0.2):
527    """ Performs `stepwise linear regression
528    <http://en.wikipedia.org/wiki/Stepwise_regression>`_:
529    on table and returns the list of remaing independent variables
530    which fit a significant linear regression model.coefficients
531
532    :param table: data instances.
533    :type table: :class:`Orange.data.Table`
534    :param weight: the weights for instances. Default: None, i.e. all data
535        instances are eqaully important in fitting the regression parameters
536    :type weight: None or list of Orange.feature.Continuous
537        which stores the weights
538    :param add_sig: lower bound of significance for which the variable
539        is included in regression model
540        default value = 0.05
541    :type add_sig: float
542    :param remove_sig: upper bound of significance for which the variable
543        is excluded from the regression model
544        default value = 0.2
545    :type remove_sig: float
546    """
547
548    inc_vars = []
549    not_inc_vars = table.domain.attributes
550
551    changed_model = True
552    while changed_model:
553        changed_model = False
554        # remove all unsignificant conditions (learn several models,
555        # where each time one variable is removed and check significance)
556        c0 = LinearRegressionLearner(table, use_vars=inc_vars)
557        reduced_model = [] # reduced model
558        for ati in range(len(inc_vars)):
559            try:
560                reduced_model.append(LinearRegressionLearner(table, weight,
561                        use_vars=inc_vars[:ati] + inc_vars[(ati + 1):]))
562            except Exception:
563                reduced_model.append(None)
564
565        sigs = [compare_models(r, c0) for r in reduced_model]
566        if sigs and max(sigs) > remove_sig:
567            # remove that variable, start again
568            crit_var = inc_vars[sigs.index(max(sigs))]
569            not_inc_vars.append(crit_var)
570            inc_vars.remove(crit_var)
571            changed_model = True
572            continue
573
574        # add all significant conditions (go through all attributes in
575        # not_inc_vars, is there one that significantly improves the model?
576        extended_model = []
577        for ati in range(len(not_inc_vars)):
578            try:
579                extended_model.append(LinearRegressionLearner(table,
580                        weight, use_vars=inc_vars + [not_inc_vars[ati]]))
581            except Exception:
582                extended_model.append(None)
583
584        sigs = [compare_models(c0, r) for r in extended_model]
585        if sigs and min(sigs) < add_sig:
586            best_var = not_inc_vars[sigs.index(min(sigs))]
587            inc_vars.append(best_var)
588            not_inc_vars.remove(best_var)
589            changed_model = True
590    return inc_vars
591
592
593if __name__ == "__main__":
594
595    import Orange
596    from Orange.regression import linear
597
598    table = Orange.data.Table("housing.tab")
599    c = LinearRegressionLearner(table)
600    print c
601
Note: See TracBrowser for help on using the repository browser.