source: orange/Orange/regression/linear.py @ 10616:0ef88f0d8128

Revision 10616:0ef88f0d8128, 21.2 KB checked in by Lan Zagar <lan.zagar@…>, 2 years ago (diff)

Do not compute whole cov matrix to get standard deviations.

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