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

Как насчет этого?

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

   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. Чтобы сравнить точечные оценки, мы также подобрали гамма-модель GLM к данным, используя glm.

   model.glm <- glm(time ~ conc, data = df, family = Gamma(link = "log"));
   

5. Оценки параметров модели Стэна:

   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
   

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

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

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

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