Logit, Probit, and Multinomial Logit models in R

1 Logit , Probit and Multinomial Logit models in R (v. ) Oscar Torres-Reyna December 2014 Intro If outcome or dependent variable is binary and in the form 0/1, then use Logit or probit models . Some examples are: Did you vote in the last election? 0 No 1 Yes Do you prefer to use public transportation or to drive a car? 0 Prefer to drive 1 Prefer public transport If outcome or dependent variable is categorical but are ordered ( low to high), then use ordered Logit or ordered probit models . Some examples are: Do you agree or disagree with the President? 1 Disagree 2 Neutral 3 Agree What is your socioeconomic status? 1 Low 2 Middle 3 High If outcome or dependent variable is categorical without any particular order, then use Multinomial Logit . Some examples are: If elections were held today, for which party would you vote? 1 Democrats 2 Independent 3 Republicans What do you like to do on the weekends? 1 Rest 2 Go to movies 3 Exercise OTR 2 Logit model # Getting sample data library(foreign) mydata <- (" ") # Running a Logit model Logit <- glm(y_bin ~ x1 + x2 + x3, family=binomial(link=" Logit "), data=mydata) summary( Logit ) Call: glm(formula = y_bin ~ x1 + x2 + x3, family = binomial(link = " Logit "), data = mydata) Deviance Residuals: Min 1Q Median 3Q Max Coefficients: Estimate Std.

2 Error z value Pr(>|z|) (Intercept) x1 x2 x3 . --- Signif. codes: 0 ** ** * . 1 (Dispersion parameter for binomial family taken to be 1) Null deviance: on 69 degrees of freedom Residual deviance: on 66 degrees of freedom AIC: Number of Fisher Scoring iterations: 5 Outcome Predictors Type of model Data source Store results The Pr(>|z|) column shows the two-tailed p-values testing the null hypothesis that the coefficient is equal to zero ( no significant effect). The usual value is , by this measure none of the coefficients have a significant effect on the log-odds ratio of the dependent variable. The coefficient for x3 is significant at 10% (< ). The z value also tests the null that the coefficient is equal to zero. For a 5% significance, the z-value should fall outside the The Estimate column shows the coefficients in log-odds form.

3 When x3 increase by one unit, the expected change in the log odds is What you get from this column is whether the effect of the predictors is positive or negative. See next page for an extended explanation. OTR 3 Logit model # The stargazer() function from the package stargazer allows a publication quality of the Logit model . # The model will be saved in the working directory under the name which you can open with Word or any other word processor. library(stargazer) stargazer( Logit , type="html", out=" ") NOTE: Use the option type = "text" if you want to see the results directly in the RStudio console. OTR 4 Logit model : odds ratio # Estimating the odds ratio by hand cbind(Estimate=round(coef( Logit ),4), OR=round(exp(coef( Logit )),4)) Estimate OR (Intercept) x1 x2 x3 The Estimate column shows the coefficients in log-odds form. When x3 increase by one unit, the expected change in the log odds is Lets hold x1 and x2 constant at their means, and vary x3 with values 1, 2, and 3, to get the predicted log-odds given each of the three values of x3: r1 <- Logit $coeff[1] + Logit $coeff[2]*mean(mydata$x1) + Logit $coeff[3]*mean(mydata$x2) + Logit $coeff[4]*1 > r1 r2 <- Logit $coeff[1] + Logit $coeff[2]*mean(mydata$x1) + Logit $coeff[3]*mean(mydata$x2) + Logit $coeff[4]*2 > r2 r3 <- Logit $coeff[1] + Logit $coeff[2]*mean(mydata$x1) + Logit $coeff[3]*mean(mydata$x2) + Logit $coeff[4]*3 > r3 When x3 increases from 1 to 2, the log-odds increases: r2-r1 When x3 increases from 2 to 3, the log-odds increases: r3-r2 Which corresponds to the estimate for x3 above.

4 The odds ratio, is the exponentiation of the difference of the log-odds > exp(r2-r1) Or, the ratio of the exponentiation of each of the log-odds. > exp(r2)/exp(r1) Which corresponds to the OR value for x3 above. Odds ratio interpretation (OR): Based on the output below, when x3 increases by one unit, the odds of y = 1 increase by 112% -( )*100-. Or, the odds of y =1 are times higher when x3 increases by one unit (keeping all other predictors constant). To get the odds ratio, you need explonentiate the Logit coefficient. # Using package - mfx-- library(mfx) logitor(y_bin ~ x1 + x2 + x3, data=mydata) Call: logitor(formula = y_bin ~ x1 + x2 + x3, data = mydata) Odds Ratio: OddsRatio Std. Err. z P>|z| x1 x2 x3 . --- Signif. codes: 0 ** ** * . 1 OTR 5 Logit model : odds ratios # Relative risk ratios allow an easier interpretation of the Logit coefficients. They are the exponentiated value of the Logit coefficients.

5 = exp(coef( Logit )) (Intercept) x1 x2 x3 library(stargazer) stargazer( Logit , type="html", coef=list( ), , out=" ") NOTE: Use the option type = "text" if you want to see the results directly in the RStudio console. OTR 6 Keeping all other variables constant, when x1 increases one unit, it is times more likely to be in the 1 category. In other words, the odds of being in the 1 category (as opposed to the 0 category) are 136% higher when x1 move one unit ( 1). The coefficient, however, is not significant. Logit model : predicted probabilities The Logit model can be written as (Gelman and Hill, 2007): Pr(yi = 1) = Logit -1(Xi ) In the example: Logit <- glm(y_bin ~ x1 + x2 + x3, family=binomial(link=" Logit "), data=mydata) coef( Logit ) (Intercept) x1 x2 x3 Pr(yi = 1) = Logit -1( + *x1 + *x2 + *x3 ) Estimating the probability at the mean point of each predictor can be done by inverting the Logit model .

6 Gelman and Hill provide a function for this (p. 81), also available in the R package arm- invlogit = function (x) {1/(1+exp(-x))} invlogit(coef( Logit )[1]+ coef( Logit )[2]*mean(mydata$x1)+ coef( Logit )[3]*mean(mydata$x2)+ coef( Logit )[4]*mean(mydata$x3)) Pr(yi = 1) = OTR 7 Logit model : predicted probabilities Adding categorical variable, the model would be: <- glm(y_bin ~ x1 + x2 + x3 + opinion, family=binomial(link=" Logit "), data=mydata) coef( ) (Intercept) x1 x2 x3 opinionAgree opinionDisag opinionStr disag Estimating the probability when opinion = Agree invlogit = function (x) {1/(1+exp(-x))} invlogit(coef( )[1]+ coef( )[2]*mean(mydata$x1)+ coef( )[3]*mean(mydata$x2)+ coef( )[4]*mean(mydata$x3)+ coef( )[5]*1) Pr(yi = 1| opinion= Agree ) = OTR 8 Logit model .

7 Predicted probabilities invlogit = function (x) {1/(1+exp(-x))} Estimating the probability when opinion = Disagree invlogit(coef( )[1]+ coef( )[2]*mean(mydata$x1)+ coef( )[3]*mean(mydata$x2)+ coef( )[4]*mean(mydata$x3)+ coef( )[6]*1) Pr(yi = 1| opinion= Disagree ) = Estimating the probability when opinion = Strongly disagree invlogit(coef( )[1]+ coef( )[2]*mean(mydata$x1)+ coef( )[3]*mean(mydata$x2)+ coef( )[4]*mean(mydata$x3)+ coef( )[7]*1) Pr(yi = 1| opinion= Strongly disagree ) = Estimating the probability when opinion = Strongly agree invlogit(coef( )[1]+ coef( )[2]*mean(mydata$x1)+ coef( )[3]*mean(mydata$x2)+ coef( )[4]*mean(mydata$x3)) Pr(yi = 1| opinion= Strongly agree ) = OTR 9 Logit model : predicted probabilities Another way to estimate the predicted probabilities is by setting initial conditions.

8 Getting predicted probabilities holding all predictors or independent variables to their means. allmean <- (x1=mean(mydata$x1), x2=mean(mydata$x2), x3=mean(mydata$x3)) allmean x1 x2 x3 1 After estimating the Logit model and creating the dataset with the mean values of the predictors, you can use the predict() function to estimate the predicted probabilities (for help/details type ? ), and add them to the allmean dataset. allmean$ <- predict( Logit , newdata=allmean, type="response") allmean x1 x2 x3 1 When all predictor values are hold to their means, the probability of y = 1 is 83%. Creating a new dataset with the mean values of the predictors The object with the Logit coefficients Dataset with the conditions Requesting predicted probabilities OTR 10 Logit model : predicted probabilities with categorical variable Logit <- glm(y_bin ~ x1+x2+x3+opinion, family=binomial(link=" Logit "), data=mydata) To estimate the predicted probabilities, we need to set the initial conditions.

9 Getting predicted probabilities holding all predictors or independent variables to their means for each category of categorical variable opinion : allmean <- (x1=rep(mean(mydata$x1),4), x2=rep(mean(mydata$x2),4), x3=rep(mean(mydata$x3),4), opinion= (c("Str agree","Agree","Disag","Str disag"))) allmean x1 x2 x3 opinion 1 Str agree 2 Agree 3 Disag 4 Str disag allmean <- cbind(allmean,predict( Logit , newdata=allmean, type="response", )) allmean (continue next page) Creating a new dataset with the mean values of the predictors for each category The object with the Logit coefficients Dataset with the conditions Requesting predicted probabilities OTR 11 Standard error of the prediction x1 x2 x3 opinion fit 1 Str agree 1 2 Agree 1 3 Disag 1 4 Str disag 1 Logit model .

10 Predicted probabilities with categorical variable # Renaming "fit" and " " columns names(allmean)[names(allmean)=="fit"] = "prob" names(allmean)[names(allmean)==" "] = " " # Estimating confidence intervals allmean$ll = allmean$prob - *allmean$ allmean$ul = allmean$prob + *allmean$ allmean x1 x2 x3 opinion prob ll ul 1 Str agree 1 2 Agree 1 3 Disag 1 4 Str disag 1 (continue next page) OTR 12 agreeAgreeDisagStr disagOpinionPr(y=1) Predicted probabilitiesadd footnote hereLogit model : predicted probabilities with categorical variable # Plotting predicted probabilities and confidence intervals using ggplot2 library(ggplot2) ggplot(allmean, aes(x=opinion, y = prob)) + geom_errorbar(aes(ymin = ll, ymax = ul), width = , lty=1, lwd=1, col="red") + geom_point(shape=18, size=5, fill="black") + scale_x_discrete(limits = c("Str agree","Agree","Disag","Str disag")) + labs(title= " Predicted probabilities", x="Opinion", y="Pr(y=1)", caption = "add footnote here") + theme( = element_text(family = "sans", face="bold", size=13, hjust= ), = element_text(family = "sans", size=9), = element_text(family = "sans", size=5)) OTR 13 Logit model : marginal effects # Using package mfx- # See ("mfx") #Do this only once library(mfx) logitmfx(y_bin ~ x1+x2+x3, data=mydata) Call: logitmfx(formula = y_bin ~ x1 + x2 + x3, data = mydata) Marginal Effects: dF/dx Std.