1# -*- coding: utf-8 -*-
2"""
3
4Created on Thu Jul 18 14:57:46 2013
5
6Author: Josef Perktold
7"""
8
9import numpy as np
10
11from statsmodels.genmod.generalized_estimating_equations import GEE, GEEMargins
12
13from statsmodels.genmod.families import Gaussian
14
15from statsmodels.genmod.tests import gee_gaussian_simulation_check as gees
16
17da,va = gees.gen_gendat_ar0(0.6)()
18ga = Gaussian()
19lhs = np.array([[0., 1, 1, 0, 0],])
20rhs = np.r_[0.,]
21
22example = []
23if 'constraint' in example:
24    md = GEE(da.endog, da.exog, da.group, da.time, ga, va,
25                     constraint=(lhs, rhs))
26    mdf = md.fit()
27    print(mdf.summary())
28
29
30md2 = GEE(da.endog, da.exog, da.group, da.time, ga, va,
31                 constraint=None)
32mdf2 = md2.fit()
33print('\n\n')
34print(mdf2.summary())
35
36
37mdf2.use_t = False
38mdf2.df_resid = np.diff(mdf2.model.exog.shape)
39tt2 = mdf2.t_test(np.eye(len(mdf2.params)))
40# need main to get wald_test
41#print mdf2.wald_test(np.eye(len(mdf2.params))[1:])
42
43mdf2.predict(da.exog.mean(0), offset=0)
44# -0.10867809062890971
45
46marg2 = GEEMargins(mdf2, ())
47print(marg2.summary())
48
49
50mdf_nc = md2.fit(cov_type='naive')
51mdf_bc = md2.fit(cov_type='bias_reduced')
52
53mdf_nc.use_t = False
54mdf_nc.df_resid = np.diff(mdf2.model.exog.shape)
55mdf_bc.use_t = False
56mdf_bc.df_resid = np.diff(mdf2.model.exog.shape)
57
58tt_nc = mdf_nc.t_test(np.eye(len(mdf2.params)))
59tt_bc = mdf_bc.t_test(np.eye(len(mdf2.params)))
60
61print('\nttest robust')
62print(tt2)
63print('\nttest naive')
64print(tt_nc)
65print('\nttest bias corrected')
66print(tt_bc)
67
68print("\nbse after fit option ")
69bse = np.column_stack((mdf2.bse, mdf2.bse, mdf_nc.bse, mdf_bc.bse))
70print(bse)
71
72print("\nimplemented `standard_errors`")
73bse2 = np.column_stack((mdf2.bse, mdf2.standard_errors(),
74                                 mdf2.standard_errors(covariance_type='naive'),
75                                 mdf2.standard_errors(covariance_type='bias_reduced')))
76print(bse2)
77print("bse and `standard_errors` agree:", np.allclose(bse, bse2))
78
79print("\nimplied standard errors in t_test")
80bse1 = np.column_stack((mdf2.bse, tt2.sd, tt_nc.sd, tt_bc.sd))
81print(bse1)
82print("t_test uses correct cov_params:", np.allclose(bse1, bse2))
83