![[datacloning0]](datacloning0.bmp)
Using Data Cloning to Calculate MLEs for the Seeds Model
model
{
for( i in 1 : N ) {
for(k in 1 : K){ # replicate data and random effects
r.rep[i, k] <- r[i]
r.rep[i, k] ~ dbin(p[i, k],n[i])
b[i, k] ~ dnorm(0.0,tau)
logit(p[i, k]) <- alpha0 + alpha1 * x1[i] + alpha2 * x2[i] +
alpha12 * x1[i] * x2[i] + b[i, k]
}
}
alpha0 ~ dnorm(0.0,1.0E-6)
alpha1 ~ dnorm(0.0,1.0E-6)
alpha2 ~ dnorm(0.0,1.0E-6)
alpha12 ~ dnorm(0.0,1.0E-6)
tau ~ dgamma(0.001,0.001)
sigma <- 1 / sqrt(tau)
}
Run simulation for K = {1, 2, 4, 8, 16, 32, 64, 128, 256}.. The MLE point estimates are given by the mean of the MCMC simulation and the MLE SE by the MCMC sd scaled by K 1/2 .Note
that Monte Carlo errors will be magnified by this scaling so accurate estimates of the MLE SE will need both large K and long MCMC runs.
Data ( click to open )
Inits for chain 1
mean sd sd*K^0.5 MC_error val2.5pc median val97.5pc start sample
K = 1
alpha0 -0.5485 0.196 0.196 0.007577 -0.9406 -0.549 -0.1511 1001 10000
alpha1 0.07439 0.3178 0.3178 0.01159 -0.5561 0.08066 0.6841 1001 10000
alpha12 -0.818 0.4501 0.4501 0.01764 -1.734 -0.8121 0.06914 1001 10000
alpha2 1.35 0.29 0.29 0.01136 0.7958 1.345 1.942 1001 10000
sigma 0.2971 0.1393 0.1393 0.006266 0.06241 0.2858 0.6034 1001 10000
deviance 101.5 6.721 0.2203 89.86 101.0 115.5 1001 10000
K = 2
alpha0 -0.5456 0.1258 0.1779 0.005102 -0.7946 -0.5452 -0.2963 1001 10000
alpha1 0.08221 0.21 0.2970 0.008157 -0.3328 0.08527 0.4923 1001 10000
alpha12 -0.7978 0.2933 0.4148 0.01172 -1.405 -0.7924 -0.2309 1001 10000
alpha2 1.335 0.1805 0.2553 0.007379 0.9932 1.332 1.71 1001 10000
sigma 0.2504 0.09355 0.1323 0.005168 0.06468 0.2501 0.4359 1001 10000
deviance 202.2 9.675 0.4351 185.4 201.5 222.4 1001 10000
K = 4
alpha0 -0.5486 0.08535 0.1707 0.002859 -0.7156 -0.5473 -0.3819 1001 10000
alpha1 0.09523 0.1402 0.2804 0.004785 -0.192 0.09921 0.3636 1001 10000
alpha12 -0.8117 0.1906 0.3812 0.006673 -1.171 -0.8153 -0.4333 1001 10000
alpha2 1.338 0.1186 0.2372 0.004158 1.112 1.339 1.575 1001 10000
sigma 0.234 0.06998 0.1400 0.004371 0.07462 0.2381 0.365 1001 10000
deviance 402.5 14.51 0.7888 378.0 401.1 436.5 1001 10000
K = 8
alpha0 -0.5516 0.06062 0.1715 0.002096 -0.6734 -0.5513 -0.4327 1001 10000
alpha1 0.1017 0.1008 0.2851 0.003262 -0.09649 0.1012 0.2982 1001 10000
alpha12 -0.8247 0.1377 0.3895 0.004714 -1.095 -0.823 -0.5542 1001 10000
alpha2 1.345 0.08519 0.2410 0.002913 1.184 1.343 1.513 1001 10000
sigma 0.2407 0.04222 0.1194 0.0018 0.1558 0.2403 0.3237 1001 10000
deviance 798.1 17.92 0.5668 765.3 797.3 834.7 1001 10000
K = 16
alpha0 -0.5493 0.04149 0.1660 0.001354 -0.6305 -0.5499 -0.4681 1001 10000
alpha1 0.09651 0.06921 0.2768 0.002323 -0.03731 0.09542 0.2342 1001 10000
alpha12 -0.8089 0.09511 0.3804 0.003167 -0.9984 -0.8081 -0.6271 1001 10000
alpha2 1.337 0.05826 0.2330 0.001839 1.224 1.336 1.452 1001 10000
sigma 0.2376 0.02841 0.1136 0.001046 0.1812 0.2378 0.2941 1001 10000
deviance 1593.0 24.57 0.6968 1546.0 1593.0 1643.0 1001 10000
K = 32
alpha0 -0.5503 0.03028 0.1713 0.001039 -0.6097 -0.5509 -0.4898 1001 10000
alpha1 0.09915 0.0504 0.2851 0.00159 9.589E-4 0.09946 0.1971 1001 10000
alpha12 -0.815 0.06983 0.3953 0.002325 -0.9544 -0.8152 -0.678 1001 10000
alpha2 1.34 0.04268 0.2414 0.001568 1.256 1.34 1.423 1001 10000
sigma 0.2371 0.01994 0.1128 7.171E-4 0.196 0.2371 0.2751 1001 10000
deviance 3184.0 34.6 0.9479 3118.0 3183.0 3254.0 1001 10000
K = 64
alpha0 -0.5488 0.02058 0.1646 6.559E-4 -0.5891 -0.5487 -0.5088 1001 10000
alpha1 0.0981 0.03465 0.2768 0.001063 0.0299 0.09791 0.1665 1001 10000
alpha12 -0.8128 0.04778 0.3829 0.001488 -0.9065 -0.8121 -0.7169 1001 10000
alpha2 1.338 0.02904 0.2323 9.41E-4 1.281 1.338 1.396 1001 10000
sigma 0.2372 0.01438 0.1150 6.108E-4 0.2088 0.2374 0.2648 1001 10000
deviance 6362.0 48.78 1.451 6270.0 6361.0 6461.0 1001 10000
K = 128
alpha0 -0.5483 0.01477 0.1671 5.063E-4 -0.5777 -0.548 -0.5198 1001 10000
alpha1 0.09684 0.02418 0.2736 8.603E-4 0.04867 0.09703 0.1443 1001 10000
alpha12 -0.8099 0.03333 0.3771 0.001116 -0.8752 -0.8102 -0.7442 1001 10000
alpha2 1.337 0.02104 0.2380 7.183E-4 1.296 1.336 1.378 1001 10000
sigma 0.2363 0.009851 0.1115 4.292E-4 0.2167 0.2363 0.2552 1001 10000
deviance 12720.0 67.86 2.166 12590.0 12720.0 12860.0 1001 10000
K = 256
alpha0 -0.5489 0.01018 0.1629 3.314E-4 -0.5692 -0.5486 -0.5293 1001 10000
alpha1 0.0972 0.01707 0.2731 5.336E-4 0.06356 0.09732 0.1304 1001 10000
alpha12 -0.8107 0.0239 0.3824 7.035E-4 -0.8577 -0.8108 -0.7642 1001 10000
alpha2 1.338 0.01457 0.2331 4.451E-4 1.31 1.337 1.367 1001 10000
sigma 0.2363 0.006876 0.1100 2.863E-4 0.2233 0.2362 0.25 1001 10000
deviance 25450.0 94.62 2.687 25260.0 25450.0 25630.0 1001 10000
We may compare simple logistic, maximum likelihood (from EGRET), penalized quasi-likelihood (PQL) Breslow and Clayton (1993) with the BUGS results
![[datacloning1]](datacloning1.bmp)