""" Test Results for discrete models from Stata """ import os import numpy as np class Namespace: pass # Discrete Model Tests # Note that there is a slight refactor of the classes, so that one dataset # might be used for more than one model cur_dir = os.path.abspath(os.path.dirname(__file__)) class Anes: def __init__(self): """r Results are from Stata 11 (checked vs R nnet package). """ self.nobs = 944 def mnlogit_basezero(): obj = Namespace() obj.nobs = 944 params = [ -.01153598, .29771435, -.024945, .08249144, .00519655, -.37340167, -.08875065, .39166864, -.02289784, .18104276, .04787398, -2.2509132, -.1059667, .57345051, -.01485121, -.00715242, .05757516, -3.6655835, -.0915567, 1.2787718, -.00868135, .19982796, .08449838, -7.6138431, -.0932846, 1.3469616, -.01790407, .21693885, .08095841, -7.0604782, -.14088069, 2.0700801, -.00943265, .3219257, .10889408, -12.105751] obj.params = np.reshape(params, (6, -1), order='F') bse = [ .0342823657, .093626795, .0065248584, .0735865799, .0176336937, .6298376313, .0391615553, .1082386919, .0079144618, .0852893563, .0222809297, .7631899491, .0570382292, .1585481337, .0113313133, .1262913234, .0336142088, 1.156541492, .0437902764, .1288965854, .0084187486, .0941250559, .0261963632, .9575809602, .0393516553, .1171860107, .0076110152, .0850070091, .0229760791, .8443638283, .042138047, .1434089089, .0081338625, .0910979921, .025300888, 1.059954821] obj.bse = np.reshape(bse, (6, -1), order='F') obj.yhat = np.loadtxt(os.path.join(cur_dir, 'yhat_mnlogit.csv')) obj.phat = np.loadtxt(os.path.join(cur_dir, 'phat_mnlogit.csv')) obj.cov_params = None obj.llf = -1461.922747312 obj.llnull = -1750.34670999 obj.llr = 576.8479253554 obj.llr_pvalue = 1.8223179e-102 obj.prsquared = .1647810465387 obj.df_model = 30 obj.df_resid = 944 - 36 obj.J = 7 obj.K = 6 obj.aic = 2995.84549462 obj.bic = 3170.45003661 z = [ -.3364988051, 3.179798597, -3.823070772, 1.121012042, .2946945327, -.5928538661, -2.266269864, 3.618564069, -2.893164162, 2.122688754, 2.148652536, -2.949348555, -1.857818873, 3.616885888, -1.310634214, -.0566342868, 1.712822091, -3.169435381, -2.090799808, 9.920912816, -1.031191864, 2.123004903, 3.225576554, -7.951122047, -2.370538224, 11.49421878, -2.352389066, 2.552011323, 3.523595639, -8.361890935, -3.34331327, 14.43480847, -1.159676452, 3.533839715, 4.303962885, -11.42100649] obj.z = np.reshape(z, (6, -1), order='F') pvalues = [ 0.7364947525, 0.0014737744, 0.0001317999, 0.2622827367, 0.7682272401, 0.5532789548, 0.0234348654, 0.0002962422, 0.0038138191, 0.0337799420, 0.0316619538, 0.0031844460, 0.0631947400, 0.0002981687, 0.1899813744, 0.9548365214, 0.0867452747, 0.0015273542, 0.0365460134, 3.37654e-23, 0.3024508550, 0.0337534410, 0.0012571921, 1.84830e-15, 0.0177622072, 1.41051e-30, 0.0186532528, 0.0107103038, 0.0004257334, 6.17209e-17, 0.0008278439, 3.12513e-47, 0.2461805610, 0.0004095694, 0.0000167770, 3.28408e-30] obj.pvalues = np.reshape(pvalues, (6, -1), order='F') conf_int = [ [[-0.0787282, 0.0556562], [0.1142092, 0.4812195], [-0.0377335, -0.0121565], [-0.0617356, 0.2267185], [-0.0293649, 0.0397580], [-1.6078610, 0.8610574]], [[-0.1655059, -0.0119954], [0.1795247, 0.6038126], [-0.0384099, -0.0073858], [0.0138787, 0.3482068], [0.0042042, 0.0915438], [-3.7467380, -0.7550884]], [[-0.2177596, 0.0058262], [0.2627019, 0.8841991], [-0.0370602, 0.0073578], [-0.2546789, 0.2403740], [-0.0083075, 0.1234578], [-5.9323630, -1.3988040]], [[-0.1773841, -0.0057293], [1.0261390, 1.5314040], [-0.0251818, 0.0078191], [0.0153462, 0.3843097], [0.0331544, 0.1358423], [-9.4906670, -5.7370190]], [[-0.1704124, -0.0161568], [1.1172810, 1.5766420], [-0.0328214, -0.0029868], [0.0503282, 0.3835495], [0.0359261, 0.1259907], [-8.7154010, -5.4055560]], [[-0.2234697, -0.0582916], [1.7890040, 2.3511560], [-0.0253747, 0.0065094], [0.1433769, 0.5004745], [0.0593053, 0.1584829], [-14.1832200, -10.0282800]]] obj.conf_int = np.asarray(conf_int) # margins, dydx(*) predict(outcome(#)) obj.margeff_dydx_overall = np.array([ [0.00868085993550, -0.09779854015456, 0.00272556969847, -0.01992376579372, -0.00603133322764], [0.00699386733148, -0.05022430802614, -0.00211003909752, -0.00536980000265, -0.00554366741814], [-0.00391040848820, -0.02824717135857, -0.00100551299310, 0.00664337806861, 0.00097987356999], [-0.00182580888015, -0.00573744730031, -0.00004249256428, -0.00546669558488, 0.00054101121854], [-0.00098558129923, 0.01985550937033, 0.00047972250012, 0.00172605778905, 0.00211291403209], [-0.00153469551647, 0.03755346502013, -0.00068531143399, 0.00472471794347, 0.00254733486106], [-0.00741820702809, 0.12459834487569, 0.00063806819375, 0.01766610701188, 0.00539385283759] ]).T obj.margeff_dydx_overall_se = np.array([ [.0038581061, .0080471125, .0007068488, .0082318967, .0020261706], [.003904378, .0073600286, .000756431, .0084381578, .0020482238], [.003137126, .0056813182, .0006601377, .0068932588, .0018481806], [.0019427783, .0031904763, .0003865411, .004361789, .0011523221], [.0029863227, .0054076092, .0005886612, .0064426365, .0018886818], [.0035806552, .0069497362, .000722511, .0078287717, .0022352393], [.0033641608, .008376629, .0006774697, .0073505286, .0021660086] ]).T obj.margeff_dydx_mean = np.array([ [0.01149887431225, -0.13784207091973, 0.00273313385873, -0.02542974260540, -0.00855346837482], [0.01114846831102, -0.09864273512889, -0.00222435063712, -0.01214617126321, -0.00903581444579], [-0.00381702868421, -0.05132297961269, -0.00116763216994, 0.00624203027060, 0.00021912081810], [-0.00233455327258, -0.00928554037343, -0.00000206561214, -0.00775415690571, 0.00060004460394], [-0.00352579921274, 0.06412187169362, 0.00073938948643, 0.00747778063206, 0.00459965010365], [-0.00574308219449, 0.11126535089794, -0.00057337915464, 0.01467424346725, 0.00641760846097], [-0.00722687818452, 0.12170608820238, 0.00049490419675, 0.01693601418978, 0.00575285798725]]).T obj.margeff_dydx_mean_se = np.array([ [.0043729758, .0110343353, .0008149907, .0092551389, .0023752071], [.004875051, .0124746358, .0009613152, .0105665812, .0026524426], [.0040718954, .0103613938, .0008554615, .0089931297, .0024374625], [.0026430804, .0070845916, .0005364369, .0057654258, .0015988838], [.0037798151, .0103849291, .0007393481, .0082021938, .0023489261], [.0045654631, .0130329403, .0009128134, .0100053262, .0028048602], [.0027682389, .0113292677, .0005325113, .0061289353, .0017330763] ]).T obj.margeff_dydx_dummy_overall = np.array([ [0.00549149574321, -0.05348235321783, 0.00298963549049, -0.01479461677951, -0.00332167981255, -0.26502967041815], [0.00345677928276, -0.00950322030929, -0.00189456107189, 0.00033893662061, -0.00314690167350, -0.21040878091828], [-0.00645089013284, 0.00401746940204, -0.00083948249351, 0.01114202556889, 0.00277069841472, -0.15967397659686], [-0.00215436802341, -0.00366545199370, -0.00000002297812, -0.00457368049644, 0.00065303026027, -0.00094772782001], [0.00058038428936, -0.00369080100124, 0.00035948233235, -0.00018863693013, 0.00079351293461, 0.12640653743480], [0.00217597030999, -0.01279456622853, -0.00091882392767, 0.00001651192759, -0.00037998290789, 0.27175070356670], [-0.00309932483642, 0.07911868907484, 0.00030378521102, 0.00805941631677, 0.00263129901425, 0.23790291475181]]).T obj.margeff_dydx_dummy_overall_se = np.array([ [.0037314453, .0094102332, .000688838, .0079744554, .0019365971, .0243914836], [.0038215262, .0095938828, .0007410885, .008259353, .0019984087, .0317628806], [.0031045718, .00785814, .0006504353, .0067892866, .0018060332, 0.0262803561], [.0019756086, .0051031194, .0003862449, .0043621673, .0011796953, .0219999601], [.0029714074, .0081732018, .0005715192, .0064742872, .0019130195, .0331694192], [.0034443743, .0097296187, .0006774867, .0075996454, .0021993881, .038600835], [.0032003518, .0098741227, .0006335772, .0070902078, .0021003227, .0255727127]]).T obj.margeff_eydx_dummy_overall = np.array([ [.03939188, -.65758371, .01750922, -.12131806, -.03613241, -3.2132513], [.02752366, -.383165, -.00830021, -.03652935, -.03286046, -1.8741853], [-.05006681, -.2719659, -.00626481, .06525323, .01012554, -2.0058029], [-.05239558, -.22549142, .00025015, -.13104416, .01114517, -.27052009], [-.00296374, .25627809, .00140513, .03358712, .02296041, 1.3302701], [.00328283, .2800168, -.0083912, .04332782, .01575863, 1.8441023], [-.03257068, .98346111, -.00122118, .10847807, .0406456, 2.9119099]]).T obj.margeff_eydx_dummy_overall_se = np.array([ [.0272085605, .0777760394, .0052427952, .0584011446, .0148618012, .5796921383], [.0262290023, .0724479385, .005174736, .0567743614, .0144447083, .3015738731], [.0321415498, .0895589422, .0067480662, .0701460193, .0190451865, .3904138447], [.0511305319, .1420904068, .0102342163, .1129912244, .0308618233, .3693799595], [.0340186217, .0991711703, .0065812158, .0737441012, .0212966336, .2346982385], [.0289250212, .0840662279, .0056743561, .0631772185, .0177278895, .2089516714], [.0318251305, .1085637405, .0062400589, .0699123044, .0201045606, .3727166284]]).T # taken from gretl obj.resid = np.loadtxt(os.path.join(cur_dir, 'mnlogit_resid.csv'), delimiter=",") return obj mnlogit_basezero = mnlogit_basezero() class DiscreteL1: def __init__(self): """ Special results for L1 models Uses the Spector data and a script to generate the baseline results """ pass def logit(): """ Results generated with: data = sm.datasets.spector.load() data.exog = sm.add_constant(data.exog, prepend=True) alpha = 3 * np.array([0, 1, 1, 1]) res2 = sm.Logit(data.endog, data.exog).fit_regularized( method="l1", alpha=alpha, disp=0, trim_mode='size', size_trim_tol=1e-5, acc=1e-10, maxiter=1000) """ obj = Namespace() nan = np.nan obj.params = [-4.10271595, 0., 0.15493781, 0.] obj.conf_int = [ [-9.15205122, 0.94661932], [nan, nan], [-0.06539482, 0.37527044], [nan, nan]] obj.bse = [2.5762388, nan, 0.11241668, nan] obj.nnz_params = 2 obj.aic = 42.091439368583671 obj.bic = 45.022911174183122 obj.cov_params = [ [6.63700638, nan, -0.28636261, nan], [nan, nan, nan, nan], [-0.28636261, nan, 0.01263751, nan], [nan, nan, nan, nan]] return obj logit = logit() def sweep(): """ Results generated with params = np.zeros((3, 4)) alphas = np.array( [[0.1, 0.1, 0.1, 0.1], [0.4, 0.4, 0.5, 0.5], [0.5, 0.5, 1, 1]]) model = sm.Logit(data.endog, data.exog) for i in range(3): alpha = alphas[i, :] res2 = model.fit_regularized(method="l1", alpha=alpha, disp=0, acc=1e-10, maxiter=1000, trim_mode='off') params[i, :] = res2.params print(params) """ obj = Namespace() obj.params = [ [-10.37593611, 2.27080968, 0.06670638, 2.05723691], [-5.32670811, 1.18216019, 0.01402395, 1.45178712], [-3.92630318, 0.90126958, -0., 1.09498178]] return obj sweep = sweep() def probit(): """ Results generated with data = sm.datasets.spector.load() data.exog = sm.add_constant(data.exog, prepend=True) alpha = np.array([0.1, 0.2, 0.3, 10]) res2 = sm.Probit(data.endog, data.exog).fit_regularized( method="l1", alpha=alpha, disp=0, trim_mode='auto', auto_trim_tol=0.02, acc=1e-10, maxiter=1000) """ obj = Namespace() nan = np.nan obj.params = [-5.40476992, 1.25018458, 0.04744558, 0.] obj.conf_int = [ [-9.44077951, -1.36876033], [0.03716721, 2.46320194], [-0.09727571, 0.19216687], [np.nan, np.nan]] obj.bse = [2.05922641, 0.61889778, 0.07383875, np.nan] obj.nnz_params = 3 obj.aic = 38.399773877542927 obj.bic = 42.796981585942106 obj.cov_params = [ [4.24041339, -0.83432592, -0.06827915, nan], [-0.83432592, 0.38303447, -0.01700249, nan], [-0.06827915, -0.01700249, 0.00545216, nan], [nan, nan, nan, nan]] return obj probit = probit() def mnlogit(): """ Results generated with anes_data = sm.datasets.anes96.load() anes_exog = anes_data.exog anes_exog = sm.add_constant(anes_exog, prepend=False) mlogit_mod = sm.MNLogit(anes_data.endog, anes_exog) alpha = 10 * np.ones((mlogit_mod.J - 1, mlogit_mod.K)) alpha[-1, :] = 0 mlogit_l1_res = mlogit_mod.fit_regularized( method='l1', alpha=alpha, trim_mode='auto', auto_trim_tol=0.02, acc=1e-10) """ obj = Namespace() obj.params = [ [0.00100163, -0.05864195, -0.06147822, -0.04769671, -0.05222987, -0.09522432], [0., 0.03186139, 0.12048999, 0.83211915, 0.92330292, 1.5680646], [-0.0218185, -0.01988066, -0.00808564, -0.00487463, -0.01400173, -0.00562079], [0., 0.03306875, 0., 0.02362861, 0.05486435, 0.14656966], [0., 0.04448213, 0.03252651, 0.07661761, 0.07265266, 0.0967758], [0.90993803, -0.50081247, -2.08285102, -5.26132955, -4.86783179, -9.31537963]] obj.conf_int = [ [[-0.0646223, 0.06662556], [np.nan, np.nan], [-0.03405931, -0.00957768], [np.nan, np.nan], [np.nan, np.nan], [0.26697895, 1.55289711]], [[-0.1337913, 0.01650741], [-0.14477255, 0.20849532], [-0.03500303, -0.00475829], [-0.11406121, 0.18019871], [0.00479741, 0.08416684], [-1.84626136, 0.84463642]], [[-0.17237962, 0.04942317], [-0.15146029, 0.39244026], [-0.02947379, 0.01330252], [np.nan, np.nan], [-0.02501483, 0.09006785], [-3.90379391, -0.26190812]], [[-0.12938296, 0.03398954], [0.62612955, 1.03810876], [-0.02046322, 0.01071395], [-0.13738534, 0.18464256], [0.03017236, 0.12306286], [-6.91227465, -3.61038444]], [[-0.12469773, 0.02023799], [0.742564, 1.10404183], [-0.02791975, -0.00008371], [-0.08491561, 0.19464431], [0.0332926, 0.11201273], [-6.29331126, -3.44235233]], [[-0.17165567, -0.01879296], [1.33994079, 1.79618841], [-0.02027503, 0.00903345], [-0.00267819, 0.29581751], [0.05343135, 0.14012026], [-11.10419107, -7.52656819]]] obj.bse = [ [0.03348221, 0.03834221, 0.05658338, 0.04167742, 0.03697408, 0.03899631], [np.nan, 0.09012101, 0.13875269, 0.10509867, 0.09221543, 0.11639184], [0.00624543, 0.00771564, 0.01091253, 0.00795351, 0.00710116, 0.00747679], [np.nan, 0.07506769, np.nan, 0.08215148, 0.07131762, 0.07614826], [np.nan, 0.02024768, 0.02935837, 0.02369699, 0.02008204, 0.02211492], [0.32804638, 0.68646613, 0.92906957, 0.84233441, 0.72729881, 0.91267567]] obj.nnz_params = 32 obj.aic = 3019.4391360294126 obj.bic = 3174.6431733460686 return obj mnlogit = mnlogit() class Spector: """ Results are from Stata 11 """ def __init__(self): self.nobs = 32 def logit(): obj = Namespace() obj.nobs = 32 obj.params = [2.82611297201, .0951576702557, 2.37868772835, -13.0213483201] obj.cov_params = [ [1.59502033639, -.036920566629, .427615725153, -4.57347950298], [-.036920566629, .0200375937069, .0149126464275, -.346255757562], [.427615725153, .0149126464275, 1.13329715236, -2.35916128427], [-4.57347950298, -.346255757562, -2.35916128427, 24.3179625937]] obj.bse = [1.26294114526, .141554207662, 1.06456430165, 4.93132462871] obj.resid_pearson = [ -.1652382, -.2515266, -.4800059, -.1630655, .8687437, -.1900454, -.165002, -.2331563, -.3535812, .6647838, -.1583799, -.4843181, -.689527, 2.043449, -.7516119, -.1764176, -.2380445, -.2003426, -1.199277, .7164842, -.255713, .3242821, -.5646816, -2.400189, .4392082, 1.038473, .75747, -.6659256, .4336657, .2404583, -1.060033, 2.829577] obj.resid_dev = [ -.2321102, -.3502712, -.6439626, -.2290982, 1.060478, -.2663844, -.2317827, -.3253788, -.4853875, .8555557, -.2225972, -.6491808, -.8819993, 1.813269, -.9463985, -.247583, -.3320177, -.2805444, -1.335131, .9103027, -.3559217, .4471892, -.744005, -1.955074, .5939538, 1.209638, .952332, -.8567857, .5870719, .335292, -1.227311, 2.096639] # from gretl obj.resid_generalized = [ -0.026578, -0.059501, -0.187260, -0.025902, 0.430107, -0.034858, -0.026504, -0.051559, -0.111127, 0.306489, -0.024470, -0.189997, -0.322240, 0.806789, -0.360990, -0.030184, -0.053626, -0.038588, -0.589872, 0.339214, -0.061376, 0.095153, -0.241772, -0.852091, 0.161709, 0.518867, 0.364579, -0.307219, 0.158296, 0.054660, -0.529117, 0.888969] obj.phat = np.array([ .02657799236476, .05950126051903, .18725991249084, .02590163610876, .56989300251007, .03485824912786, .02650404907763, .05155897513032, .11112663894892, .69351142644882, .02447037212551, .18999740481377, .32223951816559, .1932111531496, .36098992824554, .03018374741077, .05362640321255, .03858831897378, .58987241983414, .66078591346741, .06137581542134, .90484726428986, .24177247285843, .85209089517593, .8382905125618, .48113295435905, .63542068004608, .30721867084503, .84170418977737, .94534027576447, .52911710739136, .1110308393836]) obj.yhat = np.array([ -3.6007342338562, -2.7604126930237, -1.4679137468338, -3.6272060871124, .28141465783119, -3.3209850788116, -3.6035962104797, -2.9120934009552, -2.0792844295502, .81658720970154, -3.6855175495148, -1.4500269889832, -.74349880218506, -1.429278254509, -.57107019424438, -3.4698030948639, -2.8705959320068, -3.2154531478882, .36343798041344, .66679841279984, -2.7273993492126, 2.2522828578949, -1.1429864168167, 1.7510952949524, 1.6455633640289, -.07550399750471, .55554306507111, -.81315463781357, 1.6709630489349, 2.8504176139832, .11660042405128, -2.0802545547485]) obj.llf = -12.8896334653335 obj.llnull = -20.5917296966173 obj.df_model = 3 obj.df_resid = 32 - 4 # TODO: is this right? not reported in stata obj.llr = 15.4041924625676 obj.prsquared = .374038332124624 obj.llr_pvalue = .00150187761112892 obj.aic = 33.779266930667 obj.bic = 39.642210541866 obj.z = [2.237723415, 0.6722348408, 2.234423721, -2.640537645] obj.conf_int = [ [.3507938, 5.301432], [-.1822835, .3725988], [.29218, 4.465195], [-22.68657, -3.35613]] obj.pvalues = [.0252390974, .5014342039, .0254552063, .0082774596] # taken from margins command obj.margeff_nodummy_dydx = [ .36258084688424, .01220841099085, .30517768382304] obj.margeff_nodummy_dydx_se = [.1094412, .0177942, .0923796] obj.margeff_nodummy_dydxmean = [ .53385885781692, .01797548988961, .44933926079386] obj.margeff_nodummy_dydxmean_se = [.237038, .0262369, .1967626] obj.margeff_nodummy_dydxmedian = [ .25009492465091, .00842091261329, .2105003352955] obj.margeff_nodummy_dydxmedian_se = [.1546708, .0134314, .0928183] obj.margeff_nodummy_dydxzero = [ 6.252993785e-06, 2.105437138e-07, 5.263030788e-06] obj.margeff_nodummy_dydxzero_se = [.0000288, 9.24e-07, .000025] obj.margeff_nodummy_dyex = [ 1.1774000792198, .27896245178384, .16960002159996] obj.margeff_nodummy_dyex_se = [.3616481, .4090679, .0635583] obj.margeff_nodummy_dyexmean = [ 1.6641381583512, .39433730945339, .19658592659731] obj.margeff_nodummy_dyexmean_se = [.7388917, .5755722, .0860836] # NOTE: PSI at median should be a NaN or 'omitted' obj.margeff_nodummy_dyexmedian = [.76654095836557, .18947053379898, 0] obj.margeff_nodummy_dyexmedian_se = [.4740659, .302207, 0] # NOTE: all should be NaN # TODO: re above Note, since the values below are *not* NaN, # should something be done about this? obj.margeff_nodummy_dyexzero = [0, 0, 0] obj.margeff_nodummy_dyexzero_se = [0, 0, 0] obj.margeff_nodummy_eydx = [ 1.8546366266779, .06244722072812, 1.5610138123033] obj.margeff_nodummy_eydx_se = [.847903, .0930901, .7146715] obj.margeff_nodummy_eydxmean = [ 2.1116143062702, .0710998816585, 1.7773072368626] obj.margeff_nodummy_eydxmean_se = [1.076109, .1081501, .9120842] obj.margeff_nodummy_eydxmedian = [ 2.5488082240624, .0858205793373, 2.1452853812126] obj.margeff_nodummy_eydxmedian_se = [1.255377, .1283771, 1.106872] obj.margeff_nodummy_eydxzero = [ 2.8261067189993, .0951574597115, 2.3786824653103] obj.margeff_nodummy_eydxzero_se = [1.262961, .1415544, 1.064574] obj.margeff_nodummy_eyex = [ 5.4747106798973, 1.3173389907576, .44600395466634] obj.margeff_nodummy_eyex_se = [2.44682, 1.943525, .1567618] obj.margeff_nodummy_eyexmean = [ 6.5822977203268, 1.5597536538833, .77757191612739] obj.margeff_nodummy_eyexmean_se = [3.354433, 2.372543, .3990368] obj.margeff_nodummy_eyexmedian = [7.8120973525952, 1.9309630350892, 0] obj.margeff_nodummy_eyexmedian_se = [3.847731951, 2.888485089, 0] obj.margeff_nodummy_eyexzero = [0, 0, 0] obj.margeff_nodummy_eyexzero_se = [0, 0, 0] # for below GPA = 2.0, psi = 1 obj.margeff_nodummy_atexog1 = [ .1456333017086, .00490359933927, .12257689308426] obj.margeff_nodummy_atexog1_se = [.145633, .0111226, .1777101] # for below GPA at mean, tuce = 21, psi = 0 obj.margeff_nodummy_atexog2 = [ .25105129214546, .00845311433473, .2113052923675] obj.margeff_nodummy_atexog2_se = [.1735778, .012017, .0971515] # must get this from older margeff or i.psi then margins obj.margeff_dummy_dydx = [ .36258084688424, .01220841099085, .35751515254729] obj.margeff_dummy_dydx_se = [.1094412, .0177942, .1420034] obj.margeff_dummy_dydxmean = [ .53385885781692, .01797548988961, .4564984096959] obj.margeff_dummy_dydxmean_se = [.237038, .0262369, .1810537] # from margeff obj.margeff_dummy_count_dydx_median = [ 0.250110487483923, 0.008426867847905, 0.441897738279663] obj.margeff_dummy_count_dydx_median_se = [ .1546736661, .0134551951, .1792363708] # estimate with i.psi for the below then use margins obj.margeff_dummy_eydx = [ 1.8546366266779, .06244722072812, 1.5549034398832] obj.margeff_dummy_eydx_se = [.847903, .0930901, .7283702] # ie # margins, eydx(*) at((mean) _all) obj.margeff_dummy_eydxmean = [ 2.1116143062702, .0710998816585, 1.6631775707188] obj.margeff_dummy_eydxmean_se = [1.076109, .1081501, .801205] # TODO: Should we remove the commented-out code below? # Factor variables not allowed in below # test raises # obj.margeff_dummy_dydxzero # obj.margeff_dummy_eydxmedian # obj.margeff_dummy_eydxzero # obj.margeff_dummy_dyex # obj.margeff_dummy_dyexmean # obj.margeff_dummy_dyexmedian # obj.margeff_dummy_dyexzero # obj.margeff_dummy_eyex # obj.margeff_count_dummy_dydx_median # obj.margeff_count_dummy_dydx_median_se # NOTE: need old version of margeff for nodisc but at option is broken # stata command is margeff, count nodisc # this can be replicated with the new results by margeff # and then using margins for the last value obj.margeff_count_dydx = [.3625767598018, .0122068569914, .3051777] obj.margeff_count_dydx_se = [.1094379569, .0177869773, .0923796] # middle value taken from margeff rest from margins obj.margeff_count_dydxmean = [ .5338588, 0.01797186545386, .4493393] obj.margeff_count_dydxmean_se = [.237038, .0262211, .1967626] # with new version of margeff this is just a call to # margeff # mat list e(margeff_b), nonames format(%17.16g) obj.margeff_count_dummy_dydxoverall = [.362576759801767, .012206856991439, .357515163621704] # AFAICT, an easy way to get se is # mata # V = st_matrix("e(margeff_V)") # se = diagonal(cholesky(diag(V))) # last SE taken from margins with i.psi, do not know how they # do not know why margeff is different, but trust official results obj.margeff_count_dummy_dydxoverall_se = [.1094379569, .0177869773, .1420034] # from new margeff obj.margeff_count_dummy_dydxmean = [0.533849340033768, 0.017971865453858, 0.456498405282412] obj.margeff_count_dummy_dydxmean_se = [.2370202503, .0262210796, .1810536852] # for below GPA = 2.0, psi = 1 obj.margeff_dummy_atexog1 = [ .1456333017086, .00490359933927, .0494715429937] obj.margeff_dummy_atexog1_se = [.145633, .0111226, .0731368] # for below GPA at mean, tuce = 21, psi = 0 obj.margeff_dummy_atexog2 = [.25105129214546, .00845311433473, .44265645632553] obj.margeff_dummy_atexog2_se = [.1735778, .012017, .1811925] # The test for the prediction table was taken from Gretl # Gretl Output matched the Stata output here for params and SE obj.pred_table = np.array([[18, 3], [3, 8]]) return obj logit = logit() def probit(): obj = Namespace() obj.nobs = 32 obj.params = [ 1.62581025407, .051728948442, 1.42633236818, -7.45232041607] obj.cov_params = [ [.481472955383, -.01891350017, .105439226234, -1.1696681354], [-.01891350017, .00703757594, .002471864882, -.101172838897], [.105439226234, .002471864882, .354070126802, -.594791776765], [-1.1696681354, -.101172838897, -.594791776765, 6.46416639958]] obj.bse = [.693882522754, .083890261293, .595037920474, 2.54247249731] obj.llf = -12.8188033249334 obj.llnull = -20.5917296966173 obj.df_model = 3 obj.df_resid = 32 - 4 obj.llr = 15.5458527433678 obj.prsquared = .377478069409622 obj.llr_pvalue = .00140489496775855 obj.aic = 33.637606649867 obj.bic = 39.500550261066 obj.z = [2.343062695, .6166263836, 2.397044489, -2.931131182] obj.conf_int = [ [.2658255, 2.985795], [-.1126929, .2161508], [.2600795, 2.592585], [-12.43547, -2.469166]] obj.pvalues = [.0191261688, .537481188, .0165279168, .0033773013] obj.phat = [ .0181707, .0530805, .1899263, .0185707, .5545748, .0272331, .0185033, .0445714, .1088081, .6631207, .0161024, .1935566, .3233282, .1951826, .3563406, .0219654, .0456943, .0308513, .5934023, .6571863, .0619288, .9045388, .2731908, .8474501, .8341947, .488726, .6424073, .3286732, .8400168, .9522446, .5399595, .123544] obj.yhat = np.array([ -2.0930860042572, -1.615691781044, -.87816804647446, -2.0842070579529, .13722851872444, -1.9231110811234, -2.0856919288635, -1.6999372243881, -1.2328916788101, .42099541425705, -2.1418602466583, -.86486464738846, -.45841211080551, -.85895526409149, -.36825761198997, -2.0147502422333, -1.6881184577942, -1.8684275150299, .23630557954311, .40479621291161, -1.538782119751, 1.3078554868698, -.60319095849991, 1.025558590889, .97087496519089, -.02826354466379, .36490100622177, -.44357979297638, .99452745914459, 1.6670187711716, .10033150017262, -1.1574513912201]) obj.resid_dev = [ -.191509, -.3302762, -.6490455, -.1936247, 1.085867, -.2349926, -.1932698, -.3019776, -.4799906, .9064196, -.1801855, -.6559291, -.8838201, 1.807661, -.9387071, -.2107617, -.3058469, -.2503485, -1.341589, .9162835, -.3575735, .447951, -.7988633, -1.939208, .6021435, 1.196623, .9407793, -.8927477, .59048, .3128364, -1.246147, 2.045071] # Stata does not have it, but I think it's just oversight obj.resid_pearson = None # generalized residuals from gretl obj.resid_generalized = [ -0.045452, -0.114220, -0.334908, -0.046321, 0.712624, -0.064538, -0.046175, -0.098447, -0.209349, 0.550593, -0.040906, -0.340339, -0.530763, 1.413373, -0.579170, -0.053593, -0.100556, -0.071855, -0.954156, 0.559294, -0.130167, 0.187523, -0.457597, -1.545643, 0.298511, 0.815964, 0.581013, -0.538579, 0.289631, 0.104405, -0.862836, 1.652638] obj.pred_table = np.array([[18, 3], [3, 8]]) return obj probit = probit() class RandHIE: """ Results obtained from Stata 11 """ def __init__(self): self.nobs = 20190 def poisson(): obj = Namespace() obj.nobs = 20190 obj.params = [ -.052535114675, -.247086797633, .035290201794, -.03457750643, .271713973711, .033941474461, -.012635035534, .054056326828, .206115121809, .700352877227] obj.cov_params = None obj.bse = [ .00288398915279, .01061725196728, .00182833684966, .00161284852954, .01223913844387, .00056476496963, .00925061122826, .01530987068312, .02627928267502, .01116266712362] predict = np.loadtxt(os.path.join(cur_dir, 'yhat_poisson.csv'), delimiter=",") obj.phat = predict[:, 0] obj.yhat = predict[:, 1] obj.llf = -62419.588535018 obj.llnull = -66647.181687959 obj.df_model = 9 obj.df_resid = obj.nobs - obj.df_model - 1 obj.llr = 8455.186305881856 obj.prsquared = .0634324369893758 obj.llr_pvalue = 0 obj.aic = 124859.17707 obj.bic = 124938.306497 obj.z = [-18.21612769, -23.27219872, 19.30180524, -21.43878101, 22.20041672, 60.09840604, -1.36585953, 3.53081538, 7.84325525, 62.74063980] obj.conf_int = [ [-.0581876, -.0468826], [-0.2678962, -0.2262774], [0.0317067, 0.0388737], [-0.0377386, -0.0314164], [0.2477257, 0.2957022], [0.0328346, 0.0350484], [-0.0307659, 0.0054958], [0.0240495, 0.0840631], [0.1546087, 0.2576216], [0.6784745, 0.7222313]] obj.pvalues = [ 3.84415e-74, 8.4800e-120, 5.18652e-83, 5.8116e-102, 3.4028e-109, 0, .1719830562, .0004142808, 4.39014e-15, 0] # from stata # use margins and put i. in front of dummies obj.margeff_dummy_overall = [ -0.15027280560599, -0.66568074771099, 0.10094500919706, -0.09890639687842, 0.77721770295360, 0.09708707452600, -0.03608195237609, 0.15804581481115, 0.65104087597053] obj.margeff_dummy_overall_se = [ .008273103, .0269856266, .0052466639, .0046317555, .0351582169, .0016652181, .0263736472, .0457480115, .0913901155] # just use margins obj.margeff_nodummy_overall = [ -0.15027280560599, -0.70677348928158, 0.10094500919705, -0.09890639687842, 0.77721770295359, 0.09708707452600, -0.03614158359367, 0.15462412033340, 0.58957704430148] obj.margeff_nodummy_overall_se = [ .008273103, .0305119343, .0052466639, .0046317555, .0351582168, .0016652181, .0264611158, .0437974779, .0752099666] # taken from gretl obj.resid = np.loadtxt(os.path.join(cur_dir, 'poisson_resid.csv'), delimiter=",") return obj poisson = poisson() def negativebinomial_nb2_bfgs(): # R 2.15.1 MASS 7.3-22 glm.nb() obj = Namespace() obj.nobs = 20190 obj.params = [ -0.0579469537244314, -0.267787718814838, 0.0412060770911646, -0.0381376804392121, 0.268915772213171, 0.0381637446219235, -0.0441338846217674, 0.0172521803400544, 0.177960787443151, 0.663556087183864, # lnalpha from stata 1.292953339909746] # alpha and stderr from stata obj.lnalpha_std_err = .0143932 obj.lnalpha = 0.256929012449 obj.bse = [ 0.00607085853920512, 0.0226125368090765, 0.00405542008282773, 0.00344455937127785, 0.0298855063286547, 0.00142421904710063, 0.0199374393307107, 0.0358416931939136, 0.0741013728607101, 0.0250354082637892, # from stata .0186098] obj.z = [ -9.54510030998327, -11.8424447940467, 10.1607419822296, -11.071860382846, 8.99820030672628, 26.7962605187844, -2.21361850384595, 0.481343898758222, 2.40158556546135, 26.5047040652267] obj.pvalues = [ 1.35975947860026e-21, 2.35486776488278e-32, 2.96808970292151e-24, 1.71796558863781e-28, 2.2944789508802e-19, 3.57231639404726e-158, 0.0268550333379416, 0.630272102021494, 0.0163241908407114, 8.55476622951356e-155] obj.fittedvalues = [ 0.892904166867786, 0.892904166867786, 0.892904166867786, 0.892904166867786, 0.892904166867786, 0.937038051489553, 0.937038051489553, 0.937038051489553, 0.937038051489553, 0.937038051489553] # obj.aic = 86789.3241530713 # This is what R reports obj.aic = 86789.32415307125484 # from Stata obj.df_resid = 20180 obj.df_model = 9 # R conf_int: 1.96 * bse, not profile likelihood via R's confint() obj.conf_int = [ # from Stata [-.0698826, -.0460113], [-.3122654, -.2233101], [.0330781, .049334], [-.0448006, -.0314748], [.2102246, .3276069], [.0352959, .0410316], [-.0834356, -.0048321], [-.0535908, .0880951], [.0324115, .3235101], [.6150055, .7121067], # from Stata [1.256989, 1.329947] ] obj.bic = 86876.36652289562335 # stata obj.llnull = -44199.27443563430279 # stata obj.llr = 1631.224718197351 # stata obj.llf = -43383.66207653563 # stata obj.df_model = 9.0 obj.llr_pvalue = 0.0 return obj negativebinomial_nb2_bfgs = negativebinomial_nb2_bfgs() def negativebinomial_nb1_bfgs(): # Unpublished implementation intended for R's COUNT package. Sent by # J.Hilbe (of Cambridge UP NBin book) and Andrew Robinson to Vincent # Arel-Bundock on 2012-12-06. # TODO: Does the from Stata comment below apply to the # commented-out obj.params here or to everything below? If the # latter, then where do the commented-out version come from? # obj.params = [-0.065309744899923, -0.296016207412261, # 0.0411724098925173, -0.0320460533573259, 0.19083354227553, # 0.0318566232844115, -0.0331972813313092, -0.0484691550721231, # 0.111971860837541, 0.757560143822609, # 3.73086958562569] # from Stata obj = Namespace() obj.nobs = 20190 obj.params = [-.065317260803762961, -.296023807893893376, .041187021258044826, -.032028789543547605, .19065933246421754, .031871625115758778, -.033250849053302826, -.04850769174426571, .111813637465757343, .757277086555503409, 3.731151380800305] # lnalpha and lnalpha_std_err are from stata obj.lnalpha = 1.316716867203 obj.lnalpha_std_err = .0168876692 obj.bse = [ 0.00536019929563678, 0.0196998350459769, 0.00335779098766272, 0.00301145915122889, 0.0237984097096245, 0.00107360844112751, 0.0167174614755359, 0.0298037989274781, 0.0546838603596457, 0.0214703279904911, 0.0630011409376052] obj.z = [ -12.1842008660173, -15.0263292419148, 12.2617548393554, -10.6413707601675, 8.0187518663633, 29.6724784046551, -1.98578482623631, -1.62627439508848, 2.04762173155154, 35.2840508145997, # From R, this is alpha/bse(alpha) 59.2190796881069 # taken from Stata even though they do not report it # lnalpha/bse(lnalpha) # 77.968995 ] obj.conf_int = [ [-0.075815736, -0.0548037543], [-0.334627884, -0.2574045307], [0.034591140, 0.0477536802], [-0.037948513, -0.0261435934], [0.144188659, 0.2374784253], [0.029752351, 0.0339608958], [-0.065963506, -0.0004310568], [-0.106884601, 0.0099462908], [0.004791495, 0.2191522271], [3.607387349, 3.8543518219], [0.715478301, 0.7996419867]] # from Stata obj.llf = -43278.75612911823 obj.llnull = -44199.2744356343 obj.llr = 1841.036613032149 obj.aic = 86579.51225823645655 obj.bic = 86666.55462806082505 obj.llr_pvalue = 0.0 obj.df_model = 9.0 obj.df_resid = 20180.0 # Smoke tests TODO: check against other stats package obj.pvalues = [ 3.65557865e-034, 5.24431864e-051, 1.42921171e-034, 2.09797259e-026, 1.15949461e-015, 1.56785415e-193, 4.71746349e-002, 1.04731854e-001, 4.07534831e-002, 1.95504975e-272, 0.00000000e+000] obj.conf_int = [ [-.0758236, -.054811], [-.3346363, -.2574113], [.0346053, .0477687], [-.0379314, -.0261261], [.1440119, .2373067], [.0297667, .0339766], [-.0660178, -.0004839], [-.1069241, .0099087], [.0046266, .2190007], [.7151889, .7993652], # from stata for alpha no lnalpha [3.609675, 3.856716]] # [1.28360034e+00, 1.34979803e+00]] obj.fittedvalues = [ 0.8487497, 0.8487497, 0.8487497, 0.8487497, 0.8487497, 0.88201746, 0.88201746, 0.88201746, 0.88201746, 0.88201746] return obj negativebinomial_nb1_bfgs = negativebinomial_nb1_bfgs() def negativebinomial_geometric_bfgs(): # Smoke tests TODO: Cross check with other stats package obj = Namespace() obj.nobs = 20190 obj.params = [ -0.05768894, -0.26646696, 0.04088528, -0.03795503, 0.26885821, 0.03802523, -0.04308456, 0.01931675, 0.18051684, 0.66469896] obj.bse = [ 0.00553867, 0.02061988, 0.00375937, 0.0030924, 0.02701658, 0.00132201, 0.01821646, 0.03271784, 0.06666231, 0.02250053] obj.pvalues = [ 2.10310916e-025, 3.34666368e-038, 1.50697768e-027, 1.25468406e-034, 2.48155744e-023, 6.18745348e-182, 1.80230194e-002, 5.54919603e-001, 6.77044178e-003, 8.44913440e-192] obj.z = [-10.41567024, -12.92281571, 10.8755779, -12.27364916, 9.95160202, 28.76323587, -2.36514487, 0.59040434, 2.70792943, 29.54148082] obj.aic = 87101.159433012392 # old value 87101.160011780419 obj.bic = 87180.288860125467 # old value 87180.289438893495 obj.df_model = 9.0 obj.df_resid = 20180.0 obj.llf = -43540.58000589021 obj.llnull = -44586.650971362695 # old value -44199.27443567125 obj.llr = 2092.1425097129977 # old value 1317.3888595620811 obj.llr_pvalue = 0 # old value 5.4288002863296022e-278 obj.fittedvalues = [ 0.89348994, 0.89348994, 0.89348994, 0.89348994, 0.89348994, 0.9365745, 0.9365745, 0.9365745, 0.9365745, 0.9365745] obj.conf_int = [ [-0.06854453, -0.04683335], [-0.30688118, -0.22605273], [0.03351706, 0.04825351], [-0.04401602, -0.03189404], [0.21590669, 0.32180972], [0.03543415, 0.04061632], [-0.07878816, -0.00738096], [-0.04480903, 0.08344253], [0.04986111, 0.31117258], [0.62059873, 0.70879919]] return obj negativebinomial_geometric_bfgs = negativebinomial_geometric_bfgs() def generalizedpoisson_gp2(): # Stata gnpoisson function obj = Namespace() obj.nobs = 20190 obj.llf = -43326.42720093228 obj.params = [-0.0604495342, -0.277717228, 0.0438136144, -0.0395811744, 0.273044906, 0.0399108677, -0.0552626543, -0.001227569488, 0.151980519, 0.651125316, 0.448085318] obj.lnalpha_std_err = 0.0125607 obj.lnalpha = -0.8027716 obj.bse = [0.00634704, 0.02381906, 0.00443871, 0.00355094, 0.0334247, 0.00166303, 0.02102142, 0.0390845, 0.087821, 0.02626823, 0.00562825] obj.df_model = 9 obj.aic = 86674.854401865 obj.conf_int = [ [-0.07288951, -0.04800956], [-0.32440173, -0.23103272], [0.03511389, 0.05251333], [-0.04654088, -0.03262147], [0.20753371, 0.33855610], [0.03665139, 0.04317034], [-0.09646387, -0.01406144], [-0.07783191, 0.07537652], [-0.02014548, 0.32410651], [0.59964053, 0.70261011], [0.43718883, 0.45925338] ] obj.bic = 86761.896771689 obj.wald_pvalue = 4.8795019354e-254 obj.wald_statistic = 1206.46339591254 return obj generalizedpoisson_gp2 = generalizedpoisson_gp2() def zero_inflated_poisson_logit(): obj = Namespace() obj.nobs = 20190 obj.params = [.1033783, -1.045983, -.0821979, .0085692, -.0267957, 1.482363] obj.llf = -57005.72199826186 obj.bse = [0.0079912, 0.02235510, .0107145, 0.0018697, 0.0014121, 0.0085915] obj.conf_int = [ [0.0877159, 0.1190408], [-1.089798, -1.002167], [-0.1031979, -0.061198], [0.0049045, 0.0122338], [-0.0295635, -0.024028], [1.465524, 1.499202]] obj.aic = 114023.444 obj.bic = 114070.9 return obj zero_inflated_poisson_logit = zero_inflated_poisson_logit() def zero_inflated_poisson_probit(): obj = Namespace() obj.nobs = 20190 obj.params = [.0622534, -.6429324, -.0821788, .0085673, -.0267952, 1.482369] obj.llf = -57006.05 obj.bse = [.0048228, .0132516, .0107142, .0018697, .0014121, .0085913] obj.conf_int = [ [0.0528009, .0717058], [-0.6689051, -.6169597], [-0.1031783, -.0611793], [0.0049027, .0122319], [-0.0295629, -.0240275], [1.46553, 1.499208]] obj.aic = 114024.1 obj.bic = 114071.6 return obj zero_inflated_poisson_probit = zero_inflated_poisson_probit() def zero_inflated_poisson_offset(): obj = Namespace() obj.nobs = 20190 obj.params = [.1052014, -1.082434, -.0922822, .0115868, -.0283842, 1.347514] obj.llf = -58207.67 obj.bse = [.0081836, .0230043, .0107788, .0018687, .0014162, .0086309] obj.conf_int = [ [.0891619, .1212409], [-1.127522, -1.037347], [-.1134082, -.0711561], [.0079242, .0152494], [-.0311599, -.0256085], [1.330598, 1.36443]] obj.aic = 116427.3 obj.bic = 116474.8 return obj zero_inflated_poisson_offset = zero_inflated_poisson_offset() def zero_inflated_generalized_poisson(): obj = Namespace() obj.nobs = 20190 obj.params = [3.57337, -17.95797, -0.21380, 0.03847, -0.05348, 1.15666, 1.36468] obj.llf = -43630.6 obj.bse = [1.66109, 7.62052, 0.02066, 0.00339, 0.00289, 0.01680, 0.01606] obj.aic = 87275 return obj zero_inflated_generalized_poisson = zero_inflated_generalized_poisson() def zero_inflated_negative_binomial(): obj = Namespace() obj.params = [1.883859, -10.280888, -0.204769, 1.137985, 1.344457] obj.llf = -44077.91 obj.bse = [0.3653, 1.6694, 0.02178, 0.01163, 0.0217496] obj.aic = 88165.81 return obj zero_inflated_negative_binomial = zero_inflated_negative_binomial() def truncated_poisson(self): self.params = [.000241, -.0002085, .0059419, .0016508, 2.289668] self.llf = -11091.36 self.bse = [.0042718, .0157994, .002791, .0023491, .0123182] self.conf_int = [[-.0081316, .0086136], [-.0311747, .0307577], [0.0004716, .0114122], [-.0029532, .0062549], [2.265525, 2.313812]] self.aic = 22192.72 self.bic = 22222.87 def zero_truncated_poisson(self): self.params = [-.0095444, -.07581, .0080299, -.022879, 1.486142] self.llf = -44570.853 self.bse = [.0029817, .0109677, .0019165, .0016393, .0085697] self.conf_int = [[-.0153884, -0.0037005], [-.0973064, -0.0543137], [0.0042736, 0.0117861], [-.0260919, -0.019666], [1.469346, 1.502939]] self.aic = 89151.71 self.bic = 89189.39 def zero_truncted_nbp(self): self.params = [-0.0397571554, -0.1269997715, 0.0002753126, 0.9294720085, 2.246603] self.llf = -31565 self.bse = [0.0074415, 0.0308599, 0.0053180, 0.0342113, 0.1030674] self.aic = 63139.99 self.bic = 63177.69