1

我是 R 2.15.2 的 Win-7 用户

有人可以帮我为什么下面的模型不能很好地接近简单的 logit 模型估计吗？

已编辑

 Mydata <- structure(list(gg = c(13.659955, 6.621436486, 3.017166776, 2.516795069, 
3.928538296, 4.211960532, 3.235445955, 5.152860411, 18.96466673, 
5.904678823, 4.987622293, 1.170687541, 3.088224149, 4.738065529, 
3.263022593, 6.050017999, 5.650762257, 2.058924721, 3.138591919, 
7.169083435, 11.30381738, 3.036991188, 4.559013218, 3.978760664, 
3.617455798, 2.430111184, 4.440319959, 2.200267742, 6.003166495, 
3.114161526, 3.812363457, 11.12816724, 15.6564348, 13.50562576, 
5.154056904, 6.26451889, 1.849669635, 1.816757851, 3.861868285, 
2.884542233, 2.993444924, 2.724235493, 2.694159089, 1.973597356, 
4.371300647, 3.559035718, 3.59124243, 6.587196681, 10.03402072, 
4.805158339, 4.491460392, 1.627936721, 1.278291553, 0.978710462, 
3.08635052, 2.58594947, 2.354973563, 1.657519171, 2.946994656, 
2.110549733, 6.095182338, 6.000660354, 6.691960157, 1.796172588, 
2.531234555, 2.992017156, 2.882403206, 6.066420081, 5.930524609, 
5.972280022, 0.915755208, 2.398369176, 2.614677323, 3.904309459, 
2.120045842, 2.643895472, 4.756728133, 4.420426549, 5.787620043, 
10.14764507, 10.17993063, 4.529582451, 2.132173713, 3.573017249, 
2.749427487, 5.79057317, 3.30235352, 2.248377801, 3.37704153, 
3.270636543, 8.070781663, 10.63482235, 11.91206431, 2.030994094, 
3.843890691, 11.08130402, 10.56164432, 2.476733196, 2.672953595, 
4.73313584, 2.577761507, 4.397950314, 16.02648073, 1.979318271, 
4.124346958, 18.95385675, 5.364182526, 3.964093091, 4.985103684, 
5.887163026, 6.302977032, 3.594056904, 2.1592159, 2.506654073, 
2.739740977, 5.840023906, 14.66431143, 1.453209431, 1.405915159, 
6.133806412, 14.46294647, 15.7475232, 13.91882816, 3.667160495, 
2.472826568, 2.208138746, 2.313856848, 6.677046217, 5.960024468, 
2.647365332, 1.960115215, 2.392076123, 2.796502016, 2.147861348, 
3.825593606, 2.427949377, 3.111598763, 5.248771445, 3.333653698, 
3.937363645, 1.98551364, 4.090444736, 1.897938127, 3.772878129, 
1.708975251, 2.149179807, 1.985344614, 3.445809412, 1.537484672, 
3.210425705, 3.635236243, 3.340757476, 2.426647042, 16.65453923, 
2.207695228, 17.05735446, 9.148235024, 3.289077716, 2.534469392, 
1.827870535, 2.952377144, 2.218054092, 3.787767414, 1.862731602, 
6.699397675, 2.572247774, 2.42783988, 3.429506609, 2.152120465, 
6.651671979, 3.389825443, 1.787049217, 2.368181588, 2.358737508, 
3.470515609, 1.693830224, 4.124049967, 2.894483735, 1.90612459, 
3.447950783, 11.73691853, 12.78330833, 3.860448673, 2.571402925, 
2.66582863, 3.206179619, 2.546483922, 3.047299616, 2.052597638, 
4.894592388, 2.869309084, 2.851843536, 1.729619574, 2.276794788, 
12.68354926, 3.318892896, 10.71333438, 3.318892896, 10.71333438, 
1.908315292, 3.111350542, 1.847059969, 1.835544656, 2.618919198, 
2.341098688, 2.920079188, 1.830521448, 4.95152276, 10.18359208, 
11.37645094, 3.44132725, 11.51865833, 14.45180625, 12.61570573, 
10.80784104, 7.814012406, 14.48656688, 1.019333021, 2.816970615, 
7.101350458, 1.764183823, 4.821036885, 2.251823823, 1.569206635, 
1.791906542, 3.114256927, 2.496611198, 2.776086885, 2.477500021, 
4.801364906, 3.642995521, 4.279315979, 4.956759854, 3.710843927, 
15.75846719, 3.032181385, 11.60129177, 3.443892792, 2.642908469, 
5.329881656, 2.303688396, 2.681312552, 2.12796449, 2.773085833, 
4.189449094, 3.521328188, 1.761418094, 2.455309406, 2.526749219, 
10.27017188, 11.70942021, 6.537176688, 3.665906281, 4.212212875, 
7.625210531, 10.74009115, 10.74009115, 10.58339073, 11.77831083, 
2.604995813, 4.744413865, 7.29975851, 1.740470719, 1.112195677, 
2.604728, 2.58798274, 1.059064521, 4.349711208, 2.392650167, 
2.699035417, 2.153652521, 3.151369396, 3.414949031, 2.581368219, 
2.234642333, 2.392818135, 6.51472726, 2.648815969, 15.47915854, 
3.258030083, 1.990451802, 2.6811855, 2.440177469, 2.925240396, 
3.668138604, 2.917453344, 1.9656045, 2.452494323, 4.299348479, 
4.765304604, 16.62252406, 12.88866406, 2.488524583, 1.916951229, 
1.714486281, 2.816999552, 1.784181771, 2.585827844, 2.720170823, 
2.186287625, 3.07640024, 4.327105385, 4.086370688, 4.592656531, 
2.571516719, 6.396392573, 18.74246573, 14.24556208), ss = c(-7.707850005, 
-15.40300834, -49.84338333, -50.96595763, -33.57927065, -29.26871379, 
-31.13757476, -38.9446836, -18.60621918, -12.25557326, -17.3340105, 
-49.04469485, -33.97249836, -24.13226024, -23.37858067, -24.07898472, 
-32.68854786, -35.16417362, -30.17502297, -38.58802006, -13.44056998, 
-45.62119153, -30.10387457, -21.92073966, -30.01274679, -30.10716134, 
-29.60556014, -33.1849011, -24.07374895, -36.78247492, -34.72067685, 
-41.96947502, -9.624331958, -10.03069504, -41.01141746, -35.72784007, 
-50.94674135, -49.91973751, -51.05524609, -50.47122434, -46.90646667, 
-42.70813537, -50.7377332, -51.01120184, -51.0932699, -35.37166589, 
-39.92974313, -29.5446461, -8.1278036, -24.74676854, -36.47346677, 
-40.92182994, -49.92750914, -49.02615356, -35.42230149, -38.96690541, 
-35.20704603, -39.31404894, -31.98014625, -38.22213743, -25.34981719, 
-14.02261648, -20.16701322, -48.09641042, -34.72517484, -50.17069654, 
-40.66006187, -44.93973938, -43.63455705, -50.96538014, -49.74130308, 
-34.21321834, -43.03926971, -26.14166307, -51.01781663, -46.38799194, 
-51.07775476, -51.09319771, -44.82805756, -40.02547202, -15.16109778, 
-23.15247117, -41.51283866, -39.06283772, -43.02607856, -38.197886, 
-45.18138933, -50.9752189, -50.55407144, -49.39313584, -33.4766776, 
-13.54483547, -12.00068998, -33.6073029, -20.23229774, -32.90792069, 
-16.32554139, -25.90690728, -43.18838286, -51.0203703, -45.17489885, 
-50.00289866, -29.7772429, -36.57372832, -51.02230524, -50.7863748, 
-24.23635605, -44.15949002, -39.39035061, -36.81914971, -26.99636543, 
-18.03405081, -49.68132933, -32.0432924, -21.71221224, -45.80810068, 
-33.70033655, -41.01746977, -44.9545008, -40.45390925, -33.26664292, 
-9.718327927, -9.277726352, -31.89539233, -50.26845505, -50.90332052, 
-50.55203806, -49.13701041, -36.73852817, -35.43596888, -38.08838755, 
-41.12237836, -40.41028312, -37.69615075, -35.55166964, -36.72707884, 
-29.55487485, -46.23696728, -21.60602981, -47.63092247, -47.00620137, 
-42.59296241, -50.71814475, -45.10917784, -51.02946189, -20.89401331, 
-48.86575795, -36.89427487, -40.16853942, -51.00840255, -51.0075457, 
-45.30865098, -46.05110903, -35.31902034, -44.80360176, -37.35104715, 
-13.84902784, -39.74186369, -31.01123746, -35.21037218, -35.0202728, 
-48.39866317, -32.74663729, -50.90494703, -41.74675073, -27.83560326, 
-37.14191338, -39.24284429, -50.24699166, -41.91108934, -44.47568388, 
-35.2452714, -34.49466485, -31.56238493, -33.86060467, -39.88969345, 
-50.92874004, -50.61210368, -33.35654261, -49.73027655, -25.53924159, 
-8.408177651, -44.11358958, -42.78740696, -38.98041905, -50.93700572, 
-50.78977219, -47.03148495, -48.05189369, -36.41525546, -41.82258554, 
-48.0012309, -50.65991657, -34.58154214, -35.59513453, -51.08074135, 
-50.85366938, -51.08074135, -50.85366938, -50.69545917, -22.14701885, 
-43.04649469, -40.17308802, -42.42595479, -39.31338521, -40.29414865, 
-46.82535656, -38.09278958, -35.71546021, -38.48812479, -40.15711031, 
-41.26210979, -25.71862146, -26.0133424, -26.7382999, -27.87618542, 
-17.20437385, -50.21773375, -50.43229958, -51.10613073, -44.76030594, 
-24.495235, -40.8543576, -49.56576271, -50.8176775, -33.32086031, 
-34.70547198, -36.75163896, -50.58161188, -35.58172323, -50.78794125, 
-40.18822771, -34.27604302, -48.19352396, -51.01194635, -43.13363927, 
-45.60329094, -27.14570313, -39.71587979, -35.66745896, -49.72861292, 
-50.98749646, -36.75274792, -51.06588781, -50.65578125, -50.87232552, 
-51.08351406, -51.10034323, -30.78426354, -51.11306302, -19.71852042, 
-18.5394999, -51.0472551, -51.0868051, -47.5098225, -17.50247563, 
-17.50247563, -17.86539229, -16.16612115, -30.19224563, -28.90203208, 
-48.59109542, -50.87648281, -47.88816083, -28.53676927, -48.9707774, 
-50.57288479, -51.07964656, -37.37621177, -49.56755563, -50.76012125, 
-47.20241073, -51.01347313, -50.82764219, -50.78897458, -50.83598677, 
-51.11415594, -49.00892, -33.5060326, -28.38547938, -48.05074729, 
-50.29370125, -50.83504219, -50.80333563, -51.03403729, -51.0421874, 
-50.76831146, -50.87043604, -51.08487479, -51.01771458, -36.0702924, 
-20.78084771, -49.44148479, -39.28489156, -50.94829448, -51.01100906, 
-50.76106125, -50.92056417, -51.04772927, -37.12336448, -50.6384801, 
-31.84211313, -49.48543917, -51.05363948, -28.35589769, -40.21823813, 
-51.0574976, -39.4211176), dd = c(154.4236993, 149.5832755, 220.5376676, 
188.1301303, 193.4423268, 180.786566, 147.6809975, 294.3175776, 
517.5095247, 106.0196438, 126.7522827, 84.12371801, 153.8300366, 
167.6862942, 111.8452798, 213.6434143, 270.8933346, 105.7571576, 
138.8542046, 405.7345177, 222.8295022, 203.1856754, 201.2825818, 
127.8993813, 159.2263148, 107.2677566, 192.7990438, 107.0794478, 
211.9558545, 167.9804162, 194.1213087, 684.9904472, 220.9984063, 
198.6784757, 310.0050436, 328.2581513, 138.2133496, 132.9782976, 
289.1817849, 213.5229493, 205.9192931, 170.6115496, 200.483163, 
147.652789, 327.5662511, 184.6250023, 210.3079779, 285.4338146, 
119.6048186, 174.3630449, 240.233224, 97.64322396, 93.5913059, 
70.31214681, 160.2304584, 147.7143808, 121.5668885, 95.49747258, 
138.205053, 118.2084279, 226.6079498, 123.4014062, 197.9286116, 
126.6151589, 128.7066092, 220.1513453, 171.8320615, 399.8383894, 
379.5150089, 446.4220118, 66.78233993, 120.2137058, 165.0139027, 
149.6845692, 158.637227, 179.8320709, 356.3473798, 331.2533515, 
380.5121027, 595.6989688, 226.3574763, 153.7853848, 129.7965126, 
204.6755039, 173.477482, 324.3927627, 218.8154401, 168.0964657, 
250.3859192, 236.9220868, 396.2428049, 211.2501547, 209.6575701, 
100.0337171, 114.0287522, 534.8363364, 252.8850567, 93.47141464, 
169.2018281, 354.1812225, 167.0160495, 322.5169401, 699.927721, 
106.0177501, 308.6366926, 1411.84738, 190.5941033, 256.7291835, 
287.9890035, 317.9030093, 249.5487016, 94.82178213, 157.3068154, 
117.7808006, 85.46573263, 392.3572701, 724.8147089, 87.31754945, 
92.59793475, 363.9298959, 705.6552358, 224.452264, 189.3891441, 
171.4922284, 182.3002062, 164.8533233, 171.5477266, 481.1937846, 
321.1358114, 137.5508765, 109.432054, 144.223193, 165.7262492, 
118.6953876, 199.4566326, 130.7380613, 134.8405831, 355.9356426, 
105.5132727, 275.0446049, 136.8457579, 255.506731, 141.1784194, 
249.5957908, 127.9098903, 65.54666729, 142.2750164, 186.4378832, 
90.52911409, 240.1828537, 271.9479516, 221.9813631, 163.8951252, 
862.7211025, 145.0432736, 934.4263336, 185.813134, 191.6960064, 
115.1740602, 94.3482938, 151.5457861, 157.4328115, 181.8996531, 
139.0691666, 410.1901284, 104.873104, 132.2095716, 197.3701322, 
158.5931377, 408.8709478, 221.1149527, 92.28489547, 119.7748008, 
109.1577013, 172.347586, 99.0082713, 308.0480454, 214.8575513, 
93.23180651, 251.3186275, 439.633449, 157.6421674, 249.7593794, 
161.3505953, 152.2448205, 239.5266898, 189.6886754, 210.1612637, 
144.5761132, 261.4069935, 175.9787194, 200.7477829, 128.5085872, 
115.4151402, 662.1582919, 248.6459146, 799.002826, 248.6459146, 
799.002826, 141.8859927, 100.0343638, 116.5021167, 108.026373, 
162.9137625, 134.8998615, 172.545199, 125.597699, 276.6302083, 
533.4431354, 642.188624, 202.6657, 697.0803198, 545.1278615, 
481.3199646, 423.8338427, 319.4693844, 365.5350729, 75.06200271, 
208.3551844, 532.2892583, 115.6034833, 173.1313906, 134.4010906, 
114.0157948, 133.5525719, 152.0830708, 126.9199573, 149.5574833, 
183.7883156, 250.5106854, 271.3588146, 252.2147146, 249.1717438, 
262.2825646, 1178.945583, 191.7794198, 775.945224, 137.0065938, 
153.8710573, 278.8008677, 167.9884, 200.5136229, 114.6238458, 
207.6997042, 311.252324, 262.7331438, 131.9810083, 184.0258792, 
112.6623184, 770.1853177, 338.6306063, 177.741876, 274.4640479, 
315.6091333, 531.3221333, 275.6941573, 275.6941573, 277.3068708, 
279.2577229, 114.9657708, 201.0928177, 520.2239594, 129.8706885, 
77.96273083, 108.235849, 185.6859594, 78.54580635, 325.8620802, 
130.7865177, 196.174249, 160.3316583, 218.100225, 255.507301, 
192.4301396, 166.4553438, 178.4062406, 488.3921406, 190.3591188, 
760.6767625, 135.469474, 140.0819813, 197.7584458, 181.9313792, 
217.9605552, 274.5590406, 218.4073167, 146.3560396, 182.9803677, 
322.1222083, 356.5695281, 879.3790313, 392.8225521, 180.428799, 
110.1808227, 128.1094563, 210.7529375, 132.8284542, 193.1165542, 
203.6622229, 118.8509177, 228.4791198, 202.0587625, 296.5748396, 
343.8937052, 102.326891, 377.2940917, 1403.520031, 823.6412771
), OutCome = c(1L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 
0L, 0L, 0L, 1L, 0L, 0L, 0L, 1L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 
0L, 0L, 0L, 1L, 1L, 1L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 
0L, 0L, 0L, 1L, 1L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 
0L, 0L, 1L, 0L, 0L, 0L, 0L, 1L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 
0L, 1L, 0L, 1L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 1L, 
1L, 0L, 0L, 1L, 1L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 1L, 0L, 0L, 
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 1L, 1L, 1L, 0L, 
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 1L, 
1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 
0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 1L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 
0L, 0L, 0L, 0L, 0L, 1L, 0L, 1L, 1L, 1L, 1L, 1L, 1L, 0L, 0L, 1L, 
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 1L, 0L, 
1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 1L, 1L, 
0L, 0L, 1L, 1L, 1L, 1L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 
0L, 0L, 0L, 0L, 0L, 0L, 1L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 
0L, 0L, 0L, 0L, 0L, 1L, 0L, 1L)), .Names = c("gg", "ss", "dd", 
"OutCome"), class = "data.frame", row.names = c(1L, 2L, 3L, 4L, 
5L, 6L, 7L, 8L, 9L, 10L, 11L, 12L, 13L, 14L, 15L, 16L, 17L, 18L, 
19L, 20L, 21L, 22L, 23L, 24L, 25L, 26L, 27L, 28L, 29L, 30L, 31L, 
32L, 33L, 34L, 35L, 36L, 37L, 38L, 39L, 40L, 41L, 42L, 43L, 44L, 
45L, 46L, 47L, 48L, 49L, 50L, 51L, 52L, 53L, 54L, 55L, 56L, 57L, 
58L, 59L, 60L, 61L, 62L, 63L, 64L, 65L, 66L, 67L, 68L, 69L, 70L, 
71L, 72L, 73L, 74L, 75L, 76L, 77L, 78L, 79L, 80L, 81L, 82L, 83L, 
84L, 85L, 86L, 87L, 88L, 89L, 90L, 91L, 92L, 93L, 94L, 95L, 96L, 
97L, 98L, 99L, 100L, 101L, 102L, 103L, 104L, 105L, 106L, 107L, 
108L, 109L, 110L, 111L, 112L, 113L, 114L, 115L, 116L, 117L, 118L, 
119L, 120L, 121L, 122L, 123L, 124L, 125L, 126L, 127L, 128L, 129L, 
130L, 131L, 132L, 133L, 134L, 135L, 136L, 137L, 138L, 139L, 140L, 
141L, 142L, 143L, 144L, 145L, 146L, 147L, 148L, 149L, 150L, 151L, 
152L, 153L, 154L, 155L, 156L, 157L, 158L, 159L, 160L, 161L, 162L, 
163L, 164L, 165L, 166L, 167L, 168L, 169L, 170L, 171L, 172L, 173L, 
174L, 175L, 176L, 177L, 178L, 179L, 180L, 181L, 182L, 183L, 184L, 
185L, 186L, 187L, 188L, 189L, 190L, 191L, 192L, 193L, 194L, 195L, 
196L, 197L, 198L, 199L, 200L, 201L, 202L, 203L, 204L, 205L, 206L, 
207L, 208L, 209L, 210L, 211L, 212L, 233L, 234L, 235L, 236L, 237L, 
238L, 239L, 240L, 241L, 242L, 243L, 244L, 245L, 246L, 247L, 248L, 
249L, 250L, 251L, 252L, 253L, 254L, 255L, 256L, 257L, 258L, 259L, 
260L, 261L, 262L, 263L, 264L, 265L, 266L, 267L, 268L, 269L, 270L, 
271L, 272L, 273L, 274L, 275L, 276L, 277L, 278L, 279L, 280L, 281L, 
282L, 283L, 284L, 285L, 286L, 287L, 288L, 289L, 290L, 291L, 292L, 
293L, 294L, 295L, 296L, 297L, 298L, 299L, 300L, 301L, 302L, 303L, 
304L, 305L, 306L, 307L, 308L, 309L, 310L, 311L, 312L, 313L, 314L, 
315L, 316L, 317L, 318L, 319L, 320L, 321L, 322L, 323L, 324L, 325L, 
326L, 349L, 350L))

似然估计的模型代码：

Simplelogit <- glm(OutCome ~ gg+ss+dd, data = Mydata, family = "binomial")

使用R2WinBUGS的模型代码：（已编辑）

model1 ="
model
{
  # likelihood
  for(i in 1:N)
    {
    Y[i] ~ dbin(p[i],N)
    logit(p[i])<- beta1[1]+beta1[2]*X[1]+beta1[3]*X[2]+beta1[4]*X[3]
    }

    #prior
    beta1[1]~dnorm(1,1.0E-02)
    beta1[2]~dnorm(1,1.0E-02)
    beta1[3]~dnorm(1,1.0E-02)
    beta1[4]~dnorm(1,1.0E-02)
}
"
writeLines(model1,con='Model.txt')

x1 <- unlist(Mydata$gg)
x2 <- unlist(Mydata$ss)
x3 <- unlist(Mydata$dd)
N=c(nrow(Mydata))
datalist <- list(N=N,Y=c(Mydata$OutCome),X=c(x1,x2,x3))
inits <- function() list(beta1=c((Simplelogit$coefficients[,1])))
MyPara <- c("beta1")

require(R2WinBUGS)
BayesianModel <- bugs(datalist,inits,MyPara,model.file='Model.txt',n.chains=1,n.iter=54000,n.burnin=4000,n.sim=50000,program=c('WinBUGS'),debug=FALSE,codaPkg=FALSE,save.history=TRUE,bugs.directory='C:/Program Files/WinBUGS14/',working.directory = getwd()) #,over.relax=TRUE

as.numeric(BayesianModel$summary[c(1:4)),1])
#results:
-48.63550   3.47384  -0.69866   0.09043

然后使用传统方法/不使用贝叶斯方法

Simplelogit <- glm(OutCome ~ gg+ss+dd, data = Mydata, family = "binomial")
c(as.matrix(Simplelogit$coefficients[c(1:4),1]))
# result is:
-20.71281   3.47408  -0.31233  -0.03906

请建议我是否需要使用不同的模型来更改先验或语法...

4

3 回答 3

5

我还没有运行代码，但我可以发现两个错误：

没有Mydata$yy，所以向量太短（只有616，应该是3*308）。应该是x3<-unlist(Mydata$dd)。

而且您没有注意到错误，因为 logit 行中的索引是错误的。应该是这样的

logit(p[i])<- beta1[1]+beta1[2]*X[i]+beta1[3]*X[i+2*N]+beta1[4]*X[i+3*N]

于 2013-01-25T07:36:44.913 回答
3

使用 stan 的另一种解决方案

load("mydata.rdata")
library(rstan)
library(ggmcmc)
library(coda)

model1 ="
data {
  int<lower=0>  N;
  int<lower=0,upper=1> Y[N];
  real gg[N];
  real ss[N];
  real dd[N];
}
parameters{
  real  beta0;
  real  betagg;
  real  betass;
  real  betadd;
}
model
{
  #prior
  beta0 ~ normal(-2,30);
  betagg ~ normal(20,30);
  betass ~ normal(-3,30);
  betadd ~ normal(-10,40);

# likelihood
  for(i in 1:N)
  {
    Y[i] ~ bernoulli(inv_logit(beta0+betagg*gg[i]+betass*ss[i]+betadd*dd[i]));
  }

}
"

MyPar = scale(Mydata[,-4])
datalist <- list(N=nrow(Mydata),
                 Y=as.numeric(Mydata$OutCome),
                 gg=MyPar[,"gg"],ss=MyPar[,"ss"],dd=MyPar[,"dd"])

m <- stan(model_code=model1,iter=20000,data= datalist,n.chains=4)

ggmcmc(ggs(m))
print(m)
于 2013-01-25T17:22:16.853 回答
3

jags 版本（我讨厌安装 RWinBugs）

# Assuming your data have been saved in mydata.rdata
load("mydata.rdata")
library("rjags")

model1 ="
model
{
  # likelihood
  for(i in 1:N)
  {
    logit(p[i])<- beta0+betagg*gg[i]+betass*ss[i]+betadd*dd[i];
    Y[i] ~ dbin(p[i],N); # Should be dbern probably
  }

  #prior
  beta0~dnorm(1,1.0E-02);
  betagg~dnorm(1,1.0E-02);
  betass~dnorm(1,1.0E-02);
  betadd~dnorm(1,1.0E-02);
}
"
writeLines(model1,con='Model.txt')

datalist <- with(Mydata, list(N=nrow(Mydata),Y=as.numeric(OutCome),gg=gg,ss=ss,dd=dd))
# A bit of cheating: initial values adapted after first run
inits <- list(beta0=-8,betagg=0.2,betass=0.05,betadd=0.002)

m <- jags.model("Model.txt",datalist,init=inits)
update(m, 1000)
x <- coda.samples(m, c("beta0","betagg","betass","betadd"), n.iter=10000)
plot(x) # Well, not prettty, but acceptable
于 2013-01-25T11:19:10.410 回答