### R code from vignette source 'lreg-s.rnw' ################################################### ### code chunk number 1: lreg-s.rnw:5-7 ################################################### library( rjags ) library( Epi ) ################################################### ### code chunk number 2: yx ################################################### x <- c(1,2,3,4,5,6) y <- c(1,3,3,3,5,7) plot( x, y, pch=16 ) summary( m0 <- lm(y ~ x) ) abline( m0 ) ################################################### ### code chunk number 3: lreg-s.rnw:25-30 ################################################### reg.dat <- list( x=x, y=y, I=6 ) reg.ini <- list( list( alpha=-1, beta=1.0, sigma=0.7 ), list( alpha=0, beta=0.5, sigma=1.0 ), list( alpha=1, beta=2, sigma=1.2 ) ) reg.par <- c("alpha","beta","sigma" ) ################################################### ### code chunk number 4: lreg-s.rnw:34-47 ################################################### cat( "model { for( i in 1:I ) { y[i] ~ dnorm(mu[i],tau) mu[i] <- alpha + beta*x[i] } tau <- 1/pow(sigma,2) alpha ~ dnorm(0, 1.0E-6) beta ~ dnorm(0, 1.0E-6) sigma ~ dunif(0,100) }", file="reg.jag" ) ################################################### ### code chunk number 5: lreg-s.rnw:52-57 ################################################### reg.mod <- jags.model( file = "reg.jag", data = reg.dat, n.chains = 3, inits = reg.ini, n.adapt = 10000 ) str(reg.mod) ################################################### ### code chunk number 6: lreg-s.rnw:61-65 ################################################### reg.res <- coda.samples( reg.mod, var = reg.par, n.iter = 10000, thin = 10 ) ################################################### ### code chunk number 7: lreg-s.rnw:70-73 ################################################### summary( reg.res ) round(ci.lin( m0 ),4) summary( m0 )$sigma # A c.i. for an sd estimate is obtained by multiplying # the estimate by # sqrt(f/chisq-f(97.5%)) and sqrt(f/chisq-f(2.5%)) names( m0 ) df <- m0$df cim <- sqrt( c( 1, m0$df/qchisq(0.975,m0$df), m0$df/qchisq(0.025,m0$df) ) ) summary( m0 )$sigma * cim summary(reg.res) # Compare: rbind( summary(reg.res)$qu["sigma",c(3,1,5)], summary( m0 )$sigma * cim ) ################################################### ### code chunk number 8: lreg-s.rnw:88-124 ################################################### data( births ) births <- subset( births, !is.na(gestwks) ) dim( births ) mb <- lm( bweight ~ I(gestwks-35), data=births ) summary( mb ) bth.dat <- list( x=births$gestwks-35, y=births$bweight, I=nrow(births) ) bth.ini <- list( list( alpha=2400, beta=200, sigma=400 ), list( alpha=2300, beta=150, sigma=450 ), list( alpha=2500, beta=250, sigma=500 ) ) bth.par <- c("alpha","beta","sigma" ) cat( "model { for( i in 1:I ) { y[i] ~ dnorm(mu[i],tau) mu[i] <- alpha + beta*x[i] } alpha ~ dnorm(0, 1.0E-6) beta ~ dnorm(0, 1.0E-6) sigma ~ dunif(0,10000) tau <- 1/pow(sigma,2) }", file="bth.jag" ) bth.mod <- jags.model( file = "bth.jag", data = bth.dat, n.chains = 3, inits = bth.ini, n.adapt = 2000 ) bth.res <- coda.samples( bth.mod, var = bth.par, n.iter = 10000, thin = 10 ) summary( mb )$sigma ################################################### ### code chunk number 9: lreg-s.rnw:134-136 ################################################### summary( bth.res )$quan[,c(3,1,5)] ci.lin( mb )[,c(1,5,6)] ################################################### ### code chunk number 10: lreg-s.rnw:143-150 ################################################### cim <- sqrt( c( 1, mb$df/qchisq(0.975,mb$df), mb$df/qchisq(0.025,mb$df) ) ) summary(mb)$sigma * cim summary( bth.res )$quan["sigma",c(3,1,5)] summary(mb)$sigma * cim / summary( bth.res )$quan["sigma",c(3,1,5)]