#############################
# 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