### R code from vignette source 'fetal1-s.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: lmer0 ################################################### library( lme4 ) m0 <- lmer( sqrt(hc) ~ tga + (tga|id), data=fetal ) summary(m0) ################################################### ### code chunk number 7: vcoc ################################################### VarCorr( m0 ) ################################################### ### code chunk number 8: fetal1-s.rnw:111-114 ################################################### m0 <- lmer( hc ~ I(tga-18) + (I(tga-18)|id), data=fetal ) summary(m0) VarCorr( m0 ) ################################################### ### code chunk number 9: qq ################################################### r0 <- residuals(m0) par( mar=c(3,3,1,1), mgp=c(3,1,0)/1.6 ) qqnorm( r0, main="", pch=16,cex=0.6 ) qqline( r0, col="blue" ) ################################################### ### code chunk number 10: fetal1-s.rnw:169-172 ################################################### library( Epi ) tga.pt <- 14:22 ci.lin( m0, ctr.mat=cbind(1,tga.pt) ) ################################################### ### code chunk number 11: pr-mn ################################################### tr <- function( ga ) ga-0.0116638*(ga^2) ga.pt <- seq(25,41,0.5) mn.hc <- ci.lin( m0, ctr.mat=cbind(1,tr(ga.pt)-18) )[,c(1,5,6)] matplot( ga.pt, mn.hc, type="l", lty=1, col="black", lwd=c(2,1,1) ) ################################################### ### code chunk number 12: fetal1-s.rnw:205-214 ################################################### Sig.u <- as.matrix( VarCorr( m0 )$id ) Sig.b <- as.matrix( vcov( m0 ) ) sig.e <- attr( VarCorr(m0), "sc" ) pr.var <- diag( cbind(1,tr(ga.pt)-18) %*% (Sig.u+Sig.b) %*% rbind(1,tr(ga.pt)-18) ) + sig.e pr.hc <- ci.lin( m0, ctr.mat=cbind(1,tr(ga.pt)-18) )[,c(1,5,6,1,1)] pr.hc[,4] <- pr.hc[,4] - 1.96*sqrt(pr.var) pr.hc[,5] <- pr.hc[,5] + 1.96*sqrt(pr.var) ################################################### ### code chunk number 13: pr-ci ################################################### matplot( ga.pt, pr.hc, type="l", lty=1, col="black", lwd=c(2,1,1,2,2) ) ################################################### ### code chunk number 14: fetal1-s.rnw:251-257 ################################################### fetal.dat <- list( id = as.integer( factor(fetal$id) ), hc = fetal$hc, tga = fetal$tga, N = nrow(fetal), F = length( unique(fetal$id) ) ) str( fetal.dat ) ################################################### ### code chunk number 15: fetal1-s.rnw:263-307 ################################################### 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( f in 1:F ) { u[f,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]-18 ) 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 16: fetal1-s.rnw:329-330 ################################################### VarCorr( m0 ) ################################################### ### code chunk number 17: fetal1-s.rnw:336-351 ################################################### ( 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 18: fetal1-s.rnw:358-366 ################################################### library( rjags ) system.time( fetal.mod <- jags.model( file = "fetal.jag", data = fetal.dat, n.chains = 4, inits = fetal.ini, n.adapt = 10000 ) ) ################################################### ### code chunk number 19: fetal1-s.rnw:374-383 ################################################### system.time( fetal.res <- coda.samples( fetal.mod, var = c("beta","sigma.e","Sigma.u"), n.iter = 5000, thin = 20 ) ) str( fetal.res ) summary( fetal.res ) dim( as.matrix(fetal.res) ) colnames( as.matrix(fetal.res) ) ################################################### ### code chunk number 20: posterior ################################################### plot( fetal.res ) ################################################### ### code chunk number 21: joint-post ################################################### fetal.post <- as.data.frame( as.matrix( fetal.res ) ) names( fetal.post ) names( fetal.post ) <- gsub( "\\[", ".", names(fetal.post) ) names( fetal.post ) <- gsub( ",", ".", names(fetal.post) ) names( fetal.post ) <- gsub( "\\]", "", names(fetal.post) ) str( fetal.post ) with( fetal.post, plot( beta.1, beta.2, pch=16,cex=0.9, col=c("black","red","green","blue")[rep(1:4,each=250)] ) ) ################################################### ### code chunk number 22: fetal1-s.rnw:420-421 ################################################### subset( fetal, id==5 ) ################################################### ### code chunk number 23: cond-id ################################################### ( xf <- subset( fetal, id==5 )[c(1,5),] ) xf[2,"hc"] <- NA xf[,"id"] <- max(fetal$id)+1 xf ################################################### ### code chunk number 24: marg-id ################################################### xf <- xf[c(1,2,2),] xf[3,"id"] <- max(fetal$id)+2 xf ################################################### ### code chunk number 25: fetal1-s.rnw:449-457 ################################################### x.same <- data.frame( id = max(fetal$id)+3, hc = NA, ga = ga.pt, tga = tr(ga.pt) ) x.diff <- data.frame( id = max(fetal$id)+3+1:length(ga.pt), hc = NA, ga = ga.pt, tga = tr(ga.pt) ) ################################################### ### code chunk number 26: fetal1-s.rnw:462-466 ################################################### fetal.x <- rbind( fetal, xf, x.same, x.diff ) fetal.x[nrow(fetal)+0:10,] tail( fetal.x ) nrow( fetal.x ) ################################################### ### code chunk number 27: fetal1-s.rnw:479-527 ################################################### 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( f in 1:F ) { u[f,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]-18 ) hc[j] ~ dnorm( mu[j], tau.e ) } for( j in n:N ) { pr[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="fetalp.jag" ) ################################################### ### code chunk number 28: fetal1-s.rnw:534-547 ################################################### fetal.xdat <- list( id = as.integer( factor(fetal.x$id) ), hc = fetal.x$hc, tga = fetal.x$tga, n = nrow(fetal)+1, N = nrow(fetal.x), F = length( unique(fetal.x$id) ) ) system.time( fetal.xmod <- jags.model( file = "fetalp.jag", data = fetal.xdat, n.chains = 4, inits = fetal.ini, n.adapt = 5000 ) ) ################################################### ### code chunk number 29: fetal1-s.rnw:550-561 ################################################### rng <- (nrow(fetal)+1):nrow(fetal.x) ( mus <- paste("pr[",paste(range(rng),collapse=":"),"]",sep="") ) system.time( fetal.xres <- coda.samples( fetal.xmod, var = c("beta","sigma.e","Sigma.u",mus), n.iter = 5000, thin = 10 ) ) fetal.qnt <- summary( fetal.xres )$quantiles pr.rows <- rownames(fetal.qnt)[grep( "pr", rownames(fetal.qnt) )] wh <- as.numeric( gsub( "\\]","", gsub("pr\\[","", pr.rows ) ) ) cbind( fetal.x[wh,c("ga","tga")], fetal.qnt[pr.rows,c(1,3,5)] ) ################################################### ### code chunk number 30: pr-mm ################################################### matplot( ga.pt, pr.hc, type="l", lty=1, col="black", lwd=c(2,1,1,2,2) ) matlines( fetal.x[3100+1:33,"ga"], fetal.qnt[paste("pr[",3100+1:33,"]",sep=""),c(1,3,5)], col="red", lty=1, lwd=2 ) matlines( fetal.x[3100+33+1:33,"ga"], fetal.qnt[paste("pr[",3100+33+1:33,"]",sep=""),c(1,3,5)], col="blue", lty=1, lwd=2 ) ################################################### ### code chunk number 31: post ################################################### fetal.post <- as.data.frame( as.matrix( fetal.xres ) ) names( fetal.post ) <- gsub( "\\[", ".", names(fetal.post) ) names( fetal.post ) <- gsub( ",", ".", names(fetal.post) ) names( fetal.post ) <- gsub( "\\]", "", names(fetal.post) ) str( fetal.post ) plot( density( fetal.post$pr.3099 ), lwd=3, col="blue", main="", bty="n", xlab="Head circumference at 38 weeks", xlim=c(250,400) ) lines( density( fetal.post$pr.3100 ), lwd=3, col="red" ) rug( quantile( fetal.post$pr.3099, probs=1:3/4 ), lwd=2, col="blue" ) rug( quantile( fetal.post$pr.3100, probs=1:3/4 ), lwd=2, col="red" ) ################################################### ### code chunk number 32: fetal1-s.rnw:607-608 ################################################### save( fetal.res, fetal.xres, file="../data/fetal.res" )