source: orange/Orange/regression/linear.py @ 10605:2ad79c723d30

Revision 10605:2ad79c723d30, 21.3 KB checked in by Lan Zagar <lan.zagar@…>, 2 years ago (diff)

Dont regularize the intercept and avoid div by zero.

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