### R code from vignette source 'fetal-lji.rnw' ################################################### ### code chunk number 1: read ################################################### # fetal <- read.csv("http://BendixCarstensen.com/Bayes/Cph-2012/data/fetal.csv",header=TRUE) fetal <- read.csv("../data/fetal.csv",header=TRUE) head( fetal, 11 ) ################################################### ### code chunk number 2: fetal-lji.rnw:20-21 ################################################### fetal$tga <- fetal$tga - 19.5 ################################################### ### code chunk number 3: lmer0 ################################################### library( lme4 ) m0 <- lmer( hc ~ tga + (tga|id), data=fetal ) summary(m0) ################################################### ### code chunk number 4: fetal-lji.rnw:46-51 ################################################### fetal.dat <- list( id = as.integer( factor(fetal$id) ), hc = fetal$hc, tga = fetal$tga, N = nrow(fetal), I = length( unique(fetal$id) ) ) ################################################### ### code chunk number 5: jags-mod ################################################### cat(" # Fixing data to be used in model definition data { zero[1] <- 0 zero[2] <- 0 R[1,1] <- 0.1 R[1,2] <- 0 R[2,1] <- 0 R[2,2] <- 0.5 } # Then define model model { # Intercept and slope for each person, including random effects for( i in 1:I ) { u[i,1:2] ~ dmnorm(zero,Omega.u) } # Define model for each observational unit for( j in 1:N ) { mu[j] <- ( beta[1] + u[id[j],1] ) + ( beta[2] + u[id[j],2] ) * ( tga[j] ) hc[j] ~ dnorm( mu[j], tau.e ) } #------------------------------------------------------------ # Priors: # Fixed intercept and slope beta[1] ~ dnorm(0.0,1.0E-5) beta[2] ~ dnorm(0.0,1.0E-5) # Residual variance tau.e <- pow(sigma.e,-2) sigma.e ~ dunif(0,100) # Define prior for the variance-covariance matrix of the random effects Sigma.u <- inverse(Omega.u) Omega.u ~ dwish( R, 2 ) }", file="fetal.jag" ) ################################################### ### code chunk number 6: inits ################################################### ( sigma.e <- attr(VarCorr(m0),"sc") ) ( Omega.u <- solve( VarCorr(m0)$id ) ) ( beta <- fixef( m0 ) ) fetal.ini <- list( list( sigma.e = sigma.e/3, Omega.u = Omega.u/3, beta = beta /3 ), list( sigma.e = sigma.e*3, Omega.u = Omega.u*3, beta = beta *3 ), list( sigma.e = sigma.e/3, Omega.u = Omega.u*3, beta = beta /3 ), list( sigma.e = sigma.e*3, Omega.u = Omega.u/3, beta = beta *3 ) ) ################################################### ### code chunk number 7: run-jags ################################################### library( rjags ) system.time( fetal.mod <- jags.model( file = "fetal.jag", data = fetal.dat, n.chains = 4, inits = fetal.ini, n.adapt = 500 ) ) system.time( fetal.res <- coda.samples( fetal.mod, var = c("beta","sigma.e","Sigma.u"), n.iter = 500, thin = 5 ) ) str( fetal.res ) summary( fetal.res ) ################################################### ### code chunk number 8: fetal-lji.rnw:145-147 ################################################### fetal.mat <- as.matrix( fetal.res ) colnames( fetal.mat ) ################################################### ### code chunk number 9: inl-trans ################################################### fetal <- transform(fetal, tgac=tga, id=as.integer(factor(id)) ) ################################################### ### code chunk number 10: INLA1 ################################################### library( INLA ) system.time( im1 <- inla( hc ~ tga + f(id) + f(tgac), data=fetal ) ) summary( im1 ) ################################################### ### code chunk number 11: INLA2 ################################################### Nb <- max(fetal$id) fetal$xid <- fetal$id + Nb system.time( im2 <- inla( hc ~ 1 + tga + f( id, model = "iid2d", n=2*Nb ) + f( xid, tga, copy="id" ), data = fetal ) ) summary( im2 ) names( im2 ) ################################################### ### code chunk number 12: