How to perform Poisson regression with known scaling factor for mean at each value of explanary variables in R?

Question

The following webpage shows how to do Poisson regression.

https://stats.idre.ucla.edu/r/dae/poisson-regression/

summary(m1 <- glm(num_awards ~ prog + math, family="poisson", data=p))

The mean of the Poisson distribution is linked with the linear transform (X) of the explanatory variable b (here prog and math) via some function f().

mu = f(X b)

I would like to have a slightly modification to the relation by multiplying a known constant that is dependent on b.

mu = c(b) f(X b)

This should be trivial to do if one had code the regression from scratch. But it is non trivial to customize an existing piece of code that was not designed to do so.

I haven't found a R package that can do this. But I'd expect that somebody might have done it before.

Could anybody let me know whether there is a R function that can do this for Poisson regression?

Answer 1

Assuming that you wish to use a logarithmic link function (which is the default for what the glm function is implementing here), your proposed linking equation is:

$$\mu_i = c(\boldsymbol{\beta}) \exp(\boldsymbol{\beta}^\text{T} \mathbf{x}_{(i)}),$$

which can also be written as:

$$\log(\mu_i) = h(\boldsymbol{\beta}) + \boldsymbol{\beta}^\text{T} \mathbf{x}_{(i)} \quad \quad \quad h \equiv \log c.$$

Since the first term in the regression equation is an unspecified function of the coefficient vector $\boldsymbol{\beta}$, this is a complicated nonlinear model form that is going to have a different estimator for the coefficient. The log-likelihood function for the model is:

$$\begin{align} \ell_{\mathbf{y}, \mathbf{x}}(\boldsymbol{\beta}) &= \sum_{i=1}^n \log \text{Pois}(y_i | \mu_i) \\[6pt] &= \sum_{i=1}^n \log \text{Pois}(y_i | \exp(h(\boldsymbol{\beta}) + \boldsymbol{\beta}^\text{T} \mathbf{x}_{(i)})) \\[6pt] &= \sum_{i=1}^n y_i [h(\boldsymbol{\beta}) + \boldsymbol{\beta}^\text{T} \mathbf{x}_{(i)}] - \sum_{i=1}^n \exp(h(\boldsymbol{\beta}) + \boldsymbol{\beta}^\text{T} \mathbf{x}_{(i)}) + \text{const}. \\[6pt] \end{align}$$

Assuming that $h$ is differentiable you then have partial derivatives:

$$\begin{align} \frac{\partial \ell_{\mathbf{y}, \mathbf{x}}}{\partial \beta_k} (\boldsymbol{\beta}) &= \sum_{i=1}^n \bigg( \frac{\partial h}{\partial \beta_k}(\boldsymbol{\beta}) + \mathbf{x}_{(k)} \bigg) \bigg( y_i - \exp(h(\boldsymbol{\beta}) + \boldsymbol{\beta}^\text{T} \mathbf{x}_{(i)}) \bigg), \\[6pt] \end{align}$$

which lead to an alternative form for the estimated coefficients in the model.

Answer 2

Poisson regression is usually using the log link (default in R, so the link used with the code example you gave.) So your function f is the exponential function. Applying that to mu = c(b) f(X b), we obtain log(mu) = log(c(b))+ X b, so that first term have the form off an offset, a known term included in the linear predictor.

The difference here is that this "known" constant depends on b, but in a known way. That cannot be estimated with the glm function, but in principle it can be done the same way. In each iteration of the IRLS algorithm used for estimation. In each iteration there is a provisional value of b, and that value can be used to adjust the offset term for the next iteration.

I cannot find an R package that does this, but you can probably write it yourself, piecing together calls to model.matrix, glm.fit and glm.control, the last have a maxit argument you can use.