#########################################################################################
#                                                                                       
#         V Ý B Ě R   P R O N Á J M U   S T U D E N T S K É H O   B Y D L E N Í         
#	                                                                                     
#	  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
########################################################

kvantif=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$kval=km$val/km$val[refnum]
  km$kse=sqrt( (1/km$val[refnum]*km$se)^2
               +(km$kval/km$val[refnum]*km$se[refnum])^2
               -2*km$kval/km$val[refnum]^2*km$cov         )
  return(km)
}
#########################################################

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

d=read.table("data_privat.dat",header = T)
d=d[d$ok==1,]
d$zariz_1=as.factor(d$zariz_1)
d$zariz_2=as.factor(d$zariz_2)
d$zariz_3=as.factor(d$zariz_3)
d$zariz_4=as.factor(d$zariz_4)
d$posob_1=as.factor(d$posob_1)
d$posob_2=as.factor(d$posob_2)
d$posob_3=as.factor(d$posob_3)
d$posob_4=as.factor(d$posob_4)
md=mlogit.data(d,"Choice",shape="wide",varying=7:22,sep="_",id.var = "id")


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

m1=mlogit(Choice ~ najem + posob + vzdal + zariz + 0, md)
summary(m1)
LL0 = log(1/4)*nrow(d)
cat("R2(m1) = ",1-as.numeric(m1$logLik)/LL0,"\n")
AIC(m1)

#######   Mixed logit    ##########

m2=mlogit(Choice ~ najem + posob + vzdal + zariz + 0, md, panel = T,rpar=c(najem="n",posob3="n",posob2="n"), R=1000,
          start = c(c(m1$coefficients),0.0,0.0,0.0))
summary(m2)
cat("R2(m2) = ",1-as.numeric(m2$logLik)/LL0,"\n")
AIC(m2)

lrtest(m1, m2) # likelihood-ratio test

######## Kvantifikace   ###########

kvantif(m1,"najem")
kvantif(m2,"najem")
