# Cauchy M estimator of regression in R

Was wondering if anyone knows of an R package to estimate the Cauchy-M estimator of regression (see for example the end of this section, but with simultaneous estimation of the scale parameter as in section 2 of (1)).

(1) Mizera, I. Müller, C. H. (2002). Breakdown points of Cauchy regression-scale estimators, Statistics & Probability Letters, Volume 57, Issue 1, Pages 79-89.

• MASS::rlm will let you specify a $\psi$ function, though I think the choice of scale estimate is limited to a few options. Apr 29 '13 at 11:23
• Outside of that, there are some functions in package robustbase that might let you get near to what you want, but I'm not certain of that. Apr 29 '13 at 11:41
• @Glen_b: thanks. Actually, I believe this problem is somewhat more popular in the context of Bayesian estimation. I was hoping that someone more familiar with Bayesian estimation could help but I notice smillig disagrees with that assessment. Apr 29 '13 at 11:54
• A very straight-forward approach is via the EM algorithm, as Cauchy-likelihood is a scale mixture of normal likelihoods. So you get a weighted least squares regression were the weights are adaptively chosen based on the "current" residual. Apr 29 '13 at 13:33
• @probabilityislogic: I think I understand you, but can you suggest a practical way to hack a EM-algorithm for regression to achieve this? Apr 29 '13 at 14:00

This arxiv paper may be of use to you:

Data augmentation for non-Gaussian regression models using variance-mean mixtures

Very simple to code up and implement, you could use the standard linear model function lm(y~formula,weights=w) for the "M" step (including estimation of $\sigma$) where the weights are taken from the "E" step.

Note that for a Cauchy RV we have that $(e_i|\sigma^2\omega_i)\sim N(0,\sigma^2\omega_i^{-1})$ and $\omega_i\sim \chi^2(1)$ implies $(e_i|\sigma^2)\sim Cauchy(0,\sigma^2)$. The derivative of the negative log likelihood is $\frac{2e_i}{\sigma^2+e_i^2}$ which is $f'()$ in proposition of the paper. This means that the "E" step is given as:

$$\hat{\omega}_i=\frac{2\sigma^2}{\sigma^2+\hat{r}_i^2}$$

where $\hat{r}_i=y_i-x_i^T\hat{\beta}$ is the current residual (easily obtained from lm()$residuals in the previous "M" step). The tricky part with this is to note that the second derivative of the negative log likelihood is given as $$\frac{2}{\sigma^2+e_i^2}-\frac{4e_i^2}{(\sigma^2+e_i^2)^2}$$ which mean the optimisation problem is not convex - so you may have multiple modes. EM will surely find one of them, but you need to do re-starts at different values to make sure you have the actual maximum. UPDATE As mentioned in my comment, I have added some R code demonstrating the multi-modal nature of Cauchy regression. This is especially true in just the cases where Cauchy regression is most useful - when there are outliers that influence the regression fit. This example shows this by using a point of high leverage which doesn't follow the rest of the data. #em algorithm for Cauchy regression cauchylm <- function(y,x,beta_init=NULL){ n <- length(y) if(is.null(beta_init)){ ols_reg <- lm(y~x) beta_old <- ols_reg$coefficients
sigma_old <- summary(ols_reg)$sigma }else{ beta_old <- beta_init ols_reg <- list(residuals=y-beta_init[1]-x*beta_init[2]) sigma_old <- sqrt(sum(ols_reg$residuals^2)/(n-2))
}
iter<- 1
prec <- 1
while((iter<1000)&(prec>0.00000001)){
weights <- 2*sigma_old^2/(sigma_old^2+ols_reg$residuals^2) ols_reg <- lm(y~x,weights=weights) beta <- ols_reg$coefficients
sigma <- summary(ols_reg)$sigma prec <- max(abs(c(beta,sigma)-c(beta_old,sigma_old))/(abs(c(beta_old,sigma_old))+0.001)) beta_old <- beta sigma_old <- sigma iter <- iter+1 } return(list(beta=beta,sigma=sigma,ols_reg=ols_reg,iter=iter,prec=prec)) } #this example creates a point of high leverage, showing the multimodality x<- runif(10) beta <- -1 alpha <- 1 y <- alpha + beta*x + 0.3*rnorm(10) x0 <- 4 y0 <- x0+0.3*rnorm(1) # Cauchy regression using OLS estimates as starting values reg1 <- cauchylm(c(y,y0),c(x,x0)) # Cauchy regression using true values as starting values reg2 <- cauchylm(c(y,y0),c(x,x0),c(alpha,beta)) # Cauchy regression using random values as starting values reg3 <- cauchylm(c(y,y0),c(x,x0),rnorm(2)) # Cauchy regression using random values as starting values reg4 <- cauchylm(c(y,y0),c(x,x0),c(0,0)) # Cauchy regression using large values as starting values reg5 <- cauchylm(c(y,y0),c(x,x0),c(10,10)) #Shows that EM algorithm converged to different modes data.frame( ols_start=sum(dcauchy(c(y,y0)-reg1$beta[1]-reg1$beta[2]*c(x,x0),0,reg1$sigma,log=TRUE)),
true_start=sum(dcauchy(c(y,y0)-reg2$beta[1]-reg2$beta[2]*c(x,x0),0,reg2$sigma,log=TRUE)), rnorm_start=sum(dcauchy(c(y,y0)-reg3$beta[1]-reg3$beta[2]*c(x,x0),0,reg3$sigma,log=TRUE)),
zero_start=sum(dcauchy(c(y,y0)-reg4$beta[1]-reg4$beta[2]*c(x,x0),0,reg4$sigma,log=TRUE)), large_start=sum(dcauchy(c(y,y0)-reg5$beta[1]-reg5$beta[2]*c(x,x0),0,reg5$sigma,log=TRUE))
)

#plot of the fitted lines
matplot(c(x,x0),cbind(
reg1$ols_reg$fitted,
reg2$ols_reg$fitted,
reg3$ols_reg$fitted,
reg4$ols_reg$fitted,
reg5$ols_reg$fitted),type='l',ylab="Y",xlab="X",main="Cauchy Regression")
points(c(x,x0),c(y,y0))
legend(x0-1,2,c("OLS","TRUE","RAND","ZERO","LARGE"),col=1:5,lty=1:5)

• @probilityislogic: which mean the optimisation problem is not convex -so you may have multiple modes. One question (this is actually something I was trying to elucidate). Reading On the uni-modality of the likelihood for the Cauchy distribution (Copas, 1975). Biometrika 62, 701-704. I was under the impression that, except for a pathological case of bi-modality, the solution to this problem is unique. I was wondering whether you could comment on this (sorry if this sounds off base but this question is what got me interested in the Cauchy M in the first place). Apr 29 '13 at 15:53
• This paper only deals with the case of no covariates. If you think in terms of IRLS the negative hessian has the form $\sum_iw_ix_ix_i^T$. For Cauchy the weights are positive for "good" data with small residuals and negative for "bad" data with large residuals. The weights are bounded though, so as long as there aren't too many outliers you should be fine. Apr 30 '13 at 8:16
• so the negative hessian is PSD....and the optimization problem is convex? I don't understand your last sentence: did you mean 'fine' in wrt the convexity or some other aspect of the problem? Apr 30 '13 at 8:22
• I'll post some R code giving an example. The negative hessian is PSD at the mode(s) but there may be more than one mode. Depending on where you start the algorithm you get different stationary points. May 1 '13 at 7:57