Edit the Initiate Function Below

# Edit the “initiate” function below.

# Then Run these four commands in succession, to run data analysis.

#1# initiate()

#2# for (i in 1:itend) iterate()

#3# summ()

#4# clean()

# Then open output file as a comma-delimited file (that you named below), through Excel.

# You may want to run a Windows search for this output file.

initiate <- function() {

# OUTPUT CONTROLS

Notes <- "NAEP"# Notes describing the analysis

MHdraws <- "NAEP-MH1xx" # name output file for MH chain

ISOPdraws <- "NAEP-DM1xx"# name output file for ISOP chain

OUTPUT <- "NAEP2xx.OUT"# name the final output file

LItest<-"Y"# Perform test of Local Independence?

#Enter “N” to improve MCMC sampling speed.

# DATA (enter)

n <- round(P*N,0)# number of observed binomial successes, {n_ij}

N <- N# number of trials, {Nij}

P <- n/N

P <- P# observed proportion, {Pij}

I<-count.rows(N)

J<-count.cols(N)

##SPECIAL RESTRICTIONS (enter)

S.min<-0

S.max<-1

##BETA PRIORS (enter)

a<-matrix(.5,I,J)

b<-matrix(.5,I,J)

# ****** MCMC ALGORITHM SPECIFICATIONS *******

itstart <- 10#choose iterations {t=1,...,500} as burn in period

itend <- 100#MCMC iterations {t=501,...,5000} interpreted for MH and ISOP posterior

T1 <- T2<- Tstart <- matrix(sort(runif(I*J,S.min,S.max)),I,J)# Starting values

iter <- 0

W1<-matrix(c(0,1),300,300)

W1<-W1[1:I,1:J]

W12<-matrix(c(0,1,0,1),301,300)

W12<-W12[1:I,1:J]

DIC.draw.MH<-DIC.draw.ISOP<-PP.IxJ.MH<-PP.IxJ.ISOP<-DIC.draw.Un <- 0

PP.ij.MH<-PP.ij.ISOP<-U<-matrix(0,I,J)

PP.i.MH<-PP.i.ISOP<-rep(0,I)

PP.j.MH<-PP.j.ISOP<-rep(0,J)

PP.LI.MH<-PP.LI.ISOP<-matrix(0,J,J)

return()}

###############################################################

iterate<-function() {

iter <- iter+1

cat("======ITERATION ", iter, " ======",fill=T)

T1<-draw.MH(T1)

T2<-draw.ISOP(T2)

write(round(c(T1),15),file=MHdraws,ncol=(count.rows(n)*count.cols(n)),append=T)

write(round(c(T2),15),file=ISOPdraws,ncol=(count.rows(n)*count.cols(n)),append=T)

if(iter>=itstart)DIC.draw.MH<-(2*sum((n*(log((n+.5)/((N*T1)+.5))))+((N-n)*(log((N-n+.5)/((N-(N*T1))+.5))))))+ DIC.draw.MH

if(iter>=itstart)DIC.draw.ISOP<-(2*sum((n*(log((n+.5)/((N*T2)+.5))))+((N-n)*(log((N-n+.5)/((N-(N*T2))+.5))))))+ DIC.draw.ISOP

U1<-rbeta(I*J,n+.5,N-n+.5)

U<-U1 + U

DIC.draw.Un<- (2*sum((n*(log((n+.5)/((N*U1)+.5))))+((N-n)*(log((N-n+.5)/((N-(N*U1))+.5)))))) + DIC.draw.Un

if(iter>=itstart)r.1<-matrix(rbinom(I*J,N,T1),I,J)

if(iter>=itstart)r.2<-matrix(rbinom(I*J,N,T2),I,J)

if(iter>=itstart)X2.1.rep<-((r.1-(N*T1))^2)/(N*T1)

if(iter>=itstart)X2.2.rep<-((r.2-(N*T2))^2)/(N*T2)

if(iter>=itstart)X2.1<-((n-(N*T1))^2)/(N*T1)

if(iter>=itstart)X2.2<-((n-(N*T2))^2)/(N*T2)

if(iter>=itstart)PP.ij.MH<-ifelse(X2.1.rep>=X2.1,1,0) + PP.ij.MH

if(iter>=itstart)PP.ij.ISOP<-ifelse(X2.2.rep>=X2.2,1,0) + PP.ij.ISOP

if(iter>=itstart)PP.i.MH<-ifelse(rowSums(X2.1.rep,na.rm=T)>=rowSums(X2.1,na.rm=T),1,0) + PP.i.MH

if(iter>=itstart)PP.i.ISOP<-ifelse(rowSums(X2.2.rep,na.rm=T)>=rowSums(X2.2,na.rm=T),1,0) + PP.i.ISOP

if(iter>=itstart)PP.j.MH<-ifelse(colSums(X2.1.rep,na.rm=T)>=colSums(X2.1,na.rm=T),1,0) + PP.j.MH

if(iter>=itstart)PP.j.ISOP<-ifelse(colSums(X2.2.rep,na.rm=T)>=colSums(X2.2,na.rm=T),1,0) + PP.j.ISOP

if(iter>=itstart)PP.IxJ.MH<-ifelse(sum(X2.1.rep,na.rm=T)>=sum(X2.1,na.rm=T),1,0) + PP.IxJ.MH

if(iter>=itstart)PP.IxJ.ISOP<-ifelse(sum(X2.2.rep,na.rm=T)>=sum(X2.2,na.rm=T),1,0) + PP.IxJ.ISOP

if(iter>=itstart & LItest=="Y")PP.LI.MH<-ifelse(abs(cor(r.1-(N*T1)))>=abs(cor(n-(N*T1))),1,0) + PP.LI.MH

if(iter>=itstart & LItest=="Y")PP.LI.ISOP<-ifelse(abs(cor(r.2-(N*T2)))>=abs(cor(n-(N*T2))),1,0) + PP.LI.ISOP

return(iter)}

###############################################################

draw.MH <- function(T1s){

# MH: STRIPE SAMPLER

min. <- rbind(matrix(S.min,1,J),if (J==1) matrix(T1[-I,],I-1,J) else T1[-I,])

max.<-rbind(if (J==1) matrix(T1[-1,],I-1,J) else T1[-1,],matrix(S.max,1,J))

draw<-matrix(qbeta(matrix(pbeta(min.,n+a,N-n+b),I,J) + (matrix(runif(I*J),I,J)*(matrix(pbeta(max.,n+a,N-n+b),I,J)- matrix(pbeta(min.,n+a,N-n+b),I,J))),n+a,N-n+b),I,J)

T1<-ifelse(W1,draw,T1)

min. <- rbind(matrix(S.min,1,J),if (J==1) matrix(T1[-I,],I-1,J) else T1[-I,])

max.<-rbind(if (J==1) matrix(T1[-1,],I-1,J) else T1[-1,],matrix(S.max,1,J))

draw<-matrix(qbeta(matrix(pbeta(min.,n+a,N-n+b),I,J) + (matrix(runif(I*J),I,J)*(matrix(pbeta(max.,n+a,N-n+b),I,J)- matrix(pbeta(min.,n+a,N-n+b),I,J))),n+a,N-n+b),I,J)

T1<-ifelse(1-W1,draw,T1)

return(T1)}

draw.ISOP <- function(T2s){

# ISOP: CHECKER ALGORITHM

min1 <- rbind(matrix(S.min,1,J),if (J==1) matrix(T2[-I,],I-1,J) else T2[-I,])#up

min2<-cbind(matrix(S.min,I,1),if (I==1) matrix(T2[,-J],I,J-1) else T2[,-J]) # left

min.<-ifelse(min1<min2,min2,min1)

max1<-rbind(if (J==1) matrix(T2[-1,],I-1,J) else T2[-1,],matrix(S.max,1,J))

max2<-cbind(if (I==1) matrix(T2[,-1],I,J-1) else T2[,-1], matrix(S.max,I,1))

max.<-ifelse(max1<max2,max1,max2)

draw<-matrix(qbeta(matrix(pbeta(min.,n+a,N-n+b),I,J) + (matrix(runif(I*J),I,J)*(matrix(pbeta(max.,n+a,N-n+b),I,J)- matrix(pbeta(min.,n+a,N-n+b),I,J))),n+a,N-n+b),I,J)

T2<-ifelse(W12,draw,T2)

min1 <- rbind(matrix(S.min,1,J),if (J==1) matrix(T2[-I,],I-1,J) else T2[-I,])#up

min2<-cbind(matrix(S.min,I,1),if (I==1) matrix(T2[,-J],I,J-1) else T2[,-J]) # left

min.<-ifelse(min1<min2,min2,min1)

max1<-rbind(if (J==1) matrix(T2[-1,],I-1,J) else T2[-1,],matrix(S.max,1,J))

max2<-cbind(if (I==1) matrix(T2[,-1],I,J-1) else T2[,-1], matrix(S.max,I,1))

max.<-ifelse(max1<max2,max1,max2)

draw<-matrix(qbeta(matrix(pbeta(min.,n+a,N-n+b),I,J) + (matrix(runif(I*J),I,J)*(matrix(pbeta(max.,n+a,N-n+b),I,J)- matrix(pbeta(min.,n+a,N-n+b),I,J))),n+a,N-n+b),I,J)

T2<-ifelse(1-W12,draw,T2)

return(T2)}

draw.ISOP3 <- function(T2s){

# ISOP: CHECKER ALGORITHM

min1 <- rbind(matrix(S.min,1,J),if (J==1) matrix(T2[-I,],I-1,J) else T2[-I,])#up

min2<-cbind(matrix(S.min,I,1),if (I==1) matrix(T2[,-J],I,J-1) else T2[,-J]) # left

min3<-

min.<-max(min1,min2,min3)

max1<-rbind(if (J==1) matrix(T2[-1,],I-1,J) else T2[-1,],matrix(S.max,1,J))

max2<-cbind(if (I==1) matrix(T2[,-1],I,J-1) else T2[,-1], matrix(S.max,I,1))

max.<-ifelse(max1<max2,max1,max2)

draw<-matrix(qbeta(matrix(pbeta(min.,n+a,N-n+b),I,J) + (matrix(runif(I*J),I,J)*(matrix(pbeta(max.,n+a,N-n+b),I,J)- matrix(pbeta(min.,n+a,N-n+b),I,J))),n+a,N-n+b),I,J)

T2<-ifelse(W12,draw,T2)

min1 <- rbind(matrix(S.min,1,J),if (J==1) matrix(T2[-I,],I-1,J) else T2[-I,])#up

min2<-cbind(matrix(S.min,I,1),if (I==1) matrix(T2[,-J],I,J-1) else T2[,-J]) # left

min.<-ifelse(min1<min2,min2,min1)

max1<-rbind(if (J==1) matrix(T2[-1,],I-1,J) else T2[-1,],matrix(S.max,1,J))

max2<-cbind(if (I==1) matrix(T2[,-1],I,J-1) else T2[,-1], matrix(S.max,I,1))

max.<-ifelse(max1<max2,max1,max2)

draw<-matrix(qbeta(matrix(pbeta(min.,n+a,N-n+b),I,J) + (matrix(runif(I*J),I,J)*(matrix(pbeta(max.,n+a,N-n+b),I,J)- matrix(pbeta(min.,n+a,N-n+b),I,J))),n+a,N-n+b),I,J)

T2<-ifelse(1-W12,draw,T2)

return(T2)}

###############################################################

summ_function(){

import.data(DataFrame="MH.Draws",FileName=MHdraws,FileType="ASCII") # IMPORT

import.data(DataFrame="ISOP.Draws",FileName=ISOPdraws,FileType="ASCII") # IMPORT

### CALCULATE POSTERIORS: MODEL PARAMETERS###

Theta.b1 <- MH.Draws[itstart:itend,]

Theta.b2 <- ISOP.Draws[itstart:itend,]

Pmn1<-round(matrix(colMeans(Theta.b1),I,J),3)

Pmn2<-round(matrix(colMeans(Theta.b2),I,J),3)

col <- 0

ptile1<-c()

percentile <- function(){

col <- col+1

ptile1 <-rbind(ptile1,quantile(Theta.b1[,col],c(.01,.025,.975,.99)))

return(ptile1)}

for (i in 1:(I*J)) percentile()

col <- 0

ptile2<-c()

percentile <- function(){

col <- col+1

ptile2 <-rbind(ptile2,quantile(Theta.b2[,col],c(.01,.025,.975,.99)))

return(ptile1)}

for (i in 1:(I*J)) percentile()

### CALCULATE MODEL FIT ###

Gl.Fit1<-round(PP.IxJ.MH/(itend-itstart+1),3)

i.Fit1<-round(PP.i.MH/(itend-itstart+1),3)

j.Fit1<-round(PP.j.MH/(itend-itstart+1),3)

ij.Fit1<-round(PP.ij.MH/(itend-itstart+1),3)

Gl.Fit2<-round(PP.IxJ.ISOP/(itend-itstart+1),3)

i.Fit2<-round(PP.i.ISOP/(itend-itstart+1),3)

j.Fit2<-round(PP.j.ISOP/(itend-itstart+1),3)

ij.Fit2<-round(PP.ij.ISOP/(itend-itstart+1),3)

### CALCULATE DIC ###

DIC.b1 <- round((DIC.draw.MH) / (itend-itstart+1),3)

DIC.b2 <- round((DIC.draw.ISOP) / (itend-itstart+1),3)

DIC.theta1<- round(2*sum((n*(log((n+.5)/((N*Pmn1)+.5))))+((N-n)*(log((N-n+.5)/((N-(N*Pmn1))+.5))))),3)

DIC.theta2<- round(2*sum((n*(log((n+.5)/((N*Pmn2)+.5))))+((N-n)*(log((N-n+.5)/((N-(N*Pmn2))+.5))))),3)

U.mn<-round(U / itend,3)

DIC.bU<-round(DIC.draw.Un/itend,3)

DIC.U<- round(2*sum((n*(log((n+.5)/((N*U.mn)+.5))))+((N-n)*(log((N-n+.5)/((N-(N*U.mn))+.5))))),3)

DIC.rep<-rbind(c("MODEL","Dev(PM)","Mean Dev","Penalty","DIC","Pred p-val"),c("MH",DIC.theta1,DIC.b1,DIC.b1-DIC.theta1,round(DIC.b1+(DIC.b1-DIC.theta1),3),Gl.Fit1),c("ISOP",DIC.theta2,DIC.b2,DIC.b2-DIC.theta2,round(DIC.b2+(DIC.b2-DIC.theta2),3),Gl.Fit2))

DIC.rep<-rbind(DIC.rep,c("Unres",DIC.U,DIC.bU,DIC.bU-DIC.U,round(DIC.bU+(DIC.bU-DIC.U),3),""))

# CREATE OUTPUT TABLES

Data.Table<-rbind(c("",1:J,"","",1:J,"","",1:J),cbind(c(1:I),round(P,3),rep("",I),c(1:I),n,rep("",I),c(1:I),N))

Post.Table1<-rbind(c("",1:J,"","",1:J,"","",1:J,"","",1:J,"","",1:J),cbind(1:I,Pmn1,rep("",I),c(1:I),matrix(round(ptile1[,1],3),I,J),rep("",I),c(1:I),matrix(round(ptile1[,2],3),I,J),rep("",I),c(1:I),matrix(round(ptile1[,3],3),I,J),rep("",I),c(1:I),matrix(round(ptile1[,4],3),I,J)))

Post.Table2<-rbind(c("",1:J,"","",1:J,"","",1:J,"","",1:J,"","",1:J),cbind(1:I,Pmn2,rep("",I),c(1:I),matrix(round(ptile2[,1],3),I,J),rep("",I),c(1:I),matrix(round(ptile2[,2],3),I,J),rep("",I),c(1:I),matrix(round(ptile2[,3],3),I,J),rep("",I),c(1:I),matrix(round(ptile2[,4],3),I,J)))

Fit.Table1<-rbind(c("",1:J,"","Row Fit"),cbind(1:I,ij.Fit1,rep("",I),matrix(i.Fit1,I,1)),rep("",J+3),c("Col Fit",j.Fit1,"",Gl.Fit1))

Fit.Table2<-rbind(c("",1:J,"","Row Fit"),cbind(1:I,ij.Fit2,rep("",I),matrix(i.Fit2,I,1)),rep("",J+3),c("Col Fit",j.Fit2,"",Gl.Fit2))

Fit.Table<- cbind(Fit.Table1,rep("",I+3),Fit.Table2)

if(LItest=="Y")LI.Table<-rbind(c("",1:J,"","",1:J),cbind(1:J,ifelse(col(PP.LI.MH)>row(PP.LI.MH),round(PP.LI.MH/(itend-itstart+1),3),""),rep("",J),1:J,ifelse(col(PP.LI.ISOP)>row(PP.LI.ISOP),round(PP.LI.ISOP/(itend-itstart+1),3),"")))

## WRITE MH OUTPUT ##

write("======",file=OUTPUT,append=T)# Write title and analysis notes

write("**************** ISOTONIC REGRESSION MODEL ESTIMATION *****************",file=OUTPUT,append=T)

write("Using The Gibbs sampler",file=OUTPUT,append=T)

write("======",file=OUTPUT,append=T)

write(" ",file=OUTPUT,append=T)

write(Notes,file=OUTPUT,append=T)

write(" ",file=OUTPUT,append=T)

write("======",file=OUTPUT,append=T)

write("ITERATION & ALGORITHM REPORT",file=OUTPUT,append=T)

write(" ",file=OUTPUT,append=T)

write.table(matrix(c("MH sample file:",MHdraws),nrow=1,ncol=2),file=OUTPUT,append=T,dimnames.write=F,quote.strings=F,sep=" ")

write.table(matrix(c("ISOP sample file:",ISOPdraws),nrow=1,ncol=2),file=OUTPUT,append=T,dimnames.write=F,quote.strings=F,sep=" ")

write.table(matrix(c("Iterations interpreted for posterior distribution= ",itstart," to ",itend," Total iterations =",itend-(itstart-1)),nrow=1,ncol=6),file=OUTPUT,append=T,dimnames.write=F,sep="")

write.table(matrix(c("Burn in Iterations = ",1," to ",itstart-1),nrow=1,ncol=4),file=OUTPUT,append=T,dimnames.write=F,sep="")

write("Starting values for Thetas",file=OUTPUT,append=T)

write.table(round(Tstart,2),file=OUTPUT,append=T,dimnames.write=F,sep=",")

write("======",file=OUTPUT,append=T)

write(matrix(" ",nrow=5,ncol=1),file=OUTPUT,append=T)

write("======",file=OUTPUT,append=T)

write("GLOBAL MODEL FIT ANALYSIS",file=OUTPUT,append=T)

write("======",file=OUTPUT,append=T)

write.table(DIC.rep,file=OUTPUT,append=T,dimnames.write=F,sep=",")

write("======",file=OUTPUT,append=T)

write.table(matrix(c(rep("======",J+3),rep("======",J+4)),nrow=1),file=OUTPUT,append=T,sep=",")

write("POSTERIOR PREDICTIVE P-VALUES (CHI-SQUARE) : MH & ISOP",file=OUTPUT,append=T)

write.table(matrix(c("MH",rep("",J+3),"ISOP",rep("",J+3)),nrow=1),file=OUTPUT,append=T,sep=",")

write.table(matrix(c(rep("======",J+3),rep("======",J+4)),nrow=1),file=OUTPUT,append=T,sep=",")

write.table(Fit.Table,file=OUTPUT,append=T,sep=",")

write.table(matrix(c(rep("======",J+3),rep("======",J+4)),nrow=1),file=OUTPUT,append=T,sep=",")

write(matrix("",nrow=5,ncol=1),file=OUTPUT,append=T)