Regression with t-distributed errors and MASS::rlm 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)
 A: 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).
