[R-pkg-devel] How to obtain intercept of intercept-only glm in Fortran?

[Wang, Zhu]

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

Best,

Zhu 

-----Original Message-----
From: Michael Weylandt <michael.weylandt using gmail.com> 
Sent: Monday, May 20, 2019 3:41 PM
To: Wang, Zhu <wangz1 using uthscsa.edu>
Cc: Ivan Krylov <krylov.r00t using gmail.com>; R-package-devel using 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 <wangz1 using 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.weylandt using gmail.com> wrote:
>
> On Sat, May 11, 2019 at 8:28 AM Wang, Zhu <wangz1 using 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