########################### # MAC fraility example # Weibull survival model ########################### # Our dataset is unbalanced, as the number of subjects differs from site to site: # nsub=c(4,2,1,10,8,12,6,6,13,4,3), i.e. there are 4 subjects in site 1, 2 in site 2, etc. # # We use "offset" to hold the position in the dataset at which each person's data starts. # Here is a snippet of R code to calculate offset from nsub: # nsub<-c(4,2,1,10,8,12,6,6,13,4,3) # temp<-1 # offset<-rep(0,(length(nsub)+1)) # for(i in 1:length(nsub)){offset[i]<-temp; temp<-temp+nsub[i]} # offset[ (length(nsub)+1)] <-sum(nsub)+1 # Here is the dataset (used for all 5 models below): list(nsites=11, offset=c(1,5,7,8,18,26,38,44,50,63,67,70), t = c( NA,248,NA,344, NA,NA, NA, NA,64,88,NA,NA,NA,NA,NA,NA,NA, NA,NA,82,NA,NA,214,NA,262, 6,NA,76,80,202,NA,NA,NA,NA,NA,NA,NA, NA,NA,102,NA,NA,NA, NA,NA,NA,NA,NA,162, 8,NA,40,NA,NA,NA,NA,276,NA,366,NA,NA,NA, NA,NA,NA,254, NA,NA,NA), t.cen= c(74,0,272,0, 4,156, 100, 20,0,0,148,162,184,188,198,382,436, 50,64,0,186,214,0,228,0, 0,16,0,0,0,258,268,368,380,424,428,436, 32,64,0,162,182,364, 22,22,74,88,148,0, 0,16,0,120,168,174,268,0,286,0,396,466,468, 18,36,160,0, 28,70,106), drug= c(1,2,1,2, 2,1, 2, 2,2,2,2,1,1,1,1,1,1, 1,2,2,1,1,1,2,2, 1,2,1,2,2,1,1,2,1,1,2,2, 2,1,1,2,2,1, 2,1,1,1,1,2, 2,2,2,1,1,2,1,2,1,1,2,2,1, 1,1,2,2, 1,1,2 )) ## Model 1: with site specific W, b, rho ## model{ for( i in 1:nsites){ for( j in offset[i] :( offset[i+1] -1) ){ t[ j ] ~ dweib ( rho[ i ], mu[ j ] ) I(t.cen[ j ],) X[ j] <- 2* drug [ j ] -3 mu[ j ] <- exp( beta0 + beta1 * X [ j] + W[ i ] + b[i] *X[j]) } } for (i in 1:nsites) { b0star[ i ] <- beta0+ W[ i ] rho[ i ] ~ dgamma( alpha, alpha) W[ i ] ~ dnorm(0, tau.W) b[ i ] ~ dnorm(0, tau.b) } alpha <- 10 tau.W ~ dgamma( 1, 1); s2.W <- 1/tau.W; tau.b~ dgamma( 1, 1); s2.b <- 1/tau.b; beta0 ~ dnorm( 0 , 0.0001) beta1 ~ dnorm( 0, 0.0001) } # Initial values: list(rho=c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1), tau.W=1.0, tau.b=1, beta1=1.0, beta0=-10.0, W=c(0,0,0,0,0,0,0,0,0,0,0)) list(rho=c(0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5), tau.W=0.5, tau.b=0.5, beta1=0, beta0=0, W=c(-1.22201719274536, -0.084377157160674, 0.757220403084888, -0.369999602714256, -0.150438189741214, 1.02971712091728, -0.619721595818924, 1.00908379817665, -1.23491682442708, 0.79373759096893, 0.0503676864407537)) ## Model 2: with site specific W, r ## ## NOTE: This is the model summarized in the short course notes ## model{ for( i in 1:nsites){ for( j in offset[i] :( offset[i+1] -1) ){ t[ j ] ~ dweib ( rho[ i ], mu[ j ] ) I(t.cen[ j ],) X[ j] <- 2* drug [ j ] -3 mu[ j ] <- exp( beta0 + beta1 * X [ j] + W[ i ]) } } for (i in 1:nsites) { b0star[ i ] <- beta0+ W[ i ] rho[ i ] ~ dgamma( alpha, alpha) W[ i ] ~ dnorm(0, tau.W) } alpha <- 10 tau.W ~ dgamma( 1, 1); sigma.W<- 1/sqrt(tau.W) # sigma.W ~ dunif(0.01, 100); tau.W <- 1/(sigma.W * sigma.W) beta0 ~ dnorm( 0 , 0.001) beta1 ~ dnorm( 0, 0.001) relhaz <- exp(2 * beta1) } # Initial values: list(rho=c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1), tau.W=1.0, beta1=1.0, beta0=-10.0, W=c(0,0,0,0,0,0,0,0,0,0,0)) list(rho=c(0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5), tau.W=0.5, beta1=0, beta0=0, W=c(-1.22201719274536, -0.084377157160674, 0.757220403084888, -0.369999602714256, -0.150438189741214, 1.02971712091728, -0.619721595818924, 1.00908379817665, -1.23491682442708, 0.79373759096893, 0.0503676864407537)) ## Model 3: with site specific W, b ## ## NOTE: Here we use rho[i]=1 for all i ## model{ for( i in 1:nsites){ for( j in offset[i] :( offset[i+1] -1) ){ t[ j ] ~ dweib ( 1, mu[ j ] ) I(t.cen[ j ],) X[ j] <- 2* drug [ j ] -3 mu[ j ] <- exp(beta0 + beta1 * X[j] + W[i] + b[i]*X[j]) } } for (i in 1:nsites) { b0star[ i ] <- beta0+ W[ i ] W[ i ] ~ dnorm(0, tau.W) b[ i ] ~ dnorm(0, tau.b) } tau.W ~ dgamma( 1, 1) ; s2.W<- 1/tau.W; tau.b~ dgamma( 1, 1) ; s2.b <- 1/tau.b; beta0 ~ dnorm( 0 , 0.0001) beta1 ~ dnorm( 0, 0.0001) relhaz <- exp(2 * beta1) } # Initial values: list(b=c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1), tau.W=1.0, tau.b=1, beta1=1.0, beta0=-10.0, W=c(0,0,0,0,0,0,0,0,0,0,0)) list(b=c(0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5), tau.W=0.5, tau.b=0.5 beta1=0, beta0=0, W=c(-1.22201719274536, -0.084377157160674, 0.757220403084888, -0.369999602714256, -0.150438189741214, 1.02971712091728, -0.619721595818924, 1.00908379817665, -1.23491682442708, 0.79373759096893, 0.0503676864407537)) ## Model 4: with site specific W ## model{ for( i in 1:nsites){ for( j in offset[i] :( offset[i+1] -1) ){ t[ j ] ~ dweib (1, mu[ j ] ) I(t.cen[ j ],) X[ j] <- 2* drug [ j ] -3 mu[ j ] <- exp( beta0 + beta1 * X [ j] + W[ i ]) } } for (i in 1:nsites) { b0star[ i ] <- beta0+ W[ i ] W[ i ] ~ dnorm(0, tau.W) } tau.W ~ dgamma( 1, 1) ; s2.W<- 1/tau.W; beta0 ~ dnorm( 0 , 0.0001) beta1 ~ dnorm( 0, 0.0001) relhaz <- exp(2 * beta1) } # Initial values: list(tau.W=1.0, beta1=1.0, beta0=-10.0, W=c(0,0,0,0,0,0,0,0,0,0,0)) list( tau.W=0.5, beta1=0, beta0=0, W=c(-1.22201719274536, -0.084377157160674, 0.757220403084888, -0.369999602714256, -0.150438189741214, 1.02971712091728, -0.619721595818924, 1.00908379817665, -1.23491682442708, 0.79373759096893, 0.0503676864407537)) ## Model 5: no site specific random effects ## model{ for( i in 1:nsites){ for( j in offset[i] :( offset[i+1] -1) ){ t[ j ] ~ dweib (1, mu[ j ] ) I(t.cen[ j ],) X[ j] <- 2* drug [ j ] -3 mu[ j ] <- exp( beta0 + beta1 * X [ j] ) } } beta0 ~ dnorm( 0 , 0.0001) beta1 ~ dnorm( 0, 0.0001) relhaz <- exp(2 * beta1) } # Initial values: list(beta1=1.0, beta0=-10.0) list( beta1=0, beta0=0) ### DIC comparison ### ## Model 1: with site specific W, b, r ## Dbar = post.mean of -2logL; Dhat = -2LogL at post.mean of stochastic nodes Dbar Dhat pD DIC t 269.482 258.482 11.001 280.483 total 269.482 258.482 11.001 280.483 ## Model 2: with W, r ## Dbar = post.mean of -2logL; Dhat = -2LogL at post.mean of stochastic nodes Dbar Dhat pD DIC t 272.229 263.688 8.541 280.770 total 272.229 263.688 8.541 280.770 ##Model 3: with W, b ## Dbar = post.mean of -2logL; Dhat = -2LogL at post.mean of stochastic nodes Dbar Dhat pD DIC t 272.101 265.565 6.537 278.638 total 272.101 265.565 6.537 278.638 ## Model 4:with W ## Dbar = post.mean of -2logL; Dhat = -2LogL at post.mean of stochastic nodes Dbar Dhat pD DIC t 270.789 265.678 5.111 275.901 total 270.789 265.678 5.111 275.901 ## Model 5:no W ## Dbar = post.mean of -2logL; Dhat = -2LogL at post.mean of stochastic nodes Dbar Dhat pD DIC t 268.081 266.128 1.954 270.035 total 268.081 266.128 1.954 270.035