### R code from vignette source 'ox-s.rnw' ################################################### ### code chunk number 1: ox-s.rnw:28-31 ################################################### library( Epi ) oxw <- read.table( "../data/ox.dat", header=TRUE ) str(oxw) ################################################### ### code chunk number 2: ox-s.rnw:33-35 ################################################### m1 <- lm( I(pulse-co) ~ 1, data=oxw ) summary( m1 ) ################################################### ### code chunk number 3: ox-s.rnw:39-40 ################################################### ci.lin( m1 ) ################################################### ### code chunk number 4: ox-s.rnw:60-61 ################################################### qchisq(c(0.025,0.975),177-1) ################################################### ### code chunk number 5: ox-s.rnw:66-67 ################################################### sqrt( (177-1) * summary(m1)$sigma^2 / qchisq(c(0.975,0.025),177-1) ) ################################################### ### code chunk number 6: ox-s.rnw:80-82 ################################################### n <- nrow( oxw ) coef(m1) + c(-1,1) * qt(0.975,n-1) * ( summary(m1)$sigma / sqrt(n) ) ################################################### ### code chunk number 7: ox-s.rnw:94-120 ################################################### library( rjags ) cat( "model { for( i in 1:I ) { d[i] ~ dnorm( delta, tausq ) } tausq <- pow( sigma, -2 ) sigma ~ dunif( 0, 1000 ) delta ~ dnorm( 0, 0.000001 ) }", file="m1.jag" ) m1.dat <- list( d=oxw$co-oxw$pulse, I=nrow(oxw) ) m1.ini <- list( list( sigma=5, delta=0 ), list( sigma=6, delta=1 ), list( sigma=4, delta=-1 ) ) m1.par <- c("sigma","delta") m1.mod <- jags.model( file = "m1.jag", data = m1.dat, n.chains = length(m1.ini), inits = m1.ini, n.adapt = 20000 ) m1.res <- coda.samples( m1.mod, var = m1.par, n.iter = 20000, thin = 10 ) ################################################### ### code chunk number 8: post-joint ################################################### str( m1.res ) m1.mat <- as.matrix(m1.res) par( mar=c(6,6,1,1)/2, mgp=c(3,1,0)/1.6 ) plot( m1.mat[,"delta"], m1.mat[,"sigma"], xlab="delta", ylab="sigma", pch=16, cex=0.4, las=1, bty="n" ) ################################################### ### code chunk number 9: summ1 ################################################### summary( m1.res ) ################################################### ### code chunk number 10: post-marg ################################################### par( mfrow=c(1,2), mar=c(3,3,1,1), mgp=c(3,1,0)/1.6, las=1 ) plot( density( m1.mat[,"delta"] ), lwd=3, xlab="delta", ylab="", bty="n", main="" ) abline( v=quantile(m1.mat[,"delta"],probs=c(2.5,25,50,75,97.5)/100), col="gray" ) plot( density( m1.mat[,"sigma"] ), lwd=3, xlab="delta", ylab="", bty="n", main="" ) abline( v=quantile(m1.mat[,"sigma"],probs=c(2.5,25,50,75,97.5)/100), col="gray" ) ################################################### ### code chunk number 11: ox-s.rnw:164-192 ################################################### cat( "model { for( i in 1:I ) { d[i] ~ dnorm( delta, tausq ) } tausq <- pow( sigma, -2 ) sigma ~ dunif( 0, 1000 ) delta ~ dnorm( 0, 0.000001 ) agree.lo <- delta - 2*sigma agree.hi <- delta + 2*sigma }", file="m2.jag" ) m2.dat <- list( d=oxw$co-oxw$pulse, I=nrow(oxw) ) m2.ini <- list( list( sigma=5, delta=0 ), list( sigma=6, delta=1 ), list( sigma=4, delta=-1 ) ) m2.par <- c("sigma","delta","agree.lo","agree.hi") m2.mod <- jags.model( file = "m2.jag", data = m2.dat, n.chains = length(m2.ini), inits = m2.ini, n.adapt = 20000 ) m2.res <- coda.samples( m2.mod, var = m2.par, n.iter = 20000, thin = 10 ) summary( m2.res ) ################################################### ### code chunk number 12: post-LoA ################################################### M1 <- as.matrix( m1.res ) a1.lo <- M1[,"delta"] - 2*M1[,"sigma"] a1.hi <- M1[,"delta"] + 2*M1[,"sigma"] M2 <- as.matrix( m2.res ) layout( cbind(1:2), heights=1:2 ) par( mar=c(3,3,1,1), mgp=c(3,1,0)/1.6, las=1, bty="n" ) plot( density( a1.hi ), type="l", xlim=c(-20,20), lwd=3, main="" ) lines( density( a1.lo ), lwd=3 ) lines( density( M2[,"agree.hi"] ), lwd=2, col="red" ) lines( density( M2[,"agree.lo"] ), lwd=2, col="red" ) plot( M2[,"agree.hi"], M2[,"agree.lo"], xlab="Upper limit", ylab="Lower limit", pch=16, cex=0.3 ) ################################################### ### code chunk number 13: ox-s.rnw:220-221 ################################################### summary( M2[,"agree.lo"] - (M2[,"delta"]-2*M2[,"sigma"]) ) ################################################### ### code chunk number 14: ox-s.rnw:239-265 ################################################### cat( "model { for( i in 1:I ) { d[i] ~ dnorm( delta, tausq ) } tausq <- pow( sigma, -2 ) sigma ~ dunif( 0, 1000 ) delta ~ dnorm( 0, 0.4444444 ) }", file="m3.jag" ) m3.dat <- list( d=oxw$co-oxw$pulse, I=nrow(oxw) ) m3.ini <- list( list( sigma=5, delta=0 ), list( sigma=6, delta=1 ), list( sigma=4, delta=-1 ) ) m3.par <- c("sigma","delta") m3.mod <- jags.model( file = "m3.jag", data = m3.dat, n.chains = length(m3.ini), inits = m3.ini, n.adapt = 20000 ) m3.res <- coda.samples( m3.mod, var = m3.par, n.iter = 20000, thin = 10 ) summary( m3.res ) ################################################### ### code chunk number 15: post-2 ################################################### M3 <- as.matrix( m3.res ) par( mar=c(3,3,1,1), mgp=c(3,1,0)/1.6, las=1 ) plot( density( M3[,"delta"]), type="l", col=gray(0.2), lwd=3, main="", bty="n", xlab="" ) lines( density( M2[,"delta"] ), lwd=3, col="red" ) xx <- seq(0,5,,200) lines( xx, dnorm(xx,mean=0,sd=1.5), lwd=3, col=gray(0.6) ) ################################################### ### code chunk number 16: ox-s.rnw:330-362 ################################################### cat( "model { for( i in 1:I ) { mu.co[i] ~ dnorm( 0, 0.000001 ) mu.pl[i] <- mu.co[i] + delta y.co[i] ~ dnorm( mu.co[i], tausq.co ) y.pl[i] ~ dnorm( mu.pl[i], tausq.pl ) } tausq.co <- pow( sigma.co, -2 ) tausq.pl <- pow( sigma.pl, -2 ) sigma.co ~ dunif( 0, 1000 ) sigma.pl ~ dunif( 0, 1000 ) delta ~ dnorm( 0, 0.000001 ) }", file="m4.jag" ) nr <- nrow(oxw) m4.dat <- list( y.co=oxw$co, y.pl=oxw$pulse, I=nr ) m4.ini <- list( list( sigma.co=5, sigma.pl=5, mu.co=rep(80,nr), delta=0 ), list( sigma.co=6, sigma.pl=6, mu.co=rep(70,nr), delta=1 ), list( sigma.co=4, sigma.pl=4, mu.co=rep(90,nr), delta=-1 ) ) m4.par <- c("sigma.pl","sigma.co","delta") m4.mod <- jags.model( file = "m4.jag", data = m4.dat, n.chains = length(m4.ini), inits = m4.ini, n.adapt = 20000 ) m4.res <- coda.samples( m4.mod, var = m4.par, n.iter = 20000, thin = 10 ) summary( m4.res ) ################################################### ### code chunk number 17: trace4 ################################################### print( xyplot( m4.res[,c("delta","sigma.co","sigma.pl")], aspect="fill", layout=c(3,1) ) ) ################################################### ### code chunk number 18: joint4 ################################################### M4 <- as.matrix( m4.res, chains=TRUE ) plot( M4[,"sigma.co"], M4[,"sigma.pl"], pch=16, cex=0.5, col=rainbow(3)[M4[,"CHAIN"]] ) ################################################### ### code chunk number 19: dens4 ################################################### print( densityplot( m4.res[,c("delta","sigma.co","sigma.pl")], aspect="fill", layout=c(3,1), lwd=3 ) ) ################################################### ### code chunk number 20: ox-s.rnw:439-477 ################################################### cat( "model { for( i in 1:I ) { mu[i] ~ dunif( 0, 100 ) mu.co[i] <- mu[i] + a[i,repl[i]] mu.pl[i] <- mu[i] + a[i,repl[i]] + delta y.co[i] ~ dnorm( mu.co[i], tausq.co ) y.pl[i] ~ dnorm( mu.pl[i], tausq.pl ) for( r in 1:3 ) { a[i,r] ~ dnorm( 0, iomegasq ) } } tausq.co <- pow( sigma.co, -2 ) tausq.pl <- pow( sigma.pl, -2 ) iomegasq <- pow( omega, -2 ) sigma.co ~ dunif( 0, 1000 ) sigma.pl ~ dunif( 0, 1000 ) omega ~ dunif( 0, 1000 ) delta ~ dnorm( 0, 0.000001 ) }", file="m5.jag" ) m5.dat <- list( y.co=oxw$co, y.pl=oxw$pulse, repl=oxw$repl, I=nr ) m5.ini <- list( list( sigma.co=5, sigma.pl=5, omega=4, mu=rep(80,nr), delta= 0 ), list( sigma.co=6, sigma.pl=6, omega=4, mu=rep(70,nr), delta= 1 ), list( sigma.co=4, sigma.pl=4, omega=4, mu=rep(90,nr), delta=-1 ) ) m5.par <- c("sigma.pl","sigma.co","omega","delta") m5.mod <- jags.model( file = "m5.jag", data = m5.dat, n.chains = length(m5.ini), inits = m5.ini, n.adapt = 20000 ) m5.res <- coda.samples( m5.mod, var = m5.par, n.iter = 20000, thin = 10 ) summary( m5.res ) ################################################### ### code chunk number 21: trace5 ################################################### print( xyplot( m5.res[,c("delta","omega","sigma.co","sigma.pl")], aspect="fill", scales="same", layout=c(4,1) ) ) ################################################### ### code chunk number 22: pairs5 ################################################### M5 <- as.matrix( m5.res ) pairs( M5[,-1], gap=0, pch=16, cex=0.3 ) ################################################### ### code chunk number 23: ox-s.rnw:518-529 ################################################### oxl <- data.frame( y = c(oxw$co,oxw$pulse), repl = factor( rep(oxw$repl,2) ) , id = factor( rep(oxw$item,2) ), meth = factor( rep(c("co","pulse"),each=177) ) ) library( nlme ) m1 <- lme( y ~ meth + id, random = list( id = pdIdent( ~ repl-1 ) ), weights = varIdent( form = ~1 | meth ), data = oxl, control = lmeControl(returnObject=TRUE) ) m1 ################################################### ### code chunk number 24: ox-s.rnw:561-600 ################################################### cat( "model { for( i in 1:I ) { mu[i] ~ dunif( 0, 100 ) mu.co[i] <- mu[i] + a[i,repl[i]] mu.pl[i] <- alpha + beta * ( mu[i] + a[i,repl[i]] ) y.co[i] ~ dnorm( mu.co[i], tausq.co ) y.pl[i] ~ dnorm( mu.pl[i], tausq.pl ) for( r in 1:3 ) { a[i,r] ~ dnorm( 0, iomegasq ) } } tausq.co <- pow( sigma.co, -2 ) tausq.pl <- pow( sigma.pl, -2 ) iomegasq <- pow( omega, -2 ) sigma.co ~ dunif( 0, 1000 ) sigma.pl ~ dunif( 0, 1000 ) omega ~ dunif( 0, 1000 ) alpha ~ dnorm( 0, 0.000001 ) beta ~ dunif( 0, 2 ) }", file="m6.jag" ) m6.dat <- list( y.co=oxw$co, y.pl=oxw$pulse, repl=oxw$repl, I=nr ) m6.ini <- list( list( sigma.co=5, sigma.pl=5, omega=4 ), list( sigma.co=6, sigma.pl=6, omega=4 ), list( sigma.co=4, sigma.pl=4, omega=4 ) ) m6.par <- c("sigma.pl","sigma.co","omega","alpha","beta") m6.mod <- jags.model( file = "m6.jag", data = m6.dat, n.chains = length(m6.ini), inits = m6.ini, n.adapt = 20000 ) m6.res <- coda.samples( m6.mod, var = m6.par, n.iter = 20000, thin = 10 ) summary( m6.res ) ################################################### ### code chunk number 25: ox-s.rnw:616-655 ################################################### cat( "model { for( i in 1:I ) { mu[i] ~ dunif( 0, 100 ) mu.co[i] <- alpha + beta * ( mu[i] + a[i,repl[i]] ) mu.pl[i] <- mu[i] + a[i,repl[i]] y.co[i] ~ dnorm( mu.co[i], tausq.co ) y.pl[i] ~ dnorm( mu.pl[i], tausq.pl ) for( r in 1:3 ) { a[i,r] ~ dnorm( 0, iomegasq ) } } tausq.co <- pow( sigma.co, -2 ) tausq.pl <- pow( sigma.pl, -2 ) iomegasq <- pow( omega, -2 ) sigma.co ~ dunif( 0, 1000 ) sigma.pl ~ dunif( 0, 1000 ) omega ~ dunif( 0, 1000 ) alpha ~ dnorm( 0, 0.000001 ) beta ~ dunif( 0, 2 ) }", file="m7.jag" ) m7.dat <- list( y.co=oxw$co, y.pl=oxw$pulse, repl=oxw$repl, I=nr ) m7.ini <- list( list( sigma.co=5, sigma.pl=5, omega=4 ), list( sigma.co=6, sigma.pl=6, omega=4 ), list( sigma.co=4, sigma.pl=4, omega=4 ) ) m7.par <- c("sigma.pl","sigma.co","omega","alpha","beta") m7.mod <- jags.model( file = "m7.jag", data = m7.dat, n.chains = length(m7.ini), inits = m7.ini, n.adapt = 20000 ) m7.res <- coda.samples( m7.mod, var = m7.par, n.iter = 20000, thin = 10 ) summary( m7.res ) ################################################### ### code chunk number 26: ox-s.rnw:674-684 ################################################### # Mean ( ab6 <- summary( m6.res )$statistics[c("alpha","beta"),"Mean"] ) ( ab7 <- summary( m7.res )$statistics[c("alpha","beta"),"Mean"] ) abt <- c( -ab6[1]/ab6[2], 1/ab6[2] ) round( cbind( ab6, ab7, abt ), 3 ) # Median ( ab6 <- summary( m6.res )$quantiles[c("alpha","beta"),"50%"] ) ( ab7 <- summary( m7.res )$quantiles[c("alpha","beta"),"50%"] ) abt <- c( -ab6[1]/ab6[2], 1/ab6[2] ) round( cbind( ab6, ab7, abt ), 3 ) ################################################### ### code chunk number 27: ox-s.rnw:691-695 ################################################### round(ci.lin(lm(pulse~co,data=oxw))[,c(1,5,6)],3) round(summary(m6.res)$quantiles[4:5,c(3,1,5)],3) round(ci.lin(lm(co~pulse,data=oxw))[,c(1,5,6)],3) round(summary(m7.res)$quantiles[4:5,c(3,1,5)],3) ################################################### ### code chunk number 28: ox-s.rnw:752-795 ################################################### cat( "model { for( i in 1:I ) { mu[i] ~ dunif( 0, 100 ) mu.co[i] <- alpha.co + beta.co * ( mu[i] + a[i,repl[i]] ) mu.pl[i] <- alpha.pl + beta.pl * ( mu[i] + a[i,repl[i]] ) y.co[i] ~ dnorm( mu.co[i], tausq.co ) y.pl[i] ~ dnorm( mu.pl[i], tausq.pl ) for( r in 1:3 ) { a[i,r] ~ dnorm( 0, iomegasq ) } } tausq.co <- pow( sigma.co, -2 ) tausq.pl <- pow( sigma.pl, -2 ) iomegasq <- pow( omega, -2 ) sigma.co ~ dunif( 0, 1000 ) sigma.pl ~ dunif( 0, 1000 ) omega ~ dunif( 0, 1000 ) alpha.co ~ dnorm( 0, 0.000001 ) alpha.pl ~ dnorm( 0, 0.000001 ) beta.co ~ dunif( 0, 2 ) beta.pl ~ dunif( 0, 2 ) }", file="m8.jag" ) m8.dat <- list( y.co=oxw$co, y.pl=oxw$pulse, repl=oxw$repl, I=nrow(oxw) ) m8.ini <- list( list( sigma.co=5, sigma.pl=5, omega=4 ), list( sigma.co=6, sigma.pl=6, omega=4 ), list( sigma.co=4, sigma.pl=4, omega=4 ) ) m8.par <- c("sigma.pl","sigma.co","omega", "alpha.pl","alpha.co", "beta.pl", "beta.co") m8.mod <- jags.model( file = "m8.jag", data = m8.dat, n.chains = length(m8.ini), inits = m8.ini, n.adapt = 20000 ) m8.res <- coda.samples( m8.mod, var = m8.par, n.iter = 20000, thin = 10 ) summary( m8.res ) ################################################### ### code chunk number 29: trace8 ################################################### print(xyplot( m8.res[,c(7,3,6,2,5,1,4)], layout=c(2,4), aspect="fill" )) ################################################### ### code chunk number 30: ox-s.rnw:818-839 ################################################### m8b.res <- m8.res m8.res <- m8b.res str( m8.res ) for( i in 1:length(m8.res) ) { att <- attributes( m8.res[[i]] ) dfr <- as.data.frame( m8.res[[i]] ) dfr$beta.co.pl <- dfr$beta.co / dfr$beta.pl dfr$alpha.co.pl <- dfr$alpha.co - dfr$alpha.pl * dfr$beta.co.pl dfr$beta.pl.co <- dfr$beta.pl / dfr$beta.co dfr$alpha.pl.co <- dfr$alpha.pl - dfr$alpha.co * dfr$beta.pl.co dfr <- as.matrix( dfr ) att$dim <- dim( dfr ) att$dimnames <- dimnames( dfr ) attributes( dfr ) <- att m8.res[[i]] <- dfr } str( m8.res ) summary( m8.res ) round( ci.lin( lm( co ~ pulse, data=oxw ) ), 3 ) round( ci.lin( lm( pulse ~ co, data=oxw ) ), 3 ) ################################################### ### code chunk number 31: trace8x ################################################### wh <- c( grep( "sigma", varnames( m8.res ) ), grep( "omega", varnames( m8.res ) ), grep( "pl.co", varnames( m8.res ) ), grep( "co.pl", varnames( m8.res ) ) ) print(xyplot( m8.res[,wh], layout=c(4,2), aspect="fill", lwd=2 )) ################################################### ### code chunk number 32: dens8x ################################################### print( densityplot(m8.res[,wh],layout=c(4,2),lwd=2,aspect="fill") ) ################################################### ### code chunk number 33: oxdat ################################################### with( oxw, plot( co ~ pulse, pch=16, xlim=c(20,100), ylim=c(20,100) ) ) abline(0,1) abline( lm( co~pulse, data=oxw), col="red", lwd=2 ) cf <- coef( lm( pulse ~ co, data=oxw) ) abline( -cf[1]/cf[2], 1/cf[2], col="red", lwd=2 ) qnt <- summary( m8.res )$quantiles qnt <- qnt[grep("co.pl",rownames(qnt)),"50%"] abline( qnt[2], qnt[1], col="blue", lwd=2 )