Stacks: robust regression

Birkes and Dodge (1993) apply different regression models to the much-analysed stack-loss data of Brownlee (1965). This features 21 daily responses of stack loss y, the amount of ammonia escaping, with covariates being air flow x1, temperature x2 and acid concentration x3. Part of the data is shown below.
[stacks1]   
We first assume a linear regression on the expectation of y, with a variety of different error structures. Specifically   

      μi = β0 + β1z1i + β2z2i + β3z3i
      
      y
i ~ Normal(μi, τ)
      
      y
i ~ Double exp(μi, τ)
      
      y
i ~ t(μi, τ, d)
   
where zij = (xij - xbarj) /sd(xj) are covariates standardised to have zero mean and unit variance. β1, β2, β3 are initially given independent "noninformative" priors.

Maximum likelihood estimates for the double expontential (Laplace) distribution are essentially equivalent to minimising the sum of absolute deviations (LAD), while the other options are alternative heavy-tailed distributions. A t on 4 degrees of freedom has been chosen, although with more data it would be possible to allow this parameter also to be unknown.

We also consider the use of 'ridge regression', intended to avoid the instability due to correlated covariates. This has been shown Lindley and Smith (1972) to be equivalent to assuming the regression coefficients of the standardised covariates to be exchangeable, so that    

      
βj ~ Normal(0, φ), j = 1, 2, 3.
      
In the following example we extend the work of Birkes and Dodge (1993) by applying this ridge technique to each of the possible error distributions.

Birkes and Dodge (1993) suggest investigating outliers by examining residuals yi - μi greater than 2.5 standard deviations. We can calculate standardised residuals for each of these distributions, and create a variable outlier[i] taking on the value 1 whenever this condition is fulfilled. Mean values of outlier[i] then show the confidence with which this definition of outlier is fulfilled.

The BUGS language for all the models is shown below, with all models except the normal linear regression commented out:
   model
   {
   # Standardise x's 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
   }
Data
list(p = 3, N = 21,
   Y = c(42, 37, 37, 28, 18, 18, 19, 20, 15, 14, 14, 13, 11, 12, 8, 7, 8, 8, 9, 15, 15),
   x = structure(.Data = c( 80, 27, 89,
                           80, 27, 88,
                           75, 25, 90,
                           62, 24, 87,
                           62, 22, 87,
                           62, 23, 87,
                           62, 24, 93,
                           62, 24, 93,
                           58, 23, 87,
                           58, 18, 80,
                           58, 18, 89,
                           58, 17, 88,
                           58, 18, 82,
                           58, 19, 93,
                           50, 18, 89,
                           50, 18, 86,
                           50, 19, 72,
                           50, 19, 79,
                           50, 20, 80,
                           56, 20, 82,
                           70, 20, 91), .Dim = c(21, 3)))
Inits for chain 1
list(beta0 = 10, beta=c(0,0, 0), tau = 0.1)
   
Inits for chain 2
list(beta0 = 1.0, beta=c(1.0,1.0, 1.0), tau = 1.0)
Resultsa) Normal error

[stacks2]
b) Double exponential error

[stacks3]
c) t4 error

[stacks4]
d) Normal eror ridge regression

[stacks5]
e) Double exponential error ridge regression

[stacks6]
f) t4 error ridge regression

[stacks7]