如何在 glm 之后 return 对向量进行评分?
How to return score vector after glm?
我正在尝试从 R 中的概率中获取得分向量(对数似然的一阶导数)。这是我的示例代码:
library(AER) # Affairs data
data(Affairs)
mydata<-Affairs
mydata$affairs<-with(mydata,ifelse(affairs>0,1,affairs)) # convert to 1 and 0
glm(affairs ~ gender+ age + yearsmarried + children +
religiousness+education + rating,
family = binomial(link = "probit"),data = mydata)
Call: glm(formula = affairs ~ gender + age + yearsmarried + children +
religiousness + education + rating, family = binomial(link = "probit"),
data = mydata)
Coefficients:
(Intercept) gendermale age yearsmarried childrenyes religiousness education rating
0.76416 0.18882 -0.02440 0.05461 0.20807 -0.18609 0.01551 -0.27271
Degrees of Freedom: 600 Total (i.e. Null); 593 Residual
Null Deviance: 675.4
Residual Deviance: 610.5 AIC: 626.5
在 Stata 中,这将是 glm 之后的 predict anything,score
。
可以使用 sandwich
包中的 estfun()
方法提取观察到的数据的分数(也称为估计函数)。所以你先拟合模型
m <- glm(affairs ~ gender + age + yearsmarried + children + religiousness + education + rating,
family = binomial(link = "probit"), data = mydata)
然后您可以提取 n x k 观察分数矩阵:
library("sandwich")
s <- estfun(m)
dim(s)
## [1] 601 8
观测值之和基本上为零:
colSums(s)
## (Intercept) gendermale age yearsmarried childrenyes
## -0.0006048446 -0.0001081596 0.0005041310 0.0098581660 0.0005940238
## religiousness education rating
## -0.0006689870 -0.0147258776 -0.0016111599
检查前六个分数得出:
head(s)
## (Intercept) gendermale age yearsmarried childrenyes
## 4 -0.3789593 -0.3789593 -14.021496 -3.7895934 0.0000000
## 5 -0.1913665 0.0000000 -5.166896 -0.7654660 0.0000000
## 11 -0.7474727 0.0000000 -23.919126 -11.2120903 -0.7474727
## 16 -0.1564706 -0.1564706 -8.918825 -2.3470591 -0.1564706
## 23 -0.5249512 -0.5249512 -11.548926 -0.3937134 0.0000000
## 29 -0.1611754 0.0000000 -5.157613 -0.2417631 0.0000000
## religiousness education rating
## 4 -1.1368780 -6.821268 -1.515837
## 5 -0.7654660 -2.679131 -0.765466
## 11 -0.7474727 -8.969672 -2.989891
## 16 -0.7823530 -2.816471 -0.782353
## 23 -1.0499023 -8.924170 -1.574853
## 29 -0.3223508 -2.739982 -0.805877
我正在尝试从 R 中的概率中获取得分向量(对数似然的一阶导数)。这是我的示例代码:
library(AER) # Affairs data
data(Affairs)
mydata<-Affairs
mydata$affairs<-with(mydata,ifelse(affairs>0,1,affairs)) # convert to 1 and 0
glm(affairs ~ gender+ age + yearsmarried + children +
religiousness+education + rating,
family = binomial(link = "probit"),data = mydata)
Call: glm(formula = affairs ~ gender + age + yearsmarried + children +
religiousness + education + rating, family = binomial(link = "probit"),
data = mydata)
Coefficients:
(Intercept) gendermale age yearsmarried childrenyes religiousness education rating
0.76416 0.18882 -0.02440 0.05461 0.20807 -0.18609 0.01551 -0.27271
Degrees of Freedom: 600 Total (i.e. Null); 593 Residual
Null Deviance: 675.4
Residual Deviance: 610.5 AIC: 626.5
在 Stata 中,这将是 glm 之后的 predict anything,score
。
可以使用 sandwich
包中的 estfun()
方法提取观察到的数据的分数(也称为估计函数)。所以你先拟合模型
m <- glm(affairs ~ gender + age + yearsmarried + children + religiousness + education + rating,
family = binomial(link = "probit"), data = mydata)
然后您可以提取 n x k 观察分数矩阵:
library("sandwich")
s <- estfun(m)
dim(s)
## [1] 601 8
观测值之和基本上为零:
colSums(s)
## (Intercept) gendermale age yearsmarried childrenyes
## -0.0006048446 -0.0001081596 0.0005041310 0.0098581660 0.0005940238
## religiousness education rating
## -0.0006689870 -0.0147258776 -0.0016111599
检查前六个分数得出:
head(s)
## (Intercept) gendermale age yearsmarried childrenyes
## 4 -0.3789593 -0.3789593 -14.021496 -3.7895934 0.0000000
## 5 -0.1913665 0.0000000 -5.166896 -0.7654660 0.0000000
## 11 -0.7474727 0.0000000 -23.919126 -11.2120903 -0.7474727
## 16 -0.1564706 -0.1564706 -8.918825 -2.3470591 -0.1564706
## 23 -0.5249512 -0.5249512 -11.548926 -0.3937134 0.0000000
## 29 -0.1611754 0.0000000 -5.157613 -0.2417631 0.0000000
## religiousness education rating
## 4 -1.1368780 -6.821268 -1.515837
## 5 -0.7654660 -2.679131 -0.765466
## 11 -0.7474727 -8.969672 -2.989891
## 16 -0.7823530 -2.816471 -0.782353
## 23 -1.0499023 -8.924170 -1.574853
## 29 -0.3223508 -2.739982 -0.805877