1# -*- coding: utf-8 -*-
2"""Examples for Regression Plots
3
4Author: Josef Perktold
5
6"""
7
8import numpy as np
9import statsmodels.api as sm
10import matplotlib.pyplot as plt
11
12from statsmodels.sandbox.regression.predstd import wls_prediction_std
13import statsmodels.graphics.regressionplots as smrp
14from statsmodels.graphics.tests.test_regressionplots import TestPlot
15
16#example from tut.ols with changes
17#fix a seed for these examples
18np.random.seed(9876789)
19
20# OLS non-linear curve but linear in parameters
21# ---------------------------------------------
22
23nsample = 100
24sig = 0.5
25x1 = np.linspace(0, 20, nsample)
26x2 = 5 + 3* np.random.randn(nsample)
27X = np.c_[x1, x2, np.sin(0.5*x1), (x2-5)**2, np.ones(nsample)]
28beta = [0.5, 0.5, 1, -0.04, 5.]
29y_true = np.dot(X, beta)
30y = y_true + sig * np.random.normal(size=nsample)
31
32#estimate only linear function, misspecified because of non-linear terms
33exog0 = sm.add_constant(np.c_[x1, x2], prepend=False)
34
35#    plt.figure()
36#    plt.plot(x1, y, 'o', x1, y_true, 'b-')
37
38res = sm.OLS(y, exog0).fit()
39#print res.params
40#print res.bse
41
42
43plot_old = 0 #True
44if plot_old:
45
46    #current bug predict requires call to model.results
47    #print res.model.predict
48    prstd, iv_l, iv_u = wls_prediction_std(res)
49    plt.plot(x1, res.fittedvalues, 'r-o')
50    plt.plot(x1, iv_u, 'r--')
51    plt.plot(x1, iv_l, 'r--')
52    plt.title('blue: true,   red: OLS')
53
54    plt.figure()
55    plt.plot(res.resid, 'o')
56    plt.title('Residuals')
57
58    fig2 = plt.figure()
59    ax = fig2.add_subplot(2,1,1)
60    #namestr = ' for %s' % self.name if self.name else ''
61    plt.plot(x1, res.resid, 'o')
62    ax.set_title('residuals versus exog')# + namestr)
63    ax = fig2.add_subplot(2,1,2)
64    plt.plot(x2, res.resid, 'o')
65
66    fig3 = plt.figure()
67    ax = fig3.add_subplot(2,1,1)
68    #namestr = ' for %s' % self.name if self.name else ''
69    plt.plot(x1, res.fittedvalues, 'o')
70    ax.set_title('Fitted values versus exog')# + namestr)
71    ax = fig3.add_subplot(2,1,2)
72    plt.plot(x2, res.fittedvalues, 'o')
73
74    fig4 = plt.figure()
75    ax = fig4.add_subplot(2,1,1)
76    #namestr = ' for %s' % self.name if self.name else ''
77    plt.plot(x1, res.fittedvalues + res.resid, 'o')
78    ax.set_title('Fitted values plus residuals versus exog')# + namestr)
79    ax = fig4.add_subplot(2,1,2)
80    plt.plot(x2, res.fittedvalues + res.resid, 'o')
81
82    # see http://www.itl.nist.gov/div898/software/dataplot/refman1/auxillar/partregr.htm
83    fig5 = plt.figure()
84    ax = fig5.add_subplot(2,1,1)
85    #namestr = ' for %s' % self.name if self.name else ''
86    res1a = sm.OLS(y, exog0[:,[0,2]]).fit()
87    res1b = sm.OLS(x1, exog0[:,[0,2]]).fit()
88    plt.plot(res1b.resid, res1a.resid, 'o')
89    res1c = sm.OLS(res1a.resid, res1b.resid).fit()
90    plt.plot(res1b.resid, res1c.fittedvalues, '-')
91    ax.set_title('Partial Regression plot')# + namestr)
92    ax = fig5.add_subplot(2,1,2)
93    #plt.plot(x2, res.fittedvalues + res.resid, 'o')
94    res2a = sm.OLS(y, exog0[:,[0,1]]).fit()
95    res2b = sm.OLS(x2, exog0[:,[0,1]]).fit()
96    plt.plot(res2b.resid, res2a.resid, 'o')
97    res2c = sm.OLS(res2a.resid, res2b.resid).fit()
98    plt.plot(res2b.resid, res2c.fittedvalues, '-')
99
100    # see http://www.itl.nist.gov/div898/software/dataplot/refman1/auxillar/ccpr.htm
101    fig6 = plt.figure()
102    ax = fig6.add_subplot(2,1,1)
103    #namestr = ' for %s' % self.name if self.name else ''
104    x1beta = x1*res.params[1]
105    x2beta = x2*res.params[2]
106    plt.plot(x1, x1beta + res.resid, 'o')
107    plt.plot(x1, x1beta, '-')
108    ax.set_title('X_i beta_i plus residuals versus exog (CCPR)')# + namestr)
109    ax = fig6.add_subplot(2,1,2)
110    plt.plot(x2, x2beta + res.resid, 'o')
111    plt.plot(x2, x2beta, '-')
112
113
114    #print res.summary()
115
116doplots = 1
117if doplots:
118    fig1 = smrp.plot_fit(res, 0, y_true=None)
119    smrp.plot_fit(res, 1, y_true=None)
120    smrp.plot_partregress_grid(res, exog_idx=[0,1])
121    smrp.plot_regress_exog(res, exog_idx=0)
122    smrp.plot_ccpr(res, exog_idx=0)
123    smrp.plot_ccpr_grid(res, exog_idx=[0,1])
124
125tp = TestPlot()
126tp.test_plot_fit()
127
128fig1 = smrp.plot_partregress_grid(res, exog_idx=[0,1])
129#add lowess
130ax = fig1.axes[0]
131y0 = ax.get_lines()[0]._y
132x0 = ax.get_lines()[0]._x
133lres = sm.nonparametric.lowess(y0, x0, frac=0.2)
134ax.plot(lres[:,0], lres[:,1], 'r', lw=1.5)
135ax = fig1.axes[1]
136y0 = ax.get_lines()[0]._y
137x0 = ax.get_lines()[0]._x
138lres = sm.nonparametric.lowess(y0, x0, frac=0.2)
139ax.plot(lres[:,0], lres[:,1], 'r', lw=1.5)
140
141#plt.show()
142