############################# # WiNBUGS code for dugongs Binary regression model # See Carlin and Louis (2008), Example 4.4 # # Compares logistic, probit, and cloglog link functions ############################# model{ for( i in 1:n) { logage[i] <- log(x[i]) useless[i] <- Y[i] # trick to avoid deleting Y from the data list z[i] ~ dbern( p[i] ) logit(p[i]) <- beta0 + beta1*(logage[i] - mean(logage[])) # logistic model # probit(p[i]) <- beta0 + beta1*(logage[i] - mean(logage[])) # probit model (numerically unstable) # p[i] <- phi(beta0 + beta1*(logage[i] - mean(logage[]))) # probit model (numerically more stable) # cloglog(p[i]) <- beta0 + beta1*(logage[i] - mean(logage[])) # cloglog model } beta0 ~ dflat() beta1 ~ dflat() } # end of WinBUGS code #Data: list(x = c( 1.0, 1.5, 1.5, 1.5, 2.5, 4.0, 5.0, 5.0, 7.0, 8.0, 8.5, 9.0, 9.5, 9.5, 10.0, 12.0, 12.0, 13.0, 13.0, 14.5, 15.5, 15.5, 16.5, 17.0, 22.5, 29.0, 31.5), z = c(0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 1, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1), Y = c(1.80, 1.85, 1.87, 1.77, 2.02, 2.27, 2.15, 2.26, 2.47, 2.19, 2.26, 2.40, 2.39, 2.41, 2.50, 2.32, 2.32, 2.43, 2.47, 2.56, 2.65, 2.47, 2.64, 2.56, 2.70, 2.72, 2.57), n = 27) # Inits: list( beta0 = -10, beta1 = -5) list( beta0 = 0, beta1= 10) list( beta0 = 10, beta1= 25) list(beta0 = -.8, beta1 = 3.3) # "easy" starting values for probit ####################### ## Now, Haijun Ma's trick for running all 3 link functions simulataneously: run all 3 independently ## (same data but different link function) in a big loop. One can now get side-by-side boxplots of the same ## parameters of different models (using the comparison tool) *and* separate DIC scores; WinBUGS will ## decompose DIC by observation variable names :) ####################### model{ for( i in 1:n) { logage[i] <- log(x[i]) useless[i] <- Y[i] # trick to avoid deleting Y from the data list z1[i] <- z[i] # replicate data z2[i] <- z[i] z3[i] <- z[i] z1[i] ~ dbern( p1[i]) z2[i] ~ dbern( p2[i]) z3[i] ~ dbern( p3[i]) logit(p1[i]) <- beta0[1] + beta1[1]*(logage[i] - mean(logage[])) # logistic model p2[i] <- phi(beta0[2] + beta1[2]*(logage[i] - mean(logage[]))) # probit model (numerically stable version) cloglog(p3[i]) <- beta0[3] + beta1[3]*(logage[i] - mean(logage[])) # cloglog model } for(k in 1:3){ beta0[k] ~ dflat() beta1[k] ~ dflat() } } # end of WinBUGS code # Data: list(x = c( 1.0, 1.5, 1.5, 1.5, 2.5, 4.0, 5.0, 5.0, 7.0, 8.0, 8.5, 9.0, 9.5, 9.5, 10.0, 12.0, 12.0, 13.0, 13.0, 14.5, 15.5, 15.5, 16.5, 17.0, 22.5, 29.0, 31.5), z = c(0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 1, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1), Y = c(1.80, 1.85, 1.87, 1.77, 2.02, 2.27, 2.15, 2.26, 2.47, 2.19, 2.26, 2.40, 2.39, 2.41, 2.50, 2.32, 2.32, 2.43, 2.47, 2.56, 2.65, 2.47, 2.64, 2.56, 2.70, 2.72, 2.57), n = 27) # Inits: list(beta0 = c(-.8,-.8,-.8), beta1 = c(3.3, 3.3, 3.3)) # "easy" starting values to help the probit again