library(maptools)
polyOHIO<-readSplus("Ohio_splusMAP.txt")
plot(polyOHIO)
#can use fillmap with this 

#getting adj and num into INLA form
library(INLA)
source("fillmap.R")
inla.geobugs2inla(adj, num, graph.file="OHIO_map")

# S-T model data 
yL=c(5, 14, 12, 15, 9, 12, 20, 12, 16, 15,
76, 46, 53, 47, 62, 69, 53, 65, 69, 58,
11, 22, 18, 20, 27, 25, 34, 19, 28, 30,
38, 58, 56, 53, 59, 56, 56, 50, 53, 65,
30, 25, 33, 24, 21, 38, 21, 15, 27, 24,
18, 17, 20, 23, 16, 35, 22, 23, 19, 22,
40, 54, 64, 47, 45, 55, 44, 50, 62, 70,
7, 21, 12, 22, 27, 24, 19, 23, 30, 23,
91, 135, 129, 138, 150, 141, 139, 164, 159, 158,
12, 7, 12, 14, 16, 17, 11, 8, 19, 10,
12, 17, 14, 20, 19, 18, 19, 19, 11, 21,
60, 99, 80, 88, 73, 87, 92, 95, 114, 105,
51, 62, 71, 59, 51, 75, 94, 80, 89, 75,
11, 11, 16, 15, 19, 17, 26, 19, 22, 24,
55, 51, 63, 64, 62, 74, 65, 68, 67, 60,
13, 22, 13, 13, 18, 15, 11, 28, 13, 19,
23, 28, 21, 37, 21, 26, 35, 30, 24, 31,
825, 917, 913, 909, 940, 923, 937, 953, 952, 993,
22, 18, 18, 23, 23, 24, 37, 28, 28, 35,
15, 12, 13, 12, 10, 16, 12, 25, 16, 22,
25, 17, 10, 24, 27, 14, 24, 23, 28, 22,
33, 41, 32, 36, 40, 42, 40, 43, 32, 54,
41, 33, 38, 48, 53, 33, 54, 51, 54, 45,
15, 14, 22, 15, 15, 14, 16, 23, 17, 18,
402, 402, 431, 451, 478, 444, 524, 499, 488, 543,
12, 9, 15, 8, 16, 18, 9, 7, 16, 11,
20, 15, 17, 12, 15, 10, 6, 16, 19, 29,
20, 22, 29, 36, 25, 20, 20, 38, 28, 30,
49, 44, 44, 36, 57, 50, 71, 50, 60, 52,
17, 22, 24, 31, 32, 19, 22, 18, 33, 27,
528, 583, 542, 554, 568, 642, 639, 572, 593, 576,
24, 29, 22, 24, 41, 27, 31, 33, 36, 29,
11, 18, 15, 18, 25, 15, 26, 19, 20, 21,
7, 8, 12, 20, 22, 7, 7, 16, 12, 14,
6, 7, 13, 12, 2, 12, 23, 16, 9, 9,
24, 25, 13, 17, 24, 24, 25, 19, 21, 24,
9, 14, 14, 15, 15, 18, 18, 10, 17, 12,
4, 3, 8, 9, 5, 6, 8, 9, 6, 8,
26, 25, 28, 23, 18, 21, 21, 33, 23, 16,
21, 22, 12, 19, 14, 14, 15, 15, 26, 18,
57, 55, 68, 65, 54, 68, 56, 71, 68, 77,
25, 15, 22, 27, 28, 26, 26, 23, 27, 37,
87, 101, 83, 94, 114, 103, 127, 111, 108, 147,
35, 39, 26, 45, 36, 41, 37, 33, 46, 49,
57, 69, 72, 57, 68, 84, 61, 76, 61, 81,
11, 19, 31, 24, 24, 20, 25, 24, 30, 25,
121, 115, 110, 136, 125, 124, 124, 137, 136, 149,
280, 243, 282, 277, 292, 243, 308, 305, 290, 305,
10, 10, 13, 14, 24, 20, 21, 12, 15, 16,
137, 142, 159, 163, 148, 147, 183, 200, 163, 187,
25, 27, 36, 41, 37, 26, 40, 50, 40, 55,
27, 38, 40, 36, 34, 41, 45, 31, 40, 42,
19, 12, 14, 12, 9, 14, 17, 15, 19, 16,
16, 12, 15, 14, 14, 9, 15, 15, 12, 12,
44, 48, 37, 49, 63, 39, 48, 47, 44, 42,
7, 4, 7, 6, 7, 7, 6, 6, 8, 11,
259, 314, 306, 325, 345, 326, 341, 374, 361, 340,
7, 8, 8, 4, 7, 6, 11, 11, 13, 11,
9, 16, 16, 6, 13, 9, 12, 10, 14, 10,
49, 54, 59, 44, 59, 62, 60, 45, 75, 66,
6, 3, 2, 5, 5, 10, 9, 7, 2, 5,
10, 15, 21, 18, 26, 25, 28, 14, 26, 26,
13, 11, 15, 8, 9, 11, 6, 11, 14, 9,
10, 21, 15, 25, 20, 23, 23, 25, 14, 14,
27, 23, 15, 23, 21, 24, 22, 23, 26, 27,
8, 14, 10, 10, 20, 17, 8, 11, 13, 23,
38, 37, 49, 66, 55, 65, 58, 44, 62, 55,
13, 24, 13, 21, 27, 21, 22, 26, 18, 18,
11, 13, 14, 9, 13, 10, 12, 16, 14, 11,
45, 62, 66, 63, 59, 66, 50, 79, 79, 67,
27, 23, 33, 28, 34, 36, 29, 25, 36, 42,
23, 27, 23, 26, 24, 25, 34, 33, 22, 32,
41, 54, 38, 56, 44, 70, 56, 71, 58, 67,
19, 30, 17, 28, 36, 39, 25, 31, 41, 40,
9, 19, 21, 19, 30, 23, 22, 21, 16, 22,
182, 197, 183, 188, 206, 198, 205, 217, 244, 224,
247, 277, 253, 282, 310, 302, 292, 293, 323, 337,
91, 114, 112, 111, 133, 114, 139, 143, 146, 146,
28, 36, 40, 54, 57, 37, 58, 49, 38, 51,
4, 17, 9, 9, 14, 11, 10, 14, 12, 6,
9, 13, 12, 17, 13, 17, 13, 15, 17, 10,
6, 5, 7, 7, 13, 3, 6, 7, 9, 4,
43, 33, 37, 40, 42, 42, 55, 40, 56, 59,
20, 26, 34, 30, 32, 39, 27, 40, 40, 40,
36, 27, 28, 30, 30, 31, 46, 38, 45, 32,
10, 17, 23, 10, 13, 17, 17, 13, 22, 17,
28, 32, 30, 30, 44, 37, 55, 40, 45, 30,
7, 7, 11, 5, 10, 16, 10, 11, 21, 16)
eL=c(11.08265, 12.43046, 12.58617, 13.07113, 13.64872, 13.44648, 14.20172, 14.38548, 14.95299, 15.44726, 
51.54766, 57.20417, 57.04018, 58.30136, 61.00702, 60.78989, 64.1098, 64.32432, 66.04699, 67.65539, 
21.18653, 23.58796, 23.70401, 24.58979, 25.74367, 25.60512, 27.00496, 26.73573, 27.48413, 28.12705, 
48.02128, 53.02829, 52.80107, 54.51591, 56.9513, 56.11595, 58.9642, 58.08966, 58.91479, 59.60639, 
25.54036, 28.85985, 29.14606, 30.34142, 31.9459, 31.69799, 33.39785, 33.45538, 34.54226, 35.47209, 
19.65611, 21.73145, 21.82331, 22.744, 23.8938, 23.95328, 25.31747, 25.35509, 26.03825, 26.72759, 
38.56623, 42.0636, 41.97363, 43.43165, 45.88792, 44.81808, 46.23249, 44.58007, 45.22339, 45.36022, 
14.60163, 16.3542, 16.67166, 17.49995, 18.45894, 18.54305, 19.89741, 19.88945, 20.67724, 21.42489, 
118.7359, 132.3991, 133.8755, 139.09, 145.8855, 146.4967, 155.7944, 157.4248, 163.6123, 169.2534, 
11.82541, 13.03703, 13.08939, 13.6795, 14.60024, 14.62257, 15.47169, 15.4304, 15.81459, 16.15364, 
15.67667, 17.15562, 17.08506, 17.82085, 18.65204, 18.47113, 19.61539, 19.92314, 20.61976, 21.28398, 
69.8358, 76.60449, 76.37443, 78.53021, 82.3222, 81.24256, 85.3237, 85.31212, 87.24273, 89.11415, 
57.71282, 65.82852, 67.17506, 70.67236, 74.95687, 75.22547, 80.34092, 81.57098, 85.39924, 88.77426, 
16.01126, 17.66246, 17.66381, 18.29873, 19.23135, 19.14769, 20.11954, 20.21182, 20.97767, 21.57185, 
52.67453, 57.85574, 57.74398, 60.01879, 62.86696, 61.75908, 64.5058, 62.88037, 64.11284, 64.84256, 
16.55965, 18.39023, 18.49928, 19.18715, 20.17, 20.04829, 20.82488, 20.73979, 21.24078, 21.6668, 
23.48562, 25.5304, 25.47703, 26.19371, 27.21013, 27.18946, 28.543, 28.23602, 28.82986, 29.35596, 
702.4902, 762.9026, 756.0005, 779.8555, 822.0158, 809.6126, 845.5759, 844.0905, 855.8276, 869.7341, 
25.70603, 28.01751, 27.7992, 28.40706, 29.93484, 29.77232, 31.3417, 31.05596, 31.78028, 32.41493, 
18.48065, 20.33881, 20.22841, 20.71682, 21.59112, 21.23599, 22.69088, 22.68849, 23.27034, 23.87908, 
24.805, 27.5521, 27.93598, 29.38989, 31.19252, 31.27092, 33.3397, 34.18607, 35.54015, 36.99492, 
36.24261, 40.57072, 40.45171, 41.55404, 43.33169, 42.92893, 44.9724, 44.73515, 45.64945, 46.47362, 
42.30873, 47.8858, 48.49051, 50.51246, 53.0288, 52.6939, 55.83695, 56.66023, 58.79983, 60.75547, 
12.7047, 13.97347, 13.9019, 14.48464, 15.41407, 15.19235, 15.95782, 15.92759, 16.30998, 16.69189, 
406.4846, 444.3345, 447.9958, 467.6531, 495.0504, 493.6799, 522.6015, 529.1328, 546.3256, 564.1479, 
17.45652, 19.26888, 19.32098, 20.00927, 21.10864, 21.07612, 22.27337, 21.93282, 22.58947, 23.08501, 
13.63535, 15.36918, 15.39472, 15.77881, 16.67008, 16.62236, 17.48306, 17.37735, 17.9514, 18.41551, 
34.03191, 37.97007, 38.03563, 39.53481, 41.63126, 41.20462, 43.48262, 44.34541, 45.61567, 47.0657, 
60.06467, 66.19062, 66.30437, 69.08276, 72.43422, 71.29946, 75.27091, 75.69876, 78.17459, 80.2487, 
19.00359, 21.44709, 21.43185, 22.10858, 23.18576, 22.81702, 23.65325, 23.04279, 23.46055, 23.6662, 
407.7868, 445.5072, 444.9867, 461.4336, 484.2658, 478.9633, 503.4001, 507.0756, 517.8945, 530.1678, 
29.92659, 32.92372, 33.0029, 34.45148, 36.3201, 36.24579, 38.31725, 37.87029, 38.78331, 39.59003, 
14.95705, 16.66925, 16.62521, 17.09685, 17.73578, 17.54398, 18.44368, 18.25615, 18.80826, 19.24466, 
8.223136, 9.199208, 9.044219, 9.307426, 9.734615, 9.357822, 9.683545, 9.43737, 9.2791, 9.210176, 
13.21699, 14.46189, 14.36684, 14.86706, 15.76502, 15.66644, 16.51372, 16.39807, 16.78344, 17.14971, 
15.48276, 17.13619, 17.35708, 17.92799, 18.85522, 18.89875, 19.95381, 20.02827, 20.82716, 21.48658, 
11.08913, 12.39415, 12.39784, 12.78047, 13.67559, 13.66001, 14.29243, 14.34889, 14.86173, 15.29123, 
13.66543, 15.03009, 15.261, 15.9618, 16.94154, 16.78611, 17.82091, 17.68577, 18.01362, 18.35322, 
25.26315, 27.883, 27.92015, 28.84994, 30.42291, 30.3825, 31.86446, 32.1276, 33.03179, 33.8531, 
14.05833, 15.58909, 15.45903, 15.98089, 16.97009, 16.77393, 17.66914, 17.47493, 17.97451, 18.35927, 
42.13704, 46.5668, 46.22499, 47.71302, 49.26917, 48.36682, 50.33897, 48.49604, 48.90682, 49.00401, 
21.60812, 23.64268, 23.77036, 24.95947, 26.3241, 26.09414, 27.57889, 27.29623, 28.12766, 28.81045, 
98.29676, 108.4951, 108.9542, 113.7492, 119.6884, 117.9353, 123.3566, 123.7011, 126.4734, 129.3168, 
29.23195, 32.54986, 32.50121, 33.32862, 35.20066, 34.97511, 36.30995, 36.13476, 36.90544, 37.5362, 
55.50258, 61.91143, 62.72108, 65.61343, 69.11733, 68.54123, 72.42336, 72.3067, 74.21561, 76.06302, 
17.98687, 19.98131, 20.0411, 20.71045, 21.93815, 21.89872, 23.17992, 23.14502, 23.98498, 24.67195, 
126.3616, 139.9924, 139.737, 144.8432, 151.641, 149.4939, 156.9894, 155.4053, 158.835, 161.6906, 
218.5626, 240.4953, 239.8109, 248.0714, 260.3825, 256.0689, 269.0686, 269.1015, 274.4122, 280.0392, 
15.02646, 16.88968, 17.11976, 18.04467, 19.3265, 18.94633, 20.13931, 20.34076, 20.98775, 21.63898, 
133.9614, 147.3142, 146.1641, 150.3313, 158.2776, 155.0762, 161.9477, 159.0297, 161.3286, 163.2364, 
31.10297, 34.60226, 34.48909, 35.66503, 37.35222, 36.75197, 38.38354, 37.69953, 38.50183, 39.11044, 
51.48056, 57.77391, 58.37173, 60.84568, 64.23379, 63.78598, 67.20449, 66.97814, 69.27415, 71.1111, 
10.62496, 12.03359, 12.02476, 12.55771, 13.17856, 13.06422, 13.79991, 13.74192, 14.07124, 14.38829, 
17.7365, 19.55017, 19.57361, 20.45003, 21.51948, 21.33668, 22.57982, 22.33302, 22.99301, 23.53497, 
41.81355, 46.1152, 46.15864, 47.72734, 49.8552, 49.26023, 51.80665, 52.06061, 53.43704, 54.81957, 
7.922329, 8.841202, 8.741061, 9.016768, 9.394306, 9.198502, 9.534684, 9.128947, 9.167696, 9.14244, 
264.8869, 291.6328, 290.8818, 300.4148, 314.7764, 312.2435, 328.7067, 329.3978, 337.411, 345.2845, 
6.43588, 7.257283, 7.257423, 7.573024, 7.930644, 7.847051, 8.22924, 8.173478, 8.375426, 8.532823, 
12.47516, 13.51164, 13.67121, 14.16481, 14.84484, 14.80125, 15.68975, 15.67493, 16.16302, 16.61811, 
38.55466, 42.55919, 42.73356, 44.77462, 46.89149, 46.67626, 48.96491, 48.30727, 49.45495, 50.26619, 
5.165549, 5.765423, 5.743163, 5.979708, 6.258207, 6.210151, 6.660381, 6.547146, 6.696667, 6.819485, 
18.76665, 20.43444, 20.47185, 21.09605, 22.19562, 21.92915, 23.11014, 23.05906, 23.47892, 23.98129, 
9.838238, 10.83785, 10.78763, 11.07418, 11.64549, 11.46604, 12.08742, 12.00175, 12.39663, 12.69068, 
14.30453, 15.85453, 16.0189, 16.63488, 17.36357, 17.39849, 18.35704, 18.34967, 18.87226, 19.35897, 
20.39795, 22.21271, 22.12289, 22.75408, 23.67607, 23.56273, 26.124, 26.47552, 27.2708, 28.31332, 
10.36488, 11.71599, 12.09417, 12.92527, 13.78977, 13.75626, 14.58608, 14.28558, 14.77699, 15.10677, 
62.13098, 69.43773, 70.19184, 72.84647, 76.47707, 75.42849, 79.18259, 79.55898, 81.5173, 83.49575, 
17.51529, 19.46834, 19.43428, 20.05223, 21.26592, 21.41523, 22.6589, 22.46312, 23.19864, 23.76841, 
15.89973, 16.85797, 17.00544, 17.59914, 18.50092, 18.49658, 19.38861, 19.35567, 19.87964, 20.35625, 
60.26829, 66.856, 66.74379, 68.72845, 71.79167, 71.5495, 75.14473, 74.25946, 75.812, 77.07059, 
30.01128, 33.21524, 33.46478, 34.83443, 36.98001, 37.29299, 39.37498, 38.87397, 40.03778, 40.8903, 
29.4962, 32.23379, 32.03269, 33.0141, 34.59225, 34.50268, 36.09945, 35.85596, 36.62397, 37.28642, 
39.07806, 43.08545, 42.90556, 44.69506, 47.42099, 46.62039, 48.37818, 47.44009, 48.39128, 48.90422, 
29.00427, 31.49937, 31.23754, 32.38823, 34.20381, 33.9683, 35.91977, 35.55103, 36.51612, 37.27735, 
19.92545, 22.00916, 22.19128, 23.07338, 24.07907, 24.00417, 25.44191, 25.45151, 26.11824, 26.75481, 
174.7026, 193.1358, 193.0245, 200.3059, 210.1053, 207.67, 217.5818, 215.3065, 219.1714, 222.7557, 
243.7286, 267.0859, 265.44, 274.3023, 287.5316, 283.302, 295.7671, 295.0601, 301.478, 307.4574, 
111.7262, 123.315, 123.1225, 127.2642, 132.7925, 130.4624, 136.5146, 134.1485, 135.9463, 137.436, 
38.89249, 43.20411, 43.4246, 45.23235, 47.70868, 47.30137, 49.71968, 48.81783, 49.89583, 50.72703, 
13.60943, 15.1017, 15.23089, 16.07265, 16.88501, 16.89564, 17.95465, 18.30668, 18.94278, 19.6166, 
14.08841, 15.50726, 15.41871, 15.88489, 16.62306, 16.54049, 17.44178, 17.14676, 17.56031, 17.88089, 
5.246998, 5.882542, 5.80798, 6.009941, 6.37127, 6.390492, 6.649333, 6.685965, 6.812811, 6.950722, 
45.4547, 50.73245, 51.18831, 53.51505, 56.3983, 56.37872, 59.76026, 60.88114, 63.63759, 66.06119, 
29.50638, 32.76467, 32.82683, 34.19052, 36.21935, 35.71362, 37.41943, 36.98278, 37.74334, 38.37079, 
45.04144, 49.70088, 50.01446, 52.3609, 55.26991, 54.94208, 58.23152, 58.24706, 60.07682, 61.70376, 
16.69849, 18.5631, 18.56716, 19.11343, 20.07093, 20.06876, 21.27903, 20.93902, 21.5981, 22.07382, 
49.07087, 54.90015, 55.5739, 57.83726, 60.39749, 60.01984, 63.66206, 63.51696, 65.07695, 66.56497, 
10.55323, 11.56767, 11.5246, 11.92865, 12.53544, 12.49332, 13.16783, 12.85789, 13.16342, 13.39404)
m=88
T=10
t=c(1,2,3,4,5,6,7,8,9,10)

year<-rep(1:10,len=880)
region<-rep(1:88,each=10)
region2<-region
ind2<-rep(1:880)
data<-data.frame(yL,eL,year,region,region2,ind2)

#### spatial only UH
formula1<-yL~1+f(region,model="iid")
result1<-inla(formula1,family="poisson",data=data,
E=eL,control.compute=list(dic=TRUE,cpo=TRUE))
result1$dic$dic;result1$dic$p.eff
UH<-result1$summary.random$region[,2]


#### UH and CH effects
formula2<-yL~1+f(region,model="iid")+f(region2,model="besag",graph="OHIO_map.txt")
result2<-inla(formula2,family="poisson",data=data,
E=eL,control.compute=list(dic=TRUE,cpo=TRUE))
result2$dic$dic;result2$dic$p.eff

### spatial + time trend   (model 1a)
formula2<-yL~1+f(region,model="iid")+f(region2,model="besag",graph="OHIO_map.txt")+year
result2<-inla(formula2,family="poisson",data=data,
E=eL,control.compute=list(dic=TRUE,cpo=TRUE))
result2$dic$dic;result2$dic$p.eff

UH<-result2$summary.random$region[,2]
CH<-result2$summary.random$region2[,2]
result2$summary.fixed



### UH + CH + year IID
formula3<-yL~1+f(region,model="iid")+f(region2,model="besag",graph="OHIO_map")+f(year,model="iid")
result3<-inla(formula3,family="poisson",data=data,
E=eL,control.compute=list(dic=TRUE,cpo=TRUE))
result3$dic$dic;result3$dic$p.eff

### UH + CH + year RW1   (model 1b)
formula4<-yL~1+f(region,model="iid",param=c(2,1))+f(region2,model="besag",graph="OHIO_map.txt")+f(year,model="rw1",param=c(1,0.01))
result4<-inla(formula4,family="poisson",data=data,
E=eL,control.compute=list(dic=TRUE,cpo=TRUE))
result4$dic$dic;result4$dic$p.eff
UH<-result4$summary.random$region[,2]
yearR<-result4$summary.random$year[,2]


###  UH +CH +year UH +CH    (model 2)
year2<-year
formula5a<-yL~1+f(region,model="iid")+f(region2,model="besag",graph="OHIO_map")+f(year,model="rw1")+f(year2,model="iid")
result5a<-inla(formula5a,family="poisson",data=data,
E=eL,control.compute=list(dic=TRUE,cpo=TRUE))
result5a$dic$dic;result5a$dic$p.eff


###  UH+ year RW1 +INT  IID  

formula6<-yL~1+f(region,model="iid")+f(year,model="rw1")+f(ind2,model="iid")
result6<-inla(formula6,family="poisson",data=data,
E=eL,control.compute=list(dic=TRUE,cpo=TRUE))
result6$dic$dic;result6$dic$p.eff

###  UH +CH + year RW1 + INT IID  (model 3)

formula7<-yL~1+f(region,model="iid")+f(region2,model="besag",graph="OHIO_map")+f(year,model="rw1")+f(ind2,model="iid")
result7<-inla(formula7,family="poisson",data=data,
E=eL,control.compute=list(dic=TRUE,cpo=TRUE))
result7$dic$dic;result7$dic$p.eff


##########################################################
## results for model 7
###########################################################
 UH<-result6$summary.random$region[,2]
yearR<-result6$summary.random$year[,2]
STint<-result6$summary.random$ind2[,2]
fillmap(polyOHIO,"UH",UH*100,n.col=6)
time<-seq(1:10)
plot(time,yearR)
lines(time,yearR)
STest<-matrix(0,nrow=88,ncol=10)
for(k in 1:880){
i=ceiling(k/10)
j=k-10*(i-1)
STest[i,j]<-STint[k]}
ST1<-STest[,1]
ST2<-STest[,2]
ST3<-STest[,3]
ST4<-STest[,4]
ST5<-STest[,5]
ST6<-STest[,6]
ST7<-STest[,7]
fillmap(polyOHIO,"ST interaction yr 1",ST1,n.col=6)
x11()
fillmap(polyOHIO,"ST interaction yr 2",ST2,n.col=6)