Внедрение экспоненциальной общей линейной модели в 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


1 ответ

Есть решение
Maurits Evers ответил: 12 мая 2018 в 12:10

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

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

    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);
    
  2. Мы определяем модель Стэна как простую модель Гамма (экспоненциальная) с лог-линией.

    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);
    }
    "
    

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

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

    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);
    
  4. Чтобы сравнить оценки точек , мы также сопоставляем Gamma GLM с данными с помощью glm.

    model.glm <- glm(time ~ conc, data = df, family = Gamma(link = "log"));
    
  5. Оценки параметров модели 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.

  6. Мы оцениваем оценки параметров как из моделей 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)

  1. Модель Стэна с повторной параметризацией 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 это значение по умолчанию равномерным (то есть равномерным) приоритетам.

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

    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);
    
  3. Оценки параметров

    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, поэтому я узнал что-то новое из вашего примера ;-)