# RLM MODEL RESULTS import numpy as np def _shift_intercept(arr): """ A convenience function to make the SAS covariance matrix compatible with statsmodels.rlm covariance """ arr = np.asarray(arr) side = int(np.sqrt(len(arr))) return np.roll(np.roll(arr.reshape(side, side), -1, axis=1), -1, axis=0) class Huber: """ """ def __init__(self): self.params = np.array([0.82937387, 0.92610818, -0.12784916, -41.02653105]) self.bse = np.array([0.11118035, 0.3034081, 0.12885259, 9.8073472]) self.scale = 2.4407137948148447 self.weights = np.array([ 1., 1., 0.7858871, 0.50494094, 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 0.36814106]) self.resid = np.array([ 3.05027584, -2.07757332, 4.17721071, 6.50163171, -1.64615192, -2.57226011, -1.73127333, -0.73127333, -2.25476463, 0.48083217, 1.63147461, 1.42973363, -2.26346951, -0.78323693, 2.26646556, 0.88291808, -0.83307835, 0.06186577, 0.26360675, 1.54306186, -8.91752986]) self.df_model = 3 self.df_resid = 17 self.bcov_unscaled = np.array([ [1.72887367e-03, -3.47079127e-03, -6.79080082e-04, 2.73387119e-02], [-3.47079127e-03, 1.28754242e-02, 9.95952051e-07, -6.19611175e-02], [-6.79080082e-04, 9.95952051e-07, 2.32216722e-03, -1.59355028e-01], [2.73387119e-02, -6.19611175e-02, -1.59355028e-01, 1.34527267e+01]]) # From R self.fittedvalues = np.array([ 38.94972416, 39.07757332, 32.82278929, 21.49836829, 19.64615192, 20.57226011, 20.73127333, 20.73127333, 17.25476463, 13.51916783, 12.36852539, 11.57026637, 13.26346951, 12.78323693, 5.73353444, 6.11708192, 8.83307835, 7.93813423, 8.73639325, 13.45693814, 23.91752986]) self.tvalues = np.array([7.45971657, 3.0523516, -0.99221261, -4.18324448]) # from R this is equivalent to # self.res1.params/np.sqrt(np.diag(self.res1.bcov_scaled)) # The below are taken from SAS huber_h1 = [ 95.8813, 0.19485, -0.44161, -1.13577, 0.1949, 0.01232, -0.02474, -0.00484, -0.4416, -0.02474, 0.09177, 0.00001, -1.1358, -0.00484, 0.00001, 0.01655] h1 = _shift_intercept(huber_h1) huber_h2 = [ 82.6191, 0.07942, -0.23915, -0.95604, 0.0794, 0.01427, -0.03013, -0.00344, -0.2392, -0.03013, 0.10391, -0.00166, -0.9560, -0.00344, -0.00166, 0.01392] h2 = _shift_intercept(huber_h2) huber_h3 = [ 70.1633, -0.04533, -0.00790, -0.78618, -0.0453, 0.01656, -0.03608, -0.00203, -0.0079, -0.03608, 0.11610, -0.00333, -0.7862, -0.00203, -0.00333, 0.01138] h3 = _shift_intercept(huber_h3) class Hampel: """ """ def __init__(self): self.params = np.array([0.74108304, 1.22507934, -0.14552506, -40.47473236]) self.bse = np.array([0.13482596, 0.36793632, 0.1562567, 11.89315426]) self.scale = 3.0882646556217064 self.weights = np.array([ 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 0.80629719]) self.resid = np.array([ 3.06267708, -2.08284798, 4.36377602, 5.78635972, -1.7634816, -2.98856094, -2.34048993, -1.34048993, -3.02422878, 1.08249252, 2.39221804, 2.47177232, -1.62645737, -0.25076107, 2.32088237, 0.88430719, -1.37812296, -0.35944755, -0.43900184, 1.40555003, -7.65988702]) self.df_model = 3 self.df_resid = 17 self.bcov_unscaled = np.array([ [1.72887367e-03, -3.47079127e-03, -6.79080082e-04, 2.73387119e-02], [-3.47079127e-03, 1.28754242e-02, 9.95952051e-07, -6.19611175e-02], [-6.79080082e-04, 9.95952051e-07, 2.32216722e-03, -1.59355028e-01], [2.73387119e-02, -6.19611175e-02, -1.59355028e-01, 1.34527267e+01]]) self.fittedvalues = np.array([ 38.93732292, 39.08284798, 32.63622398, 22.21364028, 19.7634816, 20.98856094, 21.34048993, 21.34048993, 18.02422878, 12.91750748, 11.60778196, 10.52822768, 12.62645737, 12.25076107, 5.67911763, 6.11569281, 9.37812296, 8.35944755, 9.43900184, 13.59444997, 22.65988702]) self.tvalues = np.array([5.49659011, 3.32959607, -0.93132046, -3.40319578]) hampel_h1 = [ 141.309, 0.28717, -0.65085, -1.67388, 0.287, 0.01816, -0.03646, -0.00713, -0.651, -0.03646, 0.13524, 0.00001, -1.674, -0.00713, 0.00001, 0.02439] h1 = _shift_intercept(hampel_h1) hampel_h2 = [ 135.248, 0.18207, -0.36884, -1.60217, 0.182, 0.02120, -0.04563, -0.00567, -0.369, -0.04563, 0.15860, -0.00290, -1.602, -0.00567, -0.00290, 0.02329] h2 = _shift_intercept(hampel_h2) hampel_h3 = [ 128.921, 0.05409, -0.02445, -1.52732, 0.054, 0.02514, -0.05732, -0.00392, -0.024, -0.05732, 0.18871, -0.00652, -1.527, -0.00392, -0.00652, 0.02212] h3 = _shift_intercept(hampel_h3) class BiSquare: def __init__(self): self.params = np.array([0.9275471, 0.65073222, -0.11233103, -42.28525369]) self.bse = np.array([0.10805398, 0.29487634, 0.12522928, 9.5315672]) self.scale = 2.2818858795649497 self.weights = np.array([ 0.89283149, 0.88496132, 0.79040651, 0.3358111, 0.94617358, 0.90040725, 0.96630596, 0.99729171, 0.94968061, 0.99900087, 0.98959903, 0.99831448, 0.84731833, 0.96455873, 0.91767906, 0.98724523, 0.99762848, 0.99694419, 0.98650731, 0.95897484, 0.00222999]) self.resid = np.array([ 2.50917802, -2.60315301, 3.56070896, 6.93256033, -1.76597524, -2.41670746, -1.39345348, -.39345348, -1.70651907, -0.23917521, 0.77180408, 0.31020526, -3.01451315, -1.42960401, 2.19218084, 0.85518774, -0.36817892, 0.4181383, 0.87973712, 1.53911661, -10.43556344]) self.df_model = 3 self.df_resid = 17 self.bcov_unscaled = np.array([ [1.72887367e-03, -3.47079127e-03, -6.79080082e-04, 2.73387119e-02], [-3.47079127e-03, 1.28754242e-02, 9.95952051e-07, -6.19611175e-02], [-6.79080082e-04, 9.95952051e-07, 2.32216722e-03, -1.59355028e-01], [2.73387119e-02, -6.19611175e-02, -1.59355028e-01, 1.34527267e+01]]) self.fittedvalues = np.array([ 39.49082198, 39.60315301, 33.43929104, 21.06743967, 19.76597524, 20.41670746, 20.39345348, 20.39345348, 16.70651907, 14.23917521, 13.22819592, 12.68979474, 14.01451315, 13.42960401, 5.80781916, 6.14481226, 8.36817892, 7.5818617, 8.12026288, 13.46088339, 25.43556344]) self.tvalues = np.array([8.58410823, 2.20679698, -0.8970029, -4.43633799]) bisquare_h1 = [ 90.3354, 0.18358, -0.41607, -1.07007, 0.1836, 0.01161, -0.02331, -0.00456, -0.4161, -0.02331, 0.08646, 0.00001, -1.0701, -0.00456, 0.00001, 0.01559] h1 = _shift_intercept(bisquare_h1) bisquare_h2 = [ 67.82521, 0.091288, -0.29038, -0.78124, 0.091288, 0.013849, -0.02914, -0.00352, -0.29038, -0.02914, 0.101088, -0.001, -0.78124, -0.00352, -0.001, 0.011766] h2 = _shift_intercept(bisquare_h2) bisquare_h3 = [ 48.8983, 0.000442, -0.15919, -0.53523, 0.000442, 0.016113, -0.03461, -0.00259, -0.15919, -0.03461, 0.112728, -0.00164, -0.53523, -0.00259, -0.00164, 0.008414] h3 = _shift_intercept(bisquare_h3) class Andrews: def __init__(self): self.params = [0.9282, 0.6492, -.1123, -42.2930] self.bse = [.1061, .2894, .1229, 9.3561] self.scale = 2.2801 self.df_model = 3. self.df_resid = 17. # bcov_unscaled not defined because it is not given as part of SAS self.resid = [ 2.503338458, -2.608934536, 3.5548678338, 6.9333705014, -1.768179527, -2.417404513, -1.392991531, -0.392991531, -1.704759385, -0.244545418, 0.7659115325, 0.3028635237, -3.019999429, -1.434221475, 2.1912017882, 0.8543828047, -0.366664104, 0.4192468573, 0.8822948661, 1.5378731634, -10.44592783] self.sresids = [ 1.0979293816, -1.144242351, 1.5591155202, 3.040879735, -0.775498914, -1.06023995, -0.610946684, -0.172360612, -0.747683723, -0.107254214, 0.3359181307, 0.1328317233, -1.324529688, -0.629029563, 0.9610305856, 0.3747203984, -0.160813769, 0.1838758324, 0.3869622398, 0.6744897502, -4.581438458] self.weights = [ 0.8916509101, 0.8826581922, 0.7888664106, 0.3367252734, 0.9450252405, 0.8987321912, 0.9656622, 0.9972406688, 0.948837669, 0.9989310017, 0.9895434667, 0.998360628, 0.8447116551, 0.9636222149, 0.916330067, 0.9869982597, 0.9975977354, 0.9968600162, 0.9861384742, 0.9582432444, 0] def conf_int(self): # method to be consistent with sm return [(0.7203, 1.1360), (.0819, 1.2165), (-.3532, .1287), (-60.6305, -23.9555)] andrews_h1 = [ 87.5357, 0.177891, -0.40318, -1.03691, 0.177891, 0.01125, -0.02258, -0.00442, -0.40318, -0.02258, 0.083779, 6.481E-6, -1.03691, -0.00442, 6.481E-6, 0.01511] h1 = _shift_intercept(andrews_h1) andrews_h2 = [ 66.50472, 0.10489, -0.3246, -0.76664, 0.10489, 0.012786, -0.02651, -0.0036, -0.3246, -0.02651, 0.09406, -0.00065, -0.76664, -0.0036, -0.00065, 0.011567] h2 = _shift_intercept(andrews_h2) andrews_h3 = [ 48.62157, 0.034949, -0.24633, -0.53394, 0.034949, 0.014088, -0.02956, -0.00287, -0.24633, -0.02956, 0.100628, -0.00104, -0.53394, -0.00287, -0.00104, 0.008441] h3 = _shift_intercept(andrews_h3) # RLM Results with Huber's Proposal 2 # Obtained from SAS class HuberHuber: def __init__(self): self.h1 = [ 114.4936, 0.232675, -0.52734, -1.35624, 0.232675, 0.014714, -0.02954, -0.00578, -0.52734, -0.02954, 0.10958, 8.476E-6, -1.35624, -0.00578, 8.476E-6, 0.019764] self.h1 = _shift_intercept(self.h1) self.h2 = [ 103.2876, 0.152602, -0.33476, -1.22084, 0.152602, 0.016904, -0.03766, -0.00434, -0.33476, -0.03766, 0.132043, -0.00214, -1.22084, -0.00434, -0.00214, 0.017739] self.h2 = _shift_intercept(self.h2) self.h3 = [ 91.7544, 0.064027, -0.11379, -1.08249, 0.064027, 0.019509, -0.04702, -0.00278, -0.11379, -0.04702, 0.157872, -0.00462, -1.08249, -0.00278, -0.00462, 0.015677] self.h3 = _shift_intercept(self.h3) self.resid = [ 2.909155172, -2.225912162, 4.134132661, 6.163172632, -1.741815737, -2.789321552, -2.02642336, -1.02642336, -2.593402734, 0.698655, 1.914261011, 1.826699492, -2.031210331, -0.592975466, 2.306098648, 0.900896645, -1.037551854, -0.092080512, -0.004518993, 1.471737448, -8.498372406] self.sresids = [ 0.883018497, -0.675633129, 1.25483702, 1.870713355, -0.528694904, -0.84664529, -0.615082113, -0.311551209, -0.787177874, 0.212063383, 0.581037374, 0.554459746, -0.616535106, -0.179986379, 0.699972205, 0.273449972, -0.314929051, -0.027949281, -0.001371654, 0.446717797, -2.579518651] self.weights = [ 1, 1, 1, 0.718977066, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0.52141511] self.params = (.7990, 1.0475, -0.1351, -41.0892) self.bse = (.1213, .3310, .1406, 10.7002) self.scale = 3.2946 self.df_model = 3 self.df_resid = 17 def conf_int(self): # method for consistency with sm return [(0.5612, 1.0367), (.3987, 1.6963), (-.4106, .1405), (-62.0611, -20.1172)] class HampelHuber: def __init__(self): self.h1 = [ 147.4727, 0.299695, -0.67924, -1.7469, 0.299695, 0.018952, -0.03805, -0.00744, -0.67924, -0.03805, 0.141144, 0.000011, -1.7469, -0.00744, 0.000011, 0.025456] self.h1 = _shift_intercept(self.h1) self.h2 = [ 141.148, 0.190007, -0.38493, -1.67206, 0.190007, 0.02213, -0.04762, -0.00592, -0.38493, -0.04762, 0.165518, -0.00303, -1.67206, -0.00592, -0.00303, 0.024301] self.h2 = _shift_intercept(self.h2) self.h3 = [ 134.5444, 0.05645, -0.02552, -1.59394, 0.05645, 0.026232, -0.05982, -0.00409, -0.02552, -0.05982, 0.196946, -0.0068, -1.59394, -0.00409, -0.0068, 0.023083] self.h3 = _shift_intercept(self.h3) self.resid = [ 3.125725599, -2.022218392, 4.434082972, 5.753880172, -1.744479058, -2.995299443, -2.358455878, -1.358455878, -3.068281354, 1.150212629, 2.481708553, 2.584584946, -1.553899388, -0.177335865, 2.335744732, 0.891912757, -1.43012351, -0.394515569, -0.497391962, 1.407968887, -7.505098501] self.sresids = [ 0.952186413, -0.616026205, 1.350749906, 1.752798302, -0.531418771, -0.912454834, -0.718453867, -0.413824947, -0.934687235, 0.350388031, 0.756000196, 0.787339321, -0.473362692, -0.054021633, 0.711535395, 0.27170242, -0.43565698, -0.120180852, -0.151519976, 0.428908041, -2.28627005] self.weights = [ 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0.874787298] self.params = (.7318, 1.2508, -0.1479, -40.2712) self.bse = (.1377, .3757, .1596, 12.1438) self.scale = 3.2827 self.df_model = 3 self.df_resid = 17 def conf_int(self): return [(0.4619, 1.0016), (.5145, 1.9872), (-.4607, .1648), (-64.0727, -16.4697)] class BisquareHuber: def __init__(self): self.h1 = [ 129.9556, 0.264097, -0.59855, -1.5394, 0.264097, 0.016701, -0.03353, -0.00656, -0.59855, -0.03353, 0.124379, 9.621E-6, -1.5394, -0.00656, 9.621E-6, 0.022433] self.h1 = _shift_intercept(self.h1) self.h2 = [ 109.7685, 0.103038, -0.25926, -1.28355, 0.103038, 0.0214, -0.04688, -0.00453, -0.25926, -0.04688, 0.158535, -0.00327, -1.28355, -0.00453, -0.00327, 0.018892] self.h2 = _shift_intercept(self.h2) self.h3 = [ 91.80527, -0.09171, 0.171716, -1.05244, -0.09171, 0.027999, -0.06493, -0.00223, 0.171716, -0.06493, 0.203254, -0.0071, -1.05244, -0.00223, -0.0071, 0.015584] self.h3 = _shift_intercept(self.h3) self.resid = [ 3.034895447, -2.09863887, 4.229870063, 6.18871385, -1.715906134, -2.763596142, -2.010080245, -1.010080245, -2.590747917, 0.712961901, 1.914770759, 1.82892645, -2.019969464, -0.598781979, 2.260467209, 0.859864256, -1.057306197, -0.122565974, -0.036721665, 1.471074632, -8.432085298] self.sresids = [ 0.918227061, -0.634956635, 1.279774287, 1.872435025, -0.519158394, -0.836143718, -0.608162656, -0.305606249, -.78384738, 0.215711191, 0.579326161, 0.553353415, -0.611154703, -0.181165324, 0.683918836, 0.26015744, -0.319894764, -0.037083121, -0.011110375, 0.445083055, -2.551181429] self.weights = [ 0.924649089, 0.963600796, 0.856330585, 0.706048833, 0.975591792, 0.937309703, 0.966582366, 0.991507994, 0.944798311, 0.995764589, 0.969652425, 0.972293856, 0.966255569, 0.997011618, 0.957833493, 0.993842376, 0.990697247, 0.9998747, 0.999988752, 0.982030803, 0.494874977] self.params = (.7932, 1.0477, -0.1335, -40.8949) self.bse = (.1292, .3527, .1498, 11.3998) self.scale = 3.3052 self.df_model = 3 self.df_resid = 17 def conf_int(self): return [(0.5399, 1.0465), (.3565, 1.7389), (-.4271, .1600), (-63.2381, -18.5517)] class AndrewsHuber: def __init__(self): self.h1 = [ 129.9124, 0.264009, -0.59836, -1.53888, 0.264009, 0.016696, -0.03352, -0.00656, -0.59836, -0.03352, 0.124337, 9.618E-6, -1.53888, -0.00656, 9.618E-6, 0.022425] self.h1 = _shift_intercept(self.h1) self.h2 = [ 109.7595, 0.105022, -0.26535, -1.28332, .105022, 0.021321, -0.04664, -0.00456, -0.26535, -0.04664, 0.157885, -0.00321, -1.28332, -0.00456, -0.00321, 0.018895] self.h2 = _shift_intercept(self.h2) self.h3 = [ 91.82518, -0.08649, 0.155965, -1.05238, -0.08649, 0.027785, -0.06427, -0.0023, 0.155965, -0.06427, 0.201544, -0.00693, -1.05238, -0.0023, -0.00693, 0.015596] self.h3 = _shift_intercept(self.h3) self.resid = [ 3.040515104, -2.093093543, 4.235081748, 6.188729166, -1.714119676, -2.762695255, -2.009618953, -1.009618953, -2.591649784, 0.715967584, 1.918445405, 1.833412337, -2.016815123, -0.595695587, 2.260536347, 0.859710406, -1.059386228, -0.1241257, -0.039092633, 1.471556455, -8.424624872] self.sresids = [ 0.919639919, -0.633081011, 1.280950793, 1.871854667, -0.518455862, -0.835610004, -0.607833129, -0.305371248, -0.783875269, 0.216552902, 0.580256606, 0.554537345, -0.610009696, -0.180175208, 0.683726076, 0.260029627, -0.320423952, -0.037543293, -0.011824031, 0.445089734, -2.548127888] self.weights = [ 0.923215335, 0.963157359, 0.854300342, 0.704674258, 0.975199805, 0.936344742, 0.9660077, 0.991354016, 0.943851708, 0.995646409, 0.968993767, 0.971658421, 0.965766352, 0.99698502, 0.957106815, 0.993726436, 0.990483134, 0.999868981, 0.999987004, 0.981686004, 0.496752113] self.params = (.7928, 1.0486, -0.1336, -40.8818) self.bse = (.1292, .3526, .1498, 11.3979) self.scale = 3.3062 self.df_model = 3 self.df_resid = 17 def conf_int(self): return [(0.5395, 1.0460), (.3575, 1.7397), (-.4271, .1599), (-63.2213, -18.5423)]