Как применить функцию gmnl в R для скрытого анализа класса эксперимента с дискретным выбором?

Question

Я пытаюсь выполнить скрытый анализ класса моих данных из эксперимента с дискретным выбором. Респондентам необходимо было выбрать один из двух вариантов в качестве атрибутов: количество детей, которых они предпочитают, и уровень образования, который они предпочитают для своих детей (указывается как сочетание количества детей). Первые строки моих данных выглядят так:

Respondent Block Choice card Chosen FNoPrimary  FPrimary FSecondary FTertiary MNoPrimary
 1          1           1      1    0.0000000  0.0000000   0.00     0.0000000  0.0000000
 1          1           1      0    0.3333333  0.6666667   0.00     0.0000000  0.0000000
 1          2          12      0    0.3333333  0.3333333   0.00     0.0000000  0.0000000
 1          2          12      1    0.1666667  0.0000000   0.00     0.3333333  0.1666667
 1          3           2      0    0.0000000  0.0000000   1.00     0.0000000  0.0000000
 1          3           2      1    0.0000000  0.0000000   0.25     0.0000000  0.0000000
  MPrimary MSecondary MTertiary NChildren Age District   Religion Indigenous Ethnic group    Sex
1        0       1.00 0.0000000         1  18        0 Protestant          0      Wolaita Female
2        0       0.00 0.0000000         3  18        0 Protestant          0      Wolaita Female
3        0       0.00 0.3333333         9  18        0 Protestant          0      Wolaita Female
4        0       0.00 0.3333333        12  18        0 Protestant          0      Wolaita Female
5        0       0.00 0.0000000         1  18        0 Protestant          0      Wolaita Female
6        0       0.25 0.5000000         4  18        0 Protestant          0      Wolaita Female
       Educational level Studentornot Farmerornot Marital status Having children Ever used contraception
1 High school - grade 10            1           0              0               0                       0
2 High school - grade 10            1           0              0               0                       0
3 High school - grade 10            1           0              0               0                       0
4 High school - grade 10            1           0              0               0                       0
5 High school - grade 10            1           0              0               0                       0
6 High school - grade 10            1           0              0               0                       0
  Alternative
1           1
2           2
3           1
4           2
5           1
6           2

Я посмотрел на все пакеты, доступные в R, и думаю, что только пакет gmnl может обрабатывать данные моего типа и может добавлять ковариаты. Однако, если я сравниваю результаты моего анализа скрытого класса для простой линейной модели только с двумя ковариатами (возраст и район) (как указано ниже), я получаю совершенно другой результат для оценок параметров и логарифмического правдоподобия (-2598,6 в Statavs -2495.5 в R), когда я выполняю тот же анализ со Stata (см. код ниже).

в R:

defining_data <- mlogit.data(final_data_alternativeadded, id.var = "Respondent", choice = "Chosen", alt.var = "Alternative", chid.var="Choice.card", group.var = "Block", varying = 7:15, shape = "long")
mnl <- gmnl(Chosen ~ 1 + x2 + x3 + x4 + x5 + x6 + x7 + x8 + z | 0 | 0 | 0 | Age + District, data = defining_data, model = 'lc', Q = 3) 
summary(mnl)

в Stata:

ssc install lclogit
ssc install fmlogit
lclogit chosen fprimary fsecondary ftertiary mnoprimary mprimary msecondary mtertiary block nchildren, group(choicecard) id(respondent) nclasses(3) membership(age district)

Я попытался сделать все мои переменные числовыми, чтобы умножить пропорции смеси на количество детей вполучить значения, которые ближе друг к другу, чтобы упорядочить мой набор данных на основе значения числа карт выбора на респондентов ... но я всегда получаю другие значения для вероятностей скрытого класса. Кто-то знает почему? Я знаю, что lclogit в Stata использует алгоритм максимизации ожидания, а gmnl в R - метод максимального правдоподобия, но я не думаю, что из-за этого оценки параметров и логарифмическое правдоподобие могут полностью отличаться.

Я также пробовал с этимс существующим открытым набором данных:

Stata:

use http://fmwww.bc.edu/repec/bocode/t/traindata.dta, clear
ssc install lclogit
ssc install fmlogit
lclogit y price contract local wknown tod seasonal, group(gid) id(pid) nclasses(3) seed(12345)

R:

getwd()
traindata <- read.xlsx("traindata.xlsx", 1, header = TRUE)
library(mlogit)
library(gmnl)
traindata[1:nrow(traindata),11] <- seq(1,4) #to add a column with the alternatives per choice numbered from 1 to 4
names(traindata) <- c('y', 'price', 'contract', 'local', 'wknown', 'tod', 'seasonal', 'gid', 'pid', 'X_xi', 'Alternative')
TM <- mlogit.data(traindata, choice = "y", id.var = "pid", alt.var = "Alternative", chid.var = "gid", shape = "long")
mnl <- gmnl(y ~ price + contract + local + wknown + tod + seasonal | 0 | 0 | 0 | 1, data = TM, model = 'lc', Q = 3)
summary(mnl)

Однако Stata заявляет, что они становятся логарифмической вероятностью -1117,9997, а R становится -1329.5, и оба имеют совершенно разные оценки параметров.

Кто-нибудь знает, почему это так?

Большое спасибо заранее

С уважением

Ева

Answer 1

Здесь есть пара вещей, о которых стоит подумать.

Можно ожидать некоторых незначительных различий между ними, потому что i) модели скрытого класса подвержены сходимости к различным локальным оптимам в зависимости от начальных значений и ii) алгоритм максимизации ожидания с большей вероятностью сравнивается с максимальной вероятностью. Тем не менее, это вряд ли имеет место в данном случае.

Вероятно, проблема заключается в том, что Stata учитывает структуру панели данных, тогда как в R вы указали модель классов с поперечным сечением для каждого из них. наблюдение как независимое. Большая разница в значении LL является сильным индикатором этого.

Попробуйте вместо этого (с panel = TRUE):

mnl <- gmnl(y ~ price + contract + local + wknown + tod + seasonal | 0 | 0 | 0 | 1,
 data = TM, model = 'lc', Q = 3, panel = TRUE)