4
$\begingroup$

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). $$

$\endgroup$
  • 1
    $\begingroup$ You can make Greek letters bold using \boldsymbol $\endgroup$ – Frans Rodenburg Nov 6 '18 at 13:10
  • 1
    $\begingroup$ 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? $\endgroup$ – Ben Bolker Nov 6 '18 at 23:15
  • $\begingroup$ @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. $\endgroup$ – Confounded Nov 6 '18 at 23:39
  • 1
    $\begingroup$ (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. $\endgroup$ – Glen_b Nov 7 '18 at 1:20
  • 1
    $\begingroup$ @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. $\endgroup$ – Confounded Nov 7 '18 at 14:41

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.