Como calcular a matriz hat para regressão logística em R?

Eu quero calcular a matriz de chapéu diretamente em R para um modelo de logit. Segundo Long (1997), a matriz hat para modelos logit é definida como:

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

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

X é o vetor de variáveis independentes e V é uma matriz diagonal com na diagonal. $\sqrt{\pi(1-\pi)}$

Eu uso a optimfunção para maximizar a probabilidade e derivar o hessian. Então, acho que minha pergunta é: como faço para calcular em R? $V$

Nota: Minha função de probabilidade é assim:

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)
}

E eu alimento isso com a função otim da seguinte maneira:

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

Onde x é uma matriz de variáveis independentes e y é um vetor com a variável dependente.

Nota: Sei que existem procedimentos fixos para fazer isso, mas preciso fazê-lo do zero

r logistic deviance Thomas Jensen
fonte

De que maneira você está usando o otim (com quais opções, com ou sem fornecer uma função de gradiente, etc)? A regressão logística é um problema convexo suave. É prontamente resolvido usando o método de Newton ou similar. De fato, para obter uma estimativa da matriz de covariância, você precisa fazer (algo próximo disso) isso.

cardeal

Eu adicionei as informações ao post

Thomas Jensen

Respostas:

Para a regressão logística é calculado usando a fórmula $\pi$

π = \frac{1}{1 + \exp (- X β)}

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

Portanto, os valores diagonais de podem ser calculados da seguinte maneira: $V$

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

Agora multiplicar pela matriz diagonal da esquerda significa que cada linha é multiplicada pelo elemento correspondente da diagonal. O que em R pode ser alcançado usando multiplicação simples:

VX <- X*v

Então Hpode ser calculado da seguinte maneira:

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

Nota Como contém desvios padrão, suspeito que a fórmula correta para seja $V$ $H$

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

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

O código de exemplo funciona para esta fórmula.

mpiktas
fonte

Obrigado mpiktas, mas estou um pouco empolgado em como calcular V. V é simplesmente a diagonal da matriz de covariância?

Thomas Jensen

@Thomas, não, é a matriz diagonal que você tenha especificado-lo em seu post inicial, mas onde os são substituídas pelas estimativas , ou seja, a probabilidade estimada de que o th resposta é 1 sob o modelo.

π_{i}

$\pi_i$

{\hat{π}}_{i}

$\hat{\pi}_i$

i

$i$

cardeal

Ok, então para cada linha dos dados eu simplesmente calculo a probabilidade prevista e multiplico a raiz quadrada desse vetor pela matriz de variáveis independentes?

10138 Thomas Jensen

@ Thomas, sim, é assim que é feito no meu código. Você pode verificar com um exemplo fictício que realmente funciona.

mpiktas

@mpiktas - você está certo sobre . Efetivamente, o que você está fazendo é "padronizar" a matriz e o vetor , depois fazer os mínimos quadrados ponderados nas variáveis padronizadas e depois retroceder para a escala original. Você precisa iterate porque a padronização depende

V^{2}

$V^2$

X

$X$

Y

$Y$

β

$\beta$

probabilityislogic