Choice of algorithm for fitting Multivariate Covariance Generalized Linear Models I am using the mcglm function in R to fit a Multivariate Covariance Generalized Linear Model.
Although the help for the fit_mcglm function called by mcglm does not give the default values (giving the default values for each parameter is the usual practice on the R help files), I assume the default algorithm is Chaser's. What worries me is that using the reciprocal likelihood algorithm gives a very different fit.
How do we chose between algorithms and what are the sensible step length values for each?
 A: The chaser is the default method.
If you have different results, probably one or both algorithms did not converged. You can easily check the algorithms convergence using
plot(model, type = "algorithm")

If the algorithm converged the quasi-score functions for both regression and dispersion parameters should be very close to zero, in general < 1e-04.
Note that, if you reach the maximum number of iterations (default max_iter = 20) you have to increase the number of iterations. 
In general when using the rc method you have to change the tunning argument for some value close to 0 (example tunning = 0.001). It is important to highlight that the rc with tunning equals 0 corresponds to the chaser algorithm.
I will improve the package documentation. If you have more points for improve the documentation or the package or bugs I invite you to report at 
https://github.com/wbonat/mcglm
A simple code to illustrate how to check convergence and compare the chaser and rc methods.
require(mcglm)
x1 <- runif(100)
mu <- exp(2 - 1*x1)
y <- rpois(100, lambda = mu)
data <- data.frame(y, x1)
fit_chaser <- mcglm(linear_pred = c(y ~ x1), 
matrix_pred = list(mc_id(data)), link = "log", variance = "tweedie", 
data = data)
summary(fit_chaser)
plot(fit_chaser, type = "algorithm")
fit_rc <- mcglm(linear_pred = c(y ~ x1), matrix_pred = list(mc_id(data)),
link = "log", variance = "tweedie", data = data,
control_algorithm = list(method = "rc", tunning = 1, max_iter = 100))
summary(fit_rc)
plot(fit_rc, type = "algorithm")

cbind(coef(fit_rc), coef(fit_chaser))

All the best.
