model {
for (i in 1:N) {
for (j in 1:T) {
Y[i,j] ~ dnorm(mu[i,j], tau)
mu[i,j] <- alpha[i,1] + alpha[i,2]*(x[j] - mean(x[]))
## cross-product statistics
XY[i,j] <- Y[i,j]*(x[j] - mean(x[]))
XX[i,j] <- (x[j] - mean(x[]))*(x[j] - mean(x[]))
## posterior predictive
## N(0,1) under null
resid[i,j] <- (Y[i,j] - mu[i,j])*sqrt(tau)
resid2[i,j] <- resid[i,j]*resid[i,j]
Y.fix[i,j] <- Y[i,j] # duplicate data
Y.fix[i,j] ~ dnorm(mu.fix[i,j], tau.fixed)
mu.fix[i,j] <- alpha.fix[i,1] +
alpha.fix[i,2]*(x[j] - mean(x[]))
}
X2[i] <- sum(resid2[i,]) # X2_T under null
chi.sqr[i] ~ dchisqr(T) # comparison under null
P.resid[i] <- step(X2[i] - chi.sqr[i])
# prior for intercept and gradient
alpha[i,1:2] ~ dmnorm(mu.alpha[],R[,])
# replicated intercept and gradient
alpha.pred[i,1:2] ~ dmnorm(mu.alpha[],R[,])
# summary statistics for intercept and gradient
alpha.est[i,1] <- mean(Y[i, ])
alpha.est[i,2] <- sum(XY[i,])/sum(XX[i,])
## precision of intercept estimates
alpha.est.prec[i,1] <- tau*T
## precision of gradient estimates
alpha.est.prec[i,2] <- tau*sum(XX[i,])
alpha.est.pred[i,1] ~ dnorm(alpha.pred[i,1],
alpha.est.prec[i,1])
alpha.est.pred[i,2] ~ dnorm(alpha.pred[i,2],
alpha.est.prec[i,2])
P.alpha[i,1] <- step(alpha.est[i,1] -
alpha.est.pred[i,1])
P.alpha[i,2] <- step(alpha.est[i,2] -
alpha.est.pred[i,2])
## priors for fixed effects parameters
alpha.fix[i,1] ~ dunif(-1000,1000)
alpha.fix[i,2] ~ dunif(-1000,1000)
P.alpha.fix[i,1] <- step(alpha.fix[i,1] - alpha.pred[i,1])
P.alpha.fix[i,2] <- step(alpha.fix[i,2] - alpha.pred[i,2])
}
mu.alpha[1] ~ dunif(-1000,1000)
mu.alpha[2] ~ dunif(-1000,1000)
R[1:2,1:2] ~ dwish(Omega[,],2)
Omega[1,1] <- 1;
Omega[1,2] <- 0;
Omega[2,1] <- 0;
Omega[2,2] <- 1
tau ~ dgamma(0.001,0.001)
sigma <- 1/tau
## prevent learning about tau from duplicate data
tau.fixed <- cut(tau)
}
Inits:
list(mu.alpha = c(0,0), tau=1,
alpha = structure(
.Data = c(100,6,100,6,100,6,100,6,100,6,
100,6,100,6,100,6,100,6,100,6,
100,6,100,6,100,6,100,6,100,6,
100,6,100,6,100,6,100,6,100,6,
100,6,100,6,100,6,100,6,100,6,
100,6,100,6,100,6,100,6,100,6),
.Dim = c(30, 2)),
R = structure(.Data = c(1,0,0,1), .Dim = c(2, 2)))
Data:
list(x = c(8.0, 15.0, 22.0, 29.0, 36.0), N = 30, T = 5,
Y = structure(
.Data = c(151, 199, 246, 283, 320,
145, 199, 249, 293, 354,
147, 214, 263, 312, 328,
155, 200, 237, 272, 297,
135, 188, 230, 280, 323,
159, 210, 252, 298, 331,
141, 189, 231, 275, 305,
159, 201, 248, 297, 338,
177, 236, 285, 350, 376,
134, 182, 220, 260, 296,
160, 208, 261, 313, 352,
143, 188, 220, 273, 314,
154, 200, 244, 289, 325,
171, 221, 270, 326, 358,
163, 216, 242, 281, 312,
160, 207, 248, 288, 324,
142, 187, 234, 280, 316,
156, 203, 243, 283, 317,
157, 212, 259, 307, 336,
152, 203, 246, 286, 321,
154, 205, 253, 298, 334,
139, 190, 225, 267, 302,
146, 191, 229, 272, 302,
157, 211, 250, 285, 323,
132, 185, 237, 286, 331,
160, 207, 257, 303, 345,
169, 216, 261, 295, 333,
157, 205, 248, 289, 316,
137, 180, 219, 258, 291,
153, 200, 244, 286, 324),
.Dim = c(30,5)))
node mean sd MC error 2.5% median 97.5% start sample
P.alpha.fix[1,1] 0.4279 0.49477428 0.0023297836 0.0 0.0 1.0 1001 50000
P.alpha.fix[1,2] 0.40044 0.48998756 0.0022886739 0.0 0.0 1.0 1001 50000
P.alpha.fix[2,1] 0.64108 0.47968368 0.0021076101 0.0 1.0 1.0 1001 50000
P.alpha.fix[2,2] 0.96004 0.19586526 8.9806044E-4 0.0 1.0 1.0 1001 50000
P.alpha.fix[3,1] 0.75224 0.43171169 0.0017236038 0.0 1.0 1.0 1001 50000
P.alpha.fix[3,2] 0.72838 0.44479498 0.0020950411 0.0 1.0 1.0 1001 50000
P.alpha.fix[4,1] 0.23836 0.4260804 0.0020512911 0.0 0.0 1.0 1001 50000
P.alpha.fix[4,2] 0.04228 0.20122724 9.717087E-4 0.0 0.0 1.0 1001 50000
P.alpha.fix[5,1] 0.21458 0.41053066 0.0016739128 0.0 0.0 1.0 1001 50000
P.alpha.fix[5,2] 0.78322 0.41205149 0.0018708164 0.0 1.0 1.0 1001 50000
P.alpha.fix[6,1] 0.69188 0.46171643 0.0020426788 0.0 1.0 1.0 1001 50000
P.alpha.fix[6,2] 0.49126 0.49992361 0.0023898704 0.0 0.0 1.0 1001 50000
P.alpha.fix[7,1] 0.1625 0.36890886 0.0015942637 0.0 0.0 1.0 1001 50000
P.alpha.fix[7,2] 0.336 0.47233886 0.0020338313 0.0 0.0 1.0 1001 50000
P.alpha.fix[8,1] 0.66176 0.47311067 0.0021155238 0.0 1.0 1.0 1001 50000
P.alpha.fix[8,2] 0.68388 0.46496037 0.0021042961 0.0 1.0 1.0 1001 50000
P.alpha.fix[9,1] 0.99676 0.056828711 2.6634864E-4 1.0 1.0 1.0 1001 50000
P.alpha.fix[9,2] 0.96078 0.19411798 8.9291618E-4 0.0 1.0 1.0 1001 50000
P.alpha.fix[10,1] 0.05238 0.22279214 9.8830595E-4 0.0 0.0 1.0 1001 50000
P.alpha.fix[10,2] 0.24082 0.42758125 0.0019448526 0.0 0.0 1.0 1001 50000
P.alpha.fix[11,1] 0.862 0.34489999 0.0016284339 0.0 1.0 1.0 1001 50000
P.alpha.fix[11,2] 0.89552 0.30588222 0.0013447908 0.0 1.0 1.0 1001 50000
P.alpha.fix[12,1] 0.15348 0.36044957 0.0015957074 0.0 0.0 1.0 1001 50000
P.alpha.fix[12,2] 0.4453 0.4969989 0.0023444691 0.0 0.0 1.0 1001 50000
P.alpha.fix[13,1] 0.49418 0.49996613 0.002340263 0.0 0.0 1.0 1001 50000
P.alpha.fix[13,2] 0.48562 0.49979317 0.0023214392 0.0 0.0 1.0 1001 50000
P.alpha.fix[14,1] 0.96162 0.19211188 9.034509E-4 0.0 1.0 1.0 1001 50000
P.alpha.fix[14,2] 0.8513 0.35579251 0.0016376123 0.0 1.0 1.0 1001 50000
P.alpha.fix[15,1] 0.50362 0.4999869 0.0021594365 0.0 1.0 1.0 1001 50000
P.alpha.fix[15,2] 0.05672 0.23130681 0.0010285412 0.0 0.0 1.0 1001 50000
P.alpha.fix[16,1] 0.57664 0.4940914 0.0020640545 0.0 1.0 1.0 1001 50000
P.alpha.fix[16,2] 0.2955 0.45626719 0.0022070764 0.0 0.0 1.0 1001 50000
P.alpha.fix[17,1] 0.22688 0.41881436 0.0018302974 0.0 0.0 1.0 1001 50000
P.alpha.fix[17,2] 0.57006 0.49506726 0.0024070036 0.0 1.0 1.0 1001 50000
P.alpha.fix[18,1] 0.44102 0.49650917 0.0022131282 0.0 0.0 1.0 1001 50000
P.alpha.fix[18,2] 0.24352 0.42920626 0.0020369611 0.0 0.0 1.0 1001 50000
P.alpha.fix[19,1] 0.78336 0.41195523 0.0019710635 0.0 1.0 1.0 1001 50000
P.alpha.fix[19,2] 0.67588 0.46804511 0.0020874974 0.0 1.0 1.0 1001 50000
P.alpha.fix[20,1] 0.4699 0.49909317 0.0022474291 0.0 0.0 1.0 1001 50000
P.alpha.fix[20,2] 0.39224 0.48824971 0.0022982423 0.0 0.0 1.0 1001 50000
P.alpha.fix[21,1] 0.66224 0.47294628 0.0022574809 0.0 1.0 1.0 1001 50000
P.alpha.fix[21,2] 0.67788 0.46728867 0.0018641504 0.0 1.0 1.0 1001 50000
P.alpha.fix[22,1] 0.11024 0.31318867 0.0013769229 0.0 0.0 1.0 1001 50000
P.alpha.fix[22,2] 0.24804 0.43187517 0.0020044085 0.0 0.0 1.0 1001 50000
P.alpha.fix[23,1] 0.15888 0.36556415 0.0016742547 0.0 0.0 1.0 1001 50000
P.alpha.fix[23,2] 0.1831 0.38674848 0.0018155901 0.0 0.0 1.0 1001 50000
P.alpha.fix[24,1] 0.57042 0.49501619 0.0020697697 0.0 1.0 1.0 1001 50000
P.alpha.fix[24,2] 0.2723 0.44514347 0.0021704834 0.0 0.0 1.0 1001 50000
P.alpha.fix[25,1] 0.28266 0.45029249 0.0019273382 0.0 0.0 1.0 1001 50000
P.alpha.fix[25,2] 0.93052 0.25426862 0.0011901138 0.0 1.0 1.0 1001 50000
P.alpha.fix[26,1] 0.78888 0.40810335 0.0019562184 0.0 1.0 1.0 1001 50000
P.alpha.fix[26,2] 0.7732 0.41876218 0.0017915495 0.0 1.0 1.0 1001 50000
P.alpha.fix[27,1] 0.7969 0.40230634 0.0019423966 0.0 1.0 1.0 1001 50000
P.alpha.fix[27,2] 0.27724 0.44763599 0.00205726 0.0 0.0 1.0 1001 50000
P.alpha.fix[28,1] 0.50866 0.499925 0.0022687249 0.0 1.0 1.0 1001 50000
P.alpha.fix[28,2] 0.24124 0.42783556 0.001995231 0.0 0.0 1.0 1001 50000
P.alpha.fix[29,1] 0.04266 0.2020894 8.9487395E-4 0.0 0.0 1.0 1001 50000
P.alpha.fix[29,2] 0.1489 0.35598987 0.001520934 0.0 0.0 1.0 1001 50000
P.alpha.fix[30,1] 0.46534 0.49879724 0.002311564 0.0 0.0 1.0 1001 50000
P.alpha.fix[30,2] 0.45498 0.49796907 0.0022120855 0.0 0.0 1.0 1001 50000