# IRLS for truncated normal GLM

I have data for which responses fall in $$y \in [0,\infty)$$ for which, it seems, the standard GLMs based on, say, gamma or inverse-Gaussian fail since they don't allow responses with values equal to 0. So, I am trying to put together a GLM based on a truncated normal distribution with pdf given by

$$f_Y(y) = \frac{\phi(\frac{y-\mu}{\sigma})}{\sigma\Phi(\frac{\mu}{\sigma})}.$$

The log-likelihood for $$i$$-th observation is

$$l_i \sim -0.5\left(\frac{y_i-\mu_i}{\sigma}\right)^2 - \ln\left(\Phi\left(\frac{\mu_i}{\sigma}\right)\right)$$

and the score is

$$\frac{\partial l_i}{\partial \beta_j} = \left(y_i - E[y_i|\boldsymbol{\beta}] \right) \frac{1}{\sigma^2} \frac{\partial \mu_i}{\partial \eta_i} x_{ij}$$

where $$\mathbf{\eta} = \mathbf{X}\boldsymbol{\beta}$$ and

$$E[y_i|\boldsymbol{\beta}] = \mu_i + \sigma \frac{\phi(\frac{\mu_i}{\sigma})}{\Phi(\frac{\mu_i}{\sigma})}.$$

The Fisher information matrix is then

$$(\mathbf{F})_{jk} = E\left[\frac{\partial l_i}{\partial \beta_j}\frac{\partial l_i}{\partial \beta_k}|\boldsymbol{\beta}\right] = E\left[\left(y_i - E[y_i|\boldsymbol{\beta}] \right)^2|\boldsymbol{\beta}\right] \left(\frac{1}{\sigma^2} \frac{\partial \mu_i}{\partial \eta_i}\right)^2 x_{ij} x_{ik} = Var[y_i|\boldsymbol{\beta}]\left(\frac{1}{\sigma^2} \frac{\partial \mu_i}{\partial \eta_i}\right)^2 x_{ij} x_{ik}.$$

If we now define

$$\mathbf{W} = diag\left\{\sqrt{Var[y_i|\boldsymbol{\beta}]} \frac{1}{\sigma^2} \frac{\partial \mu_i}{\partial \eta_i} : i=1,...,n\right\}$$

and

$$\bar{y_i} = \frac{y_i - E[y_i|\boldsymbol{\beta}]}{\sqrt{Var[y_i|\boldsymbol{\beta}]}}$$

then

$$\mathbf{F} = (\mathbf{WX})^T\mathbf{WX} = \mathbf{X}^T\mathbf{W}\mathbf{WX}$$

and

$$\frac{\partial L}{\partial \boldsymbol{\beta}} = (\mathbf{WX})^T \bar{\mathbf{y}} = \mathbf{X}^T \mathbf{W} \bar{\mathbf{y}}$$

and Fisher scoring method becomes

$$\boldsymbol{\beta}^{m+1} = \boldsymbol{\beta}^{m} + \left(\mathbf{X}^T\mathbf{WWX}\right)^{-1}\mathbf{X}^T \mathbf{W} \bar{\mathbf{y}}$$

This, however, doesn't look like the "regular" expression for IRLS since it contains $$\mathbf{WW}$$ in the augmented projection matrix instead of the regular $$\mathbf{W}$$ that I see in various resources.

Have I got things wrong somewhere? How do I express this as an IRLS formula? And do I have to keep $$\sigma$$ or does it "drop out" somehow?

On a side note: how do you make Greek letters appear in a bold face in equations?

Add 1

If the above expression for Fisher scoring is correct, then I think the IRLS is obtained by setting

$$\mathbf{z} = \bar{\mathbf{y}} + \mathbf{WX}\boldsymbol{\beta}^{m}$$

which would give

$$\boldsymbol{\beta}^{m+1} = \left(\mathbf{X}^T\mathbf{WWX}\right)^{-1}\mathbf{X}^T \mathbf{Wz}$$

but it still leaves me with two weights matrices in the projection matrix.

Add 2

Given a link function $$E[y_i | \boldsymbol{\beta}] = G(\eta_i)$$, and the expression for the expectation of a truncated normal rv

$$E[y_i | \boldsymbol{\beta}] = \mu_i + \sigma \frac{\phi(\frac{\mu_i}{\sigma})}{\Phi(\frac{\mu_i}{\sigma})}$$

I get

$$G^{'}(\eta) = \frac{\partial \mu_i}{\partial \eta_i} \left( 1 - \frac{\mu_i}{\sigma} \frac{\phi(\frac{\mu_i}{\sigma})}{\Phi(\frac{\mu_i}{\sigma})} - \left( \frac{\phi(\frac{\mu_i}{\sigma})}{\Phi(\frac{\mu_i}{\sigma})}\right)^2\right)$$

I think the expression in brackets is proportional to the variance of a truncated normal RV when the lower boundary is 0, so we get

$$\frac{\partial \mu_i}{\partial \eta_i} = \sigma^2 \frac{G^{'}(\eta_i)}{Var[y_i|\beta]}$$

I think on substitution of this expression into the expressions for the score and Fisher information matrix should give the usual expression, just haven't worked through it yet.

Add 3

The pdf of truncated normal can be expressed as

$$f_Y(y) = \exp \left( \boldsymbol{\eta(\theta)} \centerdot \mathbf{T} (y) + A(\boldsymbol{\eta}) \right)$$

where

$$\boldsymbol{\eta (\theta)} = \left(\frac{-1}{2\sigma^2}, \frac{\mu}{\sigma^2}\right)^T$$

and

$$\mathbf{T}(y) = \left(y^2, y \right)^T$$

and $$A(\boldsymbol{\eta})$$ containing in addition to the terms found in normal distribution also another term given by

$$\ln \left(\Phi\left(\mu / \sigma \right)\right) = \ln \left( \Phi \left( \frac{-0.5 \eta_2 / \eta_1}{\sqrt{-0.5/\eta_1}} \right) \right).$$

• You can make Greek letters bold using \boldsymbol Commented Nov 6, 2018 at 13:10
• VERY short/ possibly uninformed comment: are things turning out ugly because the distribution is not in the exponential family, which is the context in which IRLS was originally developed/found to be useful? Commented Nov 6, 2018 at 23:15
• @BenBolker Thank you for taking time to read my post and for your comment. I believe, however, that the truncated normal is also in the exponential family as the only difference is a change to the so called "log-partition" function $A(\boldsymbol{\theta})$. It seems the "problem" was due to expressing the score and Fisher information matrix in terms of the location parameter $\mu$ rather than the mean. Commented Nov 6, 2018 at 23:39
• (i) Don't confuse probability with density. If your model doesn't have non-zero probability at $y_i=0$ then you shouldn't expect to see any $0$ values at all. (ii) A normal is exponential family but a truncated normal is not exponential family; you can't write it in exponential family form. Commented Nov 7, 2018 at 1:20
• @Glen_b It seems that usual GLM framework places some restrictions on the form of the log-partition function $A(\boldsymbol{\theta})$ in the exponential family. In particular, it requires that $A$ depends on the "dispersion" parameter in a particular way (multiplicative) while there is no restriction (as far as I can tell) on the form of $A$ in the general exponential family. With this restriction, the truncated normal is indeed not in this restricted definition. Commented Nov 7, 2018 at 14:41