#######################################################################################
#                                                                                     
#       V O L B A   D O P R A V N Í H O   P R O S T Ř E D K U   N A   C E S T U       
#       P O   B R N Ě N S K É   M E T R O P O L I T N Í   O B L A S T I              
#	                                                                                 
#       Skript ke knize                                                               
#	                                                                                  
#	VÝZKUM DISKRÉTNÍCH VOLEB:                                                   
#	Průvodce designem experimentu,                                               
#       specifikací modelu a jeho odhadem                                             
#	se skripty pro R                                                             
#	                                                                                    
#	2018   Michal Šimeček                                                        
#	                                                                                    
#######################################################################################

library("mlogit")
rm(list=ls())     # smaze promene

###############################################################
#         Funkce pro výpočet kvantifikace
###############################################################
quantif=function(model,refname) {
  km=data.frame(val=coef(model),
                se=sqrt(diag(vcov(model,"coefficient",type=c("cov")))),
                cov=vcov(model,"coefficient",type=c("cov"))[refname,])
  refnum=which(names(coef(model))==refname)
  km$qval=km$val/km$val[refnum]
  km$qse=sqrt( (1/km$val[refnum]*km$se)^2
               +(km$qval/km$val[refnum]*km$se[refnum])^2
               -2*km$qval/km$val[refnum]^2*km$cov         )
  return(km)
}
###############################################################

#######   Načtení a úprava dat   #########

d  = read.table("data_modechoice.dat",header = T)
md = mlogit.data(d,"Choice",shape="wide",varying=6:9,sep="_")
md = md[(md$cavl==1)|(md$alt!="c"),]
md = md[(md$bavl==1)|(md$alt!="b"),]
md$b = ifelse(md$alt=="b",1,0)
md$c = ifelse(md$alt=="c",1,0)
md$p = ifelse(md$alt=="p",1,0)
md$w = ifelse(md$alt=="w",1,0)
md$last = ifelse(md$alt==md$last,1,0)

#######   Multinomial logit    ##########

mn=mlogit(Choice ~ I(cost*p)+I(pcost*c) + I(trans*p) + I(parkb*b) + I(infraw*w) + I(infrab*b) + tt| 1 | last , md)
summary(mn)

#######   Nested logit   ########

nl=mlogit(Choice ~ I(cost*p)+I(pcost*c) + I(trans*p) + I(parkb*b) + I(infraw*w) + I(infrab*b) + tt | 1 | last, md, nests=list(MT=c("c","p"),NMT=c("b","w")),un.nest.el=T)
summary(nl)

(coef(nl)['iv'] - 1) / sqrt(vcov(nl)['iv', 'iv']) # test vyznamnosti nested struktury

lrtest(mn, nl) # likelihood-ratio test

AIC(mn)
AIC(nl)

quantif(nl,"I(cost * p)") # kvantifikace parametru

#######    Model pro Brno a zázemí Brna    #########

mdb=md[md$Brno==1,]
mdz=md[md$Brno!=1,]

nlB=mlogit(Choice ~ I(cost*p)+I(pcost*c) + I(trans*p) + I(parkb*b) + I(infraw*w) + I(infrab*b) + tt | 1 | last , mdb, nests=list(MT=c("c","p"),NMT=c("b","w")),un.nest.el=T)
summary(nlB)

nlZ=mlogit(Choice ~ I(cost*p)+I(pcost*c) + I(trans*p) + I(parkb*b) + I(infraw*w) + I(infrab*b) + tt| 1 | last , mdz, nests=list(MT=c("c","p"),NMT=c("b","w")),un.nest.el=T)
summary(nlZ)

#######    kvantifikace    #########

quantif(nlB,"I(cost * p)")
quantif(nlZ,"I(cost * p)")

#################################################
#            Inkrementalni model
#################################################
incremental = function(share,alt,beta,dx) {
  inc = share
  inc[alt]=inc[alt]*exp(beta*dx)
  inc=inc/sum(inc)
  return (inc)
}
#################################################

nlB2=mlogit(Choice ~ I(cost*p)+I(pcost*c) + I(trans*p) + I(parkb*b) + I(infraw*w) + I(infrab*b) + tt | 1 | last , mdb, nests=list(MT=c("c","p"),NMT=c("b","w")),un.nest.el=T,unscaled=T)
summary(nlB2)

nlZ2=mlogit(Choice ~ I(cost*p)+I(pcost*c) + I(trans*p) + I(parkb*b) + I(infraw*w) + I(infrab*b) + tt| 1 | last , mdz, nests=list(MT=c("c","p"),NMT=c("b","w")),un.nest.el=T,unscaled=T)
summary(nlZ2)

incremental(c(0.38,0.43,0.01,0.18),2,-0.0648375,-10) # Brno   - prac. den
incremental(c(0.51,0.29,0.0,0.20), 2,-0.0648375,-10) # Brno   - víkend
incremental(c(0.61,0.21,0.01,0.17),2,-0.0477251,-10) # zázemí - prac. den
incremental(c(0.66,0.12,0.01,0.21),2,-0.0477251,-10) # zázemí - víkend

################################################
#   Různá hodnota času podle dopravních módů
################################################

nl2=mlogit(Choice ~ I(cost*p)+I(pcost*c) + I(trans*p) + I(parkb*b) + I(infraw*w) + I(infrab*b) | 1 | last + tt, md, nests=list(MT=c("c","p"),NMT=c("b","w")),un.nest.el=T)
summary(nl2)
quantif(nl2,"I(cost * p)")

