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

8

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=VX(XVX)1XV

X é o vetor de variáveis ​​independentes e V é uma matriz diagonal com na diagonal.π(1π)

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

Thomas Jensen
fonte
3
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:

13

Para a regressão logística é calculado usando a fórmulaπ

π=11+exp(Xβ)

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 sejaHVH

H=VX(XV2X)1XV

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 iπiπ^ii
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
1
@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 X Y βV2XYβ
probabilityislogic