I have some data that I'm fitting a multiple regression to, with the twist that the error distribution is t (with user-defined degrees of freedom) instead of Gaussian. I've been coding up my own function to do this, but recently I saw some posts that the rlm function in the MASS package has this ability already. Is this true? I have the accompanying book* and it mentions a number of M-estimators with Huber's as the default. I don't see how any of them correspond to the t distribution.

(Background: the t is because I'm interested in simulating the distribution of the response, not just estimating its conditional mean. From looking at the residuals, it appears that a suitably scaled t distribution with 4-5 df is a reasonable fit, and certainly much better than the Gaussian. It's only tangentially related to robustness considerations.)

* Modern Applied Statistics with S 4th Ed, Venables & Ripley (2002)


1 Answer 1


It looks to me like you can supply your own psi function.

The already-supplied psi-functions are just simple functions (try MASS:::psi.huber and MASS:::psi.hampel for example). They're rather neatly set up in that the same function supplies both the psi function and its derivative, depending on whether they're called with deriv=0 (the default) or not. Note that each function takes different parameters.

I think you could just supply your own psi function (called psi.t say, with df as one of the parameters). Or if you think some particular value - say 5 df - is always going to be adequate, you could write an even simpler function with df hard-coded.

Note that with the t-distribution you have a redescending psi-function, so the issues mentioned in the help in relation to the Hampel and bisquare psi functions also apply to yours (i.e. multiple local minima). You might want to try starting at the solution for a Huber, perhaps after choosing a value of $k$ that describes the central part of the psi function for a $t_4$ or $t_5$ reasonably well.

For a parametric model, $\psi(x) = -\log(f(x))=-f'(x)/f(x)$. The $\psi-$function for a $t_\nu$ is reasonably simple; Wikipedia gives it as

$$\psi(x) = \frac{x}{x^2+\nu}$$

Consequently, assuming I made no errors, $\psi'(x) = \frac{x^2+\nu - x\cdot 2x}{(x^2+\nu)^2}= \frac{\nu - x^2}{(x^2+\nu)^2}$.

I don't think you get very much choice on the variance parameter though; I think you would have to settle for a robust scale (not that this is necessarily a bad thing; see the information relating to scale.est and method="MM" for what options there are). Alternatively, if you really want to use the t-assumption for finding the scale as well, you could instead use the various functions for maximizing likelihood to attempt an ML solution to the whole problem (perhaps after using rlm to get essentially all the way there for the mean-function).

  • $\begingroup$ Eww, looks like I'll have to do maths and stuff to define psi.t. So far I've been getting away with just throwing dt into optim and letting it do my thinking for me. Thanks though. $\endgroup$
    – Hong Ooi
    Oct 6, 2014 at 2:09
  • $\begingroup$ It's pretty simple though, it's just $-f'(x)/f(x)$. Indeed it looks like Wikipedia has it for you already. The first derivative of that is pretty straightforward. I'll update my answer. $\endgroup$
    – Glen_b
    Oct 6, 2014 at 2:25
  • $\begingroup$ Done. Your problem now is at worst one of checking the algebra (though even that can be avoided by using numerical computation to check a few values of $\psi$ and $\psi'$ against a numerical computation of $-\log(f(x))$ and a numerical derivative. Or you might use the $D$ function in R to deal with taking the derivative). $\endgroup$
    – Glen_b
    Oct 6, 2014 at 2:36
  • $\begingroup$ Sweet, thanks! I just tried it and got a result very similar to directly maximising the likelihood, but ~10x as fast. $\endgroup$
    – Hong Ooi
    Oct 6, 2014 at 5:36
  • 1
    $\begingroup$ That's true, but I've forgotten my multivariate calculus... optim expects a vector-valued function which is the derivative of the log-likelihood wrt each of the betas. With numerical derivatives it still takes only a couple of seconds though (and rlm takes 0.2 seconds). I might ask another question about that. $\endgroup$
    – Hong Ooi
    Oct 6, 2014 at 7:39

Your Answer

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

Not the answer you're looking for? Browse other questions tagged or ask your own question.