source: orange/Orange/regression/linear.py @ 10575:2718a801a48e

Revision 10575:2718a801a48e, 21.4 KB checked in by Lan Zagar <lan.zagar@…>, 2 years ago (diff)

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