0
$\begingroup$

According to Agresti(2013) pg 364-365, iterative methods such as Newton-Raphson methods, $ \begin{aligned} \beta^\text{new} &= \beta^\text{old} + (X^{T}WX)^{-1}X^{T}(V) \end{aligned} $ help to estimate the MLE of log-linear models.

These are my $Y, X , n$ matrices. The first column in $X$ matrix represents the intercept ($\beta_\text{0}$)

$$ Y= \begin{bmatrix} 1 \\ 16 \\ 30 \\ 71 \\ 102 \\ 130 \\ 133 \\ 40 \\ 4 \\ 38 \\ 119 \\ 221 \\ 259 \\ 310 \\ 226 \\ 65 \\ \end{bmatrix} ;X= \begin{bmatrix} 1 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 1 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 \\ 1 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 \\ 1 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 \\ 1 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 \\ 1 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 \\ 1 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 \\ 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 \\ 1 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 1 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 \\ 1 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 \\ 1 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 \\ 1 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 \\ 1 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 \\ 1 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 \\ 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 \\ \end{bmatrix};n= \begin{bmatrix} 172675 \\ 123065 \\ 96216 \\ 92051 \\ 72159 \\ 54722 \\ 32185 \\ 8328 \\ 181343 \\ 146207 \\ 121374 \\ [111353 \\ 83004 \\ 55932 \\ 29007 \\ 7583 \\ \end{bmatrix} $$

I am trying to fit a rate poisson model, $\log \lambda_i = \log n_i + \theta^T x$, for this data where the likelihood is $$\prod_{i=1}^N \dfrac{(n_ie^{\theta^Tx})^{y_i}}{y_i!}e^{-n_i e^{\theta^Tx}}$$

When i execute the Newton-Raphson procedure based on Agresti(2013) suggestion I am getting an error saying

Error in solve.default(t(X) %*% (V) %*% X) : 
  system is computationally singular: reciprocal condition number = 1.10707e-17

I am not sure what I am missing here for the inverse of the Hessian $H$ to throw that error. Need help. Please let me know if you need more details regarding my work.

$\endgroup$
8
  • 2
    $\begingroup$ Isn't it simply saying that in essence the Hessian is singular ? $\endgroup$
    – meh
    Jul 16, 2018 at 1:10
  • $\begingroup$ @aginensky, agree but what is the solution I am not sure. Does it mean I should not use the newton Raphson approach ? $\endgroup$ Jul 16, 2018 at 1:28
  • $\begingroup$ I'm not really an expert. I'm sure the answer is yes, but... You should google 'can I use Newton Raphson if the Hessian is singular". Unless you data set is really large or you are unlucky, you can try something else like gradient descent or stochastic gradient descent. $\endgroup$
    – meh
    Jul 16, 2018 at 1:31
  • $\begingroup$ @aginensky, I disagree. Gradient descent or Stochastic gradient descent might be relevant if my X are correlated. Here my X is not correlated. They are design matrices. $\endgroup$ Jul 16, 2018 at 1:41
  • 3
    $\begingroup$ Your X matrix is clearly linearly dependent. Add all of the columns after the first. It's all 1's, just like the first column. Hence X'X is singular. $\endgroup$
    – Glen_b
    Jul 16, 2018 at 1:55

1 Answer 1

1
$\begingroup$

Your $X$ matrix is clearly linearly dependent. Add all of the columns after the first - it's all $1$'s, just like the first column. Hence $X^\top X$ is singular.

The usual solution is to omit a single factor level, typically either the first or the last. There are other approaches (which boil done to some form of linear constraint on the coefficients); which you would use will depend on your needs (though it's possible to translate between them, so it's not crucial which approach you use).

My suggestion is simply to omit the second column of your X-matrix; the level of that is incorporated into the constant term and the other coefficients will then represent differences (on the scale of the linear predictor) from that baseline.

$\endgroup$
0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.