# Multinomial Probit and Logit Models in R

# Copyright 2013 by Ani Katchova

# install.packages("mlogit")
library(mlogit)

mydata<- read.csv("C:/Econometrics/Data/multinomial_fishing1.csv")
attach(mydata)

# Descriptive statistics
table(mode)

# Reshaping the data from wide to long format

#mydata$mode<-as.factor(mydata$mode)
mldata<-mlogit.data(mydata, varying=4:15, choice="mode", shape="wide")
mldata[1:20,]

# Multinomial logit model coefficients

mlogit.model1 <- mlogit(mode ~ 1 | income, data=mldata, reflevel="charter")
summary(mlogit.model1)

# Multinomial logit model coefficients (with different base outcome)

mlogit.model2 <- mlogit(mode ~ 1 | income, data = mldata, reflevel="pier")
summary(mlogit.model2)

# Multinomial logit model odds ratios

exp(coef(mlogit.model1))

# Conditional logit model

clogit.model1 <- mlogit(mode ~ price+catch |income, data = mldata, reflevel="charter")
summary(clogit.model1)

clogit.model2 <- mlogit(mode ~ price+catch | income, data = mldata, reflevel="pier")

summary(clogit.model2)

# Setting mean values for variables to use for marginal effects

m <- mlogit(mode ~ price+catch |income, data = mldata, reflevel="charter")
z <- with(mldata, data.frame(price = tapply(price, index(m)$alt, mean),
catch = tapply(catch, index(m)$alt, mean), income =
mean(income)))

# Multinomial logit model marginal effects

effects(mlogit.model1, covariate = "income", data = z)

# Conditional logit model marginal effects

effects(clogit.model1, covariate = "income", data = z)
effects(clogit.model1, covariate = "price", data = z)
effects(clogit.model1, covariate = "catch", data = z)

# Multinomial probit model coefficients

#mprobit.model1 <- mlogit(mode ~ 1 | income, data = mldata, reflevel="charter",
probit=TRUE)
#summary(mprobit.model1)
# Hauseman-McFadden test of independence of irrelevant alternatives
m1<- mlogit(mode ~ 1 | income, data = mldata, reflevel="beach")
m2<- mlogit(mode ~ 1 | income, data = mldata, reflevel="beach", alt.subset=c("beach",
"pier", "private"))
hmftest(m1, m2)
> # Multinomial Probit and Logit Models in R
> # Copyright 2013 by Ani Katchova
>
> # install.packages("mlogit")
> library(mlogit)
Loading required package: Formula
Loading required package: statmod
Loading required package: lmtest
Loading required package: zoo

Attaching package: ‘zoo’

The following object(s) are masked from ‘package:base’:

as.Date, as.Date.numeric

Loading required package: maxLik

Loading required package: miscTools
Loading required package: MASS
>
> mydata<- read.csv("C:/Econometrics/Data/multinomial_fishing1.csv")
> attach(mydata)
>
> # Descriptive statistics
> table(mode)
mode
beach charter pier private
134 452 178 418
>
> # Reshaping the data from wide to long format
> #mydata$mode<-as.factor(mydata$mode)
> mldata<-mlogit.data(mydata, varying=4:15, choice="mode", shape="wide")
> mldata[1:20,]
mode price catchrate income mode1 alt d catch chid
1.beach FALSE 157.930 0.5391 7.083332 4 beach 0 0.0678 1
1.charter TRUE 182.930 0.5391 7.083332 4 charter 1 0.5391 1
1.pier FALSE 157.930 0.5391 7.083332 4 pier 0 0.0503 1
1.private FALSE 157.930 0.5391 7.083332 4 private 0 0.2601 1
2.beach FALSE 15.114 0.4671 1.250000 4 beach 0 0.1049 2
2.charter TRUE 34.534 0.4671 1.250000 4 charter 1 0.4671 2
2.pier FALSE 15.114 0.4671 1.250000 4 pier 0 0.0451 2
2.private FALSE 10.534 0.4671 1.250000 4 private 0 0.1574 2
3.beach FALSE 161.874 0.2413 3.750000 3 beach 0 0.5333 3
3.charter FALSE 59.334 0.2413 3.750000 3 charter 0 1.0266 3
3.pier FALSE 161.874 0.2413 3.750000 3 pier 0 0.4522 3
3.private TRUE 24.334 0.2413 3.750000 3 private 1 0.2413 3
4.beach FALSE 15.134 0.0789 2.083333 2 beach 0 0.0678 4
4.charter FALSE 84.930 0.0789 2.083333 2 charter 0 0.5391 4
4.pier TRUE 15.134 0.0789 2.083333 2 pier 1 0.0789 4
4.private FALSE 55.930 0.0789 2.083333 2 private 0 0.1643 4
5.beach FALSE 106.930 0.1082 4.583332 3 beach 0 0.0678 5
5.charter FALSE 71.014 0.1082 4.583332 3 charter 0 0.3240 5
5.pier FALSE 106.930 0.1082 4.583332 3 pier 0 0.0503 5
5.private TRUE 41.514 0.1082 4.583332 3 private 1 0.1082 5
>
> # Multinomial logit model coefficients
> mlogit.model1 <- mlogit(mode ~ 1 | income, data=mldata, reflevel="charter")
> summary(mlogit.model1)

Call:
mlogit(formula = mode ~ 1 | income, data = mldata, reflevel = "charter",
method = "nr", print.level = 0)

Frequencies of alternatives:
charter beach pier private
0.38240 0.11337 0.15059 0.35364

nr method
4 iterations, 0h:0m:0s
g'(-H)^-1g = 8.32E-07
gradient close to zero

Coefficients :
Estimate Std. Error t-value Pr(>|t|)
beach:(intercept) -1.341291 0.194517 -6.8955 5.367e-12 ***
pier:(intercept) -0.527141 0.177784 -2.9651 0.003026 **
private:(intercept) -0.602371 0.136096 -4.4261 9.597e-06 ***
beach:income 0.031640 0.041846 0.7561 0.449591
pier:income -0.111763 0.043979 -2.5413 0.011046 *
private:income 0.123546 0.027911 4.4265 9.577e-06 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Log-Likelihood: -1477.2
McFadden R^2: 0.013736
Likelihood ratio test : chisq = 41.145 (p.value = 6.0931e-09)
>
> # Multinomial logit model coefficients (with different base outcome)
> mlogit.model2 <- mlogit(mode ~ 1 | income, data = mldata, reflevel="pier")
> summary(mlogit.model2)

Call:
mlogit(formula = mode ~ 1 | income, data = mldata, reflevel = "pier",
method = "nr", print.level = 0)

Frequencies of alternatives:
pier beach charter private
0.15059 0.11337 0.38240 0.35364

nr method
4 iterations, 0h:0m:0s
g'(-H)^-1g = 8.32E-07
gradient close to zero

Coefficients :
Estimate Std. Error t-value Pr(>|t|)
beach:(intercept) -0.814150 0.228632 -3.5610 0.0003695 ***
charter:(intercept) 0.527141 0.177784 2.9651 0.0030262 **
private:(intercept) -0.075229 0.183240 -0.4106 0.6814007
beach:income 0.143403 0.053288 2.6911 0.0071223 **
charter:income 0.111763 0.043979 2.5413 0.0110455 *
private:income 0.235309 0.043668 5.3886 7.101e-08 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Log-Likelihood: -1477.2
McFadden R^2: 0.013736
Likelihood ratio test : chisq = 41.145 (p.value = 6.0931e-09)
>
> # Multinomial logit model odds ratios
> exp(coef(mlogit.model1))
beach:(intercept) pier:(intercept) private:(intercept)
0.2615077 0.5902901 0.5475121
beach:income pier:income private:income
1.0321457 0.8942561 1.1315023
attr(,"fixed")
beach:(intercept) pier:(intercept) private:(intercept)
FALSE FALSE FALSE
beach:income pier:income private:income
FALSE FALSE FALSE
>
>
> # Conditional logit model
> clogit.model1 <- mlogit(mode ~ price+catch |income, data = mldata,
reflevel="charter")
> summary(clogit.model1)

Call:
mlogit(formula = mode ~ price + catch | income, data = mldata,
reflevel = "charter", method = "nr", print.level = 0)

Frequencies of alternatives:
charter beach pier private
0.38240 0.11337 0.15059 0.35364

nr method
7 iterations, 0h:0m:0s
g'(-H)^-1g = 1.37E-05
successive function values within tolerance limits

Coefficients :
Estimate Std. Error t-value Pr(>|t|)
beach:(intercept) -1.6943657 0.2240506 -7.5624 3.952e-14 ***
pier:(intercept) -0.9164063 0.2072648 -4.4214 9.805e-06 ***
private:(intercept) -1.1670869 0.1590475 -7.3380 2.169e-13 ***
price -0.0251166 0.0017317 -14.5042 < 2.2e-16 ***
catch 0.3577820 0.1097733 3.2593 0.001117 **
beach:income 0.0332917 0.0503409 0.6613 0.508403
pier:income -0.0942854 0.0500600 -1.8834 0.059640 .
private:income 0.1227315 0.0286306 4.2867 1.813e-05 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Log-Likelihood: -1215.1
McFadden R^2: 0.18868
Likelihood ratio test : chisq = 565.17 (p.value = < 2.22e-16)
>
> clogit.model2 <- mlogit(mode ~ price+catch | income, data = mldata, reflevel="pier")
> summary(clogit.model2)

Call:
mlogit(formula = mode ~ price + catch | income, data = mldata,
reflevel = "pier", method = "nr", print.level = 0)

Frequencies of alternatives:
pier beach charter private
0.15059 0.11337 0.38240 0.35364

nr method
7 iterations, 0h:0m:0s
g'(-H)^-1g = 1.37E-05
successive function values within tolerance limits

Coefficients :
Estimate Std. Error t-value Pr(>|t|)
beach:(intercept) -0.7779594 0.2204939 -3.5283 0.0004183 ***
charter:(intercept) 0.9164063 0.2072648 4.4214 9.805e-06 ***
private:(intercept) -0.2506806 0.2039395 -1.2292 0.2190004
price -0.0251166 0.0017317 -14.5042 < 2.2e-16 ***
catch 0.3577820 0.1097733 3.2593 0.0011170 **
beach:income 0.1275771 0.0506395 2.5193 0.0117582 *
charter:income 0.0942854 0.0500600 1.8834 0.0596396 .
private:income 0.2170169 0.0500582 4.3353 1.456e-05 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Log-Likelihood: -1215.1
McFadden R^2: 0.18868
Likelihood ratio test : chisq = 565.17 (p.value = < 2.22e-16)
>
>
> # Setting mean values for variables to use for marginal effects
> m <- mlogit(mode ~ price+catch |income, data = mldata, reflevel="charter")
> z <- with(mldata, data.frame(price = tapply(price, index(m)$alt, mean),
+ catch = tapply(catch, index(m)$alt, mean), income =
mean(income)))
>
> # Multinomial logit model marginal effects
> effects(mlogit.model1, covariate = "income", data = z)
charter beach pier private
-1.201367e-02 7.495845e-05 -2.065980e-02 3.259851e-02
>
> # Conditional logit model marginal effects
> effects(clogit.model1, covariate = "income", data = z)
charter beach pier private
-0.0217339246 -0.0007214181 -0.0093059734 0.0317613161
> effects(clogit.model1, covariate = "price", data = z)
charter beach pier private
charter -0.0062430047 6.091542e-04 7.642235e-04 0.0048696270
beach 0.0006091541 -1.249124e-03 8.681094e-05 0.0005531588
pier 0.0007642234 8.681094e-05 -1.545008e-03 0.0006939736
private 0.0048696270 5.531588e-04 6.939737e-04 -0.0061167595
> effects(clogit.model1, covariate = "catch", data = z)
charter beach pier private
charter 0.088930726 -0.008677316 -0.010886256 -0.069367154
beach -0.008677329 0.017793621 -0.001236612 -0.007879681
pier -0.010886272 -0.001236612 0.022008455 -0.009885571
private -0.069367164 -0.007879671 -0.009885559 0.087132394
>
> # Multinomial probit model coefficients
> #mprobit.model1 <- mlogit(mode ~ 1 | income, data = mldata, reflevel="charter",
probit=TRUE)
> #summary(mprobit.model1)
>
>
> # Hauseman-McFadden test of independence of irrelevant alternatives
> m1<- mlogit(mode ~ 1 | income, data = mldata, reflevel="beach")
> m2<- mlogit(mode ~ 1 | income, data = mldata, reflevel="beach",
alt.subset=c("beach", "pier", "private"))
> hmftest(m1, m2)

Hausman-McFadden test

data: mldata
chisq = 14.701, df = 4, p-value = 0.005363
alternative hypothesis: IIA is rejected