Внедрение экспоненциальной общей линейной модели в Stan / RStan

The Pointer спросил: 12 мая 2018 в 03:57 в: r

У меня есть следующая модель:

Буду признателен, если кто-нибудь сможет продемонстрировать, как реализовать эту модель.

Исследование и попытка

В качестве примера я рассмотрел следующий пример Poisson GLM:

```{stan output.var="PoissonGLMQR"}
  data{
    int<lower=1> n;     // number of observations
    int<lower=1> p;     // number of predictors
    matrix[n, p] x;     // design matrix
    int<lower=0> y[n];  // responses
  }  transformed data{
    matrix [n, p] Q_ast;
    matrix [p, p] R_ast;
    matrix [p, p] R_ast_inverse;
    Q_ast = qr_Q(x)[,1:p] * sqrt(n - 1.0);
    R_ast = qr_R(x)[1:p,] / sqrt(n - 1.0);
    R_ast_inverse = inverse(R_ast);
  }  parameters{
    vector[p] theta;  // regression coefficients for predictors
  }  transformed parameters{
    vector[p] beta;
    vector[n] mu;
    mu = exp(Q_ast*theta);
    beta = R_ast_inverse * theta;
  }  model{
     y  ~ poisson(mu);
  }
```

Я понимаю, что в этом примере была использована репараметризация QR (см. раздел 9.2 справочника Stan руководство). Однако, поскольку я новичок в Стэне, я нахожу это довольно запугивающим, и я не уверен, как реализовать экспоненциальный GLM таким же образом. Нужно ли даже использовать репараметризацию QR для экспоненциального случая?

Во всяком случае, моя наивная попытка экспоненциального GLM - это следующая адаптация кода Poisson GLM:

```{stan output.var="ExponentialGLM"}
  data{
    int<lower=1> n;     // number of observations
    int<lower=1> p;     // number of predictors
    matrix[n, p] x;     // design matrix
    real<lower=0> y[n];  // responses
  }  transformed data{
    matrix [n, p] Q_ast;
    matrix [p, p] R_ast;
    matrix [p, p] R_ast_inverse;
    Q_ast = qr_Q(x)[,1:p] * sqrt(n - 1.0);
    R_ast = qr_R(x)[1:p,] / sqrt(n - 1.0);
    R_ast_inverse = inverse(R_ast);
  }  parameters{
    vector[p] theta;  // regression coefficients for predictors
  }  transformed parameters{
    vector[p] beta;
    vector[n] mu;
    mu = exp(Q_ast*theta);
    beta = (R_ast_inverse * theta);
  }  model{
     y  ~ exponential(mu);
  }
```

Но, как я уже сказал, я совершенно не уверен, что это даже правильный способ кодирования такой модели в Стэне. Все, что я сделал, это просто попытаться адаптировать пример poisson GLM к вышеупомянутому экспоненциальному GLM.

Я был бы признателен, если бы кто-то мог продемонстрировать экспоненциальный GLM и, если это не слишком много проблем , пожалуйста, найдите время, чтобы разъяснить мои недоразумения выше.

Примеры данных

    conc  time
 1   6.1   0.8
 2   4.2   3.5
 3   0.5  12.4
 4   8.8   1.1
 5   1.5   8.9
 6   9.2   2.4
 7   8.5   0.1
 8   8.7   0.4
 9   6.7   3.5
10   6.5   8.3
11   6.3   2.6
12   6.7   1.5
13   0.2  16.6
14   8.7   0.1
15   7.5   1.3

Maurits Evers ответил: 12 мая 2018 в 09:49

Можете ли вы предоставить некоторые примеры данных? Меня тоже смущает общая модель. Вы говорите о модели Пуассона или экспоненциальной (т. Е. Гамма) модели? Это две разные вещи. Эквивалентами glm являются glm(..., family = poisson) и glm(..., family = Gamma). На первом этапе вы должны сравнить вывод вашей модели Stan с параметром glm. Согласны ли оценки баллов?

The Pointer ответил: 12 мая 2018 в 10:09

@MauritsEvers Спасибо за ответ. Я говорю об экспоненциальной модели. Модель Пуассона, которую я опубликовал, была совершенно отдельным примером, который я пытался использовать, чтобы определить, как сделать экспоненциальную модель. Я редактировал свой пост с некоторыми примерами данных. Если что-то еще неясно, пожалуйста, не стесняйтесь спрашивать. Я совершенно не уверен, как реализовать эту модель, поэтому помощь очень ценится. Спасибо за помощь.

1 ответ

Есть решение

Maurits Evers ответил: 12 мая 2018 в 12:10

Как насчет чего-то подобного?

Начнем с чтения данных образца.

df <- read.table(text =
    "    conc  time
 1   6.1   0.8
 2   4.2   3.5
 3   0.5  12.4
 4   8.8   1.1
 5   1.5   8.9
 6   9.2   2.4
 7   8.5   0.1
 8   8.7   0.4
 9   6.7   3.5
10   6.5   8.3
11   6.3   2.6
12   6.7   1.5
13   0.2  16.6
14   8.7   0.1
15   7.5   1.3", header = T);

Мы определяем модель Стэна как простую модель Гамма (экспоненциальная) с лог-линией.

model <- "
data {
    int N;                                       // Number of observations
    int K;                                       // Number of parameters
    real y[N];                                   // Response
    matrix[N,K] X;                               // Model matrix
}parameters {
    vector[K] beta;                              // Model parameters
}transformed parameters {
    vector[N] lambda;
    lambda = exp(-X * beta);                     // log-link
}model {
    // Priors
    beta[1] ~ cauchy(0, 10);
    for (i in 2:K)
        beta[i] ~ cauchy(0, 2.5);    y ~ exponential(lambda);
}
"

Здесь я предполагаю (стандартное) слабо информативное поведение Коши на бета-модели параметры.

Теперь мы сопоставим модель с данными.

library(rstan);
options(mc.cores = parallel::detectCores());
rstan_options(auto_write = TRUE);X <- model.matrix(time ~ conc, data = df);
model.stan <- stan(
    model_code = model,
    data = list(
        N = nrow(df),
        K = ncol(X),
        y = df$time,
        X = X),
    iter = 4000);

Чтобы сравнить оценки точек , мы также сопоставляем Gamma GLM с данными с помощью glm.
```
model.glm <- glm(time ~ conc, data = df, family = Gamma(link = "log"));
```

Оценки параметров модели Stan

summary(model.stan, "beta")$summary;
#              mean     se_mean         sd       2.5%        25%        50%
#beta[1]  2.9371457 0.016460000 0.62017934  1.8671652  2.5000356  2.8871936
#beta[2] -0.3099799 0.002420744 0.09252423 -0.5045937 -0.3708111 -0.3046333
#               75%      97.5%    n_eff     Rhat
#beta[1]  3.3193334  4.3078478 1419.629 1.002440
#beta[2] -0.2461567 -0.1412095 1460.876 1.002236

Оценки параметров модели GLM:

summary(model.glm)$coef;
#              Estimate Std. Error   t value     Pr(>|t|)
#(Intercept)  2.8211487 0.54432115  5.182875 0.0001762685
#conc        -0.3013355 0.08137986 -3.702827 0.0026556952

Мы видим хорошее совпадение оценок Стэна для средств и оценок параметров Gamma-GLM.

Мы оцениваем оценки параметров как из моделей Stan, так и glm.

library(tidyverse);
as.data.frame(summary(model.stan, "beta")$summary) %>%
    rownames_to_column("Parameter") %>%
    ggplot(aes(x = `50%`, y = Parameter)) +
    geom_point(size = 3) +
    geom_segment(aes(x = `2.5%`, xend = `97.5%`, yend = Parameter), size = 1) +
    geom_segment(aes(x = `25%`, xend = `75%`, yend = Parameter), size = 2) +
    labs(x = "Median (plus/minus 95% and 50% CIs)") +
    geom_point(
        data = data.frame(summary(model.glm)$coef, row.names = c("beta[1]", "beta[2]")) %>%
            rownames_to_column("Parameter"),
        aes(x = Estimate, y = Parameter),
        colour = "red", size = 4, shape = 1)

glm показаны красным цветом.

Обновить (модель Fit Stan с использованием повторной параметризации QR)

Модель Стэна с повторной параметризацией QR.

model.QR <- "
data {
    int N;                                       // Number of observations
    int K;                                       // Number of parameters
    real y[N];                                   // Response
    matrix[N,K] X;                               // Model matrix
}transformed data {
    matrix[N, K] Q;
    matrix[K, K] R;
    matrix[K, K] R_inverse;
    Q = qr_Q(X)[, 1:K];
    R = qr_R(X)[1:K, ];
    R_inverse = inverse(R);
}parameters {
    vector[K] theta;
}transformed parameters {
    vector[N] lambda;
    lambda = exp(-Q * theta);                     // log-link
}model {
    y ~ exponential(lambda);
}generated quantities {
    vector[K] beta;                              // Model parameters
    beta = R_inverse * theta;
}
"

В QR-декомпозиции мы имеем lambda = exp(- X * beta) = exp(- Q * R * beta) = exp(- Q * theta), где theta = R * beta и, следовательно, beta = R^-1 * theta.

Обратите внимание, что мы не указываем никаких приуров на тете; в Stan это значение по умолчанию равномерным (то есть равномерным) приоритетам.

Примените модель модели к данным.

library(rstan);
options(mc.cores = parallel::detectCores());
rstan_options(auto_write = TRUE);X <- model.matrix(time ~ conc, data = df);
model.stan.QR <- stan(
    model_code = model.QR,
    data = list(
        N = nrow(df),
        K = ncol(X),
        y = df$time,
        X = X),
    iter = 4000);

Оценки параметров
```
summary(model.stan.QR, "beta")$summary;
#              mean     se_mean         sd       2.5%        25%        50%
#beta[1]  2.9637547 0.009129250 0.64383609  1.8396681  2.5174800  2.9194682
#beta[2] -0.3134252 0.001321584 0.09477156 -0.5126905 -0.3740475 -0.3093937
#               75%      97.5%    n_eff      Rhat
#beta[1]  3.3498585  4.3593912 4973.710 0.9998545
#beta[2] -0.2478029 -0.1395686 5142.408 1.0003236
```
Вы можете видеть, что оценки параметров между двумя моделями Stan (с и без повторной параметризации QR) очень хорошо согласуются.

Если бы вы спросили меня, изменит ли QR re-parametrisation (большую / любую) разницу, я бы сказал, "возможно, не в этом случае"; Эндрю Гельман и другие часто подчеркивали, что использование даже очень слабых информативных пригородов будет способствовать сближению и должно быть предпочтительным по сравнению с плоскими (равномерными) призваниями; Я всегда старался использовать слабоинформативные приставки по всем параметрам и начинать с модели без повторной параметризации QR; если сходимость плохой, я бы попытался оптимизировать свою модель на следующем шаге.

The Pointer ответил: 12 мая 2018 в 10:55

Ах, это освещает! Еще раз спасибо за помощь. У меня есть пара вопросов, я не против. Таким образом, нам не нужно использовать репараметризацию QR для Exponential GLM, как это было сделано в примере Poisson GLM, который я опубликовал? Кроме того, для журнала-ссылки у вас есть lambda = exp(-X * beta). Если я не ошибаюсь, эта бета соответствует бета_1 в описании модели. Что относительно бета_0?

Maurits Evers ответил: 12 мая 2018 в 11:09

@ThePointer (1) beta[1] соответствует смещению (ваш бета_0); Стэн использует индексирование на основе 1, а не на основе 0 (с хорошим побочным эффектом, более совместимым с R). Я добавил график, который показывает сравнение оценок параметров как из моделей Stan, так и glm. (2) Что касается репараметризации QR. Я не эксперт в этом, но из того, что я понимаю, репараметризация QR может помочь ускорить процесс выборки. Это стратегия оптимизации, которая помогает, в частности, если у вас есть плоские приоритеты. [...]

Maurits Evers ответил: 12 мая 2018 в 11:13

[...] Я никогда не использовал QR-репараметризацию, потому что большую часть времени я могу обеспечить слабо информативные приоритеты, и конвергенция решений прекрасна (вы можете проверить, посмотрев на значения Rhat). Поэтому мне жаль, что я не могу дать более полный ответ по этому вопросу. Возможно, кто-то из команды Rstan / Stan может дать представление; Я знаю, что они контролируют SO для вопросов, связанных с Стэном; или, возможно, задать отдельный вопрос на форумах Stan Google.

Maurits Evers ответил: 13 мая 2018 в 06:07

Вы очень приветствуете @ ThePointer; Мне нравятся эти вопросы моделирования Stan / Rstan, так как (R) stan - фантастическая рамочная ИМО. Вы правильно интерпретируете 95% достоверный интервал. Что касается, priors, да, в модели повторного параметрирования QR вы могли бы снабжать слабоинформативные приуровни (например, приставки Коши) по параметрам тета. Но в этом случае вы можете использовать первую модель. Насколько я понимаю, QR-повторная параметризация - это стратегия оптимизации, которая подходит, когда решения не сходятся (что может случиться, когда у вас есть плоские приоритеты и / или K). [...]

Maurits Evers ответил: 13 мая 2018 в 06:08

[...] В вашем случае конвергенция на самом деле не проблема, даже с плоскими первичными в первой модели. Таким образом, дополнительное преимущество повторной параметризации QR в этом случае кажется небольшим. Но все же хорошо знать, как его реализовать; Я никогда не пользовался повторной параметризацией QR, поэтому я узнал что-то новое из вашего примера ;-)