9
$\begingroup$

I want to calculate the hat matrix directly in R for a logit model. According to Long (1997) the hat matrix for logit models is defined as:

$$H = VX(X'VX)^{-1} X'V$$

X is the vector of independent variables, and V is a diagonal matrix with $\sqrt{\pi(1-\pi)}$ on the diagonal.

I use the optim function to maximize the likelihood and derive the hessian. So I guess my question is: how do i calculuate $V$ in R?

Note: My likelihood function looks like this:

loglik <-  function(theta,x,y){
y <- y
x <- as.matrix(x)
beta <- theta[1:ncol(x)]
loglik <- sum(-y*log(1 + exp(-(x%*%beta))) - (1-y)*log(1 + exp(x%*%beta)))
return(-loglik)
}

And i feed this to the optim function as follows:

logit <- optim(c(1,1),loglik, y = y, x = x, hessian = T)

Where x is a matrix of independent variables, and y is a vector with the dependent variable.

Note: I know that there are canned procedures for doing this, but I need to do it from scratch

$\endgroup$
2
  • 3
    $\begingroup$ In what way are you using optim (with what options, with or without supplying a gradient function, etc)?? Logistic regression is a smooth convex problem. It's readily solved using Newton's method or similar. In fact, to get an estimate of the covariance matrix, you need to do (something close to) this. $\endgroup$
    – cardinal
    Mar 10 '11 at 14:41
  • $\begingroup$ I have added the info to the post $\endgroup$ Mar 10 '11 at 14:51
13
$\begingroup$

For logistic regression $\pi$ is calculated using formula

$$\pi=\frac{1}{1+\exp(-X\beta)}$$

So diagonal values of $V$ can be calculated in the following manner:

pi <- 1/(1+exp(-X%*%beta))
v <- sqrt(pi*(1-pi))

Now multiplying by diagonal matrix from left means that each row is multiplied by corresponding element from diagonal. Which in R can be achieved using simple multiplication:

VX <- X*v 

Then H can be calculated in the following way:

H <- VX%*%solve(crossprod(VX,VX),t(VX))

Note Since $V$ contains standard deviations I suspect that the correct formula for $H$ is

$$H=VX(X'V^2X)^{-1}X'V$$

The example code works for this formula.

$\endgroup$
5
  • $\begingroup$ Thanks mpiktas, but I am somewhat stuck at how to calculate V. Is V simply the diagonal of the covariance matrix? $\endgroup$ Mar 10 '11 at 14:38
  • $\begingroup$ @Thomas, no, it's the diagonal matrix as you've specified it in your initial post, but where the $\pi_i$ are replaced by the estimates $\hat{\pi}_i$, i.e., the estimated probability that the $i$th response is 1 under the model. $\endgroup$
    – cardinal
    Mar 10 '11 at 14:57
  • $\begingroup$ Ok, so for each row in the data I simply calculate the predicted probability, and multiply the square root of this vector with the matrix of independent variables? $\endgroup$ Mar 10 '11 at 15:03
  • $\begingroup$ @Thomas, yes, that is how it is done in my code. You can check with dummy example that it really works. $\endgroup$
    – mpiktas
    Mar 10 '11 at 15:05
  • 2
    $\begingroup$ @mpiktas - you are right about $V^2$. Effectively what you are doing is "standardising" the $X$ matrix, and $Y$ vector, then doing weighted least squares on the standardised variables, then backtransforming to original scale. You need to iterate because the standardisation depends on $\beta$ $\endgroup$ Mar 10 '11 at 23:14

Your Answer

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

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