model { # Standardise xs and coefficients for (j in 1 : p) { b[j] <- beta[j] / sd(x[ , j ]) for (i in 1 : N) { z[i, j] <- (x[i, j] - mean(x[, j])) / sd(x[ , j]) } } b0 <- beta0 - b[1] * mean(x[, 1]) - b[2] * mean(x[, 2]) - b[3] * mean(x[, 3]) # Model d <- 4; # degrees of freedom for t for (i in 1 : N) { Y[i] ~ dnorm(mu[i], tau) # Y[i] ~ ddexp(mu[i], tau) # Y[i] ~ dt(mu[i], tau, d) mu[i] <- beta0 + beta[1] * z[i, 1] + beta[2] * z[i, 2] + beta[3] * z[i, 3] stres[i] <- (Y[i] - mu[i]) / sigma outlier[i] <- step(stres[i] - 2.5) + step(-(stres[i] + 2.5) ) } # Priors beta0 ~ dnorm(0, 0.00001) for (j in 1 : p) { beta[j] ~ dnorm(0, 0.00001) # coeffs independent # beta[j] ~ dnorm(0, phi) # coeffs exchangeable (ridge regression) } tau ~ dgamma(1.0E-3, 1.0E-3) phi ~ dgamma(1.0E-2,1.0E-2) # standard deviation of error distribution sigma <- sqrt(1 / tau) # normal errors # sigma <- sqrt(2) / tau # double exponential errors # sigma <- sqrt(d / (tau * (d - 2))); # t errors on d degrees of freedom }