I am modeling some spectroscopic data where the response of the instrument to the size of the input is strictly positive and non-linear. Gamma regression seems like a good choice to explain the data, since the values are naturally heteroscedastic with the spread increasing as the size of the input increases. Another notable feature is that the response value will eventually saturate and stop increasing with the input. This is an important detail, since this forces me to use a specific link function:

$\mu_j = 1 + \dfrac{\beta_0}{1 + e^{\beta_1 \times (x_j - \beta_2)}} \ j=1,\dots,m$

The model is the specified as follows:

$\beta_i \sim N(\mu=0, \sigma=1000) \ i=0,1,2$

$\epsilon \sim Exp(\lambda=100)$

$\alpha = {\left(\dfrac{\mu}{\epsilon}\right)}^{2}$ $\theta=\dfrac{\mu}{\epsilon^2}$

$y_j \sim \Gamma\left(\alpha=\alpha, \beta=\dfrac{1}{\theta}\right)$

What I am interested in is fitting such models into the data for multiple spectral lines and then taking the ratios of the lines. The reason I am doing this as the Bayesian inference is because the width of the posterior distribution of the line ratios is extremely important to me due to future use of the data in the classification analysis. This is why Gamma was chosen, since this was the only way I could force the non-negativity constraint onto the ratio value.

Now, to my issue. My data has occasional outliers in it. The value of the input is actually as estimate in and of itself. In other words, when my data states that $(x,y) = (5, 3000)$, in reality $x$ can be anywhere from $4.5$ to $5.5$. It seems to me that the model should be able to handle such uncertainty by being robust.

If the residuals were normally distributed, I would replace the Gamma distribution with a Student-t or Cauchy and run a robust inference. However, I am not sure how to proceed with the Gamma. It seems to me that there should be a way to make the model more tolerant, and I would greatly appreciate any suggestions.


  • 1
    $\begingroup$ Generally speaking, using the lowest sum of (absolute value of error) is more robust to outliers that lowest sum of (error squared). Ann outlier with an error value of 5 becomes 25 when squared. $\endgroup$ Apr 10 '19 at 22:13
  • $\begingroup$ @JamesPhillips Could you please elaborate a little bit on that? To be honest, I don't quite follow. Should I modify the link function to have a square term in it? $\endgroup$
    – udushu
    Apr 11 '19 at 0:50
  • 1
    $\begingroup$ My comment was not regarding link functions but rather fitting targets, such as in ordinary least squares regression which minimizes the sum of squared error. Because errors are squared, and outliers have large errors, a regression of this type will be dominated by the outlier's even larger squared error. If the regression minimizes the sum of absolute value of error, the outlier's large error is not squared, so it will have less of a dominating effect - that is, the regression will be more robust in the presence of outliers. $\endgroup$ Apr 11 '19 at 1:27
  • 1
    $\begingroup$ Could you write a model which makes explicit the link between the real value of $x$ and its observed value? $\endgroup$ Apr 11 '19 at 6:25
  • $\begingroup$ @RobinRyder I modified the regression to have d_j ~ N(mu=x)j, sigma=0.1*x_j), and then replace all x's in the model with d's. It didn't seem to make much of a difference, unfortunately. $\endgroup$
    – udushu
    Apr 11 '19 at 15:13

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.