source("https://teaching.sociology.ul.ie/so5032/tabxw.R")
library(foreign)
hten <- read.dta("https://teaching.sociology.ul.ie/so5032/hten.dta")
hten$univ <- hten$qfedhi %in% c("higher degree", "first degree")
hten$ownocc <- hten$tenure == "owned outright"

summary(modage <- glm(data=hten, ownocc ~ age, family="binomial"))

tabxw(hten, univ, ownocc)

tab <- tabxw(hten, univ, ownocc)
sprintf("OR = (%d/%d)/(%d/%d) = %7.4f",
        tab[2,2],tab[2,1],tab[1,2],tab[1,1],
        (tab[2,2]/tab[2,1])/(tab[1,2]/tab[1,1]))

(mod1 <- glm(data=hten, ownocc ~ univ, family=binomial))
sprintf("OR = exp(%7.4f) = %7.4f",
        mod1$coefficients["univTRUE"],
        exp(mod1$coefficients["univTRUE"]))

summary(mod2 <- glm(data=hten, ownocc ~ univ + age, family=binomial))

library(lmtest)
modage <- glm(formula = ownocc ~ age, family = binomial, data = hten)
modage2 <- glm(formula = ownocc ~ age + educ, family = binomial, data = hten)
lrtest(modage2, modage)

## Install if necessary
## install.packages("margins")
library(margins)
margins(modage)

p = predict(modage, type="response")
mean(p*(1-p)*modage$coefficients["age"])

margins(modage)

modage2 <- glm(formula = ownocc ~ age + educ, family = binomial, data = hten)
p = predict(modage2, type="response")
mean(p*(1-p)*modage2$coefficients["age"])

margins(modage2)

dpdfagg <- data.frame(
    yes = c( 53, 0, 11, 4),
    no = c( 414, 16, 37, 139),
    defendant = factor(c( 1, 1, 2, 2), labels=c("White", "Black")),
    victim    = factor(c( 1, 2, 1, 2), labels=c("White", "Black")))
(mod3 <- glm(data=dpdfagg, cbind(yes, no) ~ defendant + victim, family=binomial))

(modsat <- glm(data=dpdfagg, cbind(yes, no) ~ defendant * victim, family=binomial))

dpdfagg$pp <- predict(modsat, type="response")
dpdfagg$E <- round((dpdfagg$yes +dpdfagg$no)*dpdfagg$pp)
dpdfagg

modsat <- glm(data=dpdfagg, cbind(yes, no) ~ defendant * victim, family=binomial)
mod1 <- glm(data=dpdfagg, cbind(yes, no) ~ defendant + victim, family=binomial)
lrtest(mod1, modsat)

mod1

cclogit <- glm(data=cc, card ~ income, family=binomial(link="logit"))
summary(cclogit)

ccprobit <- glm(data=cc, card ~ income, family=binomial(link="probit"))
summary(ccprobit)

cclogit$coefficients

exp(-3.5179469  + 0.1054089*50) / (1 + exp(-3.5179469  + 0.1054089*50))

ccprobit$coefficients

pnorm(-2.08933427 + 0.06222824*50)