source: orange/Orange/regression/linear.py @ 9776:600ac31393ee

Revision 9776:600ac31393ee, 21.5 KB checked in by lanumek, 2 years ago (diff)

documentation for regression: lasso, linear, pls

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