# Uploaded: April 09, 2007 ####################### avgprob <- function (n, r, measure, samp.prior, param1, param2, delta = 0.80, tdi.p0 = 0.80, alpha = 0.05, AB = 0.001, n.rep = 100, n.iter = 1000, n.burn = 250) { if (identical(measure, "ccc")) { t.nr <- replicate(n.rep, ccc.prob (n, r, samp.prior, param1, param2, delta, alpha, AB, n.iter, n.burn)) res <- apply(t.nr, 1, mean) names(res) <- c("t11star", "t22star", "t12star") return(res) } else if (identical(measure, "tdi")) { t.nr <- replicate(n.rep, tdi.prob (n, r, samp.prior, param1, param2, delta, tdi.p0, alpha, AB, n.iter, n.burn)) res <- apply(t.nr, 1, mean) names(res) <- c("t11star", "t22star", "t12star") return(res) } else (stop (" 'measure' must be either 'ccc' or 'tdi' ")) } #################### ccc.prob <- function(n, r, samp.prior, param1, param2, delta, alpha, AB, n.iter, n.burn) { # Computes posterior probabilities for concordance correlation # as the measure of agreement if ( identical(samp.prior, "normal")) { mu <- rnorm(2, mean=param1[1:2], sd=param2[1:2]) theta <- exp(rnorm(4, mean=param1[3:6], sd=param2[3:6])) } else if ( identical(samp.prior, "uniform")) { mu <- runif(2, min=param1[1:2], max=param2[1:2]) theta <- exp(runif(4, min=param1[3:6], max=param2[3:6])) } else (stop (" 'samp.prior' must be either 'uniform' or 'normal' ")) if( identical(r, 1)) theta[4] <- 0 Y <- simData(mu, theta, n, r) sim.post <- simPost(AB, theta, Y, n.iter, n.burn) mu.post <- sim.post\$mu.post theta.post <- sim.post\$theta.post ccc.post <- ccc.fun(mu.post[,1], mu.post[,2], theta.post[,1], theta.post[,2], theta.post[,3], theta.post[,4]) cc11 <- ccc.post\$intra1 cc22 <- ccc.post\$intra2 cc12 <- ccc.post\$inter t11 <- -log(cc11) t22 <- -log(cc22) t12 <- -log(cc12) t.ubd <- apply(cbind(t11, t22, t12),2,quantile,prob=(1-alpha),type=4) ubd.delta<-delta*t.ubd F11 <- ecdf(t11) F22 <- ecdf(t22) F12 <- ecdf(t12) prob <- c(F11(ubd.delta[1]), F22(ubd.delta[2]), F12(ubd.delta[3])) return(prob) } ## ccc.fun <- function(mu1, mu2, theta1, theta2, theta3, theta4) { var1 <- theta3 + theta4 + theta1 var2 <- theta3 + theta4 + theta2 intra1 <- (theta3 + theta4)/var1 intra2 <- (theta3 + theta4)/var2 inter <- (2*theta3)/( var1 + var2 + (mu1-mu2)^2) return(list(intra1=intra1, intra2=intra2, inter=inter)) } #################### tdi.prob <- function(n, r, samp.prior, param1, param2, delta, tdi.p0, alpha, AB, n.iter, n.burn) { # Computes posterior probabilities for total deviation index # as the measure of agreement # mu = (beta1, beta2) # theta = (sigma1^2, sigma2^2, psi^2, psi_I^2) if ( identical(samp.prior, "normal")) { mu <- rnorm(2, mean=param1[1:2], sd=param2[1:2]) theta <- exp(rnorm(4, mean=param1[3:6], sd=param2[3:6])) } else if ( identical(samp.prior, "uniform")) { mu <- runif(2, min=param1[1:2], max=param2[1:2]) theta <- exp(runif(4, min=param1[3:6], max=param2[3:6])) } else (stop (" 'samp.prior' must be either 'uniform' or 'normal' ")) if( identical(r, 1)) theta[4] <- 0 Y <- simData(mu, theta, n, r) sim.post <- simPost(AB, theta, Y, n.iter, n.burn) mu.post <- sim.post\$mu.post theta.post <- sim.post\$theta.post mu.p <- mu.post[,1]-mu.post[,2] sigma.p <- sqrt(theta.post[,1]+theta.post[,2]+2*theta.post[,4]) q11 <- sqrt(theta.post[,1])*sqrt(2*qchisq(tdi.p0,1)) q22 <- sqrt(theta.post[,2])*sqrt(2*qchisq(tdi.p0,1)) q12 <- q.p(tdi.p0, mu.p, sigma.p) ubd.q <- apply(cbind(q11, q22, q12),2,quantile,prob=(1-alpha),type=4) ubd.delta<-delta*ubd.q F11 <- ecdf(q11) F22 <- ecdf(q22) F12 <- ecdf(q12) prob <- c(F11(ubd.delta[1]), F22(ubd.delta[2]), F12(ubd.delta[3])) return(prob) } ## q.p <- function(p, du,sdev) { #calculates the p_th qunatile of |diff| #du and sdev can be vector, but must have the same length ans <- mapply(q.scalar, m=du, s=sdev, p=p) return(ans) } ## q.scalar <- function(p, m, s) { # m=mean,s=sigma, p=prob #p must in (0.5,1) a<-m/s p<-1+p llim <- abs(a) ulim <- llim + 1.e-6 + qnorm(p/2) f <- function(q) { pnorm(q-a)+pnorm(q+a)- p} ans <- uniroot(f, c(llim, ulim), tol=1.e-6, maxiter=1e+4)\$root return(s*ans) } ############################# simData <- function(mu,theta,n,r) { #simulates data for given parameter values. sigma<-sqrt(theta) if(identical(r, 1)) sigma[4] <- 0 bI <- rnorm(n)*sigma[3] Y1 <- mu[1]+bI+rnorm(n)*sigma[4]+rnorm(n*r)*sigma[1] Y2 <- mu[2]+bI+rnorm(n)*sigma[4]+rnorm(n*r)*sigma[2] return(list(Y1=matrix(Y1,ncol=r),Y2=matrix(Y2,ncol=r))) } ################################ simPost <- function(AB, initTheta, data, n.iter, n.burn) { #simulates from posterior using Gibbs sampler. theta<-initTheta stat<-lapply(data,data.stat) theta.par<-list(rMeans=rbind(stat\$Y1\$rMeans,stat\$Y2\$rMeans), sv=c(stat\$Y1\$sv,stat\$Y2\$sv),dim=stat\$Y1\$dim) dY<-rbind(stat\$Y1\$rMeans.c, stat\$Y2\$rMeans.c, (stat\$Y1\$rMeans.c-stat\$Y2\$rMeans.c)) gamma.par<-list(mean=c(stat\$Y1\$mean,stat\$Y2\$mean),dY=dY,dim=stat\$Y1\$dim) for(burn in seq(n.burn)){ gamma<-gamma.p(theta,gamma.par) theta<-theta.p(gamma,theta.par,AB) } post.draw<-n.iter-n.burn q.emf<-seq(post.draw) mu.post<-matrix(0,nrow=post.draw,ncol=2) theta.post<-matrix(0,nrow=post.draw,ncol=4) for( size in seq(post.draw)){ gamma<-gamma.p(theta,gamma.par) theta<-theta.p(gamma,theta.par,AB) mu.post[size,]<-gamma\$mu theta.post[size,]<-theta } return(list(mu.post=mu.post,theta.post=theta.post)) } ################################# data.stat <- function(data) { # computes summary statistics for the data Y1<-data Y1.dim<-dim(Y1) #n*r Y1.rm<-rowMeans(Y1) Y1.sv<-sum((Y1-Y1.rm)^2) Y1.m<-mean(Y1.rm) Y1.c<-Y1.rm-Y1.m return(list(mean=Y1.m,sv=Y1.sv,dim=Y1.dim,rMeans=Y1.rm,rMeans.c=Y1.c)) } ############################### theta.p <- function(gamma,parY,AB) { #computes full conditional of theta given gamma(=(mu,b)) and csY A<-rep(AB,4) B<-A n<-parY\$dim[1] r<-parY\$dim[2] mu<-gamma\$mu bs<-gamma\$bs dYb1<-parY\$rMeans[1,]-mu[1]-bs[1,]-bs[2,] dYb2<-parY\$rMeans[2,]-mu[2]-bs[1,]-bs[3,] d1<-r*sum(dYb1^2)+parY\$sv[1] d2<-r*sum(dYb2^2)+parY\$sv[2] pA1=A[1]+n*r/2 pA2=A[2]+n*r/2 pB1=B[1]+d1/2 pB2=B[2]+d2/2 pA3=A[3]+n/2 pB3=B[3]+sum((bs[1,])^2)/2 c1<-pB1/rgamma(1,shape=pA1) c2<-pB2/rgamma(1,shape=pA2) c3<-pB3/rgamma(1,shape=pA3) if( identical (r, 1)) { c4 <- 0 } else{ pA4 <- A[4]+n pB4 <- B[4]+(sum(bs[2,]^2)+sum(bs[3,]^2))/2 c4 <- 1/rgamma(1,shape=pA4,rate=pB4) } return(c(c1,c2,c3,c4)) } ################################ gamma.p<-function(theta,parY) { # computes full conditional of gamma(=mu,b) given theta and csY n<-parY\$dim[1] r<-parY\$dim[2] theta1r<-theta[1]/r;theta2r<-theta[2]/r l1<-sqrt(theta1r+theta[3]+theta[4]) l3<-theta[3]/l1 l2<-sqrt(theta2r+theta[3]+theta[4]-l3^2) N0<-rnorm(2)/sqrt(n) mu<-c(l1*N0[1],l3*N0[1]+l2*N0[2])+parY\$mean lN1<-(N0[1]-N0[2]*l3/l2)/l1 lN2<-N0[2]/l2 t1<-theta1r+theta[4];t2<-theta2r+theta[4] g<-theta[3]/(t1*t2+theta[3]*(t1+t2)) b<-sqrt(g*t1*t2)*rnorm(n) s1<-sqrt(theta1r*(theta[4]/t1))*rnorm(n)-(theta[4]/t1)*b s2<-sqrt(theta2r*(theta[4]/t2))*rnorm(n)-(theta[4]/t2)*b d1<-(t2/theta[3])*parY\$dY[1,] d2<-(t1/theta[3])*parY\$dY[2,] d12<-parY\$dY[3,] b<-b+(g*theta[3])*(d1+d2)-theta[3]*(lN1+lN2) s1<-s1+(g*theta[4])*(d1+d12)-theta[4]*lN1 s2<-s2+(g*theta[4])*(d2-d12)-theta[4]*lN2 return(list(mu=mu,bs=rbind(b,s1,s2))) } ############################################################################