# 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?

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.

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.

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 – Frans Rodenburg Nov 6 '18 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? – Ben Bolker Nov 6 '18 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. – Confounded Nov 6 '18 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. – Glen_b Nov 7 '18 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. – Confounded Nov 7 '18 at 14:41