Эффекты от полиномиальной логистической модели в млогите

Question

Я получил некоторую хорошую помощь в правильном форматировании моих данных для создания многочленной логистической модели с mlogit здесь ( Форматирование данных для mlogit)

Однако сейчас я пытаюсь проанализировать влияние ковариат в моей модели. Я нахожу файл справки в mlogit.effects() быть не очень информативным. Одна из проблем заключается в том, что модель, по-видимому, производит много строк NA (см. Ниже, index(mod1)).

1. Кто-нибудь может уточнить, почему мои данные дают эти NA?
2. Может ли кто-нибудь помочь мне получить mlogit.effects работать с данными ниже?
3. Я хотел бы рассмотреть смещение анализа на 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)

Вот моя попытка настроить фрейм данных, аналогичный тому, который изображен в файле справки. Это не работает. Признаюсь, хотя я довольно хорошо знаю семью, к которой я обращаюсь, это непристойно для меня.

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)

Answer 1

Вы работаете с ранжированными данными, а не только с данными с многочленным выбором. Структура для ранжированных данных в mlogit состоит в том, что первым набором записей для человека являются все параметры, затем вторым - все параметры, кроме того, который занимает первое место и т. Д. Но индекс предполагает одинаковое количество вариантов каждый раз. Так что куча АН. Нам просто нужно избавиться от них.

> 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

Answer 2

Я попробую вариант 3 и переключусь на multinom(), Этот код будет моделировать лог-шансы ранжирования элемента как 1-го по сравнению со ссылочным элементом (например, "Дебаты" в приведенном ниже коде). При K = 7 элементах, если мы вызываем контрольный элемент Item_K, то мы моделируем

log[ Pr(элемент_k 1-й) / Pr(элемент_K 1-й)] = α_k + x^Tβ_k

для k = 1,...,K-1, где Item_k является одним из других (то есть не ссылочных) элементов. Выбор эталонного уровня повлияет на коэффициенты и их интерпретацию, но не повлияет на прогнозируемые вероятности. (Та же история для контрольных уровней для категориальных переменных предикторов.)

Я также упомяну, что здесь я обрабатываю пропущенные данные немного иначе, чем в исходном коде. Поскольку моей модели нужно только знать, какой элемент занимает первое место, мне нужно только выбросить записи, где эта информация отсутствует. (Например, в исходной записи № 43 набора данных "Информация" заняла 1-е место, поэтому мы можем использовать эту запись, даже если 3 других элемента - NA.)

# 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

Я не проводил никакой диагностики, поэтому я не утверждаю, что используемая здесь модель хорошо подходит.