CODEWITHSUNDEEP
rglm
user227710
user227710

Reputation: 3194

How to return score vector after glm?

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.

Upvotes: 3

Views: 1381

Answers (1)

Achim Zeileis
Achim Zeileis

Reputation: 17193

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

Upvotes: 5

Related Questions