library(haven)
## Warning: package 'haven' was built under R version 3.5.1
glm_logit <- read_dta("glm-logit.dta")
View(glm_logit)
summary(glm_logit)
##      grade             gpa             tuce            psi        
##  Min.   :0.0000   Min.   :2.060   Min.   :12.00   Min.   :0.0000  
##  1st Qu.:0.0000   1st Qu.:2.812   1st Qu.:19.75   1st Qu.:0.0000  
##  Median :0.0000   Median :3.065   Median :22.50   Median :0.0000  
##  Mean   :0.3438   Mean   :3.117   Mean   :21.94   Mean   :0.4375  
##  3rd Qu.:1.0000   3rd Qu.:3.515   3rd Qu.:25.00   3rd Qu.:1.0000  
##  Max.   :1.0000   Max.   :4.000   Max.   :29.00   Max.   :1.0000
table(glm_logit$grade, glm_logit$psi)
##    
##      0  1
##   0 15  6
##   1  3  8
table(glm_logit$gpa,glm_logit$psi)
##                   
##                    0 1
##   2.05999994277954 0 1
##   2.39000010490417 0 1
##   2.63000011444092 1 0
##   2.66000008583069 1 0
##   2.67000007629395 0 1
##   2.74000000953674 1 0
##   2.75             1 0
##   2.75999999046326 1 0
##   2.82999992370605 1 1
##   2.85999989509583 1 0
##   2.86999988555908 1 0
##   2.89000010490417 1 1
##   2.92000007629395 1 0
##   3.02999997138977 1 0
##   3.09999990463257 0 1
##   3.11999988555908 0 1
##   3.16000008583069 0 1
##   3.25999999046326 1 0
##   3.27999997138977 1 0
##   3.3199999332428  1 0
##   3.39000010490417 0 1
##   3.50999999046326 0 1
##   3.52999997138977 1 0
##   3.53999996185303 0 1
##   3.5699999332428  1 0
##   3.61999988555908 0 1
##   3.65000009536743 0 1
##   3.92000007629395 1 0
##   4                1 1
boxplot(gpa ~ psi, data=glm_logit, horizontal=TRUE, yaxt="n",
        ylab="psi", xlab="gpa",
        main="Comparison of gpa of those who gave psi and Those who didn't")
axis(side=2, at=c(1,2), labels=c("gave", "notgave"))

logitMod <- glm(grade ~ psi, data=glm_logit, family=binomial(link="logit"))
summary((logitMod))
## 
## Call:
## glm(formula = grade ~ psi, family = binomial(link = "logit"), 
##     data = glm_logit)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.3018  -0.6039  -0.6039   1.0579   1.8930  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)  
## (Intercept)  -1.6094     0.6324  -2.545   0.0109 *
## psi           1.8971     0.8317   2.281   0.0225 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 41.183  on 31  degrees of freedom
## Residual deviance: 35.342  on 30  degrees of freedom
## AIC: 39.342
## 
## Number of Fisher Scoring iterations: 3

Just by using logistic regression the interpretation of the model is difficult than in oLS model. To make more meaningful interpetatin need to take exponential of the coefficients and then find the logit of odds ratio.

require(MASS)
## Loading required package: MASS
exp(cbind(coef(logitMod), confint(logitMod)))
## Waiting for profiling to be done...
##                           2.5 %     97.5 %
## (Intercept) 0.200000 0.04632242  0.6060404
## psi         6.666667 1.41559405 39.4941550

This result can be interpreted as those who are taking psi are 6.67 times more likely to get an A grade.

predict(logitMod)
##          1          2          3          4          5          6 
##  0.2876821  0.2876821 -1.6094379 -1.6094379 -1.6094379 -1.6094379 
##          7          8          9         10         11         12 
##  0.2876821 -1.6094379 -1.6094379 -1.6094379  0.2876821 -1.6094379 
##         13         14         15         16         17         18 
## -1.6094379 -1.6094379  0.2876821  0.2876821 -1.6094379  0.2876821 
##         19         20         21         22         23         24 
##  0.2876821  0.2876821 -1.6094379 -1.6094379 -1.6094379 -1.6094379 
##         25         26         27         28         29         30 
##  0.2876821  0.2876821  0.2876821 -1.6094379  0.2876821 -1.6094379 
##         31         32 
##  0.2876821 -1.6094379
library(margins)
## Warning: package 'margins' was built under R version 3.5.1
x <- glm(grade ~ psi, data = glm_logit)
(m <- margins(logitMod))
## Average marginal effects
## glm(formula = grade ~ psi, family = binomial(link = "logit"),     data = glm_logit)
##     psi
##  0.3515
summary(m)
##  factor    AME     SE      z      p  lower  upper
##     psi 0.3515 0.1005 3.4959 0.0005 0.1544 0.5485
plot(m)