rats.data <- list(x = c(8.0, 15.0, 22.0, 29.0, 36.0), xbar = 22, 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(5,30)))
plot(rats.data$x,rats.data$Y[,1],ylim=c(140,370))
for(i in 1:30){
points(rats.data$x,rats.data$Y[,i],type="l")}
Y.t <- t(rats.data$Y)
Model with Fixed Effects Only
model.1 <- "
model
{
for( i in 1 : N ) {
for( j in 1 : T ) {
Y[i , j] ~ dnorm(mu[i , j],tau.c)
mu[i , j] <- alpha.c + beta.c * (x[j] - xbar)
}
}
tau.c ~ dgamma(1,1)
alpha.c ~ dnorm(0.0,1.0E-6)
beta.c ~ dnorm(0.0,1.0E-6)
alpha0 <- alpha.c +beta.c*(8- xbar )
sigma.c <- 1/tau.c
}"
rats.data <- list(x = c(8.0, 15.0, 22.0, 29.0, 36.0), xbar = 22, N = 30, T = 5,
Y = Y.t)
jags.1 <- jags.model(textConnection(model.1),data=rats.data, n.adapt=1500)
## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 150
## Unobserved stochastic nodes: 3
## Total graph size: 181
##
## Initializing model
test.1 <- coda.samples(jags.1, c('alpha0','beta.c','sigma.c'), n.adapt=1500, n.iter=1000)
summary(test.1)
##
## Iterations = 1:1000
## Thinning interval = 1
## Number of chains = 1
## Sample size per chain = 1000
##
## 1. Empirical mean and standard deviation for each variable,
## plus standard error of the mean:
##
## Mean SD Naive SE Time-series SE
## alpha0 156.005 2.9255 0.09251 0.092511
## beta.c 6.191 0.1821 0.00576 0.006052
## sigma.c 328.996 2120.5209 67.05676 67.056758
##
## 2. Quantiles for each variable:
##
## 2.5% 25% 50% 75% 97.5%
## alpha0 151.595 154.6 156.112 157.543 160.466
## beta.c 5.937 6.1 6.183 6.276 6.447
## sigma.c 206.381 237.8 258.282 280.474 332.643
plot(test.1)
autocorr.plot(test.1)
Model with Random Intercept Only
model.1 <- "
model
{
for( i in 1 : N ) {
for( j in 1 : T ) {
Y[i , j] ~ dnorm(mu[i , j],tau.c)
mu[i , j] <- alpha[i] + beta.c * (x[j] - xbar)
}
alpha[i] ~ dnorm(alpha.c,alpha.tau)
}
tau.c ~ dgamma(1,1)
alpha.c ~ dnorm(0.0,1.0E-6)
alpha.tau ~ dgamma(1,1)
beta.c ~ dnorm(0.0,1.0E-6)
alpha.0 <- alpha.c +beta.c*(8- xbar )
sigma.c <- 1/tau.c
sigma.alpha <- 1/alpha.tau
}"
rats.data <- list(x = c(8.0, 15.0, 22.0, 29.0, 36.0), xbar = 22, N = 30, T = 5,
Y = Y.t)
jags.1 <- jags.model(textConnection(model.1),data=rats.data, n.adapt=1500)
## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 150
## Unobserved stochastic nodes: 34
## Total graph size: 359
##
## Initializing model
test.1 <- coda.samples(jags.1, c('alpha.0','beta.c','alpha','sigma.c','sigma.alpha'), n.adapt=1500, n.iter=10000)
summary(test.1)
##
## Iterations = 1:10000
## Thinning interval = 1
## Number of chains = 1
## Sample size per chain = 10000
##
## 1. Empirical mean and standard deviation for each variable,
## plus standard error of the mean:
##
## Mean SD Naive SE Time-series SE
## alpha[1] 239.961 3.6612 0.036612 0.0389688
## alpha[2] 247.604 3.6560 0.036560 0.0393994
## alpha[3] 252.068 3.6115 0.036115 0.0385306
## alpha[4] 232.931 3.6244 0.036244 0.0381590
## alpha[5] 231.948 3.5598 0.035598 0.0361369
## alpha[6] 249.507 3.6677 0.036677 0.0409042
## alpha[7] 229.192 3.5818 0.035818 0.0369421
## alpha[8] 248.227 3.6135 0.036135 0.0384150
## alpha[9] 281.840 3.8366 0.038366 0.0488233
## alpha[10] 220.089 3.6260 0.036260 0.0382045
## alpha[11] 257.718 3.6296 0.036296 0.0424008
## alpha[12] 228.763 3.5947 0.035947 0.0373002
## alpha[13] 242.420 3.6066 0.036066 0.0390420
## alpha[14] 267.357 3.7402 0.037402 0.0443182
## alpha[15] 242.711 3.6165 0.036165 0.0373019
## alpha[16] 245.217 3.5871 0.035871 0.0369102
## alpha[17] 232.580 3.6348 0.036348 0.0384003
## alpha[18] 240.594 3.6154 0.036154 0.0367471
## alpha[19] 253.340 3.6521 0.036521 0.0410524
## alpha[20] 241.651 3.6133 0.036133 0.0379880
## alpha[21] 248.355 3.6248 0.036248 0.0395676
## alpha[22] 225.905 3.5655 0.035655 0.0363322
## alpha[23] 228.992 3.5922 0.035922 0.0366401
## alpha[24] 244.986 3.5901 0.035901 0.0369641
## alpha[25] 234.769 3.6564 0.036564 0.0386407
## alpha[26] 253.526 3.6265 0.036265 0.0391739
## alpha[27] 253.973 3.6401 0.036401 0.0377402
## alpha[28] 242.919 3.6281 0.036281 0.0389186
## alpha[29] 218.822 3.6158 0.036158 0.0377267
## alpha[30] 241.502 3.5722 0.035722 0.0387127
## alpha.0 156.117 2.8620 0.028620 0.0332475
## beta.c 6.184 0.0684 0.000684 0.0006998
## sigma.alpha 194.524 56.9183 0.569183 0.6613603
## sigma.c 68.400 65.2107 0.652107 1.0162861
##
## 2. Quantiles for each variable:
##
## 2.5% 25% 50% 75% 97.5%
## alpha[1] 232.863 237.538 239.973 242.39 247.087
## alpha[2] 240.606 245.250 247.596 250.01 254.695
## alpha[3] 245.121 249.719 252.042 254.45 259.058
## alpha[4] 225.927 230.488 232.964 235.32 239.899
## alpha[5] 225.018 229.556 232.021 234.31 238.721
## alpha[6] 242.534 247.153 249.529 251.89 256.435
## alpha[7] 222.247 226.809 229.168 231.54 236.345
## alpha[8] 241.245 245.816 248.231 250.61 255.261
## alpha[9] 274.570 279.443 281.887 284.33 288.974
## alpha[10] 213.062 217.652 220.102 222.50 227.223
## alpha[11] 250.776 255.360 257.741 260.10 264.826
## alpha[12] 221.747 226.402 228.780 231.15 235.676
## alpha[13] 235.440 240.043 242.429 244.77 249.342
## alpha[14] 260.189 264.993 267.414 269.81 274.363
## alpha[15] 235.757 240.324 242.716 245.16 249.571
## alpha[16] 238.282 242.869 245.210 247.60 252.194
## alpha[17] 225.723 230.108 232.548 235.01 239.668
## alpha[18] 233.515 238.258 240.628 242.97 247.646
## alpha[19] 246.339 250.992 253.382 255.76 260.314
## alpha[20] 234.773 239.261 241.676 244.06 248.610
## alpha[21] 241.409 245.984 248.400 250.76 255.214
## alpha[22] 218.943 223.541 225.917 228.30 232.870
## alpha[23] 222.098 226.625 229.043 231.37 235.908
## alpha[24] 238.011 242.648 244.987 247.33 252.035
## alpha[25] 227.787 232.339 234.761 237.17 241.950
## alpha[26] 246.686 251.153 253.507 255.93 260.638
## alpha[27] 246.979 251.560 254.029 256.32 261.091
## alpha[28] 236.041 240.478 242.936 245.33 249.895
## alpha[29] 211.932 216.395 218.814 221.23 225.975
## alpha[30] 234.622 239.180 241.538 243.83 248.348
## alpha.0 150.595 154.263 156.135 157.97 161.602
## beta.c 6.047 6.138 6.183 6.23 6.317
## sigma.alpha 110.676 154.056 185.459 224.59 332.566
## sigma.c 52.532 61.189 66.570 72.86 87.246
plot(test.1)
autocorr.plot(test.1)
summary(test.1[,31:34])
##
## Iterations = 1:10000
## Thinning interval = 1
## Number of chains = 1
## Sample size per chain = 10000
##
## 1. Empirical mean and standard deviation for each variable,
## plus standard error of the mean:
##
## Mean SD Naive SE Time-series SE
## alpha.0 156.117 2.8620 0.028620 0.0332475
## beta.c 6.184 0.0684 0.000684 0.0006998
## sigma.alpha 194.524 56.9183 0.569183 0.6613603
## sigma.c 68.400 65.2107 0.652107 1.0162861
##
## 2. Quantiles for each variable:
##
## 2.5% 25% 50% 75% 97.5%
## alpha.0 150.595 154.263 156.135 157.97 161.602
## beta.c 6.047 6.138 6.183 6.23 6.317
## sigma.alpha 110.676 154.056 185.459 224.59 332.566
## sigma.c 52.532 61.189 66.570 72.86 87.246
out <- as.matrix(test.1)
boxplot(out[,1:30], las=2)
jags.1 <- jags.model(textConnection(model.1),data=rats.data, n.adapt=1500, n.chains=2)
## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 150
## Unobserved stochastic nodes: 34
## Total graph size: 359
##
## Initializing model
dic.1 <- dic.samples(jags.1,n.iter=1000)
Model with Random Intercept and Slope
model.1 <- "
model
{
for( i in 1 : N ) {
for( j in 1 : T ) {
Y[i , j] ~ dnorm(mu[i , j],tau.c)
mu[i , j] <- alpha[i] + beta[i] * (x[j] - xbar)
}
alpha[i] ~ dnorm(alpha.c,alpha.tau)
beta[i] ~ dnorm(beta.c,beta.tau)
}
tau.c ~ dgamma(1,1)
alpha.c ~ dnorm(0.0,1.0E-6)
alpha.tau ~ dgamma(1,1)
beta.c ~ dnorm(0.0,1.0E-6)
beta.tau ~ dgamma(1,1)
alpha.0 <- alpha.c +beta.c*(8- xbar )
sigma.c <- 1/tau.c
sigma.alpha <- 1/alpha.tau
sigma.beta <- 1/beta.tau
}"
jags.1 <- jags.model(textConnection(model.1),data=rats.data, n.adapt=1500)
## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 150
## Unobserved stochastic nodes: 65
## Total graph size: 537
##
## Initializing model
test.1 <- coda.samples(jags.1, c('alpha.0','beta.c','alpha','beta','sigma.c','sigma.alpha','sigma.beta'), n.adapt=1500, n.iter=10000)
summary(test.1)
##
## Iterations = 1:10000
## Thinning interval = 1
## Number of chains = 1
## Sample size per chain = 10000
##
## 1. Empirical mean and standard deviation for each variable,
## plus standard error of the mean:
##
## Mean SD Naive SE Time-series SE
## alpha[1] 239.8866 2.6875 0.026875 0.028058
## alpha[2] 247.7837 2.6818 0.026818 0.030910
## alpha[3] 252.3672 2.7051 0.027051 0.030225
## alpha[4] 232.5256 2.7065 0.027065 0.028179
## alpha[5] 231.6609 2.6723 0.026723 0.026723
## alpha[6] 249.7254 2.7024 0.027024 0.029075
## alpha[7] 228.7422 2.6799 0.026799 0.027359
## alpha[8] 248.3846 2.7146 0.027146 0.029563
## alpha[9] 283.2671 2.8915 0.028915 0.036844
## alpha[10] 219.2992 2.6610 0.026610 0.027216
## alpha[11] 258.2336 2.7376 0.027376 0.030502
## alpha[12] 228.1215 2.6656 0.026656 0.026370
## alpha[13] 242.4024 2.6994 0.026994 0.028920
## alpha[14] 268.1995 2.7536 0.027536 0.032301
## alpha[15] 242.7884 2.6820 0.026820 0.028889
## alpha[16] 245.2910 2.6524 0.026524 0.028613
## alpha[17] 232.1810 2.7062 0.027062 0.027920
## alpha[18] 240.4636 2.6973 0.026973 0.028792
## alpha[19] 253.7723 2.7144 0.027144 0.029840
## alpha[20] 241.6779 2.7030 0.027030 0.028297
## alpha[21] 248.5555 2.7395 0.027395 0.028475
## alpha[22] 225.2230 2.6854 0.026854 0.027524
## alpha[23] 228.5783 2.6995 0.026995 0.029292
## alpha[24] 245.0907 2.6803 0.026803 0.027628
## alpha[25] 234.5021 2.6714 0.026714 0.026714
## alpha[26] 253.9892 2.7474 0.027474 0.032002
## alpha[27] 254.3537 2.7531 0.027531 0.030069
## alpha[28] 243.0199 2.7285 0.027285 0.026647
## alpha[29] 217.9602 2.6689 0.026689 0.028009
## alpha[30] 241.4257 2.6686 0.026686 0.026686
## alpha.0 156.1013 3.1464 0.031464 0.034229
## beta[1] 6.0534 0.2471 0.002471 0.002471
## beta[2] 7.1063 0.2509 0.002509 0.002763
## beta[3] 6.5008 0.2447 0.002447 0.002486
## beta[4] 5.2859 0.2568 0.002568 0.002692
## beta[5] 6.5971 0.2467 0.002467 0.002559
## beta[6] 6.1730 0.2438 0.002438 0.002438
## beta[7] 5.9664 0.2429 0.002429 0.002429
## beta[8] 6.4332 0.2501 0.002501 0.002510
## beta[9] 7.1063 0.2545 0.002545 0.002708
## beta[10] 5.8255 0.2457 0.002457 0.002457
## beta[11] 6.8365 0.2493 0.002493 0.002606
## beta[12] 6.1195 0.2426 0.002426 0.002426
## beta[13] 6.1605 0.2463 0.002463 0.002463
## beta[14] 6.7220 0.2493 0.002493 0.002645
## beta[15] 5.3672 0.2510 0.002510 0.002616
## beta[16] 5.9044 0.2460 0.002460 0.002460
## beta[17] 6.2768 0.2470 0.002470 0.002470
## beta[18] 5.8193 0.2462 0.002462 0.002462
## beta[19] 6.4202 0.2457 0.002457 0.002457
## beta[20] 6.0406 0.2465 0.002465 0.002465
## beta[21] 6.4183 0.2450 0.002450 0.002450
## beta[22] 5.8312 0.2471 0.002471 0.002471
## beta[23] 5.7188 0.2473 0.002473 0.002546
## beta[24] 5.8692 0.2454 0.002454 0.002454
## beta[25] 6.9546 0.2498 0.002498 0.002593
## beta[26] 6.5732 0.2440 0.002440 0.002440
## beta[27] 5.8759 0.2456 0.002456 0.002518
## beta[28] 5.8190 0.2460 0.002460 0.002460
## beta[29] 5.6367 0.2484 0.002484 0.002569
## beta[30] 6.1230 0.2473 0.002473 0.002519
## beta.c 6.1843 0.1203 0.001203 0.001348
## sigma.alpha 203.1896 57.4159 0.574159 0.634742
## sigma.beta 0.3552 0.1121 0.001121 0.001493
## sigma.c 36.3713 22.8537 0.228537 0.410276
##
## 2. Quantiles for each variable:
##
## 2.5% 25% 50% 75% 97.5%
## alpha[1] 234.6347 238.1076 239.8816 241.6736 245.1479
## alpha[2] 242.6810 246.0131 247.7714 249.5333 252.9557
## alpha[3] 247.2100 250.6356 252.3735 254.1506 257.5855
## alpha[4] 227.2313 230.8042 232.5432 234.3098 237.7223
## alpha[5] 226.5473 229.9038 231.6631 233.4364 236.8658
## alpha[6] 244.6322 247.9754 249.7451 251.4778 254.9018
## alpha[7] 223.5822 226.9360 228.7678 230.5130 233.9122
## alpha[8] 243.0937 246.6370 248.3880 250.1253 253.6200
## alpha[9] 277.9569 281.5317 283.3546 285.1175 288.3453
## alpha[10] 214.1414 217.5431 219.2844 221.0398 224.5162
## alpha[11] 252.9425 256.4861 258.2681 260.0044 263.4716
## alpha[12] 222.9748 226.3641 228.1094 229.8749 233.3113
## alpha[13] 237.2123 240.6345 242.4421 244.1392 247.6284
## alpha[14] 262.9217 266.4717 268.2221 269.9641 273.2944
## alpha[15] 237.6208 241.0565 242.8004 244.5587 248.0070
## alpha[16] 240.2629 243.5614 245.2911 247.0634 250.3627
## alpha[17] 226.9913 230.3781 232.2076 233.9835 237.4504
## alpha[18] 235.3239 238.7002 240.4946 242.2248 245.6264
## alpha[19] 248.5562 252.0516 253.7918 255.5799 258.8685
## alpha[20] 236.5701 239.9194 241.6843 243.4551 246.7753
## alpha[21] 243.4290 246.7619 248.5943 250.3626 253.7396
## alpha[22] 220.1420 223.4251 225.1937 227.0424 230.4085
## alpha[23] 223.4440 226.7896 228.5346 230.3608 233.8842
## alpha[24] 239.9197 243.3475 245.0891 246.8784 250.2287
## alpha[25] 229.3097 232.7362 234.5009 236.2538 239.7416
## alpha[26] 248.7932 252.2191 253.9794 255.7931 259.1378
## alpha[27] 249.0730 252.5731 254.3651 256.1707 259.6461
## alpha[28] 237.8499 241.2614 243.0299 244.7842 248.2807
## alpha[29] 212.7666 216.2048 217.9556 219.7217 223.1279
## alpha[30] 236.3662 239.6315 241.4019 243.2235 246.5962
## alpha.0 150.0088 154.0286 156.1048 158.1517 162.2941
## beta[1] 5.5755 5.8840 6.0501 6.2187 6.5436
## beta[2] 6.6084 6.9406 7.1031 7.2755 7.5916
## beta[3] 6.0192 6.3360 6.5030 6.6618 6.9759
## beta[4] 4.7879 5.1107 5.2877 5.4595 5.7893
## beta[5] 6.1165 6.4313 6.5978 6.7620 7.0792
## beta[6] 5.6890 6.0090 6.1740 6.3357 6.6520
## beta[7] 5.4878 5.8047 5.9648 6.1293 6.4481
## beta[8] 5.9401 6.2687 6.4340 6.5971 6.9271
## beta[9] 6.6064 6.9400 7.1081 7.2750 7.6117
## beta[10] 5.3415 5.6600 5.8267 5.9901 6.3090
## beta[11] 6.3463 6.6697 6.8366 7.0024 7.3264
## beta[12] 5.6430 5.9567 6.1191 6.2807 6.5968
## beta[13] 5.6809 5.9945 6.1603 6.3274 6.6407
## beta[14] 6.2300 6.5575 6.7248 6.8905 7.2099
## beta[15] 4.8760 5.2005 5.3636 5.5324 5.8603
## beta[16] 5.4140 5.7383 5.9058 6.0672 6.3905
## beta[17] 5.7854 6.1120 6.2791 6.4385 6.7677
## beta[18] 5.3283 5.6538 5.8236 5.9848 6.2956
## beta[19] 5.9417 6.2545 6.4199 6.5881 6.8964
## beta[20] 5.5635 5.8789 6.0407 6.2051 6.5223
## beta[21] 5.9393 6.2548 6.4180 6.5815 6.9019
## beta[22] 5.3477 5.6647 5.8322 5.9975 6.3135
## beta[23] 5.2360 5.5542 5.7218 5.8777 6.2109
## beta[24] 5.3861 5.7037 5.8679 6.0323 6.3563
## beta[25] 6.4695 6.7863 6.9538 7.1214 7.4487
## beta[26] 6.0998 6.4087 6.5767 6.7396 7.0476
## beta[27] 5.3947 5.7112 5.8764 6.0387 6.3533
## beta[28] 5.3336 5.6558 5.8192 5.9829 6.3004
## beta[29] 5.1554 5.4697 5.6357 5.8024 6.1295
## beta[30] 5.6449 5.9560 6.1222 6.2890 6.6128
## beta.c 5.9481 6.1045 6.1847 6.2646 6.4211
## sigma.alpha 118.0646 162.6187 194.3456 233.9347 342.4947
## sigma.beta 0.1907 0.2751 0.3376 0.4147 0.6224
## sigma.c 27.0167 32.1527 35.4453 39.2058 47.5467
plot(test.1)
autocorr.plot(test.1)
summary(test.1[,c(31, 62:65)])
##
## Iterations = 1:10000
## Thinning interval = 1
## Number of chains = 1
## Sample size per chain = 10000
##
## 1. Empirical mean and standard deviation for each variable,
## plus standard error of the mean:
##
## Mean SD Naive SE Time-series SE
## alpha.0 156.1013 3.1464 0.031464 0.034229
## beta.c 6.1843 0.1203 0.001203 0.001348
## sigma.alpha 203.1896 57.4159 0.574159 0.634742
## sigma.beta 0.3552 0.1121 0.001121 0.001493
## sigma.c 36.3713 22.8537 0.228537 0.410276
##
## 2. Quantiles for each variable:
##
## 2.5% 25% 50% 75% 97.5%
## alpha.0 150.0088 154.0286 156.1048 158.1517 162.2941
## beta.c 5.9481 6.1045 6.1847 6.2646 6.4211
## sigma.alpha 118.0646 162.6187 194.3456 233.9347 342.4947
## sigma.beta 0.1907 0.2751 0.3376 0.4147 0.6224
## sigma.c 27.0167 32.1527 35.4453 39.2058 47.5467
out <- as.matrix(test.1)
boxplot(out[,1:30], las=2)
boxplot(out[,32:61], las=2)
jags.1 <- jags.model(textConnection(model.1),data=rats.data, n.adapt=1500, n.chains=2)
## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 150
## Unobserved stochastic nodes: 65
## Total graph size: 537
##
## Initializing model
dic.2 <- dic.samples(jags.1,n.iter=1000)
diffdic(dic.1,dic.2)
## Difference: 68.93575
## Sample standard error: 12.27836