### R code from vignette source 'fetal0.rnw' ################################################### ### code chunk number 1: read ################################################### fetal <- read.csv("http://BendixCarstensen.com/Bayes/Cph-2012/data/fetal.csv",header=TRUE) str( fetal ) head( fetal, 10 ) ################################################### ### code chunk number 2: tab-tab ################################################### with( fetal, addmargins( table( table(id) ) ) ) ################################################### ### code chunk number 3: trsf ################################################### par( mfrow=c(1,2), mar=c(3,3,1,1), mgp=c(3,1,0)/1.6 ) with( fetal, plot( tga, ga-0.0116638*(ga^2), pch=16, cex=0.5 ) ) abline(0,1,col="red") with( fetal, plot( ga, tga, pch=16, cex=0.5, xlab="Gestational age (GA)", ylab="Transformed GA" ) ) abline(0,1,col="red") ################################################### ### code chunk number 4: transfig ################################################### par( mfrow=c(1,3), mar=c(3,3,1,1), mgp=c(3,1,0)/1.6 ) id.sub <- sample( unique(fetal$id), 50 ) with( fetal, plot( ga, hc, type="n" ) ) for( i in id.sub ) with( subset(fetal,id==i), lines(ga,hc) ) with( fetal, plot( tga, hc, type="n" ) ) for( i in id.sub ) with( subset(fetal,id==i), lines(tga,hc) ) with( fetal, plot( tga, sqrt(hc), type="n" ) ) for( i in id.sub ) with( subset(fetal,id==i), lines(tga,sqrt(hc)) ) ################################################### ### code chunk number 5: lmer0 ################################################### library( lme4 ) m0 <- lmer( sqrt(hc) ~ tga + (1|id), data=fetal ) summary(m0) ################################################### ### code chunk number 6: fetal0.rnw:76-78 ################################################### fixef( m0 ) VarCorr( m0 ) ################################################### ### code chunk number 7: fetal0.rnw:82-84 ################################################### attr( VarCorr(m0)$id, "stddev" ) attr( VarCorr(m0), "sc" ) ################################################### ### code chunk number 8: fetal0.rnw:92-124 ################################################### cat(" # Fixing data to be used in model definition model { # The model for each observational unit for( j in 1:N ) { mu[j] <- beta[1] + beta[2] * ( tga[j]-18 ) + u[id[j]] hc[j] ~ dnorm( mu[j], tau.e ) } # Random effects for each person for( i in 1:I ) { u[i] ~ dnorm(0,tau.u) } # 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) # Between-person variation tau.u <- pow(sigma.u,-2) sigma.u ~ dunif(0,100) }", file="fetal0.jag" ) ################################################### ### code chunk number 9: fetal0.rnw:132-133 ################################################### length( unique(fetal$id) ) ################################################### ### code chunk number 10: fetal0.rnw:138-143 ################################################### 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 11: fetal0.rnw:147-162 ################################################### ( sigma.e <- attr(VarCorr(m0),"sc") ) ( sigma.u <- attr(VarCorr(m0)$id,"stddev") ) ( beta <- fixef( m0 ) ) fetal.ini <- list( list( sigma.e = sigma.e/3, sigma.u = sigma.u/3, beta = beta /3 ), list( sigma.e = sigma.e*3, sigma.u = sigma.u*3, beta = beta *3 ), list( sigma.e = sigma.e/3, sigma.u = sigma.u*3, beta = beta /3 ), list( sigma.e = sigma.e*3, sigma.u = sigma.u/3, beta = beta *3 ) ) ################################################### ### code chunk number 12: fetal0.rnw:167-175 ################################################### library( rjags ) system.time( fetal.mod <- jags.model( file = "fetal0.jag", data = fetal.dat, n.chains = 4, inits = fetal.ini, n.adapt = 100 ) ) ################################################### ### code chunk number 13: fetal0.rnw:181-190 ################################################### system.time( fetal.res <- coda.samples( fetal.mod, var = c("beta","sigma.e","sigma.u"), n.iter = 500, thin = 20 ) ) str( fetal.res ) summary( fetal.res ) dim( as.matrix(fetal.res) ) colnames( as.matrix(fetal.res) ) ################################################### ### code chunk number 14: fetal0.rnw:202-207 ################################################### fetal <- transform(fetal, tgac=tga) library(INLA) im0 <- inla( hc ~ tga + f(id), data=fetal ) summary( imo ) names( im0 )