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.
#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")
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:
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.