机器算法验证 - logit 的方差-协方差矩阵与矩阵计算 - 吾爱随笔录

logit 的方差-协方差矩阵与矩阵计算

机器算法验证 r 回归协方差矩阵

2022-03-27 10:56:01

我正在尝试获取逻辑回归的方差-协方差矩阵：

mydata <- read.csv("https://stats.idre.ucla.edu/stat/data/binary.csv")
mylogit <- glm(admit ~ gre + gpa, data = mydata, family = "binomial")

通过矩阵计算。我一直在关注这里发布的基本线性回归示例

X <- as.matrix(cbind(1, mydata[,c('gre','gpa')]))
beta.hat <- as.matrix(coef(mylogit))
Y <- as.matrix(mydata$admit)
y.hat <- X %*% beta.hat

n <- nrow(X)
p <- ncol(X)

sigma2 <- sum((Y - y.hat)^2)/(n - p)        
v <- solve(t(X) %*% X) * sigma2

但是我的 var/cov 矩阵不等于计算的矩阵vcov()

v == vcov(mylogit)

1   gre   gpa
1   FALSE FALSE FALSE
gre FALSE FALSE FALSE
gpa FALSE FALSE FALSE

我错过了一些日志转换吗？

3个回答

@Deep North：你是对的，不应该有'n'

逻辑回归的协方差矩阵不同于线性回归的协方差矩阵。

线性回归： $\sigma^2(X'X)^{-1}$
逻辑回归： $(X'WX)^{-1}$

其中 W 是对角矩阵 $w_{ii}=\p_i(1-\p_i)$

$\p_i$ 是观察级别上事件 = 1 的概率

来自 subra 的逻辑回归的协方差是正确的。但是。不应该有。参考。David W. Hosmer 应用逻辑回归（第 2 版）p35 和 p41 公式（2.8） $w_{ii}=\hat{\pi_i}(1-\hat{\pi_i})$ $n_i$

我修改了您的程序并与方差估计进行比较，它们很接近但不一样。

 library(Matrix)
 library(sandwich)
 mydata <- read.csv("http://www.ats.ucla.edu/stat/data/binary.csv")
 mylogit <- glm(admit ~ gre + gpa, data = mydata, family = "binomial")
 X <- as.matrix(cbind(1, mydata[,c('gre','gpa')]))
 n <- nrow(X)

 pi<-mylogit$fit

 w<-pi*(1-pi)

 v<-Diagonal(n, x = w)
 var_b<-solve(t(X)%*%v%*%X)
 var_b

         x 3 Matrix of class "dgeMatrix"
          [,1]          [,2]          [,3]
[1,]  1.1558251135 -2.818944e-04 -0.2825632388
[2,] -0.0002818944  1.118288e-06 -0.0001144821
[3,] -0.2825632388 -1.144821e-04  0.1021349767


vcov(mylogit)

         (Intercept)           gre           gpa
(Intercept)  1.1558247051 -2.818942e-04 -0.2825631552
gre         -0.0002818942  1.118287e-06 -0.0001144821
gpa         -0.2825631552 -1.144821e-04  0.1021349526

前五位数字相同

mydata <- read.csv("http://www.ats.ucla.edu/stat/data/binary.csv")
mylogit <- glm(admit ~ gre + gpa, data = mydata, family = "binomial")

X <- as.matrix(cbind(1, mydata[,c('gre','gpa')]))
beta.hat <- as.matrix(coef(mylogit))

require(slam)
p <- 1/(1+exp(-X %*% beta.hat))
V <- simple_triplet_zero_matrix(dim(X)[1])
diag(V) <- p*(1-p)
IB <- matprod_simple_triplet_matrix(t(X), V) %*% X
varcov_mat <- solve(IB)

round(solve(IB),4) == round(vcov(mylogit),4)

# 1  gre  gpa
# 1   TRUE TRUE TRUE
# gre TRUE TRUE TRUE
# gpa TRUE TRUE TRUE

其它你可能感兴趣的问题

上一篇检查硬币是否公平下一篇条件模型与边际模型