source: orange/orange/regression/linear.py @ 9669:165371b04b4a

Revision 9669:165371b04b4a, 19.3 KB checked in by anze <anze.staric@…>, 2 years ago (diff)

Moved content of Orange dir to package dir

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
10Example ::
11
12    >>> from Orange.regression import linear
13    >>> table = Orange.data.Table("housing")
14    >>> c = linear.LinearRegressionLearner(table)
15    >>> print c
16   
17      Variable  Coeff Est  Std Error    t-value          p
18     Intercept     36.459      5.103      7.144      0.000   ***
19          CRIM     -0.108      0.033     -3.287      0.001    **
20            ZN      0.046      0.014      3.382      0.001   ***
21         INDUS      0.021      0.061      0.334      0.738     
22          CHAS      2.687      0.862      3.118      0.002    **
23           NOX    -17.767      3.820     -4.651      0.000   ***
24            RM      3.810      0.418      9.116      0.000   ***
25           AGE      0.001      0.013      0.052      0.958     
26           DIS     -1.476      0.199     -7.398      0.000   ***
27           RAD      0.306      0.066      4.613      0.000   ***
28           TAX     -0.012      0.004     -3.280      0.001    **
29       PTRATIO     -0.953      0.131     -7.283      0.000   ***
30             B      0.009      0.003      3.467      0.001   ***
31         LSTAT     -0.525      0.051    -10.347      0.000   ***
32    Signif. codes:  0 *** 0.001 ** 0.01 * 0.05 . 0.1 empty 1
33       
34    >>>
35
36
37.. autoclass:: LinearRegressionLearner
38    :members:
39
40.. autoclass:: LinearRegression
41    :members:
42
43Utility functions
44-----------------
45
46.. autofunction:: stepwise
47
48
49"""
50
51
52import Orange
53from Orange.regression import base
54import numpy
55
56try:
57    from scipy import stats
58except ImportError:
59    import statc as stats
60
61from numpy import dot, sqrt
62from numpy.linalg import inv, pinv
63
64
65from Orange.misc import deprecated_members, deprecated_keywords
66
67class LinearRegressionLearner(base.BaseRegressionLearner):
68
69    """Fits the linear regression model, i.e. learns the regression parameters
70    The class is derived from
71    :class:`Orange.regression.base.BaseRegressionLearner`
72    which is used for preprocessing the data (continuization and imputation)
73    before fitting the regression parameters.
74
75    """   
76
77    def __init__(self, name='linear regression', intercept=True, \
78                 compute_stats=True, ridge_lambda=None,\
79                 imputer=None, continuizer=None, \
80                 use_vars=None, stepwise=False, add_sig=0.05,
81                 remove_sig=0.2, **kwds):
82        """
83        :param name: name of the linear model, default 'linear regression'
84        :type name: string
85        :param intercept: if True, the intercept beta0 is included
86            in the model
87        :type intercept: boolean
88        :param compute_stats: if True, statistical properties of
89            the estimators (standard error, t-scores, significances)
90            and statistical properties of the model
91            (sum of squares, R2, adjusted R2) are computed
92        :type compute_stats: boolean
93        :param ridge_lambda: if not None, the lambda parameter
94            in ridge regression
95        :type ridge_lambda: integer or None
96        :param use_vars: the list of independent varaiables included in
97            regression model. If None (default) all variables are used
98        :type use_vars: list of Orange.data.variable or None
99        :param stepwise: if True, _`stepwise regression`:
100            http://en.wikipedia.org/wiki/Stepwise_regression
101            based on F-test is performed. The significance parameters are
102            add_sig and remove_sig
103        :type stepwise: boolean
104        :param add_sig: lower bound of significance for which the variable
105            is included in regression model
106            default value = 0.05
107        :type add_sig: float
108        :param remove_sig: upper bound of significance for which
109            the variable is excluded from the regression model
110            default value = 0.2
111        :type remove_sig: float
112        """
113        self.name = name
114        self.intercept = intercept
115        self.compute_stats = compute_stats
116        self.ridge_lambda = ridge_lambda
117        self.set_imputer(imputer=imputer)
118        self.set_continuizer(continuizer=continuizer)
119        self.stepwise = stepwise
120        self.add_sig = add_sig
121        self.remove_sig = remove_sig
122        self.use_vars = use_vars
123        self.__dict__.update(kwds)
124       
125    def __call__(self, table, weight=None, verbose=0):
126        """
127        :param table: data instances.
128        :type table: :class:`Orange.data.Table`
129        :param weight: the weights for instances. Default: None, i.e.
130            all data instances are eqaully important in fitting
131            the regression parameters
132        :type weight: None or list of Orange.data.variable.Continuous
133            which stores weights for instances
134        """       
135        if not self.use_vars is None:
136            new_domain = Orange.data.Domain(self.use_vars,
137                                            table.domain.class_var)
138            new_domain.addmetas(table.domain.getmetas())
139            table = Orange.data.Table(new_domain, table)
140
141        # dicrete values are continuized       
142        table = self.continuize_table(table)
143         
144        # missing values are imputed
145        table = self.impute_table(table)
146
147        if self.stepwise:
148            use_vars = stepwise(table, weight, add_sig=self.add_sig,
149                                      remove_sig=self.remove_sig)
150            new_domain = Orange.data.Domain(use_vars, table.domain.class_var)
151            new_domain.addmetas(table.domain.getmetas())
152            table = Orange.data.Table(new_domain, table)
153
154        # convertion to numpy
155        A, y, w = table.to_numpy()
156        if A is None:
157            n, m = len(table), 0
158        else:
159            n, m = numpy.shape(A)
160     
161        if self.intercept:
162             if A is None:
163                 X = numpy.ones([n,1])
164             else:
165                 X = numpy.insert(A, 0, 1, axis=1) # adds a column of ones
166        else:
167             X = A
168             
169        domain = table.domain
170       
171        if numpy.std(y) < 10e-6: # almost constant variable
172            return Orange.regression.mean.MeanLearner(table)
173     
174        # set weights to the instances
175        W = numpy.identity(n)
176        if weight:
177            for i, ins in enumerate(table):
178                W[i, i] = float(ins[weight])
179
180        compute_stats = self.compute_stats
181        # adds some robustness by computing the pseudo inverse;
182        # normal inverse could fail due to singularity of the X.T * W * X
183        if self.ridge_lambda is None:
184            cov = pinv(dot(dot(X.T, W), X))
185        else:
186            cov = pinv(dot(dot(X.T, W), X) - self.ridge_lambda*numpy.eye(m+1))
187            compute_stats = False # TO DO: find inferential properties of the estimators
188        D = dot(dot(cov, X.T), W)
189        coefficients = dot(D, y)
190
191        mu_y, sigma_y = numpy.mean(y), numpy.std(y)
192        if A is not None:
193            cov_x = numpy.cov(X, rowvar=0)
194
195            # standardized coefficients
196            std_coefficients = (sqrt(cov_x.diagonal()) / sigma_y) \
197                                * coefficients
198        else:
199            std_coefficients = None
200
201        if compute_stats is False:
202            return LinearRegression(domain.class_var, domain, coefficients=coefficients,
203                                    std_coefficients=std_coefficients, intercept=self.intercept)
204           
205
206        fitted = dot(X, coefficients)
207        residuals = [ins.get_class() - fitted[i] \
208                     for i, ins in enumerate(table)]
209
210        # model summary       
211        # total sum of squares (total variance)
212        sst = numpy.sum((y - mu_y) ** 2)
213        # sum of squares due to regression (explained variance)
214        ssr = numpy.sum((fitted - mu_y)**2)
215        # error sum of squares (unexplaied variance)
216        sse = sst - ssr
217        # coefficient of determination
218        r2 = ssr / sst
219        r2adj = 1-(1-r2)*(n-1)/(n-m-1)
220        F = (ssr/m)/(sst-ssr/(n-m-1))
221        df = n-2 
222        sigma_square = sse/(n-m-1)
223        # standard error of the regression estimator, t-scores and p-values
224        std_error = sqrt(sigma_square*pinv(dot(X.T, X)).diagonal())
225        t_scores = coefficients/std_error
226        p_vals = [stats.betai(df*0.5,0.5,df/(df + t*t)) \
227                  for t in t_scores]
228
229        # dictionary of regression coefficients with standard errors
230        # and p-values
231        dict_model = {}
232        if self.intercept:
233            dict_model["Intercept"] = (coefficients[0],\
234                                      std_error[0], \
235                                      t_scores[0], \
236                                      p_vals[0])
237        for i, var in enumerate(domain.attributes):
238            j = i + 1 if self.intercept else i
239            dict_model[var.name] = (coefficients[j], \
240                                   std_error[j],\
241                                   t_scores[j],\
242                                   p_vals[j])
243       
244        return LinearRegression(domain.class_var, domain, coefficients, F,
245                 std_error=std_error, t_scores=t_scores, p_vals=p_vals, dict_model=dict_model,
246                 fitted=fitted, residuals=residuals, m=m, n=n, mu_y=mu_y,
247                 r2=r2, r2adj=r2adj, sst=sst, sse=sse, ssr=ssr,
248                 std_coefficients=std_coefficients, intercept=self.intercept)
249
250deprecated_members({"ridgeLambda": "ridge_lambda",
251                    "computeStats": "compute_stats",
252                    "useVars": "use_vars",
253                    "addSig": "add_sig",
254                    "removeSig": "remove_sig",
255                    }
256                   , ["__init__"],
257                   in_place=True)(LinearRegressionLearner)
258
259class LinearRegression(Orange.classification.Classifier):
260
261    """Linear regression predicts value of the response variable
262    based on the values of independent variables.
263
264    .. attribute:: F
265   
266        F-statistics of the model.
267
268    .. attribute:: coefficients
269
270        list of regression coefficients. If the intercept is included
271        the first item corresponds to the estimated intercept
272
273    .. attribute:: std_error
274
275        list of standard errors of the coefficient estimator.   
276
277    .. attribute:: t_scores
278
279        list of t-scores for the estimated regression coefficients   
280
281    .. attribute:: p_vals
282
283        list of p-values for the null hypothesis that the regression
284        coefficients equal 0 based on t-scores and two sided
285        alternative hypothesis   
286
287    .. attribute:: dict_model
288
289        dictionary of statistical properties of the model.
290        Keys - names of the independent variables (or "Intercept")
291        Values - tuples (coefficient, standard error,
292        t-value, p-value)
293
294    .. attribute:: fitted
295
296        estimated values of the dependent variable for all instances
297        from the training table
298
299    .. attribute:: residuals
300
301        differences between estimated and actual values of the
302        dependent variable for all instances from the training table
303
304    .. attribute:: m
305
306        number of independent variables   
307
308    .. attribute:: n
309
310        number of instances   
311
312    .. attribute:: mu_y
313
314        the sample mean of the dependent variable   
315
316    .. attribute:: r2
317
318        _`coefficient of determination`:
319        http://en.wikipedia.org/wiki/Coefficient_of_determination
320
321    .. attribute:: r2adj
322
323        adjusted coefficient of determination
324
325    .. attribute:: sst, sse, ssr
326
327        total sum of squares, explained sum of squares and
328        residual sum of squares respectively
329
330    .. attribute:: std_coefficients
331
332        standardized regression coefficients
333
334    """   
335
336
337   
338    def __init__(self, class_var=None, domain=None, coefficients=None, F=None,
339                 std_error=None, t_scores=None, p_vals=None, dict_model=None,
340                 fitted=None, residuals=None, m = None, n=None, mu_y=None,
341                 r2=None, r2adj=None, sst=None, sse=None, ssr=None,
342                 std_coefficients=None, intercept=None):
343        """
344        :param model: fitted linear regression model
345        :type model: :class:`LinearRegressionLearner`
346        """
347        self.class_var = class_var
348        self.domain = domain
349        self.coefficients = coefficients
350        self.F = F
351        self.std_error = std_error
352        self.t_scores = t_scores
353        self.p_vals = p_vals
354        self.dict_model = dict_model
355        self.fitted = fitted
356        self.residuals = residuals
357        self.m = m
358        self.n = n
359        self.mu_y = mu_y
360        self.r2 = r2
361        self.r2adj = r2adj
362        self.sst = sst
363        self.sse = sse
364        self.ssr = ssr
365        self.std_coefficients = std_coefficients
366        self.intercept = intercept
367
368    def __call__(self, instance, \
369                 resultType=Orange.classification.Classifier.GetValue):
370        """
371        :param instance: data instance for which the value of the response
372            variable will be predicted
373        :type instance:
374        """       
375        ins = Orange.data.Instance(self.domain, instance)
376        ins = numpy.array(ins.native())
377        if "?" in ins: # missing value -> corresponding coefficient omitted
378            def miss_2_0(x): return x if x != "?" else 0
379            ins = map(miss_2_0, ins)
380
381        if self.intercept:
382            if len(self.coefficients) > 1:
383                y_hat = self.coefficients[0] + \
384                       dot(self.coefficients[1:], ins[:-1])
385            else:
386                if len(ins) == 1:
387                    print ins
388                    y_hat = self.mu_y
389                else:
390                    y_hat = dot(self.coefficients, ins[:-1])
391        else:
392            y_hat = dot(self.coefficients, ins[:-1])
393#        y_hat = Orange.data.Value(y_hat)
394        y_hat = self.class_var(y_hat)
395        dist = Orange.statistics.distribution.Continuous(self.class_var)
396        dist[y_hat] = 1.0
397        if resultType == Orange.classification.Classifier.GetValue:
398            return y_hat
399        if resultType == Orange.classification.Classifier.GetProbabilities:
400            return dist
401        return (y_hat, dist)
402
403
404    def to_string(self):
405        """Pretty-prints linear regression model,
406        i.e. estimated regression coefficients with standard errors, t-scores
407        and significances.
408
409        """
410        from string import join
411        labels = ('Variable', 'Coeff Est', 'Std Error', 't-value', 'p')
412        lines = [join(['%10s' % l for l in labels], ' ')]
413
414        fmt = "%10s " + join(["%10.3f"]*4, " ") + " %5s"
415        if not self.p_vals:
416            raise ValueError("Model does not contain model statistics.")
417        def get_star(p):
418            if p < 0.001: return  "*"*3
419            elif p < 0.01: return "*"*2
420            elif p < 0.05: return "*"
421            elif p < 0.1: return  "."
422            else: return " "
423       
424        if self.intercept == True:
425            stars =  get_star(self.p_vals[0])
426            lines.append(fmt % ('Intercept', self.coefficients[0], 
427                         self.std_error[0], self.t_scores[0], self.p_vals[0], stars))
428            for i in range(len(self.domain.attributes)):
429                stars = get_star(self.p_vals[i+1])
430                lines.append(fmt % (self.domain.attributes[i].name,
431                             self.coefficients[i+1], self.std_error[i+1],
432                             self.t_scores[i+1], self.p_vals[i+1], stars))
433        else:
434            for i in range(len(self.domain.attributes)):
435                stars = get_star(self.p_vals[i])
436                lines.append(fmt % (self.domain.attributes[i].name,
437                             self.coefficients[i], self.std_error[i],
438                             self.t_scores[i], self.p_vals[i], stars))
439        lines.append("Signif. codes:  0 *** 0.001 ** 0.01 * 0.05 . 0.1 empty 1")
440        return "\n".join(lines)
441
442    def __str__(self):
443        return self.to_string()
444       
445
446
447
448def compare_models(c1, c2):
449    """ Compares if classifiaction model c1 is significantly better
450    than model c2. The comparison is based on F-test, the p-value
451    is returned.
452
453    :param c1, c2: linear regression model objects.
454    :type lr: :class:`LinearRegression`     
455
456    """
457    if c1 == None or c2 == None:
458        return 1.0
459    p1, p2, n = c1.m, c2.m, c1.n
460    RSS1, RSS2 = c1.sse, c2.sse
461    if RSS1 <= RSS2 or p2 <= p1 or n <= p2 or RSS2 <= 0:
462        return 1.0
463    F = ((RSS1-RSS2)/(p2-p1))/(RSS2/(n-p2))
464    return stats.fprob(int(p2-p1), int(n-p2), F)
465
466
467@deprecated_keywords({"addSig": "add_sig", "removeSig": "remove_sig"})
468def stepwise(table, weight, add_sig=0.05, remove_sig=0.2):
469    """ Performs _`stepwise linear regression`:
470    http://en.wikipedia.org/wiki/Stepwise_regression
471    on table and returns the list of remaing independent variables
472    which fit a significant linear regression model.coefficients
473
474    :param table: data instances.
475    :type table: :class:`Orange.data.Table`
476    :param weight: the weights for instances. Default: None, i.e. all data
477        instances are eqaully important in fitting the regression parameters
478    :type weight: None or list of Orange.data.variable.Continuous
479        which stores the weights
480    :param add_sig: lower bound of significance for which the variable
481        is included in regression model
482        default value = 0.05
483    :type add_sig: float
484    :param remove_sig: upper bound of significance for which the variable
485        is excluded from the regression model
486        default value = 0.2
487    :type remove_sig: float
488    """
489
490   
491    inc_vars = []
492    not_inc_vars = table.domain.attributes
493
494    changed_model = True
495    while changed_model:
496        changed_model = False
497        # remove all unsignificant conditions (learn several models,
498        # where each time one variable is removed and check significance)
499        c0 = LinearRegressionLearner(table, use_vars=inc_vars)
500        reduced_model = [] # reduced model
501        for ati in range(len(inc_vars)):
502            try:
503                reduced_model.append(LinearRegressionLearner(table, weight,
504                        use_vars=inc_vars[:ati] + inc_vars[(ati + 1):]))
505            except Exception:
506                reduced_model.append(None)
507       
508        sigs = [compare_models(r, c0) for r in reduced_model]
509        if sigs and max(sigs) > remove_sig:
510            # remove that variable, start again
511            crit_var = inc_vars[sigs.index(max(sigs))]
512            not_inc_vars.append(crit_var)
513            inc_vars.remove(crit_var)
514            changed_model = True
515            continue
516
517        # add all significant conditions (go through all attributes in
518        # not_inc_vars, is there one that significantly improves the model?
519        extended_model = []
520        for ati in range(len(not_inc_vars)):
521            try:
522                extended_model.append(LinearRegressionLearner(table,
523                        weight, use_vars=inc_vars + [not_inc_vars[ati]]))
524            except Exception:
525                extended_model.append(None)
526             
527        sigs = [compare_models(c0, r) for r in extended_model]
528        if sigs and min(sigs) < add_sig:
529            best_var = not_inc_vars[sigs.index(min(sigs))]
530            inc_vars.append(best_var)
531            not_inc_vars.remove(best_var)
532            changed_model = True
533    return inc_vars
534
535
536if __name__ == "__main__":
537
538    import Orange
539    from Orange.regression import linear
540
541    table = Orange.data.Table("housing.tab")
542    c = LinearRegressionLearner(table)
543    print c
Note: See TracBrowser for help on using the repository browser.