嵌套模型假设因变量的分类是分层次进行的,即一层嵌套一层中。
在条件Logit的设定中,假设Y_i=j,j=0,1,2,3,…,M的各种可能的选择之间是相互独立的。但实际中,可能出现这样的情况,即各种可能的选择之间可以分成若干组,组与组之间是相互独立的,但组内的选择之间却是相互关联的。这时,则需要采用嵌套模型来估计各种选择的概率。
研究了关于从悉尼到墨尔本的旅行交通方式选择的问题。市民出行时交通选择方式(mode)首先面临选择是公共交通还是私人交通方式,其中私人交通方式有air和car,公共交通方式有train和bus。该数据共有840个观测值。什么因素会影响市民出行时交通方式的选择呢?
又比如,一个高中毕业生首先面临两种选择:不上大学和上大学。在上大学的选择中又存在着上公立学校和私立学校的选择。也就是说,它面临的是三种选择:不上大学、上公立学校或上私立大学。后两种选择与第一种选择之间是相互独立的,但是后两者之间却是相关联的。
从悉尼到墨尔本的旅行交通方式选择的数据在R的AER包中有对应的数据TravelMode。wait是等待交通工具的时间;vcost是交通工具的成本;travel是旅行的时间;gcost是总的成本;income是家庭收入;size是人数规模。
mdata <- read.csv('http://data.galaxystatistics.com/blog_data/regression/TravelMode.csv',header=T)
head(mdata)
## individual mode choice wait vcost travel gcost income size
## 1 1 air no 69 59 100 70 35 1
## 2 1 train no 34 31 372 71 35 1
## 3 1 bus no 35 25 417 70 35 1
## 4 1 car yes 0 10 180 30 35 1
## 5 2 air no 64 58 68 68 30 2
## 6 2 train no 44 31 354 84 30 2
TravelMode <- mdata[-1]
mdata <- read.csv('http://data.galaxystatistics.com/blog_data/regression/TravelMode.csv',header=T)
# head(mdata)
TravelMode <- mdata[-1]
library(AER)
library(mlogit)
#####分析前,需要对数据进行一些转换,转换过程中使用with()函数,使用方法是with(data,expr,…),即对data进行表达式运算。
###数据转换###
TravelMode$avincome <- with(TravelMode, income*(mode=='air'))
TravelMode$time <- with(TravelMode, travel + wait)/60
TravelMode$timeair <- with(TravelMode, time*I(mode=='air'))
TravelMode$income <- with(TravelMode, income/10)
TravelMode$incomeother <- with(TravelMode, ifelse(mode %in% c('air','car'), income, 0))
dim(TravelMode)
## [1] 840 12
head(TravelMode)
## mode choice wait vcost travel gcost income size avincome time
## 1 air no 69 59 100 70 3.5 1 35 2.816667
## 2 train no 34 31 372 71 3.5 1 0 6.766667
## 3 bus no 35 25 417 70 3.5 1 0 7.533333
## 4 car yes 0 10 180 30 3.5 1 0 3.000000
## 5 air no 64 58 68 68 3.0 2 30 2.200000
## 6 train no 44 31 354 84 3.0 2 0 6.633333
## timeair incomeother
## 1 2.816666.... 3.5
## 2 0 0.0
## 3 0 0.0
## 4 0 3.5
## 5 2.2 3.0
## 6 0 0.0
class(TravelMode)
## [1] "data.frame"
typeof(TravelMode)
## [1] "list"
mode(TravelMode)
## [1] "list"
#####利用嵌套Logit模型进行分析,在R里仍然使用mlogit()函数分析,但是在nests里需要进行设置。
###拟合估算:方式方法1###
nl <- mlogit(choice ~ gcost + wait + incomeother, TravelMode, shape='long', alt.var='mode', nests=list(public=c('train','bus'), other=c('car','air')))
summary(nl)
##
## Call:
## mlogit(formula = choice ~ gcost + wait + incomeother, data = TravelMode,
## nests = list(public = c("train", "bus"), other = c("car",
## "air")), shape = "long", alt.var = "mode")
##
## Frequencies of alternatives:
## air bus car train
## 0.27619 0.14286 0.28095 0.30000
##
## bfgs method
## 12 iterations, 0h:0m:0s
## g'(-H)^-1g = 1.14E-07
## gradient close to zero
##
## Coefficients :
## Estimate Std. Error t-value Pr(>|t|)
## bus:(intercept) -0.7735521 0.8089799 -0.9562 0.3389677
## car:(intercept) -6.1536851 1.1783823 -5.2221 1.769e-07 ***
## train:(intercept) 0.0056402 0.6363849 0.0089 0.9929285
## gcost -0.0195488 0.0060331 -3.2402 0.0011943 **
## wait -0.1064623 0.0205635 -5.1773 2.252e-07 ***
## incomeother 0.4256608 0.1131906 3.7606 0.0001695 ***
## iv.public 0.9694958 0.3047966 3.1808 0.0014687 **
## iv.other 1.7244022 0.5186970 3.3245 0.0008858 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Log-Likelihood: -188.43
## McFadden R^2: 0.33594
## Likelihood ratio test : chisq = 190.65 (p.value = < 2.22e-16)