Setting up everything
data(titanic_train, package = "titanic")
dat <- titanic_train
rm(titanic_train)table(dat$Pclass)
1 2 3
216 184 491 reg <- glm(Survived ~ Pclass, data = dat, family = "binomial")
summary(reg)
Call:
glm(formula = Survived ~ Pclass, family = "binomial", data = dat)
Deviance Residuals:
Min 1Q Median 3Q Max
-1.4390 -0.7569 -0.7569 0.9367 1.6673
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 1.44679 0.20743 6.975 3.06e-12 ***
Pclass -0.85011 0.08715 -9.755 < 2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 1186.7 on 890 degrees of freedom
Residual deviance: 1084.4 on 889 degrees of freedom
AIC: 1088.4
Number of Fisher Scoring iterations: 4The estimated coefficients are in logits
paste("the estimated cofficients in LOGITs are", coef(reg)[1], "for the intercept. And", coef(reg)[2], "for the slope for Pclass")[1] "the estimated cofficients in LOGITs are 1.44678953058191 for the intercept. And -0.850106746726404 for the slope for Pclass" if the logit is +ve, the qualitative effect of this predicator on the binary outcome is also +ve
if the logit is -ve, the qualitative effect of this predicator on the binary outcome is also -ve
Suppose logit = 1, then the odds ratio = exp(1) approx.= 2.71 then the probability = exp(1) / (1 + exp(1)) approx.= 2.71 / (1 + 2.71) approx.= 0.75
Suppose logit = 0, then the odds ratio = exp(0) = 1 then the probability = exp(0) / (1 + exp(0)) = 1 / (1 + 1) = 0.5
if logit < 0, the associated probability < 0.5
if logit > 0, the associated probability > 0.5
From the logistic regression output, we have the following estimated equation
logit_survival = 1.44679 + (-0.85011) * Pclass; where Pclass = 1,2,3
logits_survival_1 <- coef(reg)[1] + coef(reg)[2] * 1
logits_survival_2 <- coef(reg)[1] + coef(reg)[2] * 2
logits_survival_3 <- coef(reg)[1] + coef(reg)[2] * 3
paste("the logits of survival for Pclass 1 is", logits_survival_1)[1] "the logits of survival for Pclass 1 is 0.596682783855506"paste("the logits of survival for Pclass 2 is", logits_survival_2)[1] "the logits of survival for Pclass 2 is -0.253423962870898"paste("the logits of survival for Pclass 3 is", logits_survival_3)[1] "the logits of survival for Pclass 3 is -1.1035307095973"Convert to logits of survival to the probabilities of survival
prob_survival_1 <- exp(logits_survival_1) / (1 + exp(logits_survival_1))
prob_survival_2 <- exp(logits_survival_2) / (1 + exp(logits_survival_2))
prob_survival_3 <- exp(logits_survival_3) / (1 + exp(logits_survival_3))
paste("the probabilities of survival for Pclass 1 is", prob_survival_1)[1] "the probabilities of survival for Pclass 1 is 0.644897013275819"paste("the probabilities of survival for Pclass 2 is", prob_survival_2)[1] "the probabilities of survival for Pclass 2 is 0.43698092535236"paste("the probabilities of survival for Pclass 3 is", prob_survival_3)[1] "the probabilities of survival for Pclass 3 is 0.24907893048446" prob = 0.5 * (1 + tanh(logit / 2))
prob_survival_1b <- 0.5 * (1 + tanh(logits_survival_1 / 2))
prob_survival_2b <- 0.5 * (1 + tanh(logits_survival_2 / 2))
prob_survival_3b <- 0.5 * (1 + tanh(logits_survival_3 / 2))
paste("the probabilities of survival for Pclass 1 is", prob_survival_1b)[1] "the probabilities of survival for Pclass 1 is 0.644897013275819"paste("the probabilities of survival for Pclass 2 is", prob_survival_2b)[1] "the probabilities of survival for Pclass 2 is 0.43698092535236"paste("the probabilities of survival for Pclass 3 is", prob_survival_3b)[1] "the probabilities of survival for Pclass 3 is 0.24907893048446"paste("the probabilities of survival for Pclass 1 is",
predict(reg, data.frame(Pclass = 1), type = "response")
)[1] "the probabilities of survival for Pclass 1 is 0.644897013275819"paste("the probabilities of survival for Pclass 2 is",
predict(reg, data.frame(Pclass = 2), type = "response")
)[1] "the probabilities of survival for Pclass 2 is 0.43698092535236"paste("the probabilities of survival for Pclass 3 is",
predict(reg, data.frame(Pclass = 3), type = "response")
)[1] "the probabilities of survival for Pclass 3 is 0.24907893048446"For the marginal effects, let’s examine the effect on survival from i) changing from Pclass 1 to Pclass 2 ii) changing from Pclass 2 to Pclass 3 in terms of logits
paste("the effect on survivial in LOGIT from changing from Pclass 1 to Pclass 2 is", logits_survival_2 - logits_survival_1)[1] "the effect on survivial in LOGIT from changing from Pclass 1 to Pclass 2 is -0.850106746726404"paste("the effect on survivial in LOGIT from changing from Pclass 2 to Pclass 3 is", logits_survival_3 - logits_survival_2)[1] "the effect on survivial in LOGIT from changing from Pclass 2 to Pclass 3 is -0.850106746726404"paste("note that the difference, independent of actual Pclass value, is exactly the estimated SLOPE COEFFICIENT", coef(reg)[2])[1] "note that the difference, independent of actual Pclass value, is exactly the estimated SLOPE COEFFICIENT -0.850106746726404"Now we try to compute the marginal effects in terms of probabilities by using a) compute probability ME from the difference in logits (i.e. the slope coefficient) b) compute the incremental difference in probabilities at each Pclass
prob_ME_a <- exp(coef(reg)[2]) / (1 + exp(coef(reg)[2]))
paste("the result in probability ME from the slope cofficient is", prob_ME_a)[1] "the result in probability ME from the slope cofficient is 0.299410465443492"paste("the effect on survivial in PROB from changing from Pclass 1 to Pclass 2 is", prob_survival_2 - prob_survival_1)[1] "the effect on survivial in PROB from changing from Pclass 1 to Pclass 2 is -0.207916087923458"paste("the effect on survivial in PROB from changing from Pclass 2 to Pclass 3 is", prob_survival_3 - prob_survival_2)[1] "the effect on survivial in PROB from changing from Pclass 2 to Pclass 3 is -0.1879019948679" logit(prob) = log(prob / (1 - prob)) = b0 + b1 * x
prob = exp(b0 + b1 * x) / (1 + exp(b0 + b1 * x))
=> prob(logit) = exp(logit(prob)) / (1 + exp(logit(prob)))
For discrete marginal effects (no difference in the continuous case),
A 1-unit change in x, changes the LOGIT by
logit(.|x = 2) - logit(.|x = 1)
= b0 + b1 * 2 - (b0 + b1 * 1)
= b1
A 1-unit change in x, changes the PROB by
prob(.| x = 2) - prob(. |x = 1)
= [exp(b0 + b1 * 2) / (1 + exp(b0 + b1 * 2))] - [exp(b0 + b1 * 1) / (1 + exp(b0 + b1 * 1))]
!= b1
furthermore, != a constant
That is, delta(logit) = constant
but, delta(prob) = f(x)Therefore, in step 3, approach a is WRONG and approach b is CORRECT
paste("the effect on survivial in PROB from changing from Pclass 1 to Pclass 2 is",
(exp(logits_survival_2) / (1 + exp(logits_survival_2))) - (exp(logits_survival_1) / (1 + exp(logits_survival_1))))[1] "the effect on survivial in PROB from changing from Pclass 1 to Pclass 2 is -0.207916087923458"paste("the effect on survivial in PROB from changing from Pclass 2 to Pclass 3 is",
(exp(logits_survival_3) / (1 + exp(logits_survival_3))) - (exp(logits_survival_2) / (1 + exp(logits_survival_2))))[1] "the effect on survivial in PROB from changing from Pclass 2 to Pclass 3 is -0.1879019948679"