Thanks Michael. I agree that it is possible to simplify the codes. Best,

Zhu -----Original Message----- From: Michael Weylandt <michael.weyla...@gmail.com> Sent: Monday, May 20, 2019 3:41 PM To: Wang, Zhu <wan...@uthscsa.edu> Cc: Ivan Krylov <krylov.r...@gmail.com>; R-package-devel@r-project.org Subject: Re: [R-pkg-devel] How to obtain intercept of intercept-only glm in Fortran? Negative binomial is a bit trickier since it's a two parameter family without a closed-form MLE. For the probability parameter, you can use the closed form MLE at https://en.wikipedia.org/wiki/Negative_binomial_distribution#Maximum_likelihood_estimation For the number of samples, you'll need to use an iterative method to solve the score equation (see the above link). It's probably easier to code this up yourself rather than calling into R, but if you do call into R, I'd look into using the `fitdistr` function instead of a full regression method, as demonstrated at https://stat.ethz.ch/pipermail/r-help/2012-October/338458.html Michael On Sat, May 11, 2019 at 11:10 AM Wang, Zhu <wan...@uthscsa.edu> wrote: > > Thanks Michael. > > I also need an intercept-only negative binomial model with unknown scale > parameter. So my thought was on borrowing some codes that already existed. I > think Ivan's solution is an excellent one and can be extended to other > scenarios. > > Best, > > Zhu > > On May 11, 2019, at 9:48 AM, Michael Weylandt <michael.weyla...@gmail.com> > wrote: > > On Sat, May 11, 2019 at 8:28 AM Wang, Zhu <wan...@uthscsa.edu> wrote: >> >> >> I am open to whatever suggestions but I am not aware a simple closed-form >> solution for my original question. >> > > It would help if you could clarify your original question a bit more, > but for at least the main three GLMs, there are closed form solutions, > based on means of y. Assuming canonical links, > > - Gaussian: intercept = mean(y) > - Logistic: intercept = logit(mean(y)) [Note that you have problems > here if your data is all 0 or all 1] > - Poisson: intercept = log(mean(y)) [You have problems here if your > data is all 0] > > (Check my math on these, but I'm pretty sure this is right.) > > Like I said above, this gets trickier if you add observation weights or > offsets, but the same ideas work. > > Stepping back to the statistical theory: GLMs predict the mean of y, > conditional on x. If x doesn't vary (intercept only model), then the GLM is > just predicting the mean of y and the MLE for the mean of y is exactly that > under standard GLM assumptions - the sample mean of y. > > We then just have to use the link function and its inverse to transform to > and from the observation space (where mean(y) lives) and the linear predictor > space (where the intercept term naturally lives). > > Michael ______________________________________________ R-package-devel@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-package-devel