Предсказания модели Стэна за пределами возможного диапазона данных

Question

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

Вот данные

df <- data.frame(id = rep(1:46, each = 3),
                 group = c(3,3,3,1,1,1,3,3,3,1,1,1,3,3,3,1,1,1,3,3,3,2,2,2,1,1,1,3,3,3,3,3,3,2,2,2,3,3,3,1,1,1,2,2,2,3,3,3,2,2,2,2,2,2,3,3,3,1,1,1,2,2,2,3,3,3,2,2,2,3,3,3,3,3,3,2,2,2,3,3,3,3,3,3,1,1,1,3,3,3,3,3,3,1,1,1,2,2,2,2,2,2,1,1,1,2,2,2,2,2,2,1,1,1,1,1,1,2,2,2,2,2,2,1,1,1,1,1,1,3,3,3,1,1,1,3,3,3),
                 bevType = rep(c(3,2,1), times = 46),
                 score = c(2.9,1.0,0.0,9.5,5.0,4.5,9.0,3.0,5.0,5.0,0.0,3.0,9.5,2.0,3.0,8.5,0.0,6.0,5.2,3.0,4.0,8.4,7.0,2.0,10.0,0.0,3.0,7.3,1.0,1.8,8.5,2.0,9.0,10.0,5.0,10.0,8.3,2.0,5.0,6.0,0.0,5.0,6.0,0.0,5.0,10.0,0.0,5.0,6.8,1.0,4.8,8.0,1.0,4.0,7.0,4.0,6.0,6.5,1.0,3.1,9.0,1.0,0.0,6.0,0.0,2.0,9.5,4.0,6.0,8.0,1.0,3.8,0.4,0.0,7.0,7.0,0.0,3.0,9.0,2.0,5.0,9.5,2.0,7.0,7.9,5.0,4.9,8.0,1.0,1.0,9.3,5.0,7.9,6.5,2.0,3.0,8.0,2.0,6.0,10.0,0.0,5.0,6.0,0.0,5.0,6.8,0.1,7.0,8.0,3.0,9.1,8.2,0.0,7.9,8.2,5.0,0.0,9.2,1.0,3.1,9.1,3.0,0.6,5.7,2.0,5.1,7.0,0.0,7.4,8.0,1.0,1.5,9.1,4.0,4.3,8.5,8.0,5.0))

Теперь для модели. Модель имеет параметр общего среднего a, категориальный предиктор, представляющий отклонения групп от общего среднего bGroup, термин для отклонений различных типов напитков от общего среднего bBev, термин для перехвата каждого субъекта bSubj и термин для группы по взаимодействию с напитком bGxB. Я также оценил отдельные параметры шума для каждого типа напитка.

Чтобы разрешить задние прогностические проверки, я нарисовал из заднего сустава, используя блок generated quantities и функцию normal_rng.

### Step 1: Put data into list
dList <- list(N = 138,
              nSubj = 46,
              nGroup = 3,
              nBev = 3,
              sIndex = df$id,
              gIndex = df$group,
              bIndex = df$bevType,
              score = df$score,
              gMean = 4.718841,
              gSD = 3.17)


#### Step 1 model
write("
      data{
      int<lower=1> N;
      int<lower=1> nSubj;
      int<lower=1> nGroup;
      int<lower=1> nBev;
      int<lower=1,upper=nSubj> sIndex[N];
      int<lower=1,upper=nGroup> gIndex[N];
      int<lower=1,upper=nBev> bIndex[N];
      real score[N];
      real gMean;
      real gSD;
      }

      parameters{
      real a;
      vector[nSubj] bSubj;
      vector[nGroup] bGroup;
      vector[nBev] bBev;
      vector[nBev] bGxB[nGroup];      // vector of vectors, stan no good with matrix
      vector[nBev] sigma;  
      real<lower=0> sigma_a;
      real<lower=0> sigma_s;
      real<lower=0> sigma_g;
      real<lower=0> sigma_b;
      real<lower=0> sigma_gb;
      }

      model{
      vector[N] mu;

      //hyper-priors
      sigma_s ~ normal(0,10);     
      sigma_g ~ normal(0,10);
      sigma_b ~ normal(0,10);
      sigma_gb ~ normal(0,10);


      //priors
      sigma ~ cauchy(0,1);
      a ~ normal(gMean, gSD);
      bSubj ~ normal(0, sigma_s);
      bGroup ~ normal(0,sigma_g);
      bBev ~ normal(0,sigma_b);
      for (i in 1:nGroup) {               //hierarchical prior on interaction
      bGxB[i] ~ normal(0, sigma_gb);
      }

      // likelihood
      for (i in 1:N){
      score[i] ~ normal(a + bGroup[gIndex[i]] + bBev[bIndex[i]] + bSubj[sIndex[i]] + bGxB[gIndex[i]][bIndex[i]], sigma[bIndex[i]]);
      }
      }

      generated quantities{
      real y_draw[N];
      for (i in 1:N) {
      y_draw[i] = normal_rng(a + bGroup[gIndex[i]] + bBev[bIndex[i]] + bSubj[sIndex[i]] + bGxB[gIndex[i]][bIndex[i]], sigma[bIndex[i]]);
      }
      }
      ", file = "temp.stan")

##### Step 3: generate the chains
mod <-  stan(file = "temp.stan",
             data = dList,
             iter = 5000,  
             warmup = 3000, 
             cores = 1,
             chains = 1)

Далее мы извлекаем данные из задней части сустава и генерируем оценки среднего, верхнего и нижнего 95% HPDI в группе. Для начала нам нужна функция для вычисления HPDI

HPDIFunct <- function (vector) { 
  sortVec <- sort(vector)
  ninetyFiveVec <- ceiling(.95*length(sortVec))
  fiveVec <- length(sortVec) - length(ninetyFiveVec)
  diffVec <- sapply(1:fiveVec, function (i) sortVec[i + ninetyFiveVec] - sortVec[i])
  minVal <- sortVec[which.min(diffVec)]
  maxVal <- sortVec[which.min(diffVec) + ninetyFiveVec]
  return(list(sortVec, minVal, maxVal))
}

Теперь, чтобы извлечь дро из задних

#### Step 5: Posterior predictive checks
y_draw <- data.frame(extract(mod, pars = "y_draw"))


And plot the mean, lower HPDI and upper HPDI draws of these draws against the actual data. 

df$drawMean <- apply(y_draw, 2, mean)
df$HPDI_Low <- apply(y_draw, 2, function(i) HPDIFunct(i)[[2]][1])
df$HPDI_Hi <- apply(y_draw, 2, function(i) HPDIFunct(i)[[3]][1])

### Step 6: plot posterior draws against actual data
ggplot(df, aes(x = factor(bevType), colour = factor(group))) +
       geom_jitter(aes(y = score), shape = 1, position = position_dodge(width=0.9)) +
       geom_point(aes(y = drawMean), position = position_dodge(width=0.9), stat = "summary", fun.y = "mean", shape = 3, size = 3, stroke = 2) +
       geom_point(aes(y = HPDI_Low), position = position_dodge(width=0.9), stat = "summary", fun.y = "mean", shape = 1, size = 3, stroke = 1) +
       geom_point(aes(y = HPDI_Hi), position = position_dodge(width=0.9), stat = "summary", fun.y = "mean", shape = 1, size = 3, stroke = 1) +
       scale_colour_manual(name = "Experimental Group", labels = c("Group 1", "Group 2", "Group 3"), values = c("#616a6b", "#00AFBB", "#E7B800")) +  
       scale_x_discrete(labels = c("Water", "Decaf", "Coffee")) + 
       labs(x = "Beverage Type", y = "Expectancy of Withdrawal Alleviation") +
       scale_y_continuous(breaks = seq(0,10,2)) +
       theme(axis.text.x = element_text(size = 12),
             axis.title.x = element_text(face = "bold"),
             axis.title.y = element_text(face = "bold"),
             axis.text.y = element_text(size = 12),
             legend.title = element_text(size = 13, face = "bold"))

Глядя на график, для ожиданий воды модель, по-видимому, достаточно хорошо представляет центр (крестики) и разброс (незакрашенные кружки) данных. Но это нарушает ожидания без кофеина и кофе. Для ожиданий без кофеина нижний HPDI находится ниже диапазона возможных значений (нижний предел = 0), а разброс ничьих сзади (представлен в каждой группе незакрашенными кружками) слишком велик. Верхний предел HPDI группы Кофе также выше , диапазон данных (верхний предел = 10) и разброс слишком велики для фактических данных.

Итак, мой вопрос:

Как я могу ограничить вытягивания из заднего сустава к фактическому диапазону данных?

Есть ли какой-нибудь способ грубой силы сдерживать ничьи сзади в Стэне? Или более эффективная оценка различий в дисперсии между тремя условиями напитка будет более эффективной (в этом случае это будет больше вопрос CV, чем вопрос SO)?

Answer 1

Стандартный способ ограничения апостериорной переменной - использовать функцию ссылки для ее преобразования. Так работают обобщенные линейные модели (GLM), такие как логистическая регрессия и пуассоновская регрессия. Например, чтобы перейти от положительного к неограниченному, мы используем лог-преобразование. Чтобы перейти от вероятности в (0, 1) к неограниченной, мы используем преобразование лог-шансов.

Если ваши результаты представляют собой порядковые значения по шкале 1-10, общий подход, учитывающий, что шкала данных равна порядковая логистическая регрессия .

Answer 2

Чтобы расширить ответ @Bob Carpenter, вот два способа решить эту проблему.(У меня была причина использовать оба из них недавно, и я изо всех сил пытался их запустить и запустить. Это может быть полезно для других новичков, таких как я.)

Метод 1: Упорядоченная логистическая регрессия

Мы предполагаем, что у каждого пользователя есть «истинное» ожидание для каждого ответа в произвольном непрерывном масштабе, и моделируем его как скрытую переменную.Если фактические ответы пользователя попадают в категории K, мы также моделируем K - 1 контрольные точки между этими категориями.Вероятность того, что пользователь выберет данную категорию ответа, равна площади под pdf логистики между соответствующими точками отсечения.

Модель Стэна выглядит следующим образомэтот.Основное отличие состоит в том, что модель соответствует дополнительному упорядоченному вектору cutpoints и использует распределение ordered_logistic.(Я также изменил приоры на sigma s на Коши, чтобы они оставались положительными, и переключился на нецентрированную параметризацию . Но эти изменения не зависят от рассматриваемого вопроса.)

data {
  int<lower=1> N;
  int<lower=1> nSubj;
  int<lower=1> nGroup;
  int<lower=1> nBev;
  int minResponse;
  int maxResponse;
  int<lower=1,upper=nSubj> sIndex[N];
  int<lower=1,upper=nGroup> gIndex[N];
  int<lower=1,upper=nBev> bIndex[N];
  int<lower=minResponse,upper=maxResponse> score[N];
}

parameters {
  real a;
  vector[nSubj] bSubj;
  vector[nGroup] bGroup;
  vector[nBev] bBev;
  vector[nBev] bGxB[nGroup];
  real<lower=0> sigma_s;
  real<lower=0> sigma_g;
  real<lower=0> sigma_b;
  real<lower=0> sigma_gb;
  ordered[maxResponse - minResponse] cutpoints;
}

model {
  //hyper-priors
  sigma_s ~ cauchy(0, 1);
  sigma_g ~ cauchy(0, 1);
  sigma_b ~ cauchy(0, 1);
  sigma_gb ~ cauchy(0, 1);
  //priors
  a ~ std_normal();
  bSubj ~ std_normal();
  bGroup ~ std_normal();
  bBev ~ std_normal();
  for (i in 1:nGroup) {
    bGxB[i] ~ std_normal();
  }
  // likelihood
  for(i in 1:N) {
    score[i] ~ ordered_logistic(a +
                                  (bGroup[gIndex[i]] * sigma_g) +
                                  (bBev[bIndex[i]] * sigma_b) +
                                  (bSubj[sIndex[i]] * sigma_s) +
                                  (bGxB[gIndex[i]][bIndex[i]] * sigma_gb),
                                cutpoints);
  }
}

generated quantities {
  real y_draw[N];
  for (i in 1:N) {
    y_draw[i] = ordered_logistic_rng(a +
                                       (bGroup[gIndex[i]] * sigma_g) +
                                       (bBev[bIndex[i]] * sigma_b) +
                                       (bSubj[sIndex[i]] * sigma_s) +
                                       (bGxB[gIndex[i]][bIndex[i]] * sigma_gb),
                                     cutpoints);
  }
}

Похоже, что ответы в вашем наборе данных записаны с точностью до десятой, что дает нам 101 возможную категорию от 0 до 10. Чтобы сохранить все как дружественные для Stan целые числа, мы можем умножить все ответы на 10. (Я также добавил один к каждому ответу, потому что у меня были проблемы с подгонкой модели, когда одна из возможных категорий была нулевой.)

dList <- list(N = 138,
              nSubj = 46,
              nGroup = 3,
              nBev = 3,
              minResponse = 1,
              maxResponse = 101,
              sIndex = df$id,
              gIndex = df$group,
              bIndex = df$bevType,
              score = (df$score * 10) + 1)

После того, как мы извлечем y_draw, мы можем преобразовать его обратно в исходный масштаб:

y_draw <- (data.frame(extract(mod, pars = "y_draw")) - 1) / 10

Все остальное такое же, как и раньше.Теперь апостериорные прогнозы правильно ограничены [0, 10].

Метод 2: Бета-регрессия

Как только мы получим до 101Ответные категории, называя эти возможности дискретными категориями, кажутся немного странными.Более естественно сказать, что, как пыталась сделать ваша оригинальная модель, мы фиксируем непрерывный результат, который оказывается ограниченным от 0 до 10. Кроме того, в упорядоченной логистической регрессии категории ответов не должны быть регулярноразнесены по отношению к скрытой переменной.(Это функция, а не ошибка; например, для ответов Лайкерта нет никакой гарантии, что разница между «Полностью согласен» и «Согласен» такая же, как разница между «Согласен» и «Ни согласен, ни согласен»).) в результате трудно что-либо сказать о «расстоянии», когда конкретный фактор вызывает изменение реакции в исходном масштабе (в отличие от масштаба скрытой переменной).Но точки отсчета, определяемые вышеприведенной моделью, довольно равномерно распределены, что снова говорит о том, что результат в вашем наборе данных уже достаточно масштабный:

# Get the sampled parameters
sampled.params.df = data.frame(as.array(mod)[,1,]) %>%
  select(-matches("y_draw")) %>%
  rownames_to_column("iteration")

# Plot selected cutpoints
sampled.params.df %>%
  select(matches("cutpoints")) %>%
  gather(cutpoint, value) %>%
  mutate(cutpoint.num = as.numeric(gsub("^cutpoints\\.([0-9]+)\\.$", "\\1", cutpoint))) %>%
  group_by(cutpoint.num) %>%
  summarize(mean.value = mean(value),
            lower.95 = quantile(value, 0.025),
            lower.50 = quantile(value, 0.25),
            upper.50 = quantile(value, 0.75),
            upper.95 = quantile(value, .975)) %>%
  ggplot(aes(x = cutpoint.num, y = mean.value)) +
  geom_point(size = 3) +
  geom_linerange(aes(ymin = lower.95, ymax = upper.95)) +
  geom_linerange(aes(ymin = lower.50, ymax = upper.50), size = 2) +
  scale_x_continuous("cutpoint", breaks = seq(0, 100, 10)) +
  scale_y_continuous("") +
  theme_bw()

(Толстые и тонкие линии представляют интервалы 50% и 95% соответственно. Я получаю удовольствие от небольшого «скачка» каждые 10 точек отсечения, что говорит о том, что испытуемые рассматривали, скажем, 5,9 против 6,0 как большую разницу, чем 5,8 против5.9 Но эффект кажется довольно слабым. Масштаб также, кажется, немного растягивается к верхнему пределу, но, опять же, он не слишком резок.)

Для непрерывного результата в ограниченном интервале мыможно использовать бета-дистрибутив ;см. здесь и здесь для дальнейшего обсуждения.

Для бета-распределения нам нужны два параметра, mu и phi, оба из которых должны быть положительными,В этом примере я разрешил неограниченное использование mu и применил inv_logit перед подачей его в бета-дистрибутив;Я заставил phi быть позитивным и дал ему Коши.Но вы можете сделать это любым количеством способов.Я также закодировал полный набор параметров mu, но только один phi;Опять же, вы можете поэкспериментировать с другими вариантами.

data {
  int<lower=1> N;
  int<lower=1> nSubj;
  int<lower=1> nGroup;
  int<lower=1> nBev;
  int<lower=1,upper=nSubj> sIndex[N];
  int<lower=1,upper=nGroup> gIndex[N];
  int<lower=1,upper=nBev> bIndex[N];
  vector<lower=0,upper=1>[N] score;
}

parameters {
  real a;
  real a_phi;
  vector[nSubj] bSubj;
  vector[nGroup] bGroup;
  vector[nBev] bBev;
  vector[nBev] bGxB[nGroup];
  real<lower=0> sigma_s;
  real<lower=0> sigma_g;
  real<lower=0> sigma_b;
  real<lower=0> sigma_gb;
}

model {
  vector[N] mu;
  //hyper-priors
  sigma_s ~ cauchy(0, 1);
  sigma_g ~ cauchy(0, 1);
  sigma_b ~ cauchy(0, 1);
  sigma_gb ~ cauchy(0, 1);
  //priors
  a ~ std_normal();
  a_phi ~ cauchy(0, 1);
  bSubj ~ std_normal();
  bGroup ~ std_normal();
  bBev ~ std_normal();
  for (i in 1:nGroup) {
    bGxB[i] ~ std_normal();
  }
  // likelihood
  for(i in 1:N) {
    mu[i] = a +
             (bGroup[gIndex[i]] * sigma_g) +
             (bBev[bIndex[i]] * sigma_b) +
             (bSubj[sIndex[i]] * sigma_s) +
             (bGxB[gIndex[i]][bIndex[i]] * sigma_gb);
    score[i] ~ beta(inv_logit(mu[i]) .* a_phi,
                    (1 - inv_logit(mu[i])) .* a_phi);
  }
}

generated quantities {
  real y_draw[N];
  real temp_mu;
  for (i in 1:N) {
    temp_mu = a +
               (bGroup[gIndex[i]] * sigma_g) +
               (bBev[bIndex[i]] * sigma_b) +
               (bSubj[sIndex[i]] * sigma_s) +
               (bGxB[gIndex[i]][bIndex[i]] * sigma_gb);
    y_draw[i] = beta_rng(inv_logit(temp_mu) .* a_phi,
                         (1 - inv_logit(temp_mu)) .* a_phi);
  }
}

Бета-распределение поддерживается на (0, 1), поэтому мы делим наблюдаемые оценки на 10. (Модель также дает сбой, если мы даем ей оценки ровно 0 или 1, поэтому я преобразовал все такие оценки в 0,01 и 0,99 соответственно .)

dList.beta <- list(N = 138,
                   nSubj = 46,
                   nGroup = 3,
                   nBev = 3,
                   sIndex = df$id,
                   gIndex = df$group,
                   bIndex = df$bevType,
                   score = ifelse(df$score == 0, 0.01,
                                  ifelse(df$score == 10, 0.99,
                                         df$score / 10)))

Отмените преобразование при извлечении y_draw, и тогда процедура будет такой же, как и раньше.

y_draw.beta <- data.frame(extract(mod.beta, pars = "y_draw")) * 10

Еще раз, задние ничьи правильно ограничены.