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