Cauchy M estimator of regression in R

Question

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.

Answer 1

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)