mlogit 中多项逻辑模型的影响
Effects from multinomial logistic model in mlogit
我得到了一些很好的帮助，让我的数据正确格式化，在此处使用 mlogit 生成多项逻辑模型 ()
但是，我现在正在尝试分析模型中协变量的影响。我发现 mlogit.effects() 中的帮助文件信息量不大。问题之一是该模型似乎产生了很多行的 NA（见下文，index(mod1)）。
- 谁能解释为什么我的数据会产生这些 NA？
- 谁能帮我
mlogit.effects 使用下面的数据？
- 我会考虑将分析转移到
multinom()。但是，我不知道如何格式化数据以适合使用 multinom() 的公式。我的数据是七个不同项目（可访问、信息、权衡、辩论、社交和响应）的一系列排名，我会只对他们选择的第一排名进行建模，而忽略他们在其他排名中的选择吗？我可以得到那个信息。
可重现的代码如下：
#Loadpackages
library(RCurl)
library(mlogit)
library(tidyr)
library(dplyr)
#URL where data is stored
dat.url <- 'https://raw.githubusercontent.com/sjkiss/Survey/master/mlogit.out.csv'
#Get data
dat <- read.csv(dat.url)
#Complete cases only as it seems mlogit cannot handle missing values or tied data which in this case you might get because of median imputation
dat <- dat[complete.cases(dat),]
#Change the choice index variable (X) to have no interruptions, as a result of removing some incomplete cases
dat$X <- seq(1,nrow(dat),1)
#Tidy data to get it into long format
dat.out <- dat %>%
gather(Open, Rank, -c(1,9:12)) %>%
arrange(X, Open, Rank)
#Create mlogit object
mlogit.out <- mlogit.data(dat.out, shape='long',alt.var='Open',choice='Rank', ranked=TRUE,chid.var='X')
#Fit Model
mod1 <- mlogit(Rank~1|gender+age+economic+Job,data=mlogit.out)
这是我尝试建立一个类似于帮助文件中描绘的数据框的尝试。它不起作用。我承认虽然我很了解apply家族，但tapply对我来说很模糊。
with(mlogit.out, data.frame(economic=tapply(economic, index(mod1)$alt, mean)))
从帮助中比较：
data("Fishing", package = "mlogit")
Fish <- mlogit.data(Fishing, varying = c(2:9), shape = "wide", choice = "mode")
m <- mlogit(mode ~ price | income | catch, data = Fish)
# compute a data.frame containing the mean value of the covariates in
# the sample data in the help file for effects
z <- with(Fish, data.frame(price = tapply(price, index(m)$alt, mean),
catch = tapply(catch, index(m)$alt, mean),
income = mean(income)))
# compute the marginal effects (the second one is an elasticity
effects(m, covariate = "income", data = z)
我会尝试选项 3 并切换到 multinom()。与参考项目（例如，下面代码中的 "Debate"）相比，此代码将对项目排名第一的对数几率进行建模。对于 K = 7 项，如果我们调用参考项 ItemK，那么我们正在建模
log[ Pr(Itemk 是第一个) / Pr(ItemK 是第一个) ] = α k + xTβk
for k = 1,...,K-1，其中 Itemk 是其他（即非参考）项目之一。参考水平的选择会影响系数及其解释，但不会影响预测概率。 （分类预测变量的参考水平也是如此。）
我还要提到，我在这里处理缺失数据的方式与您的原始代码略有不同。由于我的模型只需要知道哪个项目排名第一，我只需要丢弃缺少该信息的记录。 （例如，在原始数据集中记录 #43 "Information" 排名第一，因此即使其他 3 个项目不适用，我们也可以使用此记录。）
# Get data
dat.url <- 'https://raw.githubusercontent.com/sjkiss/Survey/master/mlogit.out.csv'
dat <- read.csv(dat.url)
# dataframe showing which item is ranked #1
ranks <- (dat[,2:8] == 1)
# for each combination of predictor variable values, count
# how many times each item was ranked #1
dat2 <- aggregate(ranks, by=dat[,9:12], sum, na.rm=TRUE)
# remove cases that didn't rank anything as #1 (due to NAs in original data)
dat3 <- dat2[rowSums(dat2[,5:11])>0,]
# (optional) set the reference levels for the categorical predictors
dat3$gender <- relevel(dat3$gender, ref="Female")
dat3$Job <- relevel(dat3$Job, ref="Government backbencher")
# response matrix in format needed for multinom()
response <- as.matrix(dat3[,5:11])
# (optional) set the reference level for the response by changing
# the column order
ref <- "Debate"
ref.index <- match(ref, colnames(response))
response <- response[,c(ref.index,(1:ncol(response))[-ref.index])]
# fit model (note that age & economic are continuous, while gender &
# Job are categorical)
library(nnet)
fit1 <- multinom(response ~ economic + gender + age + Job, data=dat3)
# print some results
summary(fit1)
coef(fit1)
cbind(dat3[,1:4], round(fitted(fit1),3)) # predicted probabilities
我没有做任何诊断，所以我没有断言这里使用的模型很合适。
您正在使用排名数据，而不仅仅是多项选择数据。 mlogit 中 Ranked 数据的结构是，第一组记录一个人的所有选项，然后第二组是除第一名以外的所有选项，依此类推。但是该指数每次都假定相同数量的选项。所以一堆 NA。我们只需要摆脱它们。
> with(mlogit.out, data.frame(economic=tapply(economic, index(mod1)$alt[complete.cases(index(mod1)$alt)], mean)))
economic
Accessible 5.13
Debate 4.97
Information 5.08
Officials 4.92
Responsive 5.09
Social 4.91
Trade.Offs 4.91
我得到了一些很好的帮助，让我的数据正确格式化，在此处使用 mlogit 生成多项逻辑模型 (
但是，我现在正在尝试分析模型中协变量的影响。我发现 mlogit.effects() 中的帮助文件信息量不大。问题之一是该模型似乎产生了很多行的 NA（见下文，index(mod1)）。
- 谁能解释为什么我的数据会产生这些 NA？
- 谁能帮我
mlogit.effects使用下面的数据？ - 我会考虑将分析转移到
multinom()。但是，我不知道如何格式化数据以适合使用multinom()的公式。我的数据是七个不同项目（可访问、信息、权衡、辩论、社交和响应）的一系列排名，我会只对他们选择的第一排名进行建模，而忽略他们在其他排名中的选择吗？我可以得到那个信息。
可重现的代码如下：
#Loadpackages
library(RCurl)
library(mlogit)
library(tidyr)
library(dplyr)
#URL where data is stored
dat.url <- 'https://raw.githubusercontent.com/sjkiss/Survey/master/mlogit.out.csv'
#Get data
dat <- read.csv(dat.url)
#Complete cases only as it seems mlogit cannot handle missing values or tied data which in this case you might get because of median imputation
dat <- dat[complete.cases(dat),]
#Change the choice index variable (X) to have no interruptions, as a result of removing some incomplete cases
dat$X <- seq(1,nrow(dat),1)
#Tidy data to get it into long format
dat.out <- dat %>%
gather(Open, Rank, -c(1,9:12)) %>%
arrange(X, Open, Rank)
#Create mlogit object
mlogit.out <- mlogit.data(dat.out, shape='long',alt.var='Open',choice='Rank', ranked=TRUE,chid.var='X')
#Fit Model
mod1 <- mlogit(Rank~1|gender+age+economic+Job,data=mlogit.out)
这是我尝试建立一个类似于帮助文件中描绘的数据框的尝试。它不起作用。我承认虽然我很了解apply家族，但tapply对我来说很模糊。
with(mlogit.out, data.frame(economic=tapply(economic, index(mod1)$alt, mean)))
从帮助中比较：
data("Fishing", package = "mlogit")
Fish <- mlogit.data(Fishing, varying = c(2:9), shape = "wide", choice = "mode")
m <- mlogit(mode ~ price | income | catch, data = Fish)
# compute a data.frame containing the mean value of the covariates in
# the sample data in the help file for effects
z <- with(Fish, data.frame(price = tapply(price, index(m)$alt, mean),
catch = tapply(catch, index(m)$alt, mean),
income = mean(income)))
# compute the marginal effects (the second one is an elasticity
effects(m, covariate = "income", data = z)
我会尝试选项 3 并切换到 multinom()。与参考项目（例如，下面代码中的 "Debate"）相比，此代码将对项目排名第一的对数几率进行建模。对于 K = 7 项，如果我们调用参考项 ItemK，那么我们正在建模
log[ Pr(Itemk 是第一个) / Pr(ItemK 是第一个) ] = α k + xTβk
for k = 1,...,K-1，其中 Itemk 是其他（即非参考）项目之一。参考水平的选择会影响系数及其解释，但不会影响预测概率。 （分类预测变量的参考水平也是如此。）
我还要提到，我在这里处理缺失数据的方式与您的原始代码略有不同。由于我的模型只需要知道哪个项目排名第一，我只需要丢弃缺少该信息的记录。 （例如，在原始数据集中记录 #43 "Information" 排名第一，因此即使其他 3 个项目不适用，我们也可以使用此记录。）
# Get data
dat.url <- 'https://raw.githubusercontent.com/sjkiss/Survey/master/mlogit.out.csv'
dat <- read.csv(dat.url)
# dataframe showing which item is ranked #1
ranks <- (dat[,2:8] == 1)
# for each combination of predictor variable values, count
# how many times each item was ranked #1
dat2 <- aggregate(ranks, by=dat[,9:12], sum, na.rm=TRUE)
# remove cases that didn't rank anything as #1 (due to NAs in original data)
dat3 <- dat2[rowSums(dat2[,5:11])>0,]
# (optional) set the reference levels for the categorical predictors
dat3$gender <- relevel(dat3$gender, ref="Female")
dat3$Job <- relevel(dat3$Job, ref="Government backbencher")
# response matrix in format needed for multinom()
response <- as.matrix(dat3[,5:11])
# (optional) set the reference level for the response by changing
# the column order
ref <- "Debate"
ref.index <- match(ref, colnames(response))
response <- response[,c(ref.index,(1:ncol(response))[-ref.index])]
# fit model (note that age & economic are continuous, while gender &
# Job are categorical)
library(nnet)
fit1 <- multinom(response ~ economic + gender + age + Job, data=dat3)
# print some results
summary(fit1)
coef(fit1)
cbind(dat3[,1:4], round(fitted(fit1),3)) # predicted probabilities
我没有做任何诊断，所以我没有断言这里使用的模型很合适。
您正在使用排名数据，而不仅仅是多项选择数据。 mlogit 中 Ranked 数据的结构是，第一组记录一个人的所有选项，然后第二组是除第一名以外的所有选项，依此类推。但是该指数每次都假定相同数量的选项。所以一堆 NA。我们只需要摆脱它们。
> with(mlogit.out, data.frame(economic=tapply(economic, index(mod1)$alt[complete.cases(index(mod1)$alt)], mean)))
economic
Accessible 5.13
Debate 4.97
Information 5.08
Officials 4.92
Responsive 5.09
Social 4.91
Trade.Offs 4.91