# # The code here can be used to reproduce the Bayes summaries in Table 9.13 # #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ # Normal model with y = mu + eps # Nothing to prevent flip-flop here # model { for( i in 1 : N ) { for( j in 1 : T ) { Y[i , j] ~ dnorm(mu[i , j],eps.tau) mu[i , j] <- Dose[i]*exp(theta[i,1] + theta[i,2] - theta[i,3]) * (exp(-exp(theta[i,1])*time[i,j]) - exp(-exp(theta[i,2])*time[i,j]) )/(exp(theta[i,2])-exp(theta[i,1])) } theta[i, 1:3] ~ dmnorm(beta[1:3], Dinv[1:3, 1:3]) ke[i] <- exp(theta[i,1]) ka[i] <- exp(theta[i,2]) Cl[i] <- exp(theta[i,3]) } eps.tau <- exp(logtau) logtau ~ dflat() sigma <- 1 / sqrt(eps.tau) beta[1:3] ~ dmnorm(mean[1:3], prec[1:3, 1:3]) kemed <- exp(beta[1]) kamed <- exp(beta[2]) Clmed <- exp(beta[3]) Dinv[1:3, 1:3] ~ dwish(R[1:3, 1:3], 3) D[1:3, 1:3] <- inverse(Dinv[1:3, 1:3]) for (i in 1 : 3) {sdD[i] <- sqrt(D[i, i]) } } # # Lognormal model with log y = log mu + epsilon # Parameterization prevents flip-flop # model { for( i in 1 : N ){ for( j in 1 : T ){ Y[i , j] ~ dlnorm(mu[i , j],eps.tau) mu[i , j] <- log(Dose[i]*exp(theta[i,1])*(exp(theta[i,1])+exp(theta[i,2]))*exp(- theta[i,3]) * (exp(-exp(theta[i,1])*time[i,j]) - exp(-(exp(theta[i,1])+exp(theta[i,2]))*time[i,j]) )/exp(theta[i,2])) fitted[i,j] <- mu[i,j] resid[i,j] <- log(Y[i,j])-fitted[i,j] } theta[i, 1:3] ~ dmnorm(beta[1:3], Dinv[1:3, 1:3]) ke[i] <- exp(theta[i,1]) ka[i] <- exp(theta[i,1])+exp(theta[i,2]) Cl[i] <- exp(theta[i,3]) } eps.tau <- exp(logtau) logtau ~ dflat() sigma <- 1 / sqrt(eps.tau) beta[1:3] ~ dmnorm(mean[1:3], prec[1:3, 1:3]) kemean <- exp(beta[1]+0.5*sdD[1]) kamean <- exp(beta[1]+0.5*sdD[1])+exp(beta[2]+0.5*sdD[2]) Clmean <- exp(beta[3]+0.5*sdD[3]) kemed <- exp(beta[1]) Clmed <- exp(beta[3]) Dinv[1:3, 1:3] ~ dwish(R[1:3, 1:3], 3) D[1:3, 1:3] <- inverse(Dinv[1:3, 1:3]) for (i in 1 : 3) {sdD[i] <- sqrt(D[i, i]) } } # Data for lognormal model # list( N = 12, T = 10,Dose=c(4.02,4.4,4.53,4.4,5.86,4,4.95,4.53,3.1,5.5,4.92,5.3), Y = structure( .Data = c(2.84, 6.57, 10.50, 9.66, 8.58, 8.36, 7.47, 6.89, 5.94, 3.28, 1.72, 7.91, 8.31 ,8.33, 6.85, 6.08, 5.40, 4.55, 3.01, 0.90, 4.40, 6.90, 8.20, 7.80, 7.50, 6.20, 5.30, 4.90, 3.70, 1.05, 1.89, 4.60,8.60, 8.38, 7.54, 6.88, 5.78, 5.33, 4.19, 1.15, 2.02, 5.63, 11.40, 9.33, 8.74, 7.56, 7.09, 5.90, 4.37, 1.57, 1.29, 3.08, 6.44, 6.32, 5.53, 4.94, 4.02, 3.46, 2.78, 0.92, 0.85, 2.35, 5.02, 6.58, 7.09, 6.66, 5.25, 4.39, 3.53, 1.15, 3.05, 3.05, 7.31, 7.56, 6.59, 5.88,4.73, 4.57, 3.00, 1.25, 7.37, 9.03, 7.14 , 6.33 , 5.66 , 5.67, 4.24,4.11, 3.16, 1.12, 2.89 , 5.22, 6.41, 7.83 ,10.21, 9.18, 8.02, 7.14, 5.68, 2.42, 4.86, 7.24, 8.00, 6.81, 5.87, 5.22 , 4.45, 3.62 , 2.69, 0.86 , 1.25 , 3.96 , 7.82 , 9.72 , 9.75 , 8.57, 6.59, 6.11, 4.57 , 1.17), .Dim = c(12,10)), time = structure( .Data = c(0.25 , 0.57 , 1.12 , 2.02, 3.82 , 5.10 , 7.03 , 9.05, 12.12 ,24.37 , 0.27, 0.52, 1.00, 1.92, 3.50, 5.02, 7.03, 9.00, 12.00, 24.30, 0.27, 0.58 , 1.02, 2.02 , 3.62 , 5.08, 7.07, 9.00 ,12.15, 24.17 , 0.35 , 0.60 , 1.07 , 2.13 , 3.50 , 5.02 , 7.02 , 9.02 ,11.98 ,24.65 ,0.30 , 0.52 , 1.00, 2.02 , 3.50, 5.02, 7.02 , 9.10, 12.00 ,24.35, 0.27, 0.58, 1.15, 2.03, 3.57, 5.00 , 7.00, 9.22, 12.10 ,23.85 , 0.25, 0.50 , 1.02 , 2.02 , 3.48, 5.00, 6.98 , 9.00 ,12.05 ,24.22 , 0.25 , 0.52 , 0.98 , 2.02 , 3.53 , 5.05, 7.15, 9.07, 12.10 ,24.12 ,0.30 , 0.63 , 1.05 , 2.02 , 3.53 , 5.02 , 7.17 , 8.80 ,11.60 ,24.43 , 0.37 , 0.77 , 1.02 , 2.05 , 3.55 , 5.05 , 7.08 , 9.38 , 12.10, 23.70, 0.25 , 0.50 , 0.98 , 1.98 , 3.60 , 5.02, 7.03, 9.03 ,12.12 , 24.08 , 0.25 , 0.50 , 1.00 , 2.00, 3.52 , 5.07 , 7.07, 9.03, 12.05 ,24.15 ), .Dim = c(12,10)),mean = c(0, 0, 0),R = structure(.Data = c(0.2, 0, 0, 0, 0.2, 0, 0, 0, 0.2), .Dim = c(3, 3)), prec = structure(.Data = c(1.0E-6, 0, 0, 0, 1.0E-6, 0, 0, 0, 1.0E-6), .Dim = c(3, 3))) # # Intial estimates for lognormal model # list(theta = structure( .Data = c(-2.2,0,3,-2.2,0,3,-2.2,0,3,-2.2,0,3,-2.2, 0,3,-2.2,0,3,-2.2,0,3,-2.2,0,3,-2.2,0,3,-2.2,0,3,-2.2,0,3,-2.2,0,3), .Dim = c(12, 3)), beta = c(-2, .1, -3), Dinv = structure(.Data = c(1, 0, 0, 0, 1, 0, 0, 0, 1), .Dim = c(3, 3)), logtau = 0) # #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ # Power model with var(Y) = sigma^2_0 + sigma^2_1 x mu^gamma # with gamma estimated # model { for( i in 1 : N ) { for( j in 1 : T ) { eps.tau[i,j] <- 1/(sigma0*sigma0 + sigma1*sigma1*pow(mu[i,j],gamma)) Y[i , j] ~ dnorm(mu[i , j],eps.tau[i,j]) mu[i , j] <- Dose[i]*exp(theta[i,1] + theta[i,2] - theta[i,3]) * (exp(-exp(theta[i,1])*time[i,j]) - exp(-exp(theta[i,2])*time[i,j]) )/(exp(theta[i,2])-exp(theta[i,1])) fitted[i,j] <- mu[i,j] resid[i,j] <- log(Y[i,j])-fitted[i,j] } theta[i, 1:3] ~ dmnorm(beta[1:3], Dinv[1:3, 1:3]) ke[i] <- exp(theta[i,1]) ka[i] <- exp(theta[i,2]) Cl[i] <- exp(theta[i,3]) } sigma0 ~ dunif(0,2) sigma1 ~ dunif(0,2) gamma ~ dunif(0,2) beta[1:3] ~ dmnorm(mean[1:3], prec[1:3, 1:3]) kemed <- exp(beta[1]) kamed <- exp(beta[2]) Clmed <- exp(beta[3]) Dinv[1:3, 1:3] ~ dwish(R[1:3, 1:3], 3) D[1:3, 1:3] <- inverse(Dinv[1:3, 1:3]) for (i in 1 : 3) {sdD[i] <- sqrt(D[i, i]) } } # # Data for power model # list( N = 12, T = 11,Dose=c(4.02,4.4,4.53,4.4,5.86,4,4.95,4.53,3.1,5.5,4.92,5.3), Y = structure( .Data = c(0.74, 2.84, 6.57, 10.50, 9.66, 8.58, 8.36, 7.47, 6.89, 5.94, 3.28, 0.00, 1.72, 7.91, 8.31 ,8.33, 6.85, 6.08, 5.40, 4.55, 3.01, 0.90, 0.00, 4.40, 6.90, 8.20, 7.80, 7.50, 6.20, 5.30, 4.90, 3.70, 1.05, 0.00, 1.89, 4.60,8.60, 8.38, 7.54, 6.88, 5.78, 5.33, 4.19, 1.15, 0.00, 2.02, 5.63, 11.40, 9.33, 8.74, 7.56, 7.09, 5.90, 4.37, 1.57, 0.00, 1.29, 3.08, 6.44, 6.32, 5.53, 4.94, 4.02, 3.46, 2.78, 0.92, 0.15, 0.85, 2.35, 5.02, 6.58, 7.09, 6.66, 5.25, 4.39, 3.53, 1.15, 0.00, 3.05, 3.05, 7.31, 7.56, 6.59, 5.88,4.73, 4.57, 3.00, 1.25, 0.00, 7.37, 9.03, 7.14 , 6.33 , 5.66 , 5.67, 4.24,4.11, 3.16, 1.12, 0.24 , 2.89 , 5.22, 6.41, 7.83 ,10.21, 9.18, 8.02, 7.14, 5.68, 2.42, 0.00, 4.86, 7.24, 8.00, 6.81, 5.87, 5.22 , 4.45, 3.62 , 2.69, 0.86 , 0.00, 1.25 , 3.96 , 7.82 , 9.72 , 9.75 , 8.57, 6.59, 6.11, 4.57 , 1.17), .Dim = c(12,11)), time = structure( .Data = c(0.00, 0.25 , 0.57 , 1.12 , 2.02, 3.82 , 5.10 , 7.03 , 9.05, 12.12 ,24.37 , 0.00 , 0.27, 0.52, 1.00, 1.92, 3.50, 5.02, 7.03, 9.00, 12.00, 24.30, 0.00, 0.27, 0.58 , 1.02, 2.02 , 3.62 , 5.08, 7.07, 9.00 ,12.15, 24.17 , 0.00 , 0.35 , 0.60 , 1.07 , 2.13 , 3.50 , 5.02 , 7.02 , 9.02 ,11.98 ,24.65 , 0.00 ,0.30 , 0.52 , 1.00, 2.02 , 3.50, 5.02, 7.02 , 9.10, 12.00 ,24.35, 0.00, 0.27, 0.58, 1.15, 2.03, 3.57, 5.00 , 7.00, 9.22, 12.10 ,23.85 , 0.00 , 0.25, 0.50 , 1.02 , 2.02 , 3.48, 5.00, 6.98 , 9.00 ,12.05 ,24.22 , 0.00 , 0.25 , 0.52 , 0.98 , 2.02 , 3.53 , 5.05, 7.15, 9.07, 12.10 ,24.12 , 0.00 ,0.30 , 0.63 , 1.05 , 2.02 , 3.53 , 5.02 , 7.17 , 8.80 ,11.60 ,24.43 , 0.00 , 0.37 , 0.77 , 1.02 , 2.05 , 3.55 , 5.05 , 7.08 , 9.38 , 12.10, 23.70, 0.00, 0.25 , 0.50 , 0.98 , 1.98 , 3.60 , 5.02, 7.03, 9.03 ,12.12 , 24.08 , 0.00 , 0.25 , 0.50 , 1.00 , 2.00, 3.52 , 5.07 , 7.07, 9.03, 12.05 ,24.15 ), .Dim = c(12,11)),mean = c(0, 0, 0),R = structure(.Data = c(0.2, 0, 0, 0, 0.2, 0, 0, 0, 0.2), .Dim = c(3, 3)), prec = structure(.Data = c(1.0E-6, 0, 0, 0, 1.0E-6, 0, 0, 0, 1.0E-6), .Dim = c(3, 3))) # # Initial points for power model # list(theta = structure( .Data = c(-2.2,0,3,-2.2,0,3,-2.2,0,3,-2.2,0,3,-2.2,0, 3,-2.2,0,3,-2.2,0,3,-2.2,0,3,-2.2,0,3,-2.2,0,3,-2.2,0,3,-2.2,0,3), .Dim = c(12, 3)), beta =c(-2, .1, -3), Dinv = structure(.Data = c(1, 0, 0, 0, 1, 0, 0, 0, 1), .Dim = c(3, 3)),sigma0 = .5,sigma1 =.2,gamma=1)