"n" <- 919 "J" <- 85 "y" <- c(0.78845736036427, 0.78845736036427, 1.06471073699243,.....) "county" <- c(1, 1, 1, 1, 2, 2,.....) "x" <- c(1, 0, 0, 0, 0, 0,.....) list(y=c(0.78845736036427, 0.78845736036427, 1.06471073699243,.....), county=c(1, 1, 1, 1, 2, 2,.....), x=c(1, 0, 0, 0, 0, 0,.....)) model { for (i in 1:n) { y[i] ~ dnorm(y.hat[i], tau.y) y.hat[i] <- a[county[i]] + b*x[i] } b ~ dnorm(0, 0.0001) tau.y <- pow(sigma.y, -2) sigma.y ~ dunif(0,100) for (j in 1:J) { a[j] ~ dnorm(mu.a, tau.a) } mu.a ~ dnorm(0, 0.0001) tau.a <- pow(sigma.a, -2) sigma.a ~ dunif(0, 100) } "sigma.y" <- 50 "a" <- c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0) "b" <- 0 "mu.a" <- 0 "sigma.a" <- 50 model in "radon.x.only.bug" data in "radon.x.only.jags.dat" compile inits in "radon.x.only.jags.init" initialize update 5000 monitor set sigma.y monitor set a monitor set b monitor set mu.a monitor set sigma.a update 10000 coda * exit Enter CODA index file name (or a blank line to exit) 1: /Users/jgill/Class.Multilevel/examples/radon/CODAindex.txt Enter CODA output file names, separated by return key (leave a blank line when you have finished) 1: /Users/jgill/Class.Multilevel/examples/radon/CODAchain1.txt 2: library(foreign); library(arm); options(signif=5) hiv.data <- read.csv("Class.Multilevel/examples/cd4/allvar.csv") attach.all(hiv.data) # ATTACH A BUGS OBJECT # just consider the "control" patients (treatmnt==1) and with # initial age between 1 and 5 years ok <- treatmnt==1 & !is.na(CD4PCT) & (baseage>1 & baseage<5) attach.all(hiv.data[ok,]) # setting up handy data structures y <- sqrt(CD4PCT) age.baseline <- baseage # kid's age (yrs) at the beginning of the study age.measurement <- visage # kids age (yrs) at the time of measurement treatment <- treatmnt time <- visage - baseage # set up unique patient id numbers from 1 to J unique.pid <- unique(newpid) n <- length(y) J <- length(unique.pid) person <- rep(NA, n) for (j in 1:J){ person[newpid==unique.pid[j]] <- j } sum(is.na(time)) [1] 3 sum(is.na(person)) [1] 3 cbind(person,time)[127:135,] person <- person[-c(130,131,132)]; time <- time[-c(130,131,132)] # fit multilevel model using lmer M1 <- lmer(y ~ time + (1 + time | person),na.action="na.omit") summary(M1) Linear mixed model fit by REML Formula: y ~ time + (1 + time | person) AIC BIC logLik deviance REMLdev 1108 1132 -548 1092 1096 Random effects: Groups Name Variance Std.Dev. Corr person (Intercept) 1.766 1.329 time 0.462 0.680 0.149 Residual 0.560 0.748 Number of obs: 369, groups: person, 83 Fixed effects: Correlation of Fixed Effects: Estimate Std. Error t value (Intr) (Intercept) 4.846 0.160 30.35 time -0.146 time -0.468 0.128 -3.67 plot(c(0,J),range(coef(M1)),type="n",xlab="",ylab="",xaxt="n",yaxt="n") points(1:(J-1),coef(M1)[[1]][,1],pch=20,col="forestgreen") for (j in 1:i(J-1)) lines(1:(J-1), coef(M1)[[1]][,1] + c(-1,1)*se.coef(M2)$county[j,], lwd=.5, col="gray10") cd4.rjags <- function() { for (i in 1:COUNTS) { mu[i] <- beta[person[i],1] + beta[person[i],2]*time[i] y[i] ~ dnorm(mu[i],tau.y) } for (j in 1:KIDS) { beta[j,1:2] ~ dmnorm(beta.mu,beta.W) } beta.mu[1] ~ dt(0,1,4) #(mu, tau, k) beta.mu[2] ~ dt(0,1,5) beta.W ~ dwish(R[,],2) R[1,1] <- 3.0 R[2,2] <- 3.0 R[1,2] <- 0.5 R[2,1] <- 0.5 tau.y ~ dgamma(1,0.1) } write.model(cd4.rjags, "Class.Multilevel/cd4.rjags") cd4.model <- jags.model(file="Class.Multilevel/cd4.rjags", data=cd4.list, n.chains=3, n.adapt=1000) update(cd4.model, n.iter=10000) cd4.mcmc <- coda.samples(model=cd4.model, n.iter=10000, variable.names=c("beta","beta.mu","beta.W","tau.y")) summary(cd4.mcmc) superdiag(as.mcmc.list(cd4.mcmc), burnin=0) ( cd4.dic <- dic.samples(cd4.model, n.iter=1000, type="pD") ) model { for (i in 1:n){ y[i] ~ dnorm(y.hat[i], tau.y) y.hat[i] <- mu + gamma[treatment[i]] + delta[airport[i]] } mu ~ dnorm(0, 0.0001) tau.y <- pow(sigma.y,-2) sigma.y ~ dunif(0, 100) for (j in 1:n.treatment){ gamma[j] ~ dnorm(0,tau.gamma) } tau.gamma <- pow(sigma.gamma,-2) sigma.gamma ~ dunif(0, 100) for (k in 1:airport){ delta[k] ~ dnorm(0, tau.delta) } tau.delta <- pow(sigma.delta,-2) sigma.delta~ dunif(0, 100) } for (i in 1:n) { e.y[i] = y[i] - y.hat[i] } s.y <- sd(e.y[]) # FINITE-POPULATION STANDARD DEVIATION s.g <- sd(gamma[]) # TREATMENT COEFFICIENT STANDARD DEVIATION s.d <- sd(delta[]) # AIRPORT COEFFICIENT STANDARD DEVIATION for (j in 1:n.treatment){ gamma[j] ~ dnorm(gamma.hat[j],tau.gamma) gamma.hat[j] <- a.0 + a.1*u.1[j] + a.2*u.2[j] e.g[j] <- gamma[j] - gamma.hat[j] } tau.gamma <- pow(sigma.gamma,-2) s.gamma <- sd(e.g[]) model star { for (j in 1:K) { phi[j] ~ dnorm(0.0,1.0E-6) I(-4,4); # T(-4,4) IN JAGS } tau ~ dgamma(1.0E-3,1.0E-3); # 1/sigma^2 sigma <- 1.0/sqrt(tau); for (i in 1:N) { b[i] ~ dnorm(0.0,tau); log(p[i]) <- (phi[1] + phi[2]*LOWINC[i] + phi[3]*PERASIAN[i] + phi[4]*PERBLACK[i] + phi[5]*PERHISP[i] + phi[6]*AVYRSEXP[i]+ phi[7]*PERPSPEN[i] + phi[8]*PCTAF[i])/exp(1 + phi[9]*permiNTE[i] + phi[10]*PTRATIO[i] + phi[11]*PCTYRRND[i]) + b[i]; READPASS[i] ~ dbin(p[i],READTOT[i]); } } graph.summary <- function (in.object, alpha = 0.05, digits=3) # DESIGNED BY NEAL, CODED BY JEFF { lo <- in.object$coefficient - qnorm(1-alpha/2) * summary(in.object)$coef[,2] hi <- in.object$coefficient + qnorm(1-alpha/2) * summary(in.object)$coef[,2] out.mat <- round(cbind(in.object$coefficient, summary(in.object)$coef[,2], lo, hi),digits) blanks <- " " dashes <- "--------------------------------------------------------------------------" bar.plot <- NULL scale.min <- floor(min(out.mat[,3])); scale.max <- ceiling(max(out.mat[,4])) for (i in 1:nrow(out.mat)) { ci.half.length <- abs(out.mat[i,1]-out.mat[i,3]) ci.start <- out.mat[i,1] - ci.half.length ci.stop <- out.mat[i,1] + ci.half.length bar <- paste("|",substr(dashes,1,ci.half.length), "o", substr(dashes,1,ci.half.length), "|", sep="", collapse="") start.buf <- substr(blanks,1,round(abs(scale.min - ci.start))) stop.buf <- substr(blanks,1,round(abs(scale.max - ci.stop))) bar.plot <- rbind( bar.plot, paste(start.buf,bar, stop.buf, sep="", collapse="") ) } out.df <- data.frame( matrix(NA,nrow=nrow(out.mat),ncol=ncol(out.mat)), bar.plot[1:length(bar.plot)] ) out.df[1:nrow(out.mat),1:ncol(out.mat)] <- out.mat CI.label <- paste( "CIs:", substr(blanks,1,abs(scale.min)-2-4),"ZE+RO", substr(blanks,1,abs(scale.max)-2), sep="", collapse="" ) dimnames(out.df)[[1]] <- dimnames(summary(in.object)$coef)[[1]] dimnames(out.df)[[2]] <- c("Coef","Std.Err.", paste(1-alpha,"Lower"),paste(1-alpha,"Upper"),CI.label) if (substr(in.object$call[1],1,2) == "gl") print(in.object$family) if (substr(in.object$call[1],1,2) == "lm") cat("\nFamily: gaussian\nLink function: identity\n\n") print(out.df) cat("\n") if (substr(in.object$call[1],1,2) == "gl") { cat( paste("N:", length(in.object$y)," log-likelihood:",round(ncol(out.mat)-in.object$aic/2,digits), " AIC:",round(in.object$aic,digits), " Dispersion Parameter:",summary(in.object)$dispersion,"\n") ) cat( paste(" Null deviance:",round(in.object$null.deviance,digits),"on",in.object$df.null,"degrees of freedom\n") ) cat( paste("Residual deviance: ",round(in.object$deviance,digits),"on",in.object$df.residual,"degrees of freedom\n") ) } if (substr(in.object$call[1],1,2) == "lm") { cat( paste("N:",length(in.object$fitted.values)," Estimate of Sigma:",round(summary(in.object)$sigma,digits)) ) } } dp.97 <- read.table("http://jgill.wustl.edu/data/dp.data",header=TRUE) dp.out <- glm(EXECUTIONS ~ INCOME + PERPOVERTY + PERBLACK + log(VC100k96) + SOUTH + PROPDEGREE, family=poisson,data=dp.97) graph.summary(dp.out) Family: poisson Link function: log Coef Std.Err. 0.95 Lower 0.95 Upper CIs: ZE+RO (Intercept) -6.307 4.182 -14.504 1.890 |--------o--------| INCOME 0.000 0.000 0.000 0.000 |o| PERPOVERTY 0.069 0.080 -0.088 0.226 |o| PERBLACK -0.095 0.023 -0.140 -0.050 |o| log(VC100k96) 0.221 0.443 -0.646 1.089 |o| SOUTH 2.310 0.429 1.469 3.151 |o| PROPDEGREE -19.702 4.468 -28.459 -10.945 |--------o--------| N: 17 log-likelihood: -34.738 AIC: 77.475 Dispersion Parameter: 1 Null deviance: 136.573 on 16 degrees of freedom Residual deviance: 18.212 on 10 degrees of freedom norr.df <- read.table("http://jgill.wustl.edu/data/norr.dat",header=TRUE) norr.out <- lm(d8892r2 ~ catholic+black90+urban90+law3689+ep4089n,data=norr.df) graph.summary(norr.out) Family: gaussian Link function: identity Coef Std.Err. 0.95 Lower 0.95 Upper CIs: ZE+RO (Intercept) 0.357 0.113 0.135 0.579 |o| catholic -0.011 0.002 -0.015 -0.006 |o| black90 -0.007 0.004 -0.015 0.000 |o| urban90 0.003 0.001 0.000 0.006 |o| law3689 0.005 0.002 0.001 0.009 |o| ep4089n -1.145 2.505 -6.054 3.764 |----o----| N: 48 Estimate of Sigma: 0.171 teach.df <- read.table("http://jgill.wustl.edu/data/spector.mazzeo.data", col.names=c("GPA","TUCE","PSI","GRADE")) teach.df$GRADE <- factor(teach.df$GRADE) levels(teach.df$GRADE) <- c("FAIL","PASS") teach.df$PSI <- factor(teach.df$PSI) levels(teach.df$PSI) <- c("OLD","NEW") summary(teach.df) GPA TUCE PSI GRADE Min. :2.060 Min. :12.00 OLD:18 FAIL:21 1st Qu.:2.812 1st Qu.:19.75 NEW:14 PASS:11 Mdian :3.065 Median :22.50 Mean :3.117 Mean :21.94 3rd Qu.:3.515 3rd Qu.:25.00 Max. :4.000 Max. :29.00 par(bg="slategray",mfrow=c(1,2),mar=c(2,2,2,2),oma=c(4,3,2,2), col.main="white",col.axis="white",col="white") boxplot(GPA~GRADE,col="green",main="GPA") boxplot(TUCE~GRADE,col="green",main="TUCE") teach.logit.fit <- glm(GRADE ~ GPA + TUCE + PSI, family=binomial(link=logit), data = teach.df) graph.summary(teach.logit.fit) Family: binomial Link function: logit Coef Std.Err. 0.95 Lower 0.95 Upper CIs: ZE+RO (Intercept) -13.021 4.931 -22.686 -3.356 |---------o---------| GPA 2.826 1.263 0.351 5.301 |--o--| TUCE 0.095 0.142 -0.182 0.373 |o| PSINEW 2.379 1.065 0.292 4.465 |--o--| N: 32 log-likelihood: -12.89 AIC: 33.779 Dispersion Parameter: 1 Null deviance: 41.183 on 31 degrees of freedom Residual deviance: 25.779 on 28 degrees of freedom attach(teach.df) par(bg="slategray",mfrow=c(1,2),cex=1.3,col.main="white",col.axis="white", col="white",mar=c(3,3,3,1)) ruler <- seq(min(GPA),max(GPA),length=500) xbeta <- 1/(1+exp(-teach.logit.fit$coef[1] - teach.logit.fit$coef[2]*ruler - teach.logit.fit$coef[3]*mean(teach.df$TUCE))) plot(ruler,xbeta,type="l",lwd=3,xlab="Levels of GPA",col="deepskyblue4", ylab="Probability of Passing",ylim=c(0,1),col.lab="white") text(3.5,0.02,"OLD Method",col="deepskyblue4") xbeta2 <- 1/(1+exp(-teach.logit.fit$coef[1] - teach.logit.fit$coef[2]*ruler - teach.logit.fit$coef[3]*mean(teach.df$TUCE) - teach.logit.fit$coef[4])) lines(ruler,xbeta2,lwd=3,col="deeppink3") text(3.0,0.9,"NEW Method",col="deeppink3") ruler <- seq(min(TUCE),max(TUCE),length=500) xbeta <- 1/(1+exp(-teach.logit.fit$coef[1] - teach.logit.fit$coef[3]*ruler - teach.logit.fit$coef[2]*mean(teach.df$GPA))) plot(ruler,xbeta,type="l",lwd=3,xlab="Levels of TUCE",col="deepskyblue4", ylab="Probability of Passing",ylim=c(0,1),col.lab="white") text(22,0.02,"OLD Method",col="deepskyblue4") xbeta2 <- 1/(1+exp(-teach.logit.fit$coef[1] - teach.logit.fit$coef[3]*ruler - teach.logit.fit$coef[2]*mean(teach.df$GPA) - teach.logit.fit$coef[4])) lines(ruler,xbeta2,lwd=3,col="deeppink3") text(19,0.75,"NEW Method",col="deeppink3") dp.1 <- dp.2 <- c(1,apply(dp.97[,c(3,4,5,6,7,15)],2,mean)) dp.1[6] <- 0; dp.2[6] <- 1 y.1 <- exp(dp.1 %*% dp.out$coef); y.2 <- exp(dp.2 %*% dp.out$coef) y.2 - y.1 X <- cbind(rep(1,nrow(dp.97)), as.matrix(dp.97[,3:5]), as.matrix(log(dp.97[,6])), as.matrix(dp.97[,7]), as.matrix(dp.97[,15])) X.0 <- cbind(X[,1:5],rep(0,length=nrow(X)),X[,7]) dimnames(X.0)[[2]] <- names(dp.out$coefficients) X.1 <- cbind(X[,1:5],rep(1,length=nrow(X)),X[,7]) dimnames(X.1)[[2]] <- names(dp.out$coefficients) par(mfrow=c(3,2),mar=c(4,3,2,2),oma=c(3,1,1,1),col.axis="white",col.lab="white", col.sub="white",col="white",bg="slategray") for (i in 2:(ncol(X.0)-1)) { if (i==6) i <- i+1 ruler <- seq(min(X.0[,i]),max(X.0[,i]),length=1000) xbeta0 <- exp(dp.out$coefficients[-i]%*%apply(X.0[,-i],2,mean) + dp.out$coefficients[i]*ruler) xbeta1 <- exp(dp.out$coefficients[-i]%*%apply(X.1[,-i],2,mean) + dp.out$coefficients[i]*ruler) plot(ruler,xbeta0,type="l",xlab="",ylab="", ylim=c(min(xbeta0,xbeta1)-2,max(xbeta0,xbeta1)) ) lines(ruler,xbeta1,type="b") mtext(outer=F,side=1,paste("Levels of",dimnames(X.0)[[2]][i]),cex=0.8,line=3) mtext(outer=F,side=2,"Expected Executions",cex=0.6,line=2) } plot(ruler[100:200],rep(ruler[400],101),bty="n",xaxt="n",yaxt="n",xlab="",ylab="", type="l",xlim=range(ruler),ylim=range(ruler)) lines(ruler[100:200],rep(ruler[600],101),type="b") text(ruler[445],ruler[400],"Non-South State",cex=1.4) text(ruler[390],ruler[700],"South State",cex=1.4)