### R code for paper ################################################### ### Setup ################################################### options(digits = 3) ################################################### ### Initial Example ################################################### library("gmnl") data("TravelMode", package = "AER") library("mlogit") TM <- mlogit.data(TravelMode, choice = "choice", shape = "long", alt.levels = c("air", "train", "bus", "car")) mixl <- gmnl(choice ~ vcost + travel + wait | 1, data = TM, model = "mixl", ranp = c(travel = "n"), R = 50) summary(mixl) ################################################### ### Travel Mode Data ################################################### data("TravelMode", package = "AER") with(TravelMode, prop.table(table(mode[choice == "yes"]))) head(TravelMode) ################################################### ### Transforming the data ################################################### library("mlogit") TM <- mlogit.data(TravelMode, choice = "choice", shape = "long", alt.levels = c("air", "train", "bus", "car")) ################################################### ### Data analysis ################################################### wide_TM <- reshape(TravelMode, idvar = c("individual", "income", "size"), timevar = "mode", direction = "wide") wide_TM$chosen_mode[wide_TM$choice.air == "yes"] <- "air" wide_TM$chosen_mode[wide_TM$choice.car == "yes"] <- "car" wide_TM$chosen_mode[wide_TM$choice.train == "yes"] <- "train" wide_TM$chosen_mode[wide_TM$choice.bus == "yes"] <- "bus" library("plyr") ddply(wide_TM, ~chosen_mode, summarize, mean.income = mean(income)) ddply(wide_TM, ~chosen_mode, summarize, mean.air = mean(vcost.air), mean.car = mean(vcost.car), mean.train = mean(vcost.train), mean.bus = mean(vcost.bus)) ################################################### ### SMNL ################################################### smnl.nh <- gmnl(choice ~ wait + vcost + travel | 1, data = TM, model = "smnl", R = 30, notscale = c(1, 1, 1, rep(0, 3))) summary(smnl.nh) ################################################### ### SMNL with observed heterogeneity ################################################### smnl.het <- gmnl(choice ~ wait + vcost + travel | 1 | 0 | 0 | income - 1, data = TM, model = "smnl", R = 30, notscale = c(1, 1, 1, 0, 0, 0), typeR = FALSE) summary(smnl.het) ################################################### ### Test ################################################### library("lmtest") waldtest(smnl.nh, smnl.het) lrtest(smnl.nh, smnl.het) ################################################### ### AIC and BIC ################################################### AIC(smnl.nh) AIC(smnl.het) BIC(smnl.nh) BIC(smnl.het) ################################################### ### MIXL ################################################### mixl.hier <- gmnl(choice ~ vcost + travel + wait | 1 | 0 | income + size - 1, data = TM, model = "mixl", ranp = c(travel = "t", wait = "n"), mvar = list(travel = c("income", "size"), wait = c("income")), R = 50, haltons = list(primes = c(2, 17), drop = rep(19, 2))) summary(mixl.hier) ################################################### ### Electricity data ################################################### data("Electricity", package = "mlogit") Electr <- mlogit.data(Electricity, id.var = "id", choice = "choice", varying = 3:26, shape = "wide", sep = "") ################################################### ### Correlated random parameters ################################################### Elec.cor <- gmnl(choice ~ pf + cl + loc + wk + tod + seas | 0, data = Electr, subset = 1:3000, model = "mixl", R = 50, panel = TRUE, ranp = c(cl = "n", loc = "n", wk = "n", tod = "n", seas = "n"), correlation = TRUE) summary(Elec.cor) vcov(Elec.cor, what = "ranp", type = "cov", se = "true") vcov(Elec.cor, what = "ranp", type = "sd", se = "true") vcov(Elec.cor, what = "ranp", type = "cor") ################################################### ### GMNL ################################################### Electr$asc2 <- as.numeric(Electr$alt == 2) Electr$asc3 <- as.numeric(Electr$alt == 3) Electr$asc4 <- as.numeric(Electr$alt == 4) Elec.gmnl <- gmnl(choice ~ pf + cl + loc + wk + tod + seas + asc2 + asc3 + asc4 | 0, data = Electr, subset = 1:3000, model = "gmnl", R = 50, panel = TRUE, notscale = c(rep(0, 6), 1, 1, 1), ranp = c(cl = "n", loc = "n", wk = "n", tod = "n", seas = "n", asc2 = "n", asc3 = "n", asc4 = "n")) summary(Elec.gmnl) ################################################### ### SMNL with random ASCs ################################################### Elec.smnl.re <- gmnl(choice ~ pf + cl + loc + wk + tod + seas + asc2 + asc3 + asc4 | 0, data = Electr, subset = 1:3000, model = "gmnl", R = 50, panel = TRUE, print.init = TRUE, notscale = c(rep(0, 6), 1, 1, 1), ranp = c(asc2 = "n", asc3 = "n", asc4 = "n"), init.gamma = 0, fixed = c(rep(FALSE, 16), TRUE), correlation = TRUE) summary(Elec.smnl.re) ################################################### ### Latent class ################################################### Elec.lc <- gmnl(choice ~ pf + cl + loc + wk + tod + seas | 0 | 0 | 0 | 1, data = Electr, subset = 1:3000, model = "lc", panel = TRUE, Q = 2) summary(Elec.lc) ################################################### ### MM-MIXL ################################################### Elec.mm <- gmnl(choice ~ pf + cl + loc + wk + tod + seas | 0 | 0 | 0 | 1, data = Electr, subset = 1:3000, model = "mm", R = 50, panel = TRUE, ranp = c(pf = "n", cl = "n", loc = "n", wk = "n", tod = "n", seas = "n"), Q = 2, iterlim = 500) summary(Elec.mm) ################################################### ### MM-MIXL with correlated parameters ################################################### Elec.mm.c <- gmnl(choice ~ pf + cl + loc + wk + tod + seas | 0 | 0 | 0 | 1, data = Electr, subset = 1:3000, model = "mm", R = 50, panel = TRUE, ranp = c(pf = "n", cl = "n", loc = "n", wk = "n", tod = "n", seas = "n"), Q = 2, iterlim = 500, correlation = TRUE) summary(Elec.mm.c) vcov(Elec.mm.c, what = "ranp", Q = 1, type = "sd", se = TRUE) vcov(Elec.mm.c, what = "ranp", Q = 2, type = "sd", se = TRUE) ################################################### ### Conditional Logit ################################################### clogit <- gmnl(choice ~ pf + cl + loc + wk + tod + seas | 0, data = Electr, subset = 1:3000) summary(clogit) ################################################### ### WTP ################################################### wtp.gmnl(clogit, wrt = "pf") ################################################### ### Electricity data ################################################### ElectrO <- mlogit.data(Electricity, id = "id", choice = "choice", varying = 3:26, shape = "wide", sep = "", opposite = c("pf")) ################################################### ### WTP-space ################################################### start <- c(1, 0, 0, 0, 0, 0, 0, 0) wtps <- gmnl(choice ~ pf + cl + loc + wk + tod + seas | 0 | 0 | 0 | 1, data = ElectrO, model = "smnl", subset = 1:3000, R = 1, fixed = c(TRUE, FALSE, FALSE, FALSE, FALSE, FALSE, TRUE, FALSE), panel = TRUE, start = start, method = "bhhh", iterlim = 500) summary(wtps) ################################################### ### Price coefficient ################################################### -exp(coef(wtps)["het.(Intercept)"]) ################################################### ### SE of price coefficient ################################################### library("msm") estmean <- coef(wtps) estvar <- vcov(wtps) se <- deltamethod(~-exp(x6), estmean, estvar, ses = TRUE) se ################################################### ### WTP-space with random parameters ################################################### start2 <- c(1, coef(wtps), rep(0.1, 5), 0.1, 0) wtps2 <- gmnl(choice ~ pf + cl + loc + wk + tod + seas | 0 | 0 | 0 | 1, data = ElectrO, subset = 1:3000, model = "gmnl", R = 50, fixed = c(TRUE, rep(FALSE, 12), TRUE), panel = TRUE, start = start2, ranp = c(cl = "n", loc = "n", wk = "n", tod = "n", seas = "n")) summary(wtps2) ################################################### ### WTP-space with correlated random parameters ################################################### n_ran <- 5 start3 <- c(1, coef(wtps), rep(0.1, 0.5 * n_ran * (n_ran + 1)), 0.1, 0) wtps3 <- gmnl(choice ~ pf + cl + loc + wk + tod + seas | 0 | 0 | 0 | 1, data = ElectrO, subset = 1:3000, model = "gmnl", R = 50, fixed = c(TRUE, rep(FALSE, 22), TRUE), panel = TRUE, start = start3, ranp = c(cl = "n", loc = "n", wk = "n", tod = "n", seas = "n"), correlation = TRUE) summary(wtps3) ################################################### ### Plot 1 ################################################### plot(Elec.cor, par = "loc", effect = "ce", type = "density", col = "grey") ################################################### ### Plot 2 ################################################### plot(Elec.cor, par = "loc", effect = "ce", ind = TRUE, id = 1:30) ################################################### ### Using effect.gmnl ################################################### bi.loc <- effect.gmnl(Elec.cor, par = "loc", effect = "ce") summary(bi.loc$mean) summary(bi.loc$sd.est) wtp.loc <- effect.gmnl(Elec.cor, par = "loc", effect = "wtp", wrt = "pf") summary(wtp.loc$mean) summary(wtp.loc$sd.est)