### R code from vignette source 'twins.rnw' mg <- read.csv("http://bendixcarstensen.com/Bayes/Cph-2012/data/mgram.csv") head( mg ) # Data as a list mg.dat <- c( mg, list( N=nrow(mg) ) ) str( mg.dat ) # The bugs-code for the model cat( "model { for (i in 1:N) { pdens1[i] ~ dnorm(mean.pdens1[i],tau.e) pdens2[i] ~ dnorm(mean.pdens2[i],tau.e) mean.pdens1[i] <- b.int + sqrt(rho)*a1[i] + sqrt(1-rho)*a2[i] mean.pdens2[i] <- b.int + sqrt(rho)*a1[i] + mz[i]*sqrt(1-rho)*a2[i] + dz[i]*sqrt(1-rho)*a3[i] a1[i] ~ dnorm(0,tau.a) a2[i] ~ dnorm(0,tau.a) a3[i] ~ dnorm(0,tau.a) } rho ~ dunif(0,1) b.int ~ dnorm(0,0.0001) tau.a <- pow(sigma.a,-2) sigma.a ~ dunif(0,1000) tau.e <- pow(sigma.e,-2) sigma.e ~ dunif(0,1000) sigma2.a <- pow(sigma.a,2) sigma2.e <- pow(sigma.e,2) }", file="mg.jag" ) # Inits as a list of lists mg.ini <- list( list( ), list( ), list( ) ) # Names of the parameters to monitir mg.par <- c("b.int","rho","sigma.a","sigma.e" ) # Model compilation and burn-in system.time( a.mod <- jags.model( file = "mg.jag", data = mg.dat, n.chains = 3, # inits = mg.ini, n.adapt = 100 ) ) # Sampling from the posterior a.res <- coda.samples( mg.mod, var = mg.par, n.iter = 0, thin = 0 )