<s> """<STR_LIT>""" <EOL> from statsmodels . compat . python import iterkeys , range , lrange , lzip , zip <EOL> import numpy as np <EOL> from scipy import stats <EOL> import pandas as pd <EOL> from statsmodels . tools . decorators import ( cache_readonly , <EOL> resettable_cache ) <EOL> import statsmodels . base . model as base <EOL> import statsmodels . regression . linear_model as lm <EOL> import statsmodels . base . wrapper as wrap <EOL> from statsmodels . genmod import families <EOL> from statsmodels . genmod import cov_struct as cov_structs <EOL> import statsmodels . genmod . families . varfuncs as varfuncs <EOL> from statsmodels . genmod . families . links import Link <EOL> from statsmodels . tools . sm_exceptions import ( ConvergenceWarning , <EOL> IterationLimitWarning ) <EOL> import warnings <EOL> from statsmodels . graphics . _regressionplots_doc import ( <EOL> _plot_added_variable_doc , <EOL> _plot_partial_residuals_doc , <EOL> _plot_ceres_residuals_doc ) <EOL> class ParameterConstraint ( object ) : <EOL> """<STR_LIT>""" <EOL> def __init__ ( self , lhs , rhs , exog ) : <EOL> """<STR_LIT>""" <EOL> rhs = np . atleast_1d ( rhs . squeeze ( ) ) <EOL> if rhs . ndim > <NUM_LIT:1> : <EOL> raise ValueError ( "<STR_LIT>" <EOL> "<STR_LIT>" ) <EOL> if len ( rhs ) != lhs . shape [ <NUM_LIT:0> ] : <EOL> raise ValueError ( "<STR_LIT>" <EOL> "<STR_LIT>" <EOL> "<STR_LIT>" <EOL> "<STR_LIT>" ) <EOL> self . lhs = lhs <EOL> self . rhs = rhs <EOL> lhs_u , lhs_s , lhs_vt = np . linalg . svd ( lhs . T , full_matrices = <NUM_LIT:1> ) <EOL> self . lhs0 = lhs_u [ : , len ( lhs_s ) : ] <EOL> self . lhs1 = lhs_u [ : , <NUM_LIT:0> : len ( lhs_s ) ] <EOL> self . lhsf = np . hstack ( ( self . lhs0 , self . lhs1 ) ) <EOL> self . param0 = np . dot ( self . lhs1 , np . dot ( lhs_vt , self . rhs ) / <EOL> lhs_s ) <EOL> self . _offset_increment = np . dot ( exog , self . param0 ) <EOL> self . orig_exog = exog <EOL> self . exog_fulltrans = np . dot ( exog , self . lhsf ) <EOL> def offset_increment ( self ) : <EOL> """<STR_LIT>""" <EOL> return self . _offset_increment <EOL> def reduced_exog ( self ) : <EOL> """<STR_LIT>""" <EOL> return self . exog_fulltrans [ : , <NUM_LIT:0> : self . lhs0 . shape [ <NUM_LIT:1> ] ] <EOL> def restore_exog ( self ) : <EOL> """<STR_LIT>""" <EOL> return self . orig_exog <EOL> def unpack_param ( self , params ) : <EOL> """<STR_LIT>""" <EOL> return self . param0 + np . dot ( self . lhs0 , params ) <EOL> def unpack_cov ( self , bcov ) : <EOL> """<STR_LIT>""" <EOL> return np . dot ( self . lhs0 , np . dot ( bcov , self . lhs0 . T ) ) <EOL> _gee_init_doc = """<STR_LIT>""" <EOL> _gee_family_doc = """<STR_LIT>""" <EOL> _gee_ordinal_family_doc = """<STR_LIT>""" <EOL> _gee_nominal_family_doc = """<STR_LIT>""" <EOL> _gee_fit_doc = """<STR_LIT>""" <EOL> _gee_results_doc = """<STR_LIT>""" <EOL> _gee_example = """<STR_LIT>""" <EOL> _gee_ordinal_example = """<STR_LIT>""" <EOL> _gee_nominal_example = """<STR_LIT>""" <EOL> class GEE ( base . Model ) : <EOL> __doc__ = ( <EOL> "<STR_LIT>" <EOL> "<STR_LIT>" + _gee_init_doc % <EOL> { '<STR_LIT>' : base . _missing_param_doc , <EOL> '<STR_LIT>' : _gee_family_doc , <EOL> '<STR_LIT>' : _gee_example } ) <EOL> cached_means = None <EOL> def __init__ ( self , endog , exog , groups , time = None , family = None , <EOL> cov_struct = None , missing = '<STR_LIT:none>' , offset = None , <EOL> exposure = None , dep_data = None , constraint = None , <EOL> update_dep = True , weights = None , ** kwargs ) : <EOL> if ( family is not None ) and not isinstance ( family . link , tuple ( family . safe_links ) ) : <EOL> import warnings <EOL> warnings . warn ( "<STR_LIT>" % <EOL> ( family . link . __class__ . __name__ , family . __class__ . __name__ ) ) <EOL> self . missing = missing <EOL> self . dep_data = dep_data <EOL> self . constraint = constraint <EOL> self . update_dep = update_dep <EOL> groups = np . array ( groups ) <EOL> super ( GEE , self ) . __init__ ( endog , exog , groups = groups , <EOL> time = time , offset = offset , <EOL> exposure = exposure , weights = weights , <EOL> dep_data = dep_data , missing = missing , <EOL> ** kwargs ) <EOL> self . _init_keys . extend ( [ "<STR_LIT>" , "<STR_LIT>" , "<STR_LIT>" , <EOL> "<STR_LIT>" ] ) <EOL> if family is None : <EOL> family = families . Gaussian ( ) <EOL> else : <EOL> if not issubclass ( family . __class__ , families . Family ) : <EOL> raise ValueError ( "<STR_LIT>" <EOL> "<STR_LIT>" ) <EOL> self . family = family <EOL> if cov_struct is None : <EOL> cov_struct = cov_structs . Independence ( ) <EOL> else : <EOL> if not issubclass ( cov_struct . __class__ , cov_structs . CovStruct ) : <EOL> raise ValueError ( "<STR_LIT>" <EOL> "<STR_LIT>" ) <EOL> self . cov_struct = cov_struct <EOL> self . _offset_exposure = None <EOL> if offset is not None : <EOL> self . _offset_exposure = self . offset . copy ( ) <EOL> self . offset = offset <EOL> if exposure is not None : <EOL> if not isinstance ( self . family . link , families . links . Log ) : <EOL> raise ValueError ( "<STR_LIT>" ) <EOL> if self . _offset_exposure is not None : <EOL> self . _offset_exposure += np . log ( exposure ) <EOL> else : <EOL> self . _offset_exposure = np . log ( exposure ) <EOL> self . exposure = exposure <EOL> self . constraint = None <EOL> if constraint is not None : <EOL> if len ( constraint ) != <NUM_LIT:2> : <EOL> raise ValueError ( "<STR_LIT>" ) <EOL> if constraint [ <NUM_LIT:0> ] . shape [ <NUM_LIT:1> ] != self . exog . shape [ <NUM_LIT:1> ] : <EOL> raise ValueError ( "<STR_LIT>" <EOL> "<STR_LIT>" <EOL> "<STR_LIT>" ) <EOL> self . constraint = ParameterConstraint ( constraint [ <NUM_LIT:0> ] , <EOL> constraint [ <NUM_LIT:1> ] , <EOL> self . exog ) <EOL> if self . _offset_exposure is not None : <EOL> self . _offset_exposure += self . constraint . offset_increment ( ) <EOL> else : <EOL> self . _offset_exposure = self . constraint . offset_increment ( ) . copy ( ) <EOL> self . exog = self . constraint . reduced_exog ( ) <EOL> group_labels , ix = np . unique ( self . groups , return_inverse = True ) <EOL> se = pd . Series ( index = np . arange ( len ( ix ) ) ) <EOL> gb = se . groupby ( ix ) . groups <EOL> dk = [ ( lb , np . asarray ( gb [ k ] ) ) for k , lb in enumerate ( group_labels ) ] <EOL> self . group_indices = dict ( dk ) <EOL> self . group_labels = group_labels <EOL> self . endog_li = self . cluster_list ( self . endog ) <EOL> self . exog_li = self . cluster_list ( self . exog ) <EOL> if self . weights is not None : <EOL> self . weights_li = self . cluster_list ( self . weights ) <EOL> self . weights_li = [ x [ <NUM_LIT:0> ] for x in self . weights_li ] <EOL> self . weights_li = np . asarray ( self . weights_li ) <EOL> self . num_group = len ( self . endog_li ) <EOL> if self . time is not None : <EOL> self . time = np . asarray ( self . time , np . float64 ) <EOL> if self . time . ndim == <NUM_LIT:1> : <EOL> self . time = self . time [ : , None ] <EOL> self . time_li = self . cluster_list ( self . time ) <EOL> else : <EOL> self . time_li = [ np . arange ( len ( y ) , dtype = np . float64 ) [ : , None ] <EOL> for y in self . endog_li ] <EOL> self . time = np . concatenate ( self . time_li ) <EOL> if self . _offset_exposure is not None : <EOL> self . offset_li = self . cluster_list ( self . _offset_exposure ) <EOL> else : <EOL> self . offset_li = None <EOL> if constraint is not None : <EOL> self . constraint . exog_fulltrans_li = self . cluster_list ( self . constraint . exog_fulltrans ) <EOL> self . family = family <EOL> self . cov_struct . initialize ( self ) <EOL> group_ns = [ len ( y ) for y in self . endog_li ] <EOL> self . nobs = sum ( group_ns ) <EOL> self . df_model = self . exog . shape [ <NUM_LIT:1> ] - <NUM_LIT:1> <EOL> self . df_resid = self . nobs - self . exog . shape [ <NUM_LIT:1> ] <EOL> maxgroup = max ( [ len ( x ) for x in self . endog_li ] ) <EOL> if maxgroup == <NUM_LIT:1> : <EOL> self . update_dep = False <EOL> @ classmethod <EOL> def from_formula ( cls , formula , groups , data , subset = None , <EOL> time = None , offset = None , exposure = None , <EOL> * args , ** kwargs ) : <EOL> """<STR_LIT>""" % { '<STR_LIT>' : base . _missing_param_doc } <EOL> if type ( groups ) == str : <EOL> groups = data [ groups ] <EOL> if type ( time ) == str : <EOL> time = data [ time ] <EOL> if type ( offset ) == str : <EOL> offset = data [ offset ] <EOL> if type ( exposure ) == str : <EOL> exposure = data [ exposure ] <EOL> model = super ( GEE , cls ) . from_formula ( formula , data , subset , <EOL> groups , time = time , <EOL> offset = offset , <EOL> exposure = exposure , <EOL> * args , ** kwargs ) <EOL> return model <EOL> def cluster_list ( self , array ) : <EOL> """<STR_LIT>""" <EOL> if array . ndim == <NUM_LIT:1> : <EOL> return [ np . array ( array [ self . group_indices [ k ] ] ) <EOL> for k in self . group_labels ] <EOL> else : <EOL> return [ np . array ( array [ self . group_indices [ k ] , : ] ) <EOL> for k in self . group_labels ] <EOL> def estimate_scale ( self ) : <EOL> """<STR_LIT>""" <EOL> if isinstance ( self . family , ( families . Binomial , families . Poisson , <EOL> _Multinomial ) ) : <EOL> return <NUM_LIT:1.> <EOL> endog = self . endog_li <EOL> exog = self . exog_li <EOL> cached_means = self . cached_means <EOL> nobs = self . nobs <EOL> varfunc = self . family . variance <EOL> scale = <NUM_LIT:0.> <EOL> fsum = <NUM_LIT:0.> <EOL> for i in range ( self . num_group ) : <EOL> if len ( endog [ i ] ) == <NUM_LIT:0> : <EOL> continue <EOL> expval , _ = cached_means [ i ] <EOL> f = self . weights_li [ i ] if self . weights is not None else <NUM_LIT:1.> <EOL> sdev = np . sqrt ( varfunc ( expval ) ) <EOL> resid = ( endog [ i ] - expval ) / sdev <EOL> scale += f * np . sum ( resid ** <NUM_LIT:2> ) <EOL> fsum += f * len ( endog [ i ] ) <EOL> scale /= ( fsum * ( nobs - self . ddof_scale ) / float ( nobs ) ) <EOL> return scale <EOL> def mean_deriv ( self , exog , lin_pred ) : <EOL> """<STR_LIT>""" <EOL> idl = self . family . link . inverse_deriv ( lin_pred ) <EOL> dmat = exog * idl [ : , None ] <EOL> return dmat <EOL> def mean_deriv_exog ( self , exog , params , offset_exposure = None ) : <EOL> """<STR_LIT>""" <EOL> lin_pred = np . dot ( exog , params ) <EOL> if offset_exposure is not None : <EOL> lin_pred += offset_exposure <EOL> idl = self . family . link . inverse_deriv ( lin_pred ) <EOL> dmat = np . outer ( idl , params ) <EOL> return dmat <EOL> def _update_mean_params ( self ) : <EOL> """<STR_LIT>""" <EOL> endog = self . endog_li <EOL> exog = self . exog_li <EOL> cached_means = self . cached_means <EOL> varfunc = self . family . variance <EOL> bmat , score = <NUM_LIT:0> , <NUM_LIT:0> <EOL> for i in range ( self . num_group ) : <EOL> expval , lpr = cached_means [ i ] <EOL> resid = endog [ i ] - expval <EOL> dmat = self . mean_deriv ( exog [ i ] , lpr ) <EOL> sdev = np . sqrt ( varfunc ( expval ) ) <EOL> rslt = self . cov_struct . covariance_matrix_solve ( expval , i , <EOL> sdev , ( dmat , resid ) ) <EOL> if rslt is None : <EOL> return None , None <EOL> vinv_d , vinv_resid = tuple ( rslt ) <EOL> f = self . weights_li [ i ] if self . weights is not None else <NUM_LIT:1.> <EOL> bmat += f * np . dot ( dmat . T , vinv_d ) <EOL> score += f * np . dot ( dmat . T , vinv_resid ) <EOL> update = np . linalg . solve ( bmat , score ) <EOL> self . _fit_history [ "<STR_LIT>" ] . append ( <EOL> self . cov_struct . cov_adjust ) <EOL> return update , score <EOL> def update_cached_means ( self , mean_params ) : <EOL> """<STR_LIT>""" <EOL> endog = self . endog_li <EOL> exog = self . exog_li <EOL> offset = self . offset_li <EOL> linkinv = self . family . link . inverse <EOL> self . cached_means = [ ] <EOL> for i in range ( self . num_group ) : <EOL> if len ( endog [ i ] ) == <NUM_LIT:0> : <EOL> continue <EOL> lpr = np . dot ( exog [ i ] , mean_params ) <EOL> if offset is not None : <EOL> lpr += offset [ i ] <EOL> expval = linkinv ( lpr ) <EOL> self . cached_means . append ( ( expval , lpr ) ) <EOL> def _covmat ( self ) : <EOL> """<STR_LIT>""" <EOL> endog = self . endog_li <EOL> exog = self . exog_li <EOL> varfunc = self . family . variance <EOL> cached_means = self . cached_means <EOL> bmat , cmat = <NUM_LIT:0> , <NUM_LIT:0> <EOL> for i in range ( self . num_group ) : <EOL> expval , lpr = cached_means [ i ] <EOL> resid = endog [ i ] - expval <EOL> dmat = self . mean_deriv ( exog [ i ] , lpr ) <EOL> sdev = np . sqrt ( varfunc ( expval ) ) <EOL> rslt = self . cov_struct . covariance_matrix_solve ( expval , i , <EOL> sdev , ( dmat , resid ) ) <EOL> if rslt is None : <EOL> return None , None , None , None <EOL> vinv_d , vinv_resid = tuple ( rslt ) <EOL> f = self . weights_li [ i ] if self . weights is not None else <NUM_LIT:1.> <EOL> bmat += f * np . dot ( dmat . T , vinv_d ) <EOL> dvinv_resid = f * np . dot ( dmat . T , vinv_resid ) <EOL> cmat += np . outer ( dvinv_resid , dvinv_resid ) <EOL> scale = self . estimate_scale ( ) <EOL> bmati = np . linalg . inv ( bmat ) <EOL> cov_naive = bmati * scale <EOL> cov_robust = np . dot ( bmati , np . dot ( cmat , bmati ) ) <EOL> bcm = <NUM_LIT:0> <EOL> for i in range ( self . num_group ) : <EOL> expval , lpr = cached_means [ i ] <EOL> resid = endog [ i ] - expval <EOL> dmat = self . mean_deriv ( exog [ i ] , lpr ) <EOL> sdev = np . sqrt ( varfunc ( expval ) ) <EOL> rslt = self . cov_struct . covariance_matrix_solve ( expval , <EOL> i , sdev , ( dmat , ) ) <EOL> if rslt is None : <EOL> return None , None , None , None <EOL> vinv_d = rslt [ <NUM_LIT:0> ] <EOL> vinv_d /= scale <EOL> hmat = np . dot ( vinv_d , cov_naive ) <EOL> hmat = np . dot ( hmat , dmat . T ) . T <EOL> f = self . weights_li [ i ] if self . weights is not None else <NUM_LIT:1.> <EOL> aresid = np . linalg . solve ( np . eye ( len ( resid ) ) - hmat , resid ) <EOL> rslt = self . cov_struct . covariance_matrix_solve ( expval , i , <EOL> sdev , ( aresid , ) ) <EOL> if rslt is None : <EOL> return None , None , None , None <EOL> srt = rslt [ <NUM_LIT:0> ] <EOL> srt = f * np . dot ( dmat . T , srt ) / scale <EOL> bcm += np . outer ( srt , srt ) <EOL> cov_robust_bc = np . dot ( cov_naive , np . dot ( bcm , cov_naive ) ) <EOL> cov_naive *= self . scaling_factor <EOL> cov_robust *= self . scaling_factor <EOL> cov_robust_bc *= self . scaling_factor <EOL> return ( cov_robust , cov_naive , cov_robust_bc , cmat ) <EOL> def predict ( self , params , exog = None , offset = None , <EOL> exposure = None , linear = False ) : <EOL> """<STR_LIT>""" <EOL> if exposure is not None and not isinstance ( self . family . link , families . links . Log ) : <EOL> raise ValueError ( "<STR_LIT>" ) <EOL> _offset = <NUM_LIT:0.> <EOL> if exog is None : <EOL> exog = self . exog <EOL> if not isinstance ( self . family . link , families . links . Log ) : <EOL> if offset is None : <EOL> if self . _offset_exposure is not None : <EOL> _offset = self . _offset_exposure . copy ( ) <EOL> else : <EOL> _offset = offset <EOL> else : <EOL> if offset is None and exposure is None : <EOL> if self . _offset_exposure is not None : <EOL> _offset = self . _offset_exposure <EOL> elif offset is None and exposure is not None : <EOL> _offset = np . log ( exposure ) <EOL> if hasattr ( self , "<STR_LIT>" ) : <EOL> _offset = _offset + self . offset <EOL> elif offset is not None and exposure is None : <EOL> _offset = offset <EOL> if hasattr ( self , "<STR_LIT>" ) : <EOL> _offset = offset + np . log ( self . exposure ) <EOL> else : <EOL> _offset = offset + np . log ( exposure ) <EOL> else : <EOL> if offset is not None : <EOL> _offset = _offset + offset <EOL> if exposure is not None : <EOL> _offset += np . log ( exposure ) <EOL> lin_pred = _offset + np . dot ( exog , params ) <EOL> if not linear : <EOL> return self . family . link . inverse ( lin_pred ) <EOL> return lin_pred <EOL> def _starting_params ( self ) : <EOL> return np . zeros ( self . exog . shape [ <NUM_LIT:1> ] ) <EOL> def fit ( self , maxiter = <NUM_LIT> , ctol = <NUM_LIT> , start_params = None , <EOL> params_niter = <NUM_LIT:1> , first_dep_update = <NUM_LIT:0> , <EOL> cov_type = '<STR_LIT>' , ddof_scale = None , scaling_factor = <NUM_LIT:1.> ) : <EOL> if ddof_scale is None : <EOL> self . ddof_scale = self . exog . shape [ <NUM_LIT:1> ] <EOL> else : <EOL> if not ddof_scale >= <NUM_LIT:0> : <EOL> raise ValueError ( "<STR_LIT>" ) <EOL> self . ddof_scale = ddof_scale <EOL> self . scaling_factor = scaling_factor <EOL> self . _fit_history = { '<STR_LIT>' : [ ] , <EOL> '<STR_LIT>' : [ ] , <EOL> '<STR_LIT>' : [ ] , <EOL> '<STR_LIT>' : [ ] } <EOL> if self . weights is not None and cov_type == '<STR_LIT>' : <EOL> raise ValueError ( "<STR_LIT>" ) <EOL> if start_params is None : <EOL> mean_params = self . _starting_params ( ) <EOL> else : <EOL> start_params = np . asarray ( start_params ) <EOL> mean_params = start_params . copy ( ) <EOL> self . update_cached_means ( mean_params ) <EOL> del_params = - <NUM_LIT:1.> <EOL> num_assoc_updates = <NUM_LIT:0> <EOL> for itr in range ( maxiter ) : <EOL> update , score = self . _update_mean_params ( ) <EOL> if update is None : <EOL> warnings . warn ( "<STR_LIT>" , <EOL> ConvergenceWarning ) <EOL> break <EOL> mean_params += update <EOL> self . update_cached_means ( mean_params ) <EOL> del_params = np . sqrt ( np . sum ( score ** <NUM_LIT:2> ) ) <EOL> self . _fit_history [ '<STR_LIT>' ] . append ( mean_params . copy ( ) ) <EOL> self . _fit_history [ '<STR_LIT>' ] . append ( score ) <EOL> self . _fit_history [ '<STR_LIT>' ] . append ( <EOL> self . cov_struct . dep_params ) <EOL> if ( del_params < ctol and <EOL> ( num_assoc_updates > <NUM_LIT:0> or self . update_dep == False ) ) : <EOL> break <EOL> if ( self . update_dep and ( itr % params_niter ) == <NUM_LIT:0> <EOL> and ( itr >= first_dep_update ) ) : <EOL> self . _update_assoc ( mean_params ) <EOL> num_assoc_updates += <NUM_LIT:1> <EOL> if del_params >= ctol : <EOL> warnings . warn ( "<STR_LIT>" , <EOL> IterationLimitWarning ) <EOL> if mean_params is None : <EOL> warnings . warn ( "<STR_LIT>" , <EOL> ConvergenceWarning ) <EOL> return None <EOL> bcov , ncov , bc_cov , _ = self . _covmat ( ) <EOL> if bcov is None : <EOL> warnings . warn ( "<STR_LIT>" <EOL> "<STR_LIT>" , ConvergenceWarning ) <EOL> return None <EOL> if self . constraint is not None : <EOL> mean_params , bcov = self . _handle_constraint ( mean_params , bcov ) <EOL> if mean_params is None : <EOL> warnings . warn ( "<STR_LIT>" <EOL> "<STR_LIT>" , ConvergenceWarning ) <EOL> return None <EOL> scale = self . estimate_scale ( ) <EOL> res_kwds = dict ( cov_type = cov_type , <EOL> cov_robust = bcov , <EOL> cov_naive = ncov , <EOL> cov_robust_bc = bc_cov ) <EOL> results = GEEResults ( self , mean_params , bcov / scale , scale , <EOL> cov_type = cov_type , use_t = False , <EOL> attr_kwds = res_kwds ) <EOL> results . fit_history = self . _fit_history <EOL> delattr ( self , "<STR_LIT>" ) <EOL> results . score_norm = del_params <EOL> results . converged = ( del_params < ctol ) <EOL> results . cov_struct = self . cov_struct <EOL> results . params_niter = params_niter <EOL> results . first_dep_update = first_dep_update <EOL> results . ctol = ctol <EOL> results . maxiter = maxiter <EOL> results . _props = [ "<STR_LIT>" , "<STR_LIT>" , <EOL> "<STR_LIT>" , "<STR_LIT>" , <EOL> "<STR_LIT>" , "<STR_LIT>" , <EOL> "<STR_LIT>" , <EOL> "<STR_LIT>" , "<STR_LIT>" , "<STR_LIT>" , <EOL> "<STR_LIT>" , "<STR_LIT>" , "<STR_LIT>" , <EOL> "<STR_LIT>" ] <EOL> return GEEResultsWrapper ( results ) <EOL> fit . __doc__ = _gee_fit_doc <EOL> def _handle_constraint ( self , mean_params , bcov ) : <EOL> """<STR_LIT>""" <EOL> red_p = len ( mean_params ) <EOL> full_p = self . constraint . lhs . shape [ <NUM_LIT:1> ] <EOL> mean_params0 = np . r_ [ mean_params , np . zeros ( full_p - red_p ) ] <EOL> save_exog_li = self . exog_li <EOL> self . exog_li = self . constraint . exog_fulltrans_li <EOL> import copy <EOL> save_cached_means = copy . deepcopy ( self . cached_means ) <EOL> self . update_cached_means ( mean_params0 ) <EOL> _ , score = self . _update_mean_params ( ) <EOL> if score is None : <EOL> warnings . warn ( "<STR_LIT>" , <EOL> ConvergenceWarning ) <EOL> return None , None <EOL> _ , ncov1 , _ , cmat = self . _covmat ( ) <EOL> scale = self . estimate_scale ( ) <EOL> cmat = cmat / scale ** <NUM_LIT:2> <EOL> score2 = score [ red_p : ] / scale <EOL> amat = np . linalg . inv ( ncov1 ) <EOL> bmat_11 = cmat [ <NUM_LIT:0> : red_p , <NUM_LIT:0> : red_p ] <EOL> bmat_22 = cmat [ red_p : , red_p : ] <EOL> bmat_12 = cmat [ <NUM_LIT:0> : red_p , red_p : ] <EOL> amat_11 = amat [ <NUM_LIT:0> : red_p , <NUM_LIT:0> : red_p ] <EOL> amat_12 = amat [ <NUM_LIT:0> : red_p , red_p : ] <EOL> score_cov = bmat_22 - np . dot ( amat_12 . T , np . linalg . solve ( amat_11 , bmat_12 ) ) <EOL> score_cov -= np . dot ( bmat_12 . T , <EOL> np . linalg . solve ( amat_11 , amat_12 ) ) <EOL> score_cov += np . dot ( amat_12 . T , <EOL> np . dot ( np . linalg . solve ( amat_11 , bmat_11 ) , <EOL> np . linalg . solve ( amat_11 , amat_12 ) ) ) <EOL> from scipy . stats . distributions import chi2 <EOL> score_statistic = np . dot ( score2 , <EOL> np . linalg . solve ( score_cov , score2 ) ) <EOL> score_df = len ( score2 ) <EOL> score_pvalue = <NUM_LIT:1> - chi2 . cdf ( score_statistic , score_df ) <EOL> self . score_test_results = { "<STR_LIT>" : score_statistic , <EOL> "<STR_LIT>" : score_df , <EOL> "<STR_LIT>" : score_pvalue } <EOL> mean_params = self . constraint . unpack_param ( mean_params ) <EOL> bcov = self . constraint . unpack_cov ( bcov ) <EOL> self . exog_li = save_exog_li <EOL> self . cached_means = save_cached_means <EOL> self . exog = self . constraint . restore_exog ( ) <EOL> return mean_params , bcov <EOL> def _update_assoc ( self , params ) : <EOL> """<STR_LIT>""" <EOL> self . cov_struct . update ( params ) <EOL> def _derivative_exog ( self , params , exog = None , transform = '<STR_LIT>' , <EOL> dummy_idx = None , count_idx = None ) : <EOL> """<STR_LIT>""" <EOL> offset_exposure = None <EOL> if exog is None : <EOL> exog = self . exog <EOL> offset_exposure = self . _offset_exposure <EOL> margeff = self . mean_deriv_exog ( exog , params , offset_exposure ) <EOL> if '<STR_LIT>' in transform : <EOL> margeff *= exog <EOL> if '<STR_LIT>' in transform : <EOL> margeff /= self . predict ( params , exog ) [ : , None ] <EOL> if count_idx is not None : <EOL> from statsmodels . discrete . discrete_margins import ( <EOL> _get_count_effects ) <EOL> margeff = _get_count_effects ( margeff , exog , count_idx , transform , <EOL> self , params ) <EOL> if dummy_idx is not None : <EOL> from statsmodels . discrete . discrete_margins import ( <EOL> _get_dummy_effects ) <EOL> margeff = _get_dummy_effects ( margeff , exog , dummy_idx , transform , <EOL> self , params ) <EOL> return margeff <EOL> class GEEResults ( base . LikelihoodModelResults ) : <EOL> __doc__ = ( <EOL> "<STR_LIT>" <EOL> + _gee_results_doc ) <EOL> def __init__ ( self , model , params , cov_params , scale , <EOL> cov_type = '<STR_LIT>' , use_t = False , ** kwds ) : <EOL> super ( GEEResults , self ) . __init__ ( model , params , <EOL> normalized_cov_params = cov_params , scale = scale ) <EOL> self . df_resid = model . df_resid <EOL> self . df_model = model . df_model <EOL> self . family = model . family <EOL> attr_kwds = kwds . pop ( '<STR_LIT>' , { } ) <EOL> self . __dict__ . update ( attr_kwds ) <EOL> if not ( hasattr ( self , '<STR_LIT>' ) and <EOL> hasattr ( self , '<STR_LIT>' ) ) : <EOL> self . cov_type = cov_type <EOL> covariance_type = self . cov_type . lower ( ) <EOL> allowed_covariances = [ "<STR_LIT>" , "<STR_LIT>" , "<STR_LIT>" ] <EOL> if covariance_type not in allowed_covariances : <EOL> msg = "<STR_LIT>" + "<STR_LIT:U+002CU+0020>" . join ( allowed_covariances ) <EOL> raise ValueError ( msg ) <EOL> if cov_type == "<STR_LIT>" : <EOL> cov = self . cov_robust <EOL> elif cov_type == "<STR_LIT>" : <EOL> cov = self . cov_naive <EOL> elif cov_type == "<STR_LIT>" : <EOL> cov = self . cov_robust_bc <EOL> self . cov_params_default = cov <EOL> else : <EOL> if self . cov_type != cov_type : <EOL> raise ValueError ( '<STR_LIT>' <EOL> '<STR_LIT>' ) <EOL> def standard_errors ( self , cov_type = "<STR_LIT>" ) : <EOL> """<STR_LIT>""" <EOL> covariance_type = cov_type . lower ( ) <EOL> allowed_covariances = [ "<STR_LIT>" , "<STR_LIT>" , "<STR_LIT>" ] <EOL> if covariance_type not in allowed_covariances : <EOL> msg = "<STR_LIT>" + "<STR_LIT:U+002CU+0020>" . join ( allowed_covariances ) <EOL> raise ValueError ( msg ) <EOL> if covariance_type == "<STR_LIT>" : <EOL> return np . sqrt ( np . diag ( self . cov_robust ) ) <EOL> elif covariance_type == "<STR_LIT>" : <EOL> return np . sqrt ( np . diag ( self . cov_naive ) ) <EOL> elif covariance_type == "<STR_LIT>" : <EOL> return np . sqrt ( np . diag ( self . cov_robust_bc ) ) <EOL> @ cache_readonly <EOL> def bse ( self ) : <EOL> return self . standard_errors ( self . cov_type ) <EOL> @ cache_readonly <EOL> def resid ( self ) : <EOL> """<STR_LIT>""" <EOL> return self . model . endog - self . fittedvalues <EOL> @ cache_readonly <EOL> def resid_split ( self ) : <EOL> """<STR_LIT>""" <EOL> sresid = [ ] <EOL> for v in self . model . group_labels : <EOL> ii = self . model . group_indices [ v ] <EOL> sresid . append ( self . resid [ ii ] ) <EOL> return sresid <EOL> @ cache_readonly <EOL> def resid_centered ( self ) : <EOL> """<STR_LIT>""" <EOL> cresid = self . resid . copy ( ) <EOL> for v in self . model . group_labels : <EOL> ii = self . model . group_indices [ v ] <EOL> cresid [ ii ] -= cresid [ ii ] . mean ( ) <EOL> return cresid <EOL> @ cache_readonly <EOL> def resid_centered_split ( self ) : <EOL> """<STR_LIT>""" <EOL> sresid = [ ] <EOL> for v in self . model . group_labels : <EOL> ii = self . model . group_indices [ v ] <EOL> sresid . append ( self . centered_resid [ ii ] ) <EOL> return sresid <EOL> split_resid = resid_split <EOL> centered_resid = resid_centered <EOL> split_centered_resid = resid_centered_split <EOL> @ cache_readonly <EOL> def resid_response ( self ) : <EOL> return self . model . endog - self . fittedvalues <EOL> @ cache_readonly <EOL> def resid_pearson ( self ) : <EOL> val = self . model . endog - self . fittedvalues <EOL> val = val / np . sqrt ( self . family . variance ( self . fittedvalues ) ) <EOL> return val <EOL> @ cache_readonly <EOL> def resid_working ( self ) : <EOL> val = self . resid_response <EOL> val = val / self . family . link . deriv ( self . fittedvalues ) <EOL> return val <EOL> @ cache_readonly <EOL> def resid_anscombe ( self ) : <EOL> return self . family . resid_anscombe ( self . model . endog , self . fittedvalues ) <EOL> @ cache_readonly <EOL> def resid_deviance ( self ) : <EOL> return self . family . resid_dev ( self . model . endog , self . fittedvalues ) <EOL> @ cache_readonly <EOL> def fittedvalues ( self ) : <EOL> """<STR_LIT>""" <EOL> return self . model . family . link . inverse ( np . dot ( self . model . exog , <EOL> self . params ) ) <EOL> def plot_added_variable ( self , focus_exog , resid_type = None , <EOL> use_glm_weights = True , fit_kwargs = None , <EOL> ax = None ) : <EOL> from statsmodels . graphics . regressionplots import plot_added_variable <EOL> fig = plot_added_variable ( self , focus_exog , <EOL> resid_type = resid_type , <EOL> use_glm_weights = use_glm_weights , <EOL> fit_kwargs = fit_kwargs , ax = ax ) <EOL> return fig <EOL> plot_added_variable . __doc__ = _plot_added_variable_doc % { <EOL> '<STR_LIT>' : '<STR_LIT>' } <EOL> def plot_partial_residuals ( self , focus_exog , ax = None ) : <EOL> from statsmodels . graphics . regressionplots import plot_partial_residuals <EOL> return plot_partial_residuals ( self , focus_exog , ax = ax ) <EOL> plot_partial_residuals . __doc__ = _plot_partial_residuals_doc % { <EOL> '<STR_LIT>' : '<STR_LIT>' } <EOL> def plot_ceres_residuals ( self , focus_exog , frac = <NUM_LIT> , cond_means = None , <EOL> ax = None ) : <EOL> from statsmodels . graphics . regressionplots import plot_ceres_residuals <EOL> return plot_ceres_residuals ( self , focus_exog , frac , <EOL> cond_means = cond_means , ax = ax ) <EOL> plot_ceres_residuals . __doc__ = _plot_ceres_residuals_doc % { <EOL> '<STR_LIT>' : '<STR_LIT>' } <EOL> def conf_int ( self , alpha = <NUM_LIT> , cols = None , cov_type = None ) : <EOL> """<STR_LIT>""" <EOL> if cov_type is None : <EOL> bse = self . bse <EOL> else : <EOL> bse = self . standard_errors ( cov_type = cov_type ) <EOL> params = self . params <EOL> dist = stats . norm <EOL> q = dist . ppf ( <NUM_LIT:1> - alpha / <NUM_LIT:2> ) <EOL> if cols is None : <EOL> lower = self . params - q * bse <EOL> upper = self . params + q * bse <EOL> else : <EOL> cols = np . asarray ( cols ) <EOL> lower = params [ cols ] - q * bse [ cols ] <EOL> upper = params [ cols ] + q * bse [ cols ] <EOL> return np . asarray ( lzip ( lower , upper ) ) <EOL> def summary ( self , yname = None , xname = None , title = None , alpha = <NUM_LIT> ) : <EOL> """<STR_LIT>""" <EOL> top_left = [ ( '<STR_LIT>' , None ) , <EOL> ( '<STR_LIT>' , None ) , <EOL> ( '<STR_LIT>' , [ '<STR_LIT>' ] ) , <EOL> ( '<STR_LIT>' , [ '<STR_LIT>' ] ) , <EOL> ( '<STR_LIT>' , [ self . model . family . __class__ . __name__ ] ) , <EOL> ( '<STR_LIT>' , <EOL> [ self . model . cov_struct . __class__ . __name__ ] ) , <EOL> ( '<STR_LIT>' , None ) , <EOL> ( '<STR_LIT>' , [ self . cov_type , ] ) <EOL> ] <EOL> NY = [ len ( y ) for y in self . model . endog_li ] <EOL> top_right = [ ( '<STR_LIT>' , [ sum ( NY ) ] ) , <EOL> ( '<STR_LIT>' , [ len ( self . model . endog_li ) ] ) , <EOL> ( '<STR_LIT>' , [ min ( NY ) ] ) , <EOL> ( '<STR_LIT>' , [ max ( NY ) ] ) , <EOL> ( '<STR_LIT>' , [ "<STR_LIT>" % np . mean ( NY ) ] ) , <EOL> ( '<STR_LIT>' , [ '<STR_LIT>' % <EOL> len ( self . fit_history [ '<STR_LIT>' ] ) ] ) , <EOL> ( '<STR_LIT>' , [ "<STR_LIT>" % self . scale ] ) , <EOL> ( '<STR_LIT>' , None ) , <EOL> ] <EOL> skew1 = stats . skew ( self . resid ) <EOL> kurt1 = stats . kurtosis ( self . resid ) <EOL> skew2 = stats . skew ( self . centered_resid ) <EOL> kurt2 = stats . kurtosis ( self . centered_resid ) <EOL> diagn_left = [ ( '<STR_LIT>' , [ "<STR_LIT>" % skew1 ] ) , <EOL> ( '<STR_LIT>' , [ "<STR_LIT>" % skew2 ] ) ] <EOL> diagn_right = [ ( '<STR_LIT>' , [ "<STR_LIT>" % kurt1 ] ) , <EOL> ( '<STR_LIT>' , [ "<STR_LIT>" % kurt2 ] ) <EOL> ] <EOL> if title is None : <EOL> title = self . model . __class__ . __name__ + '<STR_LIT:U+0020>' + "<STR_LIT>" <EOL> if xname is not None : <EOL> xna = xname <EOL> else : <EOL> xna = self . model . exog_names <EOL> from statsmodels . iolib . summary import Summary <EOL> smry = Summary ( ) <EOL> smry . add_table_2cols ( self , gleft = top_left , gright = top_right , <EOL> yname = self . model . endog_names , xname = xna , <EOL> title = title ) <EOL> smry . add_table_params ( self , yname = yname , xname = xna , <EOL> alpha = alpha , use_t = False ) <EOL> smry . add_table_2cols ( self , gleft = diagn_left , <EOL> gright = diagn_right , yname = yname , <EOL> xname = xna , title = "<STR_LIT>" ) <EOL> return smry <EOL> def get_margeff ( self , at = '<STR_LIT>' , method = '<STR_LIT>' , atexog = None , <EOL> dummy = False , count = False ) : <EOL> """<STR_LIT>""" <EOL> if self . model . constraint is not None : <EOL> warnings . warn ( "<STR_LIT>" ) <EOL> return GEEMargins ( self , ( at , method , atexog , dummy , count ) ) <EOL> def plot_isotropic_dependence ( self , ax = None , xpoints = <NUM_LIT:10> , <EOL> min_n = <NUM_LIT:50> ) : <EOL> """<STR_LIT>""" <EOL> from statsmodels . graphics import utils as gutils <EOL> resid = self . model . cluster_list ( self . resid ) <EOL> time = self . model . cluster_list ( self . model . time ) <EOL> xre , xdt = [ ] , [ ] <EOL> for re , ti in zip ( resid , time ) : <EOL> ix = np . tril_indices ( re . shape [ <NUM_LIT:0> ] , <NUM_LIT:0> ) <EOL> re = re [ ix [ <NUM_LIT:0> ] ] * re [ ix [ <NUM_LIT:1> ] ] / self . scale ** <NUM_LIT:2> <EOL> xre . append ( re ) <EOL> dists = np . sqrt ( ( ( ti [ ix [ <NUM_LIT:0> ] , : ] - ti [ ix [ <NUM_LIT:1> ] , : ] ) ** <NUM_LIT:2> ) . sum ( <NUM_LIT:1> ) ) <EOL> xdt . append ( dists ) <EOL> xre = np . concatenate ( xre ) <EOL> xdt = np . concatenate ( xdt ) <EOL> if ax is None : <EOL> fig , ax = gutils . create_mpl_ax ( ax ) <EOL> else : <EOL> fig = ax . get_figure ( ) <EOL> ii = np . flatnonzero ( xdt == <NUM_LIT:0> ) <EOL> v0 = np . mean ( xre [ ii ] ) <EOL> xre /= v0 <EOL> if np . isscalar ( xpoints ) : <EOL> xpoints = np . linspace ( <NUM_LIT:0> , max ( xdt ) , xpoints ) <EOL> dg = np . digitize ( xdt , xpoints ) <EOL> dgu = np . unique ( dg ) <EOL> hist = np . asarray ( [ np . sum ( dg == k ) for k in dgu ] ) <EOL> ii = np . flatnonzero ( hist >= min_n ) <EOL> dgu = dgu [ ii ] <EOL> dgy = np . asarray ( [ np . mean ( xre [ dg == k ] ) for k in dgu ] ) <EOL> dgx = np . asarray ( [ np . mean ( xdt [ dg == k ] ) for k in dgu ] ) <EOL> ax . plot ( dgx , dgy , '<STR_LIT:->' , color = '<STR_LIT>' , lw = <NUM_LIT:5> ) <EOL> ax . set_xlabel ( "<STR_LIT>" ) <EOL> ax . set_ylabel ( "<STR_LIT>" ) <EOL> return fig <EOL> def sensitivity_params ( self , dep_params_first , <EOL> dep_params_last , num_steps ) : <EOL> """<STR_LIT>""" <EOL> model = self . model <EOL> import copy <EOL> cov_struct = copy . deepcopy ( self . model . cov_struct ) <EOL> update_dep = model . update_dep <EOL> model . update_dep = False <EOL> dep_params = [ ] <EOL> results = [ ] <EOL> for x in np . linspace ( <NUM_LIT:0> , <NUM_LIT:1> , num_steps ) : <EOL> dp = x * dep_params_last + ( <NUM_LIT:1> - x ) * dep_params_first <EOL> dep_params . append ( dp ) <EOL> model . cov_struct = copy . deepcopy ( cov_struct ) <EOL> model . cov_struct . dep_params = dp <EOL> rslt = model . fit ( start_params = self . params , <EOL> ctol = self . ctol , <EOL> params_niter = self . params_niter , <EOL> first_dep_update = self . first_dep_update , <EOL> cov_type = self . cov_type ) <EOL> results . append ( rslt ) <EOL> model . update_dep = update_dep <EOL> return results <EOL> params_sensitivity = sensitivity_params <EOL> class GEEResultsWrapper ( lm . RegressionResultsWrapper ) : <EOL> _attrs = { <EOL> '<STR_LIT>' : '<STR_LIT>' , <EOL> } <EOL> _wrap_attrs = wrap . union_dicts ( lm . RegressionResultsWrapper . _wrap_attrs , <EOL> _attrs ) <EOL> wrap . populate_wrapper ( GEEResultsWrapper , GEEResults ) <EOL> class OrdinalGEE ( GEE ) : <EOL> __doc__ = ( <EOL> "<STR_LIT>" <EOL> "<STR_LIT>" + <EOL> _gee_init_doc % { '<STR_LIT>' : base . _missing_param_doc , <EOL> '<STR_LIT>' : _gee_ordinal_family_doc , <EOL> '<STR_LIT>' : _gee_ordinal_example } ) <EOL> def __init__ ( self , endog , exog , groups , time = None , family = None , <EOL> cov_struct = None , missing = '<STR_LIT:none>' , offset = None , <EOL> dep_data = None , constraint = None , ** kwargs ) : <EOL> if family is None : <EOL> family = families . Binomial ( ) <EOL> else : <EOL> if not isinstance ( family , families . Binomial ) : <EOL> raise ValueError ( "<STR_LIT>" ) <EOL> if cov_struct is None : <EOL> cov_struct = cov_structs . OrdinalIndependence ( ) <EOL> endog , exog , groups , time , offset = self . setup_ordinal ( endog , <EOL> exog , groups , time , offset ) <EOL> super ( OrdinalGEE , self ) . __init__ ( endog , exog , groups , time , <EOL> family , cov_struct , missing , <EOL> offset , dep_data , constraint ) <EOL> def setup_ordinal ( self , endog , exog , groups , time , offset ) : <EOL> """<STR_LIT>""" <EOL> self . endog_orig = endog . copy ( ) <EOL> self . exog_orig = exog . copy ( ) <EOL> self . groups_orig = groups . copy ( ) <EOL> if offset is not None : <EOL> self . offset_orig = offset . copy ( ) <EOL> else : <EOL> self . offset_orig = None <EOL> offset = np . zeros ( len ( endog ) ) <EOL> if time is not None : <EOL> self . time_orig = time . copy ( ) <EOL> else : <EOL> self . time_orig = None <EOL> time = np . zeros ( ( len ( endog ) , <NUM_LIT:1> ) ) <EOL> exog = np . asarray ( exog ) <EOL> endog = np . asarray ( endog ) <EOL> groups = np . asarray ( groups ) <EOL> time = np . asarray ( time ) <EOL> offset = np . asarray ( offset ) <EOL> self . endog_values = np . unique ( endog ) <EOL> endog_cuts = self . endog_values [ <NUM_LIT:0> : - <NUM_LIT:1> ] <EOL> ncut = len ( endog_cuts ) <EOL> nrows = ncut * len ( endog ) <EOL> exog_out = np . zeros ( ( nrows , exog . shape [ <NUM_LIT:1> ] ) , <EOL> dtype = np . float64 ) <EOL> endog_out = np . zeros ( nrows , dtype = np . float64 ) <EOL> intercepts = np . zeros ( ( nrows , ncut ) , dtype = np . float64 ) <EOL> groups_out = np . zeros ( nrows , dtype = groups . dtype ) <EOL> time_out = np . zeros ( ( nrows , time . shape [ <NUM_LIT:1> ] ) , <EOL> dtype = np . float64 ) <EOL> offset_out = np . zeros ( nrows , dtype = np . float64 ) <EOL> jrow = <NUM_LIT:0> <EOL> zipper = zip ( exog , endog , groups , time , offset ) <EOL> for ( exog_row , endog_value , group_value , time_value , <EOL> offset_value ) in zipper : <EOL> for thresh_ix , thresh in enumerate ( endog_cuts ) : <EOL> exog_out [ jrow , : ] = exog_row <EOL> endog_out [ jrow ] = ( int ( endog_value > thresh ) ) <EOL> intercepts [ jrow , thresh_ix ] = <NUM_LIT:1> <EOL> groups_out [ jrow ] = group_value <EOL> time_out [ jrow ] = time_value <EOL> offset_out [ jrow ] = offset_value <EOL> jrow += <NUM_LIT:1> <EOL> exog_out = np . concatenate ( ( intercepts , exog_out ) , axis = <NUM_LIT:1> ) <EOL> xnames = [ "<STR_LIT>" % v for v in endog_cuts ] <EOL> if type ( self . exog_orig ) == pd . DataFrame : <EOL> xnames . extend ( self . exog_orig . columns ) <EOL> else : <EOL> xnames . extend ( [ "<STR_LIT>" % k for k in range ( <NUM_LIT:1> , exog . shape [ <NUM_LIT:1> ] + <NUM_LIT:1> ) ] ) <EOL> exog_out = pd . DataFrame ( exog_out , columns = xnames ) <EOL> if type ( self . endog_orig ) == pd . Series : <EOL> endog_out = pd . Series ( endog_out , name = self . endog_orig . name ) <EOL> return endog_out , exog_out , groups_out , time_out , offset_out <EOL> def _starting_params ( self ) : <EOL> model = GEE ( self . endog , self . exog , self . groups , <EOL> time = self . time , family = families . Binomial ( ) , <EOL> offset = self . offset , exposure = self . exposure ) <EOL> result = model . fit ( ) <EOL> return result . params <EOL> def fit ( self , maxiter = <NUM_LIT> , ctol = <NUM_LIT> , start_params = None , <EOL> params_niter = <NUM_LIT:1> , first_dep_update = <NUM_LIT:0> , <EOL> cov_type = '<STR_LIT>' ) : <EOL> rslt = super ( OrdinalGEE , self ) . fit ( maxiter , ctol , start_params , <EOL> params_niter , first_dep_update , <EOL> cov_type = cov_type ) <EOL> rslt = rslt . _results <EOL> res_kwds = dict ( ( ( k , getattr ( rslt , k ) ) for k in rslt . _props ) ) <EOL> ord_rslt = OrdinalGEEResults ( self , rslt . params , <EOL> rslt . cov_params ( ) / rslt . scale , <EOL> rslt . scale , <EOL> cov_type = cov_type , <EOL> attr_kwds = res_kwds ) <EOL> return OrdinalGEEResultsWrapper ( ord_rslt ) <EOL> fit . __doc__ = _gee_fit_doc <EOL> class OrdinalGEEResults ( GEEResults ) : <EOL> __doc__ = ( <EOL> "<STR_LIT>" <EOL> "<STR_LIT>" <EOL> + _gee_results_doc ) <EOL> def plot_distribution ( self , ax = None , exog_values = None ) : <EOL> """<STR_LIT>""" <EOL> from statsmodels . graphics import utils as gutils <EOL> if ax is None : <EOL> fig , ax = gutils . create_mpl_ax ( ax ) <EOL> else : <EOL> fig = ax . get_figure ( ) <EOL> if exog_values is None : <EOL> exog_values = [ { } , ] <EOL> exog_means = self . model . exog . mean ( <NUM_LIT:0> ) <EOL> ix_icept = [ i for i , x in enumerate ( self . model . exog_names ) if <EOL> x . startswith ( "<STR_LIT>" ) ] <EOL> for ev in exog_values : <EOL> for k in ev . keys ( ) : <EOL> if k not in self . model . exog_names : <EOL> raise ValueError ( "<STR_LIT>" <EOL> % k ) <EOL> pr = [ ] <EOL> for j in ix_icept : <EOL> xp = np . zeros_like ( self . params ) <EOL> xp [ j ] = <NUM_LIT:1.> <EOL> for i , vn in enumerate ( self . model . exog_names ) : <EOL> if i in ix_icept : <EOL> continue <EOL> if vn in ev : <EOL> xp [ i ] = ev [ vn ] <EOL> else : <EOL> xp [ i ] = exog_means [ i ] <EOL> p = <NUM_LIT:1> / ( <NUM_LIT:1> + np . exp ( - np . dot ( xp , self . params ) ) ) <EOL> pr . append ( p ) <EOL> pr . insert ( <NUM_LIT:0> , <NUM_LIT:1> ) <EOL> pr . append ( <NUM_LIT:0> ) <EOL> pr = np . asarray ( pr ) <EOL> prd = - np . diff ( pr ) <EOL> ax . plot ( self . model . endog_values , prd , '<STR_LIT>' ) <EOL> ax . set_xlabel ( "<STR_LIT>" ) <EOL> ax . set_ylabel ( "<STR_LIT>" ) <EOL> ax . set_ylim ( <NUM_LIT:0> , <NUM_LIT:1> ) <EOL> return fig <EOL> class OrdinalGEEResultsWrapper ( GEEResultsWrapper ) : <EOL> pass <EOL> wrap . populate_wrapper ( OrdinalGEEResultsWrapper , OrdinalGEEResults ) <EOL> class NominalGEE ( GEE ) : <EOL> __doc__ = ( <EOL> "<STR_LIT>" <EOL> "<STR_LIT>" + <EOL> _gee_init_doc % { '<STR_LIT>' : base . _missing_param_doc , <EOL> '<STR_LIT>' : _gee_nominal_family_doc , <EOL> '<STR_LIT>' : _gee_nominal_example } ) <EOL> def __init__ ( self , endog , exog , groups , time = None , family = None , <EOL> cov_struct = None , missing = '<STR_LIT:none>' , offset = None , <EOL> dep_data = None , constraint = None , ** kwargs ) : <EOL> endog , exog , groups , time , offset = self . setup_nominal ( endog , <EOL> exog , groups , time , offset ) <EOL> if family is None : <EOL> family = _Multinomial ( self . ncut + <NUM_LIT:1> ) <EOL> if cov_struct is None : <EOL> cov_struct = cov_structs . NominalIndependence ( ) <EOL> super ( NominalGEE , self ) . __init__ ( endog , exog , groups , <EOL> time , family , cov_struct , missing , offset , dep_data , <EOL> constraint ) <EOL> def _starting_params ( self ) : <EOL> model = GEE ( self . endog , self . exog , self . groups , <EOL> time = self . time , family = families . Binomial ( ) , <EOL> offset = self . offset , exposure = self . exposure ) <EOL> result = model . fit ( ) <EOL> return result . params <EOL> def setup_nominal ( self , endog , exog , groups , time , offset ) : <EOL> """<STR_LIT>""" <EOL> self . endog_orig = endog . copy ( ) <EOL> self . exog_orig = exog . copy ( ) <EOL> self . groups_orig = groups . copy ( ) <EOL> if offset is not None : <EOL> self . offset_orig = offset . copy ( ) <EOL> else : <EOL> self . offset_orig = None <EOL> offset = np . zeros ( len ( endog ) ) <EOL> if time is not None : <EOL> self . time_orig = time . copy ( ) <EOL> else : <EOL> self . time_orig = None <EOL> time = np . zeros ( ( len ( endog ) , <NUM_LIT:1> ) ) <EOL> exog = np . asarray ( exog ) <EOL> endog = np . asarray ( endog ) <EOL> groups = np . asarray ( groups ) <EOL> time = np . asarray ( time ) <EOL> offset = np . asarray ( offset ) <EOL> self . endog_values = np . unique ( endog ) <EOL> endog_cuts = self . endog_values [ <NUM_LIT:0> : - <NUM_LIT:1> ] <EOL> ncut = len ( endog_cuts ) <EOL> self . ncut = ncut <EOL> nrows = len ( endog_cuts ) * exog . shape [ <NUM_LIT:0> ] <EOL> ncols = len ( endog_cuts ) * exog . shape [ <NUM_LIT:1> ] <EOL> exog_out = np . zeros ( ( nrows , ncols ) , dtype = np . float64 ) <EOL> endog_out = np . zeros ( nrows , dtype = np . float64 ) <EOL> groups_out = np . zeros ( nrows , dtype = np . float64 ) <EOL> time_out = np . zeros ( ( nrows , time . shape [ <NUM_LIT:1> ] ) , <EOL> dtype = np . float64 ) <EOL> offset_out = np . zeros ( nrows , dtype = np . float64 ) <EOL> jrow = <NUM_LIT:0> <EOL> zipper = zip ( exog , endog , groups , time , offset ) <EOL> for ( exog_row , endog_value , group_value , time_value , <EOL> offset_value ) in zipper : <EOL> for thresh_ix , thresh in enumerate ( endog_cuts ) : <EOL> u = np . zeros ( len ( endog_cuts ) , dtype = np . float64 ) <EOL> u [ thresh_ix ] = <NUM_LIT:1> <EOL> exog_out [ jrow , : ] = np . kron ( u , exog_row ) <EOL> endog_out [ jrow ] = ( int ( endog_value == thresh ) ) <EOL> groups_out [ jrow ] = group_value <EOL> time_out [ jrow ] = time_value <EOL> offset_out [ jrow ] = offset_value <EOL> jrow += <NUM_LIT:1> <EOL> if type ( self . exog_orig ) == pd . DataFrame : <EOL> xnames_in = self . exog_orig . columns <EOL> else : <EOL> xnames_in = [ "<STR_LIT>" % k for k in range ( <NUM_LIT:1> , exog . shape [ <NUM_LIT:1> ] + <NUM_LIT:1> ) ] <EOL> xnames = [ ] <EOL> for tr in endog_cuts : <EOL> xnames . extend ( [ "<STR_LIT>" % ( v , tr ) for v in xnames_in ] ) <EOL> exog_out = pd . DataFrame ( exog_out , columns = xnames ) <EOL> exog_out = pd . DataFrame ( exog_out , columns = xnames ) <EOL> if type ( self . endog_orig ) == pd . Series : <EOL> endog_out = pd . Series ( endog_out , name = self . endog_orig . name ) <EOL> return endog_out , exog_out , groups_out , time_out , offset_out <EOL> def mean_deriv ( self , exog , lin_pred ) : <EOL> """<STR_LIT>""" <EOL> expval = np . exp ( lin_pred ) <EOL> expval_m = np . reshape ( expval , ( len ( expval ) // self . ncut , <EOL> self . ncut ) ) <EOL> denom = <NUM_LIT:1> + expval_m . sum ( <NUM_LIT:1> ) <EOL> denom = np . kron ( denom , np . ones ( self . ncut , dtype = np . float64 ) ) <EOL> mprob = expval / denom <EOL> dmat = mprob [ : , None ] * exog <EOL> ddenom = expval [ : , None ] * exog <EOL> dmat -= mprob [ : , None ] * ddenom / denom [ : , None ] <EOL> return dmat <EOL> def mean_deriv_exog ( self , exog , params , offset_exposure = None ) : <EOL> """<STR_LIT>""" <EOL> if offset_exposure is not None : <EOL> warnings . warn ( "<STR_LIT>" ) <EOL> lpr = np . dot ( exog , params ) <EOL> expval = np . exp ( lpr ) <EOL> expval_m = np . reshape ( expval , ( len ( expval ) // self . ncut , <EOL> self . ncut ) ) <EOL> denom = <NUM_LIT:1> + expval_m . sum ( <NUM_LIT:1> ) <EOL> denom = np . kron ( denom , np . ones ( self . ncut , dtype = np . float64 ) ) <EOL> bmat0 = np . outer ( np . ones ( exog . shape [ <NUM_LIT:0> ] ) , params ) <EOL> qmat = [ ] <EOL> for j in range ( self . ncut ) : <EOL> ee = np . zeros ( self . ncut , dtype = np . float64 ) <EOL> ee [ j ] = <NUM_LIT:1> <EOL> qmat . append ( np . kron ( ee , np . ones ( len ( params ) // self . ncut ) ) ) <EOL> qmat = np . array ( qmat ) <EOL> qmat = np . kron ( np . ones ( ( exog . shape [ <NUM_LIT:0> ] // self . ncut , <NUM_LIT:1> ) ) , qmat ) <EOL> bmat = bmat0 * qmat <EOL> dmat = expval [ : , None ] * bmat / denom [ : , None ] <EOL> expval_mb = np . kron ( expval_m , np . ones ( ( self . ncut , <NUM_LIT:1> ) ) ) <EOL> expval_mb = np . kron ( expval_mb , np . ones ( ( <NUM_LIT:1> , self . ncut ) ) ) <EOL> dmat -= expval [ : , None ] * ( bmat * expval_mb ) / denom [ : , None ] ** <NUM_LIT:2> <EOL> return dmat <EOL> def fit ( self , maxiter = <NUM_LIT> , ctol = <NUM_LIT> , start_params = None , <EOL> params_niter = <NUM_LIT:1> , first_dep_update = <NUM_LIT:0> , <EOL> cov_type = '<STR_LIT>' ) : <EOL> rslt = super ( NominalGEE , self ) . fit ( maxiter , ctol , start_params , <EOL> params_niter , first_dep_update , <EOL> cov_type = cov_type ) <EOL> if rslt is None : <EOL> warnings . warn ( "<STR_LIT>" , <EOL> ConvergenceWarning ) <EOL> return None <EOL> rslt = rslt . _results <EOL> res_kwds = dict ( ( ( k , getattr ( rslt , k ) ) for k in rslt . _props ) ) <EOL> nom_rslt = NominalGEEResults ( self , rslt . params , <EOL> rslt . cov_params ( ) / rslt . scale , <EOL> rslt . scale , <EOL> cov_type = cov_type , <EOL> attr_kwds = res_kwds ) <EOL> return NominalGEEResultsWrapper ( nom_rslt ) <EOL> fit . __doc__ = _gee_fit_doc <EOL> class NominalGEEResults ( GEEResults ) : <EOL> __doc__ = ( <EOL> "<STR_LIT>" <EOL> "<STR_LIT>" <EOL> + _gee_results_doc ) <EOL> def plot_distribution ( self , ax = None , exog_values = None ) : <EOL> """<STR_LIT>""" <EOL> from statsmodels . graphics import utils as gutils <EOL> if ax is None : <EOL> fig , ax = gutils . create_mpl_ax ( ax ) <EOL> else : <EOL> fig = ax . get_figure ( ) <EOL> if exog_values is None : <EOL> exog_values = [ { } , ] <EOL> link = self . model . family . link . inverse <EOL> ncut = self . model . family . ncut <EOL> k = self . model . exog . shape [ <NUM_LIT:1> ] / ncut <EOL> exog_means = self . model . exog . mean ( <NUM_LIT:0> ) [ <NUM_LIT:0> : k ] <EOL> exog_names = self . model . exog_names [ <NUM_LIT:0> : k ] <EOL> exog_names = [ x . split ( "<STR_LIT:[>" ) [ <NUM_LIT:0> ] for x in exog_names ] <EOL> params = np . reshape ( self . params , <EOL> ( ncut , len ( self . params ) / ncut ) ) <EOL> for ev in exog_values : <EOL> exog = exog_means . copy ( ) <EOL> for k in ev . keys ( ) : <EOL> if k not in exog_names : <EOL> raise ValueError ( "<STR_LIT>" <EOL> % k ) <EOL> ii = exog_names . index ( k ) <EOL> exog [ ii ] = ev [ k ] <EOL> lpr = np . dot ( params , exog ) <EOL> pr = link ( lpr ) <EOL> pr = np . r_ [ pr , <NUM_LIT:1> - pr . sum ( ) ] <EOL> ax . plot ( self . model . endog_values , pr , '<STR_LIT>' ) <EOL> ax . set_xlabel ( "<STR_LIT>" ) <EOL> ax . set_ylabel ( "<STR_LIT>" ) <EOL> ax . set_xticks ( self . model . endog_values ) <EOL> ax . set_xticklabels ( self . model . endog_values ) <EOL> ax . set_ylim ( <NUM_LIT:0> , <NUM_LIT:1> ) <EOL> return fig <EOL> class NominalGEEResultsWrapper ( GEEResultsWrapper ) : <EOL> pass <EOL> wrap . populate_wrapper ( NominalGEEResultsWrapper , NominalGEEResults ) <EOL> class _MultinomialLogit ( Link ) : <EOL> """<STR_LIT>""" <EOL> def __init__ ( self , ncut ) : <EOL> self . ncut = ncut <EOL> def inverse ( self , lpr ) : <EOL> """<STR_LIT>""" <EOL> expval = np . exp ( lpr ) <EOL> denom = <NUM_LIT:1> + np . reshape ( expval , ( len ( expval ) // self . ncut , <EOL> self . ncut ) ) . sum ( <NUM_LIT:1> ) <EOL> denom = np . kron ( denom , np . ones ( self . ncut , dtype = np . float64 ) ) <EOL> prob = expval / denom <EOL> return prob <EOL> class _Multinomial ( families . Family ) : <EOL> """<STR_LIT>""" <EOL> links = [ _MultinomialLogit , ] <EOL> variance = varfuncs . binary <EOL> safe_links = [ _MultinomialLogit , ] <EOL> def __init__ ( self , nlevels ) : <EOL> """<STR_LIT>""" <EOL> self . initialize ( nlevels ) <EOL> def initialize ( self , nlevels ) : <EOL> self . ncut = nlevels - <NUM_LIT:1> <EOL> self . link = _MultinomialLogit ( self . ncut ) <EOL> from statsmodels . discrete . discrete_margins import ( _get_margeff_exog , <EOL> _get_const_index , _check_margeff_args , _effects_at , <EOL> margeff_cov_with_se , _check_at_is_all , _transform_names , <EOL> _check_discrete_args , _get_dummy_index , _get_count_index ) <EOL> class GEEMargins ( object ) : <EOL> """<STR_LIT>""" <EOL> def __init__ ( self , results , args , kwargs = { } ) : <EOL> self . _cache = resettable_cache ( ) <EOL> self . results = results <EOL> self . get_margeff ( * args , ** kwargs ) <EOL> def _reset ( self ) : <EOL> self . _cache = resettable_cache ( ) <EOL> @ cache_readonly <EOL> def tvalues ( self ) : <EOL> _check_at_is_all ( self . margeff_options ) <EOL> return self . margeff / self . margeff_se <EOL> def summary_frame ( self , alpha = <NUM_LIT> ) : <EOL> """<STR_LIT>""" <EOL> _check_at_is_all ( self . margeff_options ) <EOL> from pandas import DataFrame <EOL> names = [ _transform_names [ self . margeff_options [ '<STR_LIT>' ] ] , <EOL> '<STR_LIT>' , '<STR_LIT:z>' , '<STR_LIT>' , <EOL> '<STR_LIT>' , '<STR_LIT>' ] <EOL> ind = self . results . model . exog . var ( <NUM_LIT:0> ) != <NUM_LIT:0> <EOL> exog_names = self . results . model . exog_names <EOL> var_names = [ name for i , name in enumerate ( exog_names ) if ind [ i ] ] <EOL> table = np . column_stack ( ( self . margeff , self . margeff_se , self . tvalues , <EOL> self . pvalues , self . conf_int ( alpha ) ) ) <EOL> return DataFrame ( table , columns = names , index = var_names ) <EOL> @ cache_readonly <EOL> def pvalues ( self ) : <EOL> _check_at_is_all ( self . margeff_options ) <EOL> return stats . norm . sf ( np . abs ( self . tvalues ) ) * <NUM_LIT:2> <EOL> def conf_int ( self , alpha = <NUM_LIT> ) : <EOL> """<STR_LIT>""" <EOL> _check_at_is_all ( self . margeff_options ) <EOL> me_se = self . margeff_se <EOL> q = stats . norm . ppf ( <NUM_LIT:1> - alpha / <NUM_LIT:2> ) <EOL> lower = self . margeff - q * me_se <EOL> upper = self . margeff + q * me_se <EOL> return np . asarray ( lzip ( lower , upper ) ) <EOL> def summary ( self , alpha = <NUM_LIT> ) : <EOL> """<STR_LIT>""" <EOL> _check_at_is_all ( self . margeff_options ) <EOL> results = self . results <EOL> model = results . model <EOL> title = model . __class__ . __name__ + "<STR_LIT>" <EOL> method = self . margeff_options [ '<STR_LIT>' ] <EOL> top_left = [ ( '<STR_LIT>' , [ model . endog_names ] ) , <EOL> ( '<STR_LIT>' , [ method ] ) , <EOL> ( '<STR_LIT>' , [ self . margeff_options [ '<STR_LIT>' ] ] ) , ] <EOL> from statsmodels . iolib . summary import ( Summary , summary_params , <EOL> table_extend ) <EOL> exog_names = model . exog_names [ : ] <EOL> smry = Summary ( ) <EOL> _ , const_idx = _get_const_index ( model . exog ) <EOL> if const_idx is not None : <EOL> exog_names . pop ( const_idx ) <EOL> J = int ( getattr ( model , "<STR_LIT>" , <NUM_LIT:1> ) ) <EOL> if J > <NUM_LIT:1> : <EOL> yname , yname_list = results . _get_endog_name ( model . endog_names , <EOL> None , all = True ) <EOL> else : <EOL> yname = model . endog_names <EOL> yname_list = [ yname ] <EOL> smry . add_table_2cols ( self , gleft = top_left , gright = [ ] , <EOL> yname = yname , xname = exog_names , title = title ) <EOL> table = [ ] <EOL> conf_int = self . conf_int ( alpha ) <EOL> margeff = self . margeff <EOL> margeff_se = self . margeff_se <EOL> tvalues = self . tvalues <EOL> pvalues = self . pvalues <EOL> if J > <NUM_LIT:1> : <EOL> for eq in range ( J ) : <EOL> restup = ( results , margeff [ : , eq ] , margeff_se [ : , eq ] , <EOL> tvalues [ : , eq ] , pvalues [ : , eq ] , conf_int [ : , : , eq ] ) <EOL> tble = summary_params ( restup , yname = yname_list [ eq ] , <EOL> xname = exog_names , alpha = alpha , use_t = False , <EOL> skip_header = True ) <EOL> tble . title = yname_list [ eq ] <EOL> header = [ '<STR_LIT>' , _transform_names [ method ] , '<STR_LIT>' , '<STR_LIT:z>' , <EOL> '<STR_LIT>' , '<STR_LIT>' % ( <NUM_LIT:100> - alpha * <NUM_LIT:100> ) ] <EOL> tble . insert_header_row ( <NUM_LIT:0> , header ) <EOL> table . append ( tble ) <EOL> table = table_extend ( table , keep_headers = True ) <EOL> else : <EOL> restup = ( results , margeff , margeff_se , tvalues , pvalues , conf_int ) <EOL> table = summary_params ( restup , yname = yname , xname = exog_names , <EOL> alpha = alpha , use_t = False , skip_header = True ) <EOL> header = [ '<STR_LIT>' , _transform_names [ method ] , '<STR_LIT>' , '<STR_LIT:z>' , <EOL> '<STR_LIT>' , '<STR_LIT>' % ( <NUM_LIT:100> - alpha * <NUM_LIT:100> ) ] <EOL> table . insert_header_row ( <NUM_LIT:0> , header ) <EOL> smry . tables . append ( table ) <EOL> return smry <EOL> def get_margeff ( self , at = '<STR_LIT>' , method = '<STR_LIT>' , atexog = None , <EOL> dummy = False , count = False ) : <EOL> self . _reset ( ) <EOL> method = method . lower ( ) <EOL> at = at . lower ( ) <EOL> _check_margeff_args ( at , method ) <EOL> self . margeff_options = dict ( method = method , at = at ) <EOL> results = self . results <EOL> model = results . model <EOL> params = results . params <EOL> exog = model . exog . copy ( ) <EOL> effects_idx , const_idx = _get_const_index ( exog ) <EOL> if dummy : <EOL> _check_discrete_args ( at , method ) <EOL> dummy_idx , dummy = _get_dummy_index ( exog , const_idx ) <EOL> else : <EOL> dummy_idx = None <EOL> if count : <EOL> _check_discrete_args ( at , method ) <EOL> count_idx , count = _get_count_index ( exog , const_idx ) <EOL> else : <EOL> count_idx = None <EOL> exog = _get_margeff_exog ( exog , at , atexog , effects_idx ) <EOL> effects = model . _derivative_exog ( params , exog , method , <EOL> dummy_idx , count_idx ) <EOL> effects = _effects_at ( effects , at ) <EOL> if at == '<STR_LIT:all>' : <EOL> self . margeff = effects [ : , effects_idx ] <EOL> else : <EOL> margeff_cov , margeff_se = margeff_cov_with_se ( model , params , exog , <EOL> results . cov_params ( ) , at , <EOL> model . _derivative_exog , <EOL> dummy_idx , count_idx , <EOL> method , <NUM_LIT:1> ) <EOL> self . margeff_cov = margeff_cov [ effects_idx ] [ : , effects_idx ] <EOL> self . margeff_se = margeff_se [ effects_idx ] <EOL> self . margeff = effects [ effects_idx ] </s>
