# FILENAME: WinBUGSmodel.txt
# joint model for linear regression with missing responses

model
	{  
		for (i in 1:N) {
			# Model of Interest
			y[i]~dnorm(mu[i],tau)
			mu[i]<-beta0+beta1*x[i]

			# Missingness model
			m[i]~dbern(p[i])
			logit(p[i])<-theta0+theta1*(y[i]-mean.y)

			linkp[i]<-theta0+theta1*(y[i]-mean.y)
			
			# log likelihood contribution of individuals
			logL.MoI[i]<- log(2*3.14159265/tau)+(pow((y[i]-(beta0+beta1*x[i])),2)*tau)
			logL.MoM[i]<-(m[i]*log(p[i]))+((1-m[i])*log(1-p[i])) 
			
			# Observed individuals only
			logL.obs[i]<-logL.MoI[i]*(1-m[i])
			
		}
		
		# priors		
		beta0~dnorm(0,0.00000001)
		beta1~dnorm(0,0.00000001)
		tau~dgamma(0.001,0.001)
		sigma<-1/sqrt(tau)
		theta0~dlogis(0,1)
		theta1~dnorm(0,1.48)  # limit to 5 fold change
		
		D.MoI<- sum(logL.MoI[])  # calculate Deviance of model of interest
		D.MoM<- -2*sum(logL.MoM[])  # calculate Bernoulli Deviance (missingness model)
		D.MoI.obs<-sum(logL.obs[])
		
	}