# WinBUGS code for JARE paper "A methodology for evaluating how
# product characteristics impact choice in retail settings with many
# zero observations: an application to restaurant wine purchase." 

# Uses "zeros trick" for modeling non-standard distributions.
# Data available on request."

# Poisson Model

model {for (i in 1:N) {
	Q[i] ~ dpois(mu[i])
	ll[i] <- -mu[i] + Q[i]*log(mu[i]) - logfact(Q[i])
	lls[i] <- -Q[i] + log(pow(Q[i],Q[i])) - logfact(Q[i])
# Poisson means
	log(mu[i]) <- beta[1]*White[i] +
			beta[2]*OrChard[i] +
			beta[43]*GlassW[i] +
			beta[44]*GlassR[i] +
			part1[i] + part2[i]
		part1[i] <- beta[3]*OrPinGris[i] +
			beta[4]*CaOthW[i] +
			beta[5]*NwOthW[i] +
			beta[6]*FrOthW[i] +
			beta[7]*PriceW[i] + 
			beta[8]*LowpriceW[i] +
			beta[9]*BodyW[i] +
			beta[10]*FinW[i] +
			beta[11]*OakW[i] +
			beta[12]*RichW[i] +
			beta[13]*SpiceW[i] +
			beta[14]*Wbuttery[i] +
			beta[15]*Wcreamy[i] +
			beta[16]*Wdry[i] +
			beta[17]*Whoney[i] +
			beta[18]*Wmelon[i] +
			beta[19]*WCitrus[i] +
			beta[20]*WTreeFruit[i] +
			beta[21]*WTropFruit[i] +
			beta[22]*Red[i]
		part2[i] <- beta[23]*CaCab[i] +
			beta[24]*CaZinf[i] +
			beta[25]*OrPinNoir[i] +
			beta[26]*CaOthR[i] +
			beta[27]*NwOthR[i] +
			beta[28]*FrOthR[i] +
			beta[29]*ItOthR[i] +
			beta[30]*PriceR[i] +
			beta[31]*LowpriceR[i] +
			beta[32]*BodyR[i] +
			beta[33]*FinR[i] +
			beta[34]*OakR[i] +
			beta[35]*RichR[i] +
			beta[36]*SpiceR[i] +
			beta[37]*TannR[i] +
			beta[38]*VanR[i] +
			beta[39]*Rbcurrant[i] +
			beta[40]*Rberry[i] +
			beta[41]*Rcherry[i] +
			beta[42]*Rchoc[i]
	GlassW[i] <- Glass[i]*White[i]
	GlassR[i] <- Glass[i]*Red[i]
	notused1[i] <- Price[i]
	notused2[i] <- Lowprice[i]
	notused3[i] <- Wine[i]}
	dev <- -2*sum(ll[])
	devr <- dev + 2*sum(lls[])
# Priors
	for (j in 1:44) {beta[j] ~ dnorm(0.0, 0.01)}}

# Zero-Inflated Poisson Model

model {for (i in 1:N) {
	z[i] <- 0
	z[i] ~ dpois(phi[i])
	phi[i] <- -ll[i]
# group membership
	logit(p[i]) <- eta[1] + eta[2]*Price[i]
# log-likelihood
	ll[i] <- zero[i] * (Glass[i]*-mu[i] + 
		(1-Glass[i])*log(1 - p[i] + p[i]*exp(-mu[i]))) +
		(1-zero[i]) * (-mu[i] + Q[i]*log(mu[i]) - logfact(Q[i]) +
		(1-Glass[i])*log(p[i]))
	zero[i] <- equals(Q[i], 0)
# Poisson means
	log(mu[i]) <- beta[1]*White[i] +
			beta[2]*OrChard[i] +
			beta[43]*GlassW[i] +
			beta[44]*GlassR[i] +
			part1[i] + part2[i]
		part1[i] <- beta[3]*OrPinGris[i] +
			beta[4]*CaOthW[i] +
			beta[5]*NwOthW[i] +
			beta[6]*FrOthW[i] +
			beta[7]*PriceW[i] + 
			beta[8]*LowpriceW[i] +
			beta[9]*BodyW[i] +
			beta[10]*FinW[i] +
			beta[11]*OakW[i] +
			beta[12]*RichW[i] +
			beta[13]*SpiceW[i] +
			beta[14]*Wbuttery[i] +
			beta[15]*Wcreamy[i] +
			beta[16]*Wdry[i] +
			beta[17]*Whoney[i] +
			beta[18]*Wmelon[i] +
			beta[19]*WCitrus[i] +
			beta[20]*WTreeFruit[i] +
			beta[21]*WTropFruit[i] +
			beta[22]*Red[i]
		part2[i] <- beta[23]*CaCab[i] +
			beta[24]*CaZinf[i] +
			beta[25]*OrPinNoir[i] +
			beta[26]*CaOthR[i] +
			beta[27]*NwOthR[i] +
			beta[28]*FrOthR[i] +
			beta[29]*ItOthR[i] +
			beta[30]*PriceR[i] +
			beta[31]*LowpriceR[i] +
			beta[32]*BodyR[i] +
			beta[33]*FinR[i] +
			beta[34]*OakR[i] +
			beta[35]*RichR[i] +
			beta[36]*SpiceR[i] +
			beta[37]*TannR[i] +
			beta[38]*VanR[i] +
			beta[39]*Rbcurrant[i] +
			beta[40]*Rberry[i] +
			beta[41]*Rcherry[i] +
			beta[42]*Rchoc[i]
	GlassW[i] <- Glass[i]*White[i]
	GlassR[i] <- Glass[i]*Red[i]
	notused1[i] <- Lowprice[i]
	notused2[i] <- Wine[i]}
	dev <- -2*sum(ll[])
# priors
	for (j in 1:44) {beta[j] ~ dnorm(0.0, 0.01)}
	for (j in 1:2) {eta[j] ~ dnorm(0.0, 0.01)}}

# Hurdle Model

model {for (i in 1:N) {
	z[i] <- 0
	z[i] ~ dpois(phi[i])
	phi[i] <- -ll[i]
# group membership
	logit(p[i]) <- eta[1] + eta[2]*Price[i]
# log-likelihood
	ll[i] <- zero[i] * (Glass[i]*-mu[i] + (1-Glass[i])*log(1 - p[i])) +
		(1-zero[i]) * (-mu[i] + Q[i]*log(mu[i]) - logfact(Q[i]) + 
		(1-Glass[i])*(log(p[i]/(1-exp(-mu[i])))))
	zero[i] <- equals(Q[i], 0)
# Poisson means
	log(mu[i]) <- beta[1]*White[i] +
			beta[2]*OrChard[i] +
			beta[43]*GlassW[i] +
			beta[44]*GlassR[i] +
			part1[i] + part2[i]
		part1[i] <- beta[3]*OrPinGris[i] +
			beta[4]*CaOthW[i] +
			beta[5]*NwOthW[i] +
			beta[6]*FrOthW[i] +
			beta[7]*PriceW[i] + 
			beta[8]*LowpriceW[i] +
			beta[9]*BodyW[i] +
			beta[10]*FinW[i] +
			beta[11]*OakW[i] +
			beta[12]*RichW[i] +
			beta[13]*SpiceW[i] +
			beta[14]*Wbuttery[i] +
			beta[15]*Wcreamy[i] +
			beta[16]*Wdry[i] +
			beta[17]*Whoney[i] +
			beta[18]*Wmelon[i] +
			beta[19]*WCitrus[i] +
			beta[20]*WTreeFruit[i] +
			beta[21]*WTropFruit[i] +
			beta[22]*Red[i]
		part2[i] <- beta[23]*CaCab[i] +
			beta[24]*CaZinf[i] +
			beta[25]*OrPinNoir[i] +
			beta[26]*CaOthR[i] +
			beta[27]*NwOthR[i] +
			beta[28]*FrOthR[i] +
			beta[29]*ItOthR[i] +
			beta[30]*PriceR[i] +
			beta[31]*LowpriceR[i] +
			beta[32]*BodyR[i] +
			beta[33]*FinR[i] +
			beta[34]*OakR[i] +
			beta[35]*RichR[i] +
			beta[36]*SpiceR[i] +
			beta[37]*TannR[i] +
			beta[38]*VanR[i] +
			beta[39]*Rbcurrant[i] +
			beta[40]*Rberry[i] +
			beta[41]*Rcherry[i] +
			beta[42]*Rchoc[i]
	GlassW[i] <- Glass[i]*White[i]
	GlassR[i] <- Glass[i]*Red[i]
	notused1[i] <- Lowprice[i]
	notused2[i] <- Wine[i]}
	dev <- -2*sum(ll[])
# priors
	for (j in 1:44) {beta[j] ~ dnorm(0.0, 0.01)I(-5, 5)}
	for (j in 1:2) {eta[j] ~ dnorm(0.0, 0.01)}}

# Negative Binomial Model

model {for (i in 1:N) {
	Q[i] ~ dnegbin(pstar[i],alphastar)
	pstar[i] <- alphastar/(alphastar+mu[i])
# log-likelihood
	ll[i] <- loggam(Q[i]+alphastar) + alphastar*log(pstar[i]) + 
		Q[i]*log(1-pstar[i]) - logfact(Q[i]) - loggam(alphastar)
# Poisson means
	log(mu[i]) <- beta[1]*White[i] +
			beta[2]*OrChard[i] +
			beta[43]*GlassW[i] +
			beta[44]*GlassR[i] +
			part1[i] + part2[i]
		part1[i] <- beta[3]*OrPinGris[i] +
			beta[4]*CaOthW[i] +
			beta[5]*NwOthW[i] +
			beta[6]*FrOthW[i] +
			beta[7]*PriceW[i] + 
			beta[8]*LowpriceW[i] +
			beta[9]*BodyW[i] +
			beta[10]*FinW[i] +
			beta[11]*OakW[i] +
			beta[12]*RichW[i] +
			beta[13]*SpiceW[i] +
			beta[14]*Wbuttery[i] +
			beta[15]*Wcreamy[i] +
			beta[16]*Wdry[i] +
			beta[17]*Whoney[i] +
			beta[18]*Wmelon[i] +
			beta[19]*WCitrus[i] +
			beta[20]*WTreeFruit[i] +
			beta[21]*WTropFruit[i] +
			beta[22]*Red[i]
		part2[i] <- beta[23]*CaCab[i] +
			beta[24]*CaZinf[i] +
			beta[25]*OrPinNoir[i] +
			beta[26]*CaOthR[i] +
			beta[27]*NwOthR[i] +
			beta[28]*FrOthR[i] +
			beta[29]*ItOthR[i] +
			beta[30]*PriceR[i] +
			beta[31]*LowpriceR[i] +
			beta[32]*BodyR[i] +
			beta[33]*FinR[i] +
			beta[34]*OakR[i] +
			beta[35]*RichR[i] +
			beta[36]*SpiceR[i] +
			beta[37]*TannR[i] +
			beta[38]*VanR[i] +
			beta[39]*Rbcurrant[i] +
			beta[40]*Rberry[i] +
			beta[41]*Rcherry[i] +
			beta[42]*Rchoc[i]
	GlassW[i] <- Glass[i]*White[i]
	GlassR[i] <- Glass[i]*Red[i]
	notused1[i] <- Price[i]
	notused2[i] <- Lowprice[i]
	notused3[i] <- Wine[i]}
	dev <- -2*sum(ll[])
# Priors
	for (j in 1:44) {beta[j] ~ dnorm(0.0, 0.01)}
	alphastar <- exp(logalpha)
	logalpha ~ dnorm(0.0, 0.01)}

# Negative Binomial Model using Random Effects

model {for (i in 1:N) {
	Q[i] ~ dpois(mustar[i])
	mustar[i] <- r[i]*mu[i]
	r[i] ~ dgamma(alphastar, alphastar)
# log-likelihood
	ll[i] <- loggam(Q[i]+alphastar) + alphastar*log(pstar[i]) + 
		Q[i]*log(1-pstar[i]) - logfact(Q[i]) - loggam(alphastar)
	pstar[i] <- alphastar/(alphastar+mu[i])
	llc[i] <- -mustar[i] + Q[i]*log(mustar[i]) - logfact(Q[i])
# Poisson means
	log(mu[i]) <- beta[1]*White[i] +
			beta[2]*OrChard[i] +
			beta[43]*GlassW[i] +
			beta[44]*GlassR[i] +
			part1[i] + part2[i]
		part1[i] <- beta[3]*OrPinGris[i] +
			beta[4]*CaOthW[i] +
			beta[5]*NwOthW[i] +
			beta[6]*FrOthW[i] +
			beta[7]*PriceW[i] + 
			beta[8]*LowpriceW[i] +
			beta[9]*BodyW[i] +
			beta[10]*FinW[i] +
			beta[11]*OakW[i] +
			beta[12]*RichW[i] +
			beta[13]*SpiceW[i] +
			beta[14]*Wbuttery[i] +
			beta[15]*Wcreamy[i] +
			beta[16]*Wdry[i] +
			beta[17]*Whoney[i] +
			beta[18]*Wmelon[i] +
			beta[19]*WCitrus[i] +
			beta[20]*WTreeFruit[i] +
			beta[21]*WTropFruit[i] +
			beta[22]*Red[i]
		part2[i] <- beta[23]*CaCab[i] +
			beta[24]*CaZinf[i] +
			beta[25]*OrPinNoir[i] +
			beta[26]*CaOthR[i] +
			beta[27]*NwOthR[i] +
			beta[28]*FrOthR[i] +
			beta[29]*ItOthR[i] +
			beta[30]*PriceR[i] +
			beta[31]*LowpriceR[i] +
			beta[32]*BodyR[i] +
			beta[33]*FinR[i] +
			beta[34]*OakR[i] +
			beta[35]*RichR[i] +
			beta[36]*SpiceR[i] +
			beta[37]*TannR[i] +
			beta[38]*VanR[i] +
			beta[39]*Rbcurrant[i] +
			beta[40]*Rberry[i] +
			beta[41]*Rcherry[i] +
			beta[42]*Rchoc[i]
	GlassW[i] <- Glass[i]*White[i]
	GlassR[i] <- Glass[i]*Red[i]
	notused1[i] <- Price[i]
	notused2[i] <- Lowprice[i]
	notused3[i] <- Wine[i]}
	dev <- -2*sum(ll[])
	devc <- -2*sum(llc[])
# Priors
	for (j in 1:44) {beta[j] ~ dnorm(0.0, 0.01)}
	alphastar <- exp(logalpha)
	logalpha ~ dnorm(0.0, 0.01)}

# Zero-Inflated Negative Binomial Model

model {for (i in 1:N) {
	z[i] <- 0
	z[i] ~ dpois(phi[i])
	phi[i] <- -ll[i]
# group membership
	logit(p[i]) <- eta[1] + eta[2]*Price[i]
# log-likelihood
	ll[i] <- zero[i] * (Glass[i]*alphastar*log(pstar[i]) + 
		(1-Glass[i])*log(1 - p[i] + p[i]*pow(pstar[i],alphastar))) +
		(1-zero[i]) * (loggam(Q[i]+alphastar) + alphastar*log(pstar[i]) +
		Q[i]*log(1-pstar[i]) - logfact(Q[i]) - loggam(alphastar) +
		(1-Glass[i])*log(p[i]))
	pstar[i] <- alphastar/(alphastar+mu[i])
	zero[i] <- equals(Q[i], 0)
# Poisson means
	log(mu[i]) <- beta[1]*White[i] +
			beta[2]*OrChard[i] +
			beta[43]*GlassW[i] +
			beta[44]*GlassR[i] +
			part1[i] + part2[i]
		part1[i] <- beta[3]*OrPinGris[i] +
			beta[4]*CaOthW[i] +
			beta[5]*NwOthW[i] +
			beta[6]*FrOthW[i] +
			beta[7]*PriceW[i] + 
			beta[8]*LowpriceW[i] +
			beta[9]*BodyW[i] +
			beta[10]*FinW[i] +
			beta[11]*OakW[i] +
			beta[12]*RichW[i] +
			beta[13]*SpiceW[i] +
			beta[14]*Wbuttery[i] +
			beta[15]*Wcreamy[i] +
			beta[16]*Wdry[i] +
			beta[17]*Whoney[i] +
			beta[18]*Wmelon[i] +
			beta[19]*WCitrus[i] +
			beta[20]*WTreeFruit[i] +
			beta[21]*WTropFruit[i] +
			beta[22]*Red[i]
		part2[i] <- beta[23]*CaCab[i] +
			beta[24]*CaZinf[i] +
			beta[25]*OrPinNoir[i] +
			beta[26]*CaOthR[i] +
			beta[27]*NwOthR[i] +
			beta[28]*FrOthR[i] +
			beta[29]*ItOthR[i] +
			beta[30]*PriceR[i] +
			beta[31]*LowpriceR[i] +
			beta[32]*BodyR[i] +
			beta[33]*FinR[i] +
			beta[34]*OakR[i] +
			beta[35]*RichR[i] +
			beta[36]*SpiceR[i] +
			beta[37]*TannR[i] +
			beta[38]*VanR[i] +
			beta[39]*Rbcurrant[i] +
			beta[40]*Rberry[i] +
			beta[41]*Rcherry[i] +
			beta[42]*Rchoc[i]
	GlassW[i] <- Glass[i]*White[i]
	GlassR[i] <- Glass[i]*Red[i]
	notused1[i] <- Lowprice[i]
	notused2[i] <- Wine[i]}
	dev <- -2*sum(ll[])
# priors
	for (j in 1:44) {beta[j] ~ dnorm(0.0, 0.01)}
	for (j in 1:2) {eta[j] ~ dnorm(0.0, 0.01)}
	alphastar <- exp(logalpha)
	logalpha ~ dnorm(0.0, 0.01)I(-5, 5)}

# ZIP White, ZINB Red

model {for (i in 1:N) {
	z[i] <- 0
	z[i] ~ dpois(phi[i])
	phi[i] <- -ll[i]
# group membership
	logit(p[i]) <- eta[1] + eta[2]*Price[i]
# log-likelihood
	ll[i] <- llw[i] + llr[i]
	llw[i] <- White[i] * (zero[i] * (Glass[i]*-mu[i] + 
		(1-Glass[i])*log(1 - p[i] + p[i]*exp(-mu[i]))) +
		(1-zero[i]) * (-mu[i] + Q[i]*log(mu[i]) - logfact(Q[i]) +
		(1-Glass[i])*log(p[i])))
	llr[i] <- Red[i] * (zero[i] * (Glass[i]*alphastar*log(pstar[i]) + 
		(1-Glass[i])*log(1 - p[i] + p[i]*pow(pstar[i],alphastar))) +
		(1-zero[i]) * (loggam(Q[i]+alphastar) + alphastar*log(pstar[i]) +
		Q[i]*log(1-pstar[i]) - logfact(Q[i]) - loggam(alphastar) +
		(1-Glass[i])*log(p[i])))
	pstar[i] <- alphastar/(alphastar+mu[i])
	zero[i] <- equals(Q[i], 0)
# Poisson means
	log(mu[i]) <- beta[1]*White[i] +
			beta[2]*OrChard[i] +
			beta[43]*GlassW[i] +
			beta[44]*GlassR[i] +
			part1[i] + part2[i]
		part1[i] <- beta[3]*OrPinGris[i] +
			beta[4]*CaOthW[i] +
			beta[5]*NwOthW[i] +
			beta[6]*FrOthW[i] +
			beta[7]*PriceW[i] + 
			beta[8]*LowpriceW[i] +
			beta[9]*BodyW[i] +
			beta[10]*FinW[i] +
			beta[11]*OakW[i] +
			beta[12]*RichW[i] +
			beta[13]*SpiceW[i] +
			beta[14]*Wbuttery[i] +
			beta[15]*Wcreamy[i] +
			beta[16]*Wdry[i] +
			beta[17]*Whoney[i] +
			beta[18]*Wmelon[i] +
			beta[19]*WCitrus[i] +
			beta[20]*WTreeFruit[i] +
			beta[21]*WTropFruit[i] +
			beta[22]*Red[i]
		part2[i] <- beta[23]*CaCab[i] +
			beta[24]*CaZinf[i] +
			beta[25]*OrPinNoir[i] +
			beta[26]*CaOthR[i] +
			beta[27]*NwOthR[i] +
			beta[28]*FrOthR[i] +
			beta[29]*ItOthR[i] +
			beta[30]*PriceR[i] +
			beta[31]*LowpriceR[i] +
			beta[32]*BodyR[i] +
			beta[33]*FinR[i] +
			beta[34]*OakR[i] +
			beta[35]*RichR[i] +
			beta[36]*SpiceR[i] +
			beta[37]*TannR[i] +
			beta[38]*VanR[i] +
			beta[39]*Rbcurrant[i] +
			beta[40]*Rberry[i] +
			beta[41]*Rcherry[i] +
			beta[42]*Rchoc[i]
	GlassW[i] <- Glass[i]*White[i]
	GlassR[i] <- Glass[i]*Red[i]
	notused1[i] <- Lowprice[i]
	notused2[i] <- Wine[i]}
	dev <- -2*sum(ll[])
	devw <- -2*sum(llw[])
	devr <- -2*sum(llr[])
# priors
	for (j in 1:44) {beta[j] ~ dnorm(0.0, 0.01)}
	for (j in 1:2) {eta[j] ~ dnorm(0.0, 0.01)}
	alphastar <- exp(logalpha)
	logalpha ~ dnorm(0.0, 0.01)I(-5, 5)}