Why Gemfury?
Push, build, and install
RubyGems
npm packages
Python packages
Maven artifacts
PHP packages
Go Modules
Debian packages
RPM packages
NuGet packages
alkaline-ml / statsmodels python
Repository URL to install this package:
|
Version:
0.10.2 ▾
|
"""
Test Results for discrete models from Stata
"""
import os
import numpy as np
class Namespace(object):
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(object):
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(object):
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(as_pandas=False)
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(as_pandas=False)
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(as_pandas=False)
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(object):
"""
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, don't know how they
# don't 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 doesn't 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(object):
"""
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 don't 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()