I am trying to get the score vector (first order derivative of log likelihood) from probit in R. Here is my sample code:
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
In Stata, this would be predict anything,score after glm.
The scores (aka estimating functions) for the observed data can be extracted with the estfun() method from the sandwich package. So you first fit the model
m <- glm(affairs ~ gender + age + yearsmarried + children + religiousness + education + rating,
family = binomial(link = "probit"), data = mydata)
Then you can extract the n x k matrix of observation-wise scores:
library("sandwich")
s <- estfun(m)
dim(s)
## [1] 601 8
The sums across observations are essentially zero:
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
And inspecting the first six scores gives:
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