Поместите биноминальную смесь во flexmix - что я делаю не так? - PullRequest
2 голосов
/ 01 октября 2019

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

set.seed(42)
data=data.frame(class=rbinom(1000,size=1,prob=0.3)) %>% # 30% in class 1, 70% in class 0
  mutate(yes=map_dbl(class,~ifelse(.x,rbinom(1,1,prob=0.8),rbinom(1,1,prob=0.4)))) # class 1 = 80%, class 2 = 40%

# Algebraic
(mean(data$yes==1)-0.4)/(0.8-0.4)
# = 0.295 this is what the model should recover, right?


library(flexmix)
mo1=FLXMRglm(offset=rep(log(.4/.6), nrow(data)),family="binomial")
mo2=FLXMRglm(offset=rep(log(.8/.2), nrow(data)),family="binomial")
flexfit <- flexmix(cbind(yes,1-yes) ~ -1, data=data, k=2, model=list(mo1, mo2))
flexfit
summary(flexfit)

Этот код возвращаетсявсе точки данных принадлежат первому классу. Я неправильно настраиваю модель? Или у меня есть более фундаментальное недоразумение о том, как работают смешанные модели?

1 Ответ

0 голосов
/ 04 октября 2019

Поскольку у вас нет регрессоров, я думаю, что вам необходимо указать перехват в вашей линейной модели, то есть cbind(yes, 1 - yes) ~ 1.

flexfit <- initFlexmix(cbind(yes, 1 - yes) ~ 1, data=data, k=2, model=list(mo1, mo2))
summary(flexfit)

разбивает набор данных на кластеры по 482 и 518 единиц. Модель имеет логарифмическое правдоподобие = -692,499, AIC = 1390,998, BIC = 1419,537.

Обратите внимание, что смеси биномиальных распределений в общем случае не идентифицируемы. Это не совпадение, что в группе 2 518 единиц: это именно то количество наблюдений, для которых yes == 1. Если мы передадим метки модели, мы, очевидно, получим точный ответ (что делает модель совершенно бесполезной):

data$label <- ifelse(data$class == 1, 2, 1)
flexfit <- initFlexmix(cbind(yes, 1 - yes) ~ 1 | label, data=data, k=2, model=list(mo1, mo2))
summary(flexfit)

Я надеюсь, что кто-то может доказать, что я неправ, но я не думаю, что есть способвосстановить скрытые группы, наблюдая только результаты 0/1 и без сильных предварительных предположений.

Подумайте об этом следующим образом: вы можете классифицировать ваши двоичные наблюдения на два кластера разными способами (общее число - это число Стирлингавторого рода), и каждое выделение является вполне разумным объяснением ваших данных. Оценка максимального правдоподобия является просто наиболее разумной: оцените вероятность скрытого класса по наблюдаемой относительной частоте, затем установите условные вероятности равными 0 или 1 в зависимости от класса.

Pr (Y= 1 | π, θ, ϕ) = Pr (C = 1) Pr (Y = 1 | C = 1) + Pr (C = 2) Pr (Y = 1 | C = 2) = πθ + (1 -π) ϕ

Установите π на наблюдаемую относительную частоту, ϕ = 1 и θ = 0 . Модель flexmix более сложна, чем эта, но я уверен, что ее можно свести к моему примеру, учитывая, что у нас есть только перехват.

...