betareg fails to converge is dependent variable has more than 8 digits

Question

I am having problems using the betareg package. When I try to run the model m1, the model does not converge. However, when I round y to 8 digits, it works well (127 iterations). See below for the code. Does anybody know why this tiny change as such a big influence on the model?
See code below.
Thanks a lot!

#install.packages("betareg", dependencies=TRUE)
library(betareg)
data <- data.frame("x" = c(194, 194, 194,  73,  73,  73, 105, 105, 105, 222, 222, 222,   0,   0,   0,   0,   0,   0),
                   "y" = c(0.9717500000, 0.9191161111, 0.9456172222, 0.0464116667, 0.0413683333, 0.034105555, 0.9178222222, 0.9661872222, 0.9619844444,
                           0.9576777778, 0.9710794444, 0.9562516667, 0.0277777778, 0.0277777778, 0.0277777778, 0.0277777778, 0.0277777778, 0.0277777778))


library(betareg)

m1 <-  betareg(formula = y ~ x, data = data, link = "cauchit", type = "ML")

m2 <-  betareg(formula = round(y,8) ~ x, data = data, link = "cauchit", type = "ML")

Answer 1

Nice MCVE. Since the point at which this infelicity occurs is exactly the point where the y values will differ from the values where it does not generate the warning by an amount comparable to the default value for fstol in betareg.control, I reasoned that there is some test of the difference of y versus an estimate of y that is not being met. I tried your code with the addition of a somewhat larger value for

fstol ::
numeric tolerance for convergence in (quasi) Fisher scoring.

... , control=betareg.control(fstol =1e-7)

And no warning appears. Since it doesn't throw an error I think you should call it an "infelicity" rather than a "bug" when you report this to the package author.

maintainer("betareg")
[1] "Achim Zeileis <Achim.Zeileis@R-project.org>"

He should be able to fix the test that generates this warning. (I'm not able to figure out which of those conditions are being met or unmet.) Here's the testing code>

if ((fsmaxit == 0 & opt$convergence > 0) | iter >= fsmaxit) {
    converged <- FALSE
    warning("optimization failed to converge")

I tried playing around with other values for fsmaxit in the betareg.control values, but wasn't able to narrow down the problem yet. I wondered if the second nested term in the logical tests should be:

opt$convergence > fstol

.... but didn't think that it was likely that R's optim would work that way.

Answer 2

I don't think that this really has to do something with with the rounding. It's just that this is a small data set where y switches from very close to zero to very close to 1 if x is less or more than 100. Hence, the cauchit link is a bit challenging to optimize here. With the default settings it almost converges but any of the following modifications leads to successful convergence:

• use a log-link for the precision parameter
• increase the allowed iterations for the Fisher-Scoring
• use better starting values

Additionally, relaxing the strictness of the Fisher-Scoring as suggested in a previous answer by IRTFM would also lead to successful convergence.

Out of the options above I would always recommend to try out a log-link first because it typically improves inference about the precision parameter (e.g., coverage of Wald confidence intervals is better etc.). The only reason that this is not the default in betareg() is for backward compatibility with the first version of the R package and for consistency with the original Ferrari & Cribari Neto (2004) paper. But the log-link is the default when explicitly specifying a two-part formula, here just with an intercept:

m1 <- betareg(y ~ x | 1, data = data, link = "cauchit")
summary(m1)
## Call:
## betareg(formula = y ~ x | 1, data = data, link = "cauchit")
## 
## Standardized weighted residuals 2:
##     Min      1Q  Median      3Q     Max 
## -1.6746 -1.1704 -0.4335  0.8269  1.9560 
## 
## Coefficients (mean model with cauchit link):
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -19.16562    3.62912  -5.281 1.28e-07 ***
## x             0.21393    0.03946   5.422 5.90e-08 ***
## 
## Phi coefficients (precision model with log link):
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   3.6936     0.3397   10.87   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
## 
## Type of estimator: ML (maximum likelihood)
## Log-likelihood: 31.22 on 3 Df
## Pseudo R-squared: 0.7941
## Number of iterations: 29 (BFGS) + 124 (Fisher scoring) 

Using the corresponding coefficients as starting values then also leads to a successful convergence with an identity link:

s <- coef(m1)
s[3] <- exp(s[3])
m2 <- betareg(y ~ x, data = data, link = "cauchit", start = s)
summary(m2)
## Call:
## betareg(formula = y ~ x, data = data, link = "cauchit", start = s)
## 
## Standardized weighted residuals 2:
##     Min      1Q  Median      3Q     Max 
## -1.6746 -1.1704 -0.4335  0.8269  1.9560 
## 
## Coefficients (mean model with cauchit link):
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -19.16562    3.62912  -5.281 1.28e-07 ***
## x             0.21393    0.03946   5.422 5.90e-08 ***
## 
## Phi coefficients (precision model with identity link):
##       Estimate Std. Error z value Pr(>|z|)   
## (phi)    40.19      13.65   2.944  0.00324 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
## 
## Type of estimator: ML (maximum likelihood)
## Log-likelihood: 31.22 on 3 Df
## Pseudo R-squared: 0.7941
## Number of iterations: 1 (BFGS) + 2 (Fisher scoring) 

As you see, the models are virtually identical with just a couple of additional Fisher-Scoring iterations, leading to the same coefficients and inference in the mean model.