#######################################
# R Examples for Hierarchical Models  #
#######################################

# Math score data (from Hoff (2009))

Y <- structure(c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 
2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 
3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 
4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 5, 
5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6, 6, 
6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 
7, 7, 7, 7, 7, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 
8, 8, 8, 8, 8, 8, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 
9, 9, 9, 9, 9, 9, 9, 9, 9, 10, 10, 10, 10, 10, 10, 10, 10, 10, 
10, 10, 10, 10, 10, 10, 10, 10, 10, 11, 11, 11, 11, 11, 11, 11, 
11, 11, 11, 11, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 
12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 
12, 12, 13, 13, 13, 13, 13, 13, 13, 13, 13, 13, 13, 13, 13, 13, 
13, 13, 13, 13, 13, 14, 14, 14, 14, 14, 14, 14, 14, 14, 14, 14, 
14, 14, 14, 14, 14, 14, 14, 15, 15, 15, 15, 15, 15, 15, 15, 15, 
15, 15, 15, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 
16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 17, 17, 17, 17, 17, 17, 
17, 18, 18, 18, 18, 18, 18, 18, 18, 18, 18, 18, 18, 18, 18, 19, 
19, 19, 19, 19, 19, 19, 19, 19, 19, 19, 19, 20, 20, 20, 20, 20, 
20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 21, 
21, 21, 21, 21, 21, 21, 21, 21, 21, 21, 21, 21, 21, 21, 22, 22, 
22, 22, 22, 22, 22, 22, 22, 22, 22, 22, 22, 22, 22, 22, 22, 22, 
22, 22, 22, 22, 22, 22, 23, 23, 23, 23, 23, 23, 23, 23, 23, 23, 
23, 23, 23, 23, 24, 24, 24, 24, 24, 24, 24, 24, 24, 24, 24, 24, 
24, 24, 24, 24, 24, 24, 24, 24, 25, 25, 25, 25, 25, 25, 25, 25, 
25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 26, 26, 26, 26, 26, 
26, 26, 26, 26, 26, 26, 26, 26, 26, 26, 26, 26, 26, 26, 26, 26, 
26, 26, 26, 26, 26, 26, 26, 26, 26, 27, 27, 27, 27, 27, 27, 27, 
27, 27, 27, 27, 27, 28, 28, 28, 28, 28, 28, 28, 28, 28, 28, 28, 
28, 28, 28, 28, 28, 28, 28, 28, 28, 28, 28, 29, 29, 29, 29, 29, 
29, 29, 29, 29, 29, 29, 29, 29, 29, 29, 29, 29, 29, 29, 29, 30, 
30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 31, 
31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 
31, 31, 31, 31, 31, 32, 32, 32, 32, 32, 32, 32, 32, 32, 32, 32, 
32, 32, 32, 32, 33, 33, 33, 33, 33, 33, 33, 33, 33, 33, 33, 33, 
33, 33, 33, 33, 33, 33, 33, 33, 33, 33, 33, 33, 33, 33, 34, 34, 
34, 34, 34, 34, 34, 34, 34, 34, 34, 34, 34, 34, 34, 34, 34, 34, 
34, 34, 34, 34, 34, 34, 35, 35, 35, 35, 35, 35, 35, 35, 35, 35, 
35, 35, 35, 35, 35, 35, 35, 36, 36, 36, 36, 36, 36, 36, 36, 36, 
36, 36, 36, 37, 37, 37, 37, 37, 37, 37, 37, 37, 37, 37, 37, 37, 
37, 37, 37, 37, 37, 37, 37, 37, 37, 37, 37, 37, 37, 37, 37, 37, 
37, 37, 37, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 
38, 38, 38, 39, 39, 39, 39, 39, 39, 39, 39, 39, 39, 39, 39, 39, 
39, 39, 39, 39, 39, 39, 39, 39, 39, 39, 39, 39, 39, 39, 40, 40, 
40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 
40, 40, 40, 40, 40, 40, 40, 41, 41, 41, 41, 41, 41, 41, 41, 41, 
41, 41, 41, 41, 41, 41, 41, 41, 41, 41, 41, 41, 41, 42, 42, 42, 
42, 42, 42, 42, 42, 42, 42, 42, 42, 42, 42, 42, 42, 42, 42, 42, 
42, 42, 42, 42, 42, 42, 42, 42, 42, 43, 43, 43, 43, 43, 43, 43, 
43, 43, 43, 43, 43, 43, 43, 43, 44, 44, 44, 44, 44, 44, 44, 44, 
44, 44, 44, 44, 44, 44, 44, 44, 44, 45, 45, 45, 45, 45, 45, 45, 
45, 45, 45, 45, 45, 45, 45, 45, 46, 46, 46, 46, 46, 46, 46, 46, 
46, 46, 46, 46, 46, 46, 46, 46, 46, 46, 46, 46, 46, 47, 47, 47, 
47, 47, 47, 47, 47, 47, 47, 47, 47, 47, 47, 47, 47, 47, 47, 47, 
47, 48, 48, 48, 48, 48, 48, 48, 48, 48, 48, 48, 48, 48, 48, 48, 
48, 48, 48, 48, 48, 48, 48, 49, 49, 49, 49, 49, 49, 49, 49, 49, 
49, 49, 49, 49, 49, 49, 50, 50, 50, 50, 50, 50, 50, 50, 50, 50, 
50, 50, 50, 50, 50, 50, 50, 50, 51, 51, 51, 51, 51, 51, 51, 51, 
51, 51, 51, 51, 51, 51, 51, 51, 51, 51, 51, 52, 52, 52, 52, 52, 
52, 52, 52, 52, 52, 52, 52, 52, 52, 52, 52, 52, 52, 52, 52, 52, 
52, 52, 52, 52, 53, 53, 53, 53, 53, 53, 53, 53, 53, 53, 53, 54, 
54, 54, 54, 54, 54, 54, 54, 54, 54, 54, 54, 54, 54, 54, 54, 54, 
54, 54, 54, 54, 54, 54, 54, 54, 54, 54, 55, 55, 55, 55, 55, 55, 
55, 55, 55, 55, 55, 55, 55, 55, 55, 55, 55, 55, 55, 55, 55, 55, 
55, 55, 56, 56, 56, 56, 56, 56, 56, 56, 56, 56, 56, 56, 56, 56, 
56, 56, 56, 56, 56, 56, 56, 56, 56, 56, 57, 57, 57, 57, 57, 57, 
57, 57, 57, 57, 57, 57, 57, 57, 57, 57, 57, 57, 57, 57, 57, 57, 
57, 57, 57, 58, 58, 58, 58, 58, 58, 58, 58, 58, 58, 58, 58, 58, 
58, 58, 58, 58, 58, 58, 58, 58, 58, 58, 58, 58, 58, 58, 58, 59, 
59, 59, 59, 59, 59, 59, 59, 59, 59, 59, 59, 59, 59, 59, 59, 59, 
59, 60, 60, 60, 60, 60, 60, 60, 60, 60, 60, 60, 60, 60, 60, 60, 
60, 60, 60, 60, 60, 60, 60, 60, 60, 60, 60, 60, 60, 61, 61, 61, 
61, 61, 61, 61, 61, 61, 61, 61, 61, 61, 61, 61, 61, 61, 61, 61, 
61, 61, 61, 61, 61, 62, 62, 62, 62, 62, 62, 62, 62, 62, 62, 62, 
62, 62, 63, 63, 63, 63, 63, 63, 63, 63, 63, 63, 63, 63, 63, 63, 
63, 63, 64, 64, 64, 64, 64, 64, 64, 64, 64, 64, 64, 64, 64, 64, 
64, 64, 64, 64, 64, 64, 64, 64, 64, 64, 64, 64, 64, 65, 65, 65, 
65, 65, 65, 65, 65, 65, 65, 65, 65, 65, 65, 65, 65, 65, 66, 66, 
66, 66, 66, 66, 66, 66, 66, 66, 66, 66, 66, 66, 66, 67, 67, 67, 
67, 68, 68, 68, 68, 68, 68, 68, 68, 68, 68, 68, 68, 68, 68, 68, 
68, 68, 68, 68, 69, 69, 69, 69, 69, 69, 69, 69, 69, 69, 69, 69, 
69, 69, 69, 69, 69, 69, 69, 69, 69, 69, 69, 69, 70, 70, 70, 70, 
70, 70, 70, 70, 70, 70, 70, 70, 70, 70, 70, 70, 70, 70, 70, 70, 
70, 70, 70, 70, 70, 70, 70, 71, 71, 71, 71, 71, 71, 71, 71, 71, 
71, 71, 71, 71, 71, 71, 71, 71, 71, 71, 71, 71, 71, 72, 72, 72, 
72, 72, 72, 72, 72, 72, 72, 72, 72, 72, 73, 73, 73, 73, 73, 73, 
73, 73, 73, 73, 73, 73, 73, 73, 73, 73, 73, 73, 73, 73, 73, 73, 
73, 73, 73, 73, 73, 74, 74, 74, 74, 74, 74, 74, 74, 74, 74, 74, 
74, 74, 74, 74, 74, 74, 75, 75, 75, 75, 75, 75, 75, 75, 75, 75, 
75, 75, 75, 75, 75, 75, 75, 75, 75, 75, 75, 75, 75, 76, 76, 76, 
76, 76, 76, 76, 76, 76, 76, 76, 76, 76, 77, 77, 77, 77, 77, 77, 
77, 77, 77, 77, 77, 77, 77, 77, 77, 77, 77, 77, 77, 77, 77, 78, 
78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 79, 79, 79, 79, 
79, 79, 79, 79, 79, 79, 79, 79, 79, 80, 80, 80, 80, 80, 80, 80, 
80, 80, 80, 80, 80, 81, 81, 81, 81, 81, 81, 81, 81, 81, 81, 81, 
81, 81, 81, 81, 81, 81, 81, 81, 81, 82, 82, 82, 82, 82, 83, 83, 
83, 83, 83, 83, 83, 83, 83, 83, 83, 83, 83, 83, 84, 84, 84, 84, 
84, 84, 84, 84, 84, 84, 84, 84, 84, 84, 84, 84, 84, 84, 84, 84, 
84, 84, 84, 84, 85, 85, 85, 85, 85, 85, 85, 85, 85, 85, 85, 85, 
85, 85, 85, 86, 86, 86, 86, 86, 86, 86, 86, 86, 86, 86, 86, 86, 
86, 86, 86, 86, 87, 87, 87, 87, 87, 87, 87, 87, 87, 87, 87, 87, 
87, 87, 87, 87, 87, 87, 87, 87, 87, 87, 88, 88, 88, 88, 88, 88, 
88, 88, 88, 88, 88, 88, 88, 88, 88, 88, 88, 88, 88, 88, 88, 89, 
89, 89, 89, 89, 89, 89, 89, 89, 89, 89, 89, 89, 89, 89, 89, 89, 
89, 89, 89, 89, 89, 89, 89, 89, 89, 90, 90, 90, 90, 90, 90, 90, 
90, 90, 90, 90, 90, 90, 90, 90, 90, 90, 90, 90, 90, 90, 90, 90, 
90, 90, 90, 90, 90, 90, 90, 91, 91, 91, 91, 91, 91, 91, 91, 91, 
91, 91, 91, 91, 91, 91, 91, 91, 91, 91, 91, 92, 92, 92, 92, 92, 
92, 92, 92, 92, 92, 92, 92, 92, 92, 92, 92, 92, 92, 92, 92, 93, 
93, 93, 93, 93, 93, 93, 93, 93, 93, 93, 93, 93, 93, 93, 93, 93, 
93, 93, 93, 93, 93, 93, 93, 93, 94, 94, 94, 94, 94, 94, 94, 94, 
95, 95, 95, 95, 95, 95, 95, 95, 95, 95, 95, 95, 95, 95, 95, 95, 
95, 95, 95, 95, 95, 95, 95, 95, 95, 96, 96, 96, 96, 96, 96, 96, 
96, 96, 96, 96, 96, 96, 96, 96, 96, 96, 96, 96, 96, 96, 97, 97, 
97, 97, 97, 97, 97, 97, 97, 97, 97, 97, 97, 97, 97, 97, 97, 97, 
97, 97, 97, 97, 97, 97, 97, 97, 97, 98, 98, 98, 98, 98, 98, 98, 
98, 98, 98, 98, 98, 98, 98, 98, 98, 98, 98, 98, 98, 98, 98, 98, 
98, 98, 98, 98, 98, 98, 98, 98, 99, 99, 99, 99, 99, 99, 99, 99, 
99, 99, 99, 99, 99, 99, 99, 99, 99, 99, 99, 99, 99, 99, 99, 99, 
100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 
100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 
100, 100, 100, 100, 52.11, 57.65, 66.44, 44.68, 40.57, 35.04, 
50.71, 66.17, 39.43, 46.17, 58.76, 47.97, 39.18, 64.63, 69.38, 
32.38, 29.98, 59.32, 43.04, 57.83, 46.07, 47.74, 48.66, 40.8, 
66.32, 53.7, 52.42, 71.38, 59.66, 47.52, 39.51, 57.57, 42.4, 
41.41, 55.22, 43.9, 53.04, 49, 62.45, 53.78, 49.08, 40.25, 43.08, 
52.43, 21.73, 53.68, 41.45, 45.47, 34.06, 33.45, 60.78, 35.92, 
52.4, 40.69, 43.22, 45.57, 46.73, 49.2, 50.27, 42.79, 61.03, 
53.03, 37.98, 50.42, 54.19, 55.95, 40.53, 55.59, 69.12, 54.18, 
41.88, 47.14, 42.86, 49.12, 52.04, 38.34, 55, 54.23, 53.4, 34.34, 
48.07, 53.65, 46.6, 60.93, 56.3, 31.21, 31.56, 31.53, 59.25, 
31.31, 51.54, 57.84, 42.3, 56.57, 43.38, 33.3, 33.24, 39.66, 
49.48, 24.88, 32.72, 51.1, 48.38, 41.14, 27.89, 35.12, 28.26, 
47.54, 24.28, 46.25, 37.99, 39.1, 27.82, 31.98, 29.64, 38.47, 
39.92, 38.33, 37.38, 34.24, 39.49, 40.33, 36.46, 38.03, 47.06, 
40.75, 38.4, 52.93, 26.14, 32.95, 36.48, 44.63, 42.71, 40.72, 
36.7, 48.91, 36.11, 36.46, 32.08, 47.58, 34.62, 39.79, 48.43, 
37, 33.62, 49.87, 35.84, 46.25, 38.48, 43.39, 61.3, 58.4, 41.75, 
41.86, 55.91, 42.07, 36.3, 42.6, 42.65, 44.02, 47.64, 66.83, 
47.41, 61.74, 54.13, 43.8, 41.02, 80, 47.67, 35.73, 53.08, 56.46, 
44.66, 51.6, 46.42, 41.11, 39.67, 48.66, 48.99, 52.12, 43.04, 
42.9, 64.95, 36.89, 53.68, 55.06, 46.62, 43.12, 44.06, 57.07, 
39.24, 64.21, 48.96, 57.18, 43.19, 39.72, 47.76, 33.5, 30.87, 
34.06, 43.95, 37.08, 49.7, 49.83, 35.36, 53.48, 34.91, 43.36, 
49.23, 27.9, 56.5, 28.83, 60.86, 60.7, 67.69, 61.74, 46.47, 56.55, 
65.89, 47.89, 63.34, 44.12, 62.18, 45.18, 29.99, 63.71, 59.97, 
65.09, 64.5, 53.72, 62.08, 46.41, 62.56, 27.04, 61.06, 39.03, 
41.68, 50.54, 41.53, 35.41, 64.69, 38.45, 32.06, 41.42, 79.01, 
47.2, 41.93, 68.93, 67.64, 53.5, 34.03, 45.83, 51.62, 49.08, 
63.32, 45.7, 50.52, 54.41, 46.11, 46.81, 55.59, 41.36, 63.6, 
51.93, 35.46, 59.07, 53.82, 37.53, 48.66, 57.47, 33.41, 45.55, 
53.66, 63.85, 50.86, 65.24, 74.73, 56.98, 47.92, 52.15, 51.38, 
72.17, 45.71, 61.74, 55.51, 62.97, 63.22, 59.2, 64.09, 55.36, 
64.33, 41.46, 64.61, 61.06, 51.38, 52.42, 53.56, 45.25, 54.74, 
61.15, 72.08, 55.13, 48.18, 56.3, 60.36, 64.65, 62.42, 50.32, 
57.32, 57.96, 52.73, 54.95, 60.49, 43.01, 60.25, 50.76, 64.58, 
47.55, 53.35, 55.9, 64.38, 62.8, 60.49, 53.81, 33.85, 47.51, 
35.25, 36.51, 37.49, 41.33, 27.09, 40.31, 50, 69.18, 57.08, 58.62, 
52.35, 43.85, 38.6, 44.71, 54.67, 56.13, 37.27, 45.1, 57.79, 
41, 47.34, 29.54, 52.91, 51.36, 33.92, 36.91, 44.53, 56.28, 50.14, 
37.72, 52.18, 28.26, 61.95, 62.39, 50.35, 42.87, 33.17, 53.31, 
31.59, 50.33, 29.94, 47.93, 46.02, 51.61, 50.06, 39.9, 47.67, 
47.89, 48.23, 45.02, 36.33, 40.69, 55.84, 44.34, 41.06, 55.06, 
56.05, 55.66, 44.44, 59.63, 60.78, 42.84, 54.13, 60.07, 50.73, 
36.07, 43.85, 56.5, 60.1, 59.46, 49.95, 46.1, 64.68, 30.95, 47.95, 
40.74, 61.53, 40.16, 40.33, 46.7, 40.64, 48.21, 59.07, 53.81, 
35.9, 32.13, 51.67, 55.1, 46.63, 50.33, 31.09, 52.33, 35.94, 
48.13, 48.31, 63.37, 53.66, 65.18, 48.21, 44.47, 44.97, 68.85, 
45.95, 55.01, 53.52, 36.85, 37.33, 45.89, 47.74, 62.68, 57.83, 
33.14, 26.73, 32.41, 49.68, 50.72, 53.03, 49.08, 51.83, 48.34, 
46.94, 56.26, 38.38, 39.84, 46.22, 42.9, 48.89, 52.65, 43.33, 
26.33, 55.05, 32.43, 54.12, 41.91, 48.47, 58.24, 59.8, 35.43, 
39.51, 55.59, 38.48, 38.27, 41.68, 38.66, 51.03, 26.07, 47.08, 
43.41, 50.69, 46.07, 50.94, 39.71, 42.96, 40.73, 54.64, 46.59, 
58.86, 42.51, 39.88, 55.37, 45.78, 42.42, 54.48, 57.12, 62.45, 
34.02, 50.25, 45.87, 29.29, 32.75, 51.69, 52.64, 50.24, 52.85, 
69.02, 44.67, 39.38, 41.27, 43.9, 44.08, 23.71, 35.95, 51.23, 
28.43, 52.74, 46.96, 52.52, 55.59, 72.27, 50.47, 54.91, 61.6, 
28.34, 47.49, 66.12, 35.97, 31.54, 58.87, 50.02, 46.9, 51.75, 
55.08, 46.69, 57.64, 57.8, 51.95, 56.83, 49.12, 58.61, 50.49, 
52.13, 66.68, 42.49, 54.6, 50.42, 40.2, 50.97, 23, 35.38, 44.29, 
41.89, 52.86, 45.74, 42.89, 41.18, 30.39, 53.76, 45.65, 38.93, 
49.44, 43.03, 48.99, 64.64, 52.17, 56.17, 40.49, 55.72, 45.73, 
57.96, 42.69, 44.5, 57.67, 46.95, 48.4, 61.35, 42.11, 47.24, 
62.73, 53.7, 51.61, 59.88, 42.29, 48.44, 28.28, 33.3, 54.71, 
65.09, 50.07, 46.23, 69.53, 42.83, 45.59, 53.91, 46.81, 47.19, 
29.78, 47.43, 54.6, 46.82, 43.87, 39.77, 56.65, 54.88, 55.58, 
53.54, 47.03, 48.43, 57.62, 54.61, 51.2, 46.28, 48.56, 61.41, 
48.18, 45.4, 46.43, 42.42, 43.13, 45.04, 50.79, 40.94, 39.68, 
40.89, 56.92, 46.66, 62.13, 32.96, 64.23, 51.18, 30.22, 45.72, 
56.13, 39.2, 69.96, 40.29, 47.09, 36.36, 53.09, 60.93, 36.74, 
53.17, 55.87, 42.1, 35.48, 51.42, 42.94, 40.95, 42.25, 37.05, 
51.34, 40.63, 56.11, 38.38, 42.16, 50.08, 55.57, 37.14, 48.86, 
40.8, 50.39, 33.39, 54.16, 71.24, 44.45, 69.88, 40.11, 51.02, 
30.98, 72.24, 64.95, 57.94, 51.11, 47.75, 56.79, 61.87, 56.46, 
60.73, 58.01, 59.28, 58.58, 58.9, 52.69, 76.17, 57.19, 48.41, 
41.89, 50.8, 53.61, 42.9, 42.45, 30.35, 42.99, 63.21, 56.56, 
54.17, 44.47, 33.19, 42.14, 41.91, 42.02, 52.83, 49.6, 44.37, 
53.14, 39.47, 45.52, 49.79, 39.66, 46.24, 62.17, 49.7, 49.25, 
45.03, 59.1, 25.86, 42.98, 46.3, 51.18, 34.29, 52.13, 36.66, 
48.32, 49.9, 49.49, 42.31, 56.69, 36.17, 60.91, 55.56, 49.56, 
38.58, 70.51, 57.09, 55.89, 51.39, 53.39, 55.84, 47.93, 62.49, 
40.71, 45, 51.98, 47.92, 33.51, 36.24, 46.31, 45.34, 37.03, 44.82, 
59.26, 66.58, 34.61, 44.79, 38.2, 52.56, 41.62, 44.18, 52.65, 
32.13, 51.08, 43.51, 44.29, 49.34, 58.64, 42.78, 36.94, 46.98, 
47.72, 56.35, 54.04, 36.68, 71.04, 45.19, 33.5, 29.68, 47.56, 
47.76, 50.26, 58.06, 45.39, 43.27, 63.01, 50.36, 53.35, 49.95, 
44.5, 52.41, 52.28, 49.08, 59.25, 59.42, 44.55, 48.82, 51.96, 
49.13, 69.36, 55.68, 46.28, 58.64, 54.75, 53.94, 72.89, 70.76, 
44.91, 73.88, 60.89, 47.69, 59.76, 57.45, 62.14, 61.47, 48.15, 
52.87, 50.03, 41.51, 37.42, 64.42, 45.44, 46.06, 46.37, 46.66, 
29.01, 35.69, 49.16, 55.9, 45.84, 35.44, 43.21, 48.36, 74.14, 
46.76, 36.97, 43.84, 43.24, 56.9, 47.64, 38.84, 42.96, 41.58, 
45.96, 80.02, 62.47, 66.48, 52.83, 44.79, 48.08, 41.06, 49.99, 
33.41, 35.25, 52.75, 47.9, 47.48, 48.89, 59.58, 49.72, 49.23, 
57, 35.73, 58.75, 42.01, 28.94, 50.73, 48.34, 51.79, 50.66, 63.51, 
46.16, 58.61, 40.67, 46.2, 59.87, 45.39, 50.05, 42.21, 40.26, 
47.12, 52.32, 56.05, 34.62, 42.59, 39.96, 37.14, 55.71, 36.06, 
44.56, 53.34, 48.97, 64.48, 30.88, 24.8, 32.03, 35.77, 27.67, 
35.16, 25.84, 31.74, 56.34, 34.66, 38.08, 46.73, 46.04, 47.03, 
37.63, 45.26, 40.39, 49.96, 44.24, 59.68, 33.7, 53.87, 53, 51.18, 
50.93, 39.05, 47.81, 37.77, 46.15, 49.05, 38.3, 47.66, 46.36, 
37.35, 41.99, 31.38, 38.33, 33.31, 41.66, 43.82, 27.5, 42.43, 
49.32, 62.03, 35.68, 49.26, 48.02, 55.13, 48.68, 43.07, 32.8, 
49.26, 45.14, 45.21, 32.34, 60.05, 50.35, 42.91, 58.76, 58.37, 
49.1, 34.89, 42.81, 43, 48.71, 24.36, 35.08, 40.33, 36.05, 43.37, 
29.88, 39.48, 30.9, 41.1, 46.97, 35.51, 41.55, 29.84, 49.14, 
57.8, 44.15, 29.21, 58.24, 47.22, 37.36, 29.06, 39.84, 47.1, 
38.98, 40.91, 46.81, 35.26, 40.45, 35.52, 66.26, 66.12, 71.22, 
54.9, 61.98, 69.42, 61.22, 62.99, 57.99, 61.33, 66.85, 67.87, 
63.94, 73.7, 70.36, 64.01, 57.35, 68.25, 57.39, 66.5, 36.71, 
50.16, 50.28, 52.18, 51.46, 35.73, 53.69, 68.71, 44.75, 43.91, 
33.41, 45.72, 44.99, 47.73, 33.51, 60.51, 44.62, 53.45, 55, 31.21, 
47.96, 66.59, 61.97, 53.8, 49.44, 31.68, 63.99, 36.71, 40.77, 
52.12, 30.88, 40.01, 38.55, 42.27, 41.62, 34.24, 43.66, 37.74, 
25.47, 54.77, 52.72, 38.93, 34.49, 45.02, 64.05, 44.41, 59.1, 
66.11, 43.47, 64.42, 39.82, 32.3, 54.9, 45.63, 48.73, 29.44, 
43.16, 54.72, 56.88, 38.84, 43.39, 55.71, 34.7, 44.97, 43.9, 
40.48, 40.45, 44.63, 38.93, 37.42, 44.24, 42.12, 43.33, 53.41, 
55.7, 44.3, 51.11, 28.57, 45.45, 32.07, 43.68, 47.25, 49.88, 
47.64, 42.71, 38.03, 44.64, 67.06, 42.99, 59.1, 36.68, 52, 45.4, 
46.9, 54.94, 52.67, 41.2, 43.28, 49.12, 43.39, 54.9, 50.91, 55.56, 
55.29, 43.53, 38.92, 41.53, 70.07, 46.2, 44.7, 33.36, 35.15, 
34.73, 48.76, 51.32, 49.63, 22.61, 42.76, 54.3, 35.11, 42.16, 
44.64, 38.93, 37.08, 40.98, 29.51, 54.33, 64.63, 38.54, 67.82, 
55.08, 35.43, 41.85, 51.39, 48.03, 48.07, 41.27, 63.51, 51.78, 
49.75, 41.73, 46.13, 58.3, 39.72, 62, 47.12, 43.94, 64.65, 60.22, 
53.22, 35.86, 38.2, 43.39, 40.64, 55.21, 53.15, 46.47, 50.94, 
43.83, 60.47, 34.19, 35.1, 48.93, 45.41, 39.53, 51.06, 27.55, 
46.58, 49.02, 61.95, 30.7, 48.79, 47.7, 79.84, 59.75, 63.72, 
39.62, 38.44, 50.18, 62.13, 53.53, 42.59, 42.51, 23.36, 39.96, 
33.67, 47.56, 51.11, 30.66, 38.06, 32.37, 27.56, 38.39, 37.28, 
32.78, 41.86, 44.23, 38.77, 41.69, 83.27, 43.06, 47.27, 50.33, 
45.37, 47.38, 53.18, 51.68, 27.88, 43.16, 41.5, 56.09, 34.78, 
46.6, 50.36, 45.68, 40.92, 34.71, 42.76, 43.4, 50.04, 54.07, 
42.22, 49.36, 51.26, 45.4, 44.92, 50.74, 39.49, 33.12, 47.15, 
40.96, 55.28, 40.45, 41.88, 57.84, 37.46, 52.7, 63.67, 42.04, 
53.76, 45.25, 57.92, 59.6, 56.59, 57.74, 44.62, 45.11, 74.87, 
38.94, 52.01, 45.17, 42.74, 57.96, 30.35, 61.15, 51.59, 49.49, 
43.96, 60.98, 35.58, 60.82, 43.77, 57.6, 44.07, 28.23, 33.66, 
36.74, 39.27, 46.27, 44.11, 31.81, 37.26, 34.8, 51.8, 40.55, 
50.17, 50.97, 37.04, 48.1, 41.99, 34.94, 34.83, 59.81, 55.26, 
46.32, 61.99, 34.37, 37.51, 45.82, 32.71, 46.37, 60.96, 33.26, 
52.83, 48.55, 30.61, 27.48, 44.92, 48.1, 39.44, 63.57, 59.24, 
45.78, 73, 50.6, 56.43, 42.65, 48.71, 46.9, 35.65, 35.05, 45.17, 
59.52, 59.19, 58.84, 30.41, 51.99, 56.08, 49.96, 50.7, 52.87, 
75.62, 55.86, 66.16, 62.43, 41.58, 38.96, 56.95, 48.14, 51.62, 
39.05, 47.17, 50.11, 49, 57.17, 52.5, 37.95, 37.61, 34.05, 46.11, 
34.13, 45.14, 35, 47.95, 59.97, 58.57, 70.4, 52.84, 28.62, 69.1, 
52.88, 56.82, 39.14, 54.97, 50.73, 67.46, 50.87, 63.87, 57.6, 
47.45, 46.44, 53.28, 58.69, 52.45, 42.02, 33.19, 39.03, 38.47, 
48.23, 34.81, 55.7, 43.17, 40.6, 53.69, 38.07, 36.24, 42.16, 
54.11, 40.03, 59.01, 32.85, 47.01, 39.07, 43.82, 40.19, 40.4, 
50.35, 40.58, 47.45, 47.09, 48.25, 35.4, 34.16, 43.34, 37.92, 
51.02, 52.66, 47.09, 35.97, 36.31, 58.43, 49, 31.41, 62.18, 37.45, 
37.66, 45.94, 55.05, 37.14, 52.02, 60.94, 39.51, 52.8, 50.77, 
61.42, 36.36, 36.37, 45.44, 40.68, 36, 41.23, 37.95, 31.78, 35.52, 
30.69, 36.61, 43.73, 33.07, 31.97, 35.68, 62.09, 49.85, 59.68, 
65.1, 41.49, 50.53, 56.85, 61.21, 58.09, 61.68, 52.6, 60.88, 
56.05, 48.99, 55.85, 61.86, 38.58, 47.76, 45.49, 51.27, 58.56, 
50.59, 52.33, 51.51, 52.04, 49.87, 52.27, 46.74, 41.73, 42.34, 
34.11, 42.09, 43.17, 41.04, 37.77, 35.89, 27.1, 39.63, 28.81, 
50.26, 36.96, 43.79, 48.46, 45.37, 46.92, 46.08, 71.31, 33.09, 
68.67, 49.75, 53.83, 55.72, 45.45, 45.37, 46.23, 43.32, 59.8, 
22.58, 47.71, 45.53, 47.37, 45.21, 50.99, 53.14, 42.25, 60.76, 
42.59, 63.86, 57.05, 39.52, 60.05, 56.02, 61.22, 57.02, 54.91, 
57.76, 53.31, 46.86, 57.36, 56.55, 36.46, 48.25, 30.78, 57.02, 
58.51, 35.82, 39.13, 47.8, 44.15, 59.43, 53.94, 50.34, 49.39, 
29.01, 35.99, 35.07, 41.96, 53.31, 47.48, 32.08, 52.27, 55.08, 
42.62, 38.36, 54.73, 34.9, 30.75, 32.9, 33.32, 26.46, 34.64, 
52.21, 34.92, 39.53, 63.51, 63.4, 49.91, 53.43, 48.83, 49.93, 
62.79, 64.09, 79.49, 72.55, 84, 57.69, 52.53, 62.48, 42.74, 64.11, 
53.75, 55.42, 53.54, 40.13, 45.13, 61.12, 53.87, 47.97, 43.9, 
48.23, 37.81, 47.35, 56.62, 44.97, 25.28, 54.35, 56.86, 41.63, 
55.61, 62.35, 50.82, 57.17, 36.79, 39.84, 32.85, 61.18, 39.68, 
41.66, 46.64, 35.15, 32.17, 52.36, 48.27, 25.87, 56.17, 37.91, 
38.19, 58.9, 57.52, 52.41, 31.98, 33.64, 47.22, 46.94, 33.09, 
47.7, 36.01, 39.1, 38.26, 55.33, 37.47, 30.15, 44.88, 46.05, 
32.28, 53.84, 46.8, 60.35, 47.24, 32.95, 43.38, 58.92, 36.55, 
49.53, 27.88, 42.73, 57.23, 44.32, 43.47, 40.07, 34.82, 41.19, 
39.95, 66.62, 54.45, 30.85, 52.96, 47.43, 48.87, 58.11, 46.53, 
50.76, 40.58, 58.26, 50.4, 46.87, 40.59, 47.81, 46.43, 64.73, 
40.79, 54.47, 43.66, 54.13, 61.37, 57.48, 63.36, 44.31, 45.53, 
51.42, 60.42, 62.56, 54.03, 60.38, 40.5, 59.45, 55.11, 56.21, 
69.23, 63.83, 45.59, 68.05, 62.95, 60.81, 60.33, 50.91, 58.43, 
55.54, 68.99, 56.11, 47.85, 46.85, 69.55, 47.48, 72.19, 56.78, 
53.64, 60.68, 47.62, 35.31, 54.37, 61.85, 56.25, 47.67, 54.13, 
49.83, 58.58, 61.38, 56.11, 67.77, 72.48, 35.62, 53, 60.69, 51.47, 
43.35, 52.37, 47.32, 49.36, 68.77, 62.32, 60.13, 48.04, 35.16, 
48.05, 56.48, 52.96, 43.68, 49.61, 68.98, 55.09, 64.44, 48.16, 
56.4, 56.35, 51.24, 50.22, 68.46, 53.58, 49.65, 65.95, 34.77, 
52.78, 38.31, 49.74, 57.46, 48.22, 48.12, 59.15, 63.94, 40.1, 
43.41, 59.86, 47.28, 28.55, 49.9, 42.25, 44.14, 46.67, 46.91, 
56.72, 49.39, 49.56, 38.97, 44.27, 66.26, 42.34, 40.88, 64.58, 
45.16, 54.19, 50.73, 40.54, 60.44, 40.63, 57.74, 40.79, 43.21, 
50.13, 37.08, 48.18, 29.71, 72.43, 43.77, 49.18, 47.63, 44.94, 
59.78, 38.93, 43.63, 51.34, 51, 46.3, 56.14, 38.46, 57.21, 50.47, 
47.34, 47.54, 47.13, 66.16, 59.96, 37.97, 30.17, 54.08, 51.48, 
48.79, 47.46, 46.31, 60.17, 47.56, 26.38, 42.77, 61.1, 36.34, 
33.71, 67.09, 48.08, 44.74, 50.43, 66.19, 64.52, 49.29, 60.45, 
50.49, 44.1, 51.03, 34.54, 62.06, 62.59, 62.27, 50.11, 39.9, 
52.78, 52.19, 42.4, 58.32, 56, 36.67, 56.94, 50.32, 51.76, 38.85, 
43.8, 45.33, 49.19, 59.4, 30.76, 48.84, 46.12, 45.57, 38.25, 
31.2, 51.25, 44.4, 47.2, 60.47, 59.81, 41.22, 37.78, 46.38, 40.2, 
32.59, 46.54, 41.39, 49.8, 44.59, 31.56, 43.61, 43.66, 62.82, 
34.08, 42.45, 51.4, 49.71, 46.13, 68.22, 54.47, 29.19, 37.22, 
47.24, 45.11, 52.29, 51.76, 47.08, 41.26, 46.85, 46.12, 47.52, 
46.23, 51.02, 52.68, 32.11, 53.89, 46.11, 51.3, 45.62, 41.36, 
48.98, 57.27, 36.79, 31.97, 53.89, 50.54, 34.89, 51.57, 35, 57.81, 
58.98, 37.4, 31.91, 56.78, 48.43, 40.65, 57.64, 49.19, 36.75, 
35.97, 52.68, 49.6, 64.88, 54.14, 49.68, 57.53, 61.81, 57.62, 
27.87, 41.18, 68.95, 50.43, 65.62, 34.81, 34.45, 48.77, 59.68, 
71.65, 57.7, 59.19, 47.8, 62.73, 59.31, 36.19, 55.36, 45.34, 
60.27, 57.31, 53.32, 43.81, 52.49, 52.61, 44.37, 58.12, 50.76, 
53.41, 59.47, 47.92, 58.59, 51.18, 62.03, 52.52, 49.41, 45.14, 
56.6, 50.14, 45.59, 64.56, 53.19, 49.29, 62.04, 46.03, 48.52, 
34.53, 40.65, 49.38, 49.45, 62.16, 50.38, 39.93, 33.23, 48.59, 
62.84, 50.1, 34.55, 59.14, 50.54, 39.09, 54.6, 39.59, 47.87, 
58.69, 52.53, 50.49, 36.77, 50.11, 62.42, 50.71, 43.24, 56.25, 
38.42, 54.95, 30.29, 47.66, 55.43, 29.73), .Dim = c(1993, 2), 
.Dimnames = list(NULL, c("school", "mathscore")))

## Prior hyperparameter specification:

nu1 <- 1; nu2 <- 100;
eta1 <- 1; eta2 <- 100;
phi0 <- 50; g2 <- 25;

# Starting values for Gibbs algorithm:
# Here we let the starting values be based on some sample statistics:

m <- length(unique(Y[,1])) # m = number of schools in the study
n <- sv <- ybar <- rep(NA,m)  # setting up placeholder vectors

for (j in 1:m){
 ybar[j] <- mean(Y[Y[,1]==j,2])  # within-school sample means for schools j=1,...,m
 sv[j] <- var(Y[Y[,1]==j,2])  # within-school sample variances for schools j=1,...,m
 n[j] <- sum(Y[,1]==j)  # within-school sample sizes for schools j=1,...,m
}

mu <- ybar; sigma2 <- mean(sv)
phi <- mean(mu); tau2 <- var(mu)

## Setting up the MCMC:

#set.seed(1)

S <- 5000  # Number of MCMC iterations

MU <- matrix(nrow=S,ncol=m)  # matrix to hold the iterated values of (mu_1,..., mu_m)
SPT <- matrix(nrow=S,ncol=3)  # matrix to hold the iterated values of (sigma^2, phi, tau^2)

## Running the MCMC algorithm:

for (s in 1:S){
 
# sampling new values of the mu's from their full conditional:
 for (j in 1:m){
  vmu <- 1/(n[j]/sigma2+1/tau2)
  emu <- vmu*(ybar[j]*n[j]/sigma2+phi/tau2)
  mu[j] <- rnorm(1,emu,sqrt(vmu))
 }

# sampling new values of sigma^2 from its full conditional:
nun <- nu1+sum(n)
ss<-nu1*nu2
 for (j in 1:m){
  ss <- ss+sum((Y[Y[,1]==j,2]-mu[j])^2)
 }
sigma2 <- 1/rgamma(1,nun/2,ss/2)

# sampling new values of phi from its full conditional:
vphi<- 1/(m/tau2+1/g2)
ephi<- vphi*(m*mean(mu)/tau2 + phi0/g2)
phi <- rnorm(1,ephi,sqrt(vphi))

# sampling new values of tau^2 from its full conditional:
etam <- eta1+m
ss <- eta1*eta2 + sum((mu-phi)^2)
tau2 <- 1/rgamma(1,etam/2,ss/2)

#  Storing iterated values of (mu_1,..., mu_m):
MU[s,] <- mu

#  Storing iterated values of (sigma^2, phi, tau^2):
SPT[s,] <- c(sigma2,phi,tau2)

}
#### END Gibbs sampler.


### Diagnostics for chains of sigma^2 values, phi values, and tau^2 values:

# Trace plots:
par(mfrow=c(3,1))  # 3-by-1 plotting window:

plot(SPT[,1],type='l',xlab='Iteration',ylab='sigma^2 values')
plot(SPT[,2],type='l',xlab='Iteration',ylab='phi values')
plot(SPT[,3],type='l',xlab='Iteration',ylab='tau^2 values')

# Autocorrelation function plots:
par(mfrow=c(3,1))  # 3-by-1 plotting window:

acf(SPT[,1], main='Autocorrelations for sigma^2')
acf(SPT[,2], main='Autocorrelations for phi')
acf(SPT[,3], main='Autocorrelations for tau^2')

# Estimated posterior density plots:
par(mfrow=c(1,1))  # regular plotting window:

plot(density(SPT[,1]),xlab='sigma^2 values',ylab='Density')
windows()
plot(density(SPT[,2]),xlab='phi values',ylab='Density')
windows()
plot(density(SPT[,3]),xlab='tau^2 values',ylab='Density')


# Point Estimate of sigma, the within-class std. deviation:

sqrt(median(SPT[,1]))

# Point Estimate of tau, the std. deviation among the class mean scores:

sqrt(median(SPT[,3]))

### Examining the posterior distributions of the school means:

# posterior medians for mu_1,...,mu_100:
post.median.mus <- apply(MU,2,median)

# posterior means for mu_1,...,mu_100:
post.mean.mus <- apply(MU,2,mean)

# Plot of 95% credible interval for each mu_j:

CIs.mus <- apply(MU,2,quantile,probs=c(0.025,0.5,0.975))
matplot(matrix(1:m,nrow=3,ncol=m,byrow=T),CIs.mus,type='l',col=1:5,lty=1,xlab='School',ylab='Score', main='Pt. Estimates & 95% Cred. Intervals for School Means')
matpoints(matrix(1:m,nrow=1,ncol=m,byrow=T),t(as.matrix(post.median.mus)),col=1:5,cex=.5,pch=19)
abline(h=median(SPT[,2])) # horizontal line at point estimate for phi

## Plots to show shrinkage:

plot(ybar, post.mean.mus, xlab="Within-school sample mean", ylab="Posterior mean of Within-school population mean")
abline(0,1)
abline(h=median(SPT[,2])) # horizontal line at point estimate for phi

# Note the shrinkage:  Estimates of a school's mu_j have been pulled toward phi, slightly away from its y-bar.

## How much shrinkage?

plot(n, abs(ybar - post.mean.mus), xlab='within-school sample size', ylab='Amount of Shrinkage')

# Ranking schools based on their ordinary sample mean scores:

(1:m)[order(ybar)]

# Ranking schools based on their posterior estimated population mean scores:

(1:m)[order(post.mean.mus)]

# Within-school sample sizes:

cbind(1:m,n)

#####################################################

## Empirical Bayes Estimation of the lambda_i's in the Italian Marriage example:

x <- c(7,9,8,7,7,6,6,5,5,7,9,10,8,8,8,7) # reading in the data

xbar <- mean(x)

alpha <- 5

beta.hat <- alpha/xbar

# Consider, say, lambda_7:

# A point estimate for lambda_7 is the mean of its posterior distribution:
pt.est.lambda.7 <- (xbar/(xbar+alpha))*(x[7] + alpha)

print(pt.est.lambda.7)

# Plot of the posterior distribution for lambda_7:

my.seq <- seq(0,15,length=2000)

plot(my.seq, dgamma(my.seq, shape= x[7] + alpha, rate= 1 + beta.hat), type='l', xlab='lambda_7 values', ylab='Density')

# 95% quantile-based credible interval for lambda_7 (with posterior median printed in the middle):

qgamma(p=c(0.025, 0.5, 0.975), shape= x[7] + alpha, rate= 1 + beta.hat)

#