Comparing two groups

Comparing two groups

This section is based on chapter 8 of Hoff, P. D. (2009).
We can use the following R code to approximate the posterior distribution
p(μ,δ,σ2y1,y2)p(\mu, \delta, \sigma^2\mid \mathbf y_1,\mathbf y_2)
.
1
#!/usr/bin/env Rscript
2
## R program for Comparing two groups
3
4
## get data
5
## you can found it on Peter Hoff's website
6
## https://www.stat.washington.edu/~pdhoff/book.php
7
8
Y.school.mathscore<-dget("Y.school.mathscore")
9
y.school1<-dget("y.school1")
10
y.school2<-dget("y.school2")
11
12
y1 = y.school1; n1 = length(y1)
13
y2 = y.school2; n2 = length(y2)
14
15
## prior parameters
16
mu0 = 50; g02 = 625
17
del0 = 0; t02 = 625
18
s20 = 100; nu0 = 1
19
20
## starting value
21
mu = (mean(y1) + mean(y2)) / 2
22
del = (mean(y1) - mean(y2)) / 2
23
24
## gibbs sampler
25
MU = DEL = S2 = NULL
26
set.seed(1)
27
for (s in 1:5000)
28
{
29
# update s2
30
s2 = 1/rgamma(1, (nu0+n1+n2)/2, (nu0*s20+sum((y1-mu-del)^2)+sum((y2-mu+del)^2))/2)
31
32
# update mu
33
var.mu = 1/(1/g02 + (n1+n2)/s2)
34
mean.mu = var.mu*(mu0/g02+sum(y1-del)/s2+sum(y2+del)/s2)
35
mu = rnorm(1, mean.mu, sqrt(var.mu))
36
37
# update del
38
var.del = 1/(1/t02 + (n1+n2)/s2)
39
mean.del = var.del*(del0/t02 + sum(y1-mu)/s2 - sum(y2 - mu)/s2)
40
del = rnorm(1, mean.del, sqrt(var.del))
41
42
# save parameter values
43
MU = c(MU, mu)
44
DEL = c(DEL, del)
45
S2 = c(S2, s2)
46
}
47
48
# plot
49
png("reproduce-fig-8-2l.png")
50
plot(density(MU), main="", xlab=expression(mu), lwd=2, col="black")
51
lines(density(rnorm(5000, 50, 25)), col="red", lwd=2)
52
legend("topleft", c("prior", "posterior"), col=c("black","red"), lwd=2)
53
dev.off()
54
55
png("reproduce-fig-8-2r.png")
56
plot(density(DEL), main="", xlab=expression(delta), lwd=2, col="black")
57
lines(density(rnorm(5000, 0, 25)), col="red", lwd=2)
58
legend("topleft", c("prior", "posterior"), col=c("black","red"), lwd=2)
59
dev.off()
Copied!
We can reproduce the figures as show in Fig. 8.2.

Comparing multiple groups

Heterogeneity across group means

posterior factorization
full conditional distributions
We can use the following R code to implement this Gibbs sampling.
1
## comparing multiple groups
2
Y = Y.school.mathscore
3
## weakly informative priors
4
nu0 = 1; s20 = 100
5
eta0 = 1; t20 = 100
6
mu0 = 50; g20 = 25
7
8
## starting values
9
m = length(unique(Y[, 1]))
10
n = sv = ybar = rep(NA, m)
11
12
for (j in 1:m)
13
{
14
ybar[j] = mean(Y[Y[, 1]==j, 2])
15
sv[j] = var(Y[Y[, 1]==j, 2])
16
n[j] = sum(Y[, 1]==j)
17
}
18
theta = ybar
19
sigma2 = mean(sv)
20
mu = mean(theta)
21
tau2 = var(theta)
22
23
## setup MCMC
24
set.seed(1)
25
S = 5000
26
THETA = matrix(nrow = S, ncol = m)
27
SMT = matrix(nrow = S, ncol = 3)
28
29
## MCMC algorithm
30
for (s in 1:S)
31
{
32
# sample new values of the thetas
33
for (j in 1:m)
34
{
35
vtheta = 1/(n[j]/sigma2+1/tau2)
36
etheta = vtheta*(ybar[j]*n[j]/sigma2+mu/tau2)
37
theta[j] = rnorm(1, etheta, sqrt(vtheta))
38
}
39
# sample new value of sigma2
40
nun = nu0 + sum(n)
41
ss = nu0*s20
42
for (j in 1:m)
43
{
44
ss = ss + sum(Y[Y[,1]==j, 2]-theta[j])^2
45
}
46
sigma2 = 1/rgamma(1, nun/2, ss/2)
47
48
# sample new value of mu
49
vmu = 1/(m/tau2+1/g20)
50
emu = vmu*(m*mean(theta)/tau2 + mu0/g20)
51
mu = rnorm(1, emu, sqrt(vmu))
52
53
# sample a new value of tau2
54
etam = eta0 + m
55
ss = eta0*t20 + sum((theta-mu)^2)
56
tau2 = 1/rgamma(1, etam/2, ss/2)
57
58
# store results
59
THETA[s, ] = theta
60
SMT[s, ] = c(sigma2, mu, tau2)
61
}
Copied!

Heterogeneity across group means and variances

Similarly, we can get the posterior factorization and full conditional distributions.
We can use the following Julia code to implement this Gibbs sampler.
1
## Julia program for Hierarchical Modeling of means and variances
2
## author: weiya <[email protected]>
3
## date: 27 August, 2018
4
5
using Distributions
6
using SpecialFunctions
7
using StatsBase
8
using DelimitedFiles
9
10
function higibbs(Y, T, mu0 = 50.0, gamma20 = 25.0, nu0 = 1.0, sigma20 = 100.0, eta0 = 1.0, tau20 = 100.0, a = 1.0, b = 1/100.0, alpha = 1.0, NUMAX = 40)
11
m = size(unique(Y[:,1]), 1)
12
# starting value
13
ybar = ones(m)
14
sv = ones(m)
15
n = ones(m)
16
for j = 1:m
17
yj = Y[ [Y[i,1] == j for i = 1:end], 2]
18
ybar[j] = mean(yj)
19
sv[j] = var(yj)
20
n[j] = size(yj, 1)
21
end
22
theta = ybar
23
sigma2 = copy(sv)
24
# sigma20 = 1 / mean(sigma2)
25
# nu0 = 2 * mean(sigma2)^2 / var(sigma2)
26
mu = mean(theta)
27
tau2 = var(theta)
28
29
THETA = ones(T, m)
30
SIGMA2 = ones(T, m)
31
# mu tau2 sigma20 nu0
32
MTSN = ones(T, 4)
33
34
for t = 1:T
35
# sample mu
36
varmu = 1 / (m / tau2 + 1 / gamma20)
37
meanmu = varmu * (m * mean(theta) / tau2 + mu0 / gamma20)
38
rnorm = Normal(meanmu, sqrt(varmu))
39
mu = rand(rnorm, 1)[1]
40
41
# sample 1/tau2
42
shapetau = (eta0 + m) / 2
43
ratetau = ( eta0 * tau20 + sum((theta .- mu).^2) ) / 2
44
rgamma = Gamma(shapetau, 1/ratetau)
45
tau2 = 1 / rand(rgamma, 1)[1]
46
47
# sample theta
48
for j = 1:m
49
vartheta = 1 / (n[j] / sigma2[j] + 1 / tau2)
50
meantheta = vartheta * ( n[j] * mean(Y[ [Y[i,1] == j for i = 1:end], 2]) / sigma2[j] + mu / tau2)
51
rnorm = Normal(meantheta, sqrt(vartheta))
52
theta[j] = rand(rnorm, 1)[1]
53
end
54
THETA[t, :] .= theta
55
56
# sample sigma2
57
for j = 1:m
58
shapesig = (nu0 + n[j])/2
59
yj = Y[ [Y[i,1] == j for i = 1:end], 2]
60
ratesig = ( nu0*sigma20 + sum( (yj .- theta[j]).^2 ) )/2
61
rgamma = Gamma(shapesig, 1/ratesig)
62
sigma2[j] = 1 / rand(rgamma, 1)[1]
63
end
64
SIGMA2[t, :] .= sigma2
65
66
# sample sigma20
67
shapesig = a + 0.5 * m * nu0
68
ratesig = b + 0.5 * nu0 * sum(1 ./ sigma2)
69
rgamma = Gamma(shapesig, 1/ratesig)
70
sigma20 = rand(rgamma, 1)[1]
71
72
# sample nu0
73
x = 1:NUMAX
74
lpnu0 = ones(NUMAX)
75
lpnu0 .= m * ( .5 * x .* log.(sigma20 * x / 2) .- lgamma.(x/2) ) .+ (x / 2 .+ 1) * sum(log.(1 ./ sigma2)) .- x .* (alpha + .5 * sigma20 * sum(1 ./ sigma2))
76
#println(lpnu0)
77
nu0 = sample(x, pweights(exp.(lpnu0 .- maximum(lpnu0))))
78
# println(pweights(exp.(lpnu0 .- maximum(lpnu0))))
79
80
# store results
81
MTSN[t, :] .= [mu, tau2, sigma20, nu0]
82
end
83
return THETA, SIGMA2, MTSN, sv, n
84
end
85
86
# run
87
Y = readdlm("math-score-Y.csv")
88
THETA, SIGMA2, MTSN, sv, n = higibbs(Y, 5000)
89
90
using PyPlot
91
# histogram
92
plt[:subplot](221)
93
plt[:hist](MTSN[:,2])
94
ylabel(L"$\muquot;)
95
plt[:subplot](222)
96
plt[:hist](MTSN[:,2])
97
ylabel(L"$\tau^2quot;)
98
plt[:subplot](223)
99
plt[:hist](MTSN[:,3])
100
ylabel(L"$\nu_0quot;)
101
plt[:subplot](224)
102
plt[:hist](MTSN[:,4])
103
ylabel(L"$\sigma_0^2quot;)
104
plt[:tight_layout]()
105
show()
106
107
# shrinkage
108
f, (ax1, ax2) = plt[:subplots](1, 2)
109
ax1[:scatter](sv, SIGMA2[end,:])
110
ax1[:plot](sv, sv)
111
ax1[:set_xlabel](L"$s^2quot;)
112
ax1[:set_ylabel](L"$\hat \sigma^2quot;)
113
ax2[:scatter](n, sv-SIGMA2[end,:])
114
ax2[:plot](n, zeros(size(n, 1)))
115
ax2[:set_xlabel]("sample size")
116
ax2[:set_ylabel](L"$s^2-\hat \sigma^2quot;)
117
plt[:tight_layout]()
118
show()
Copied!
And the results are:
Last modified 3yr ago