komparasi linear dan bayes dengan broom
Copy Right @ Prof Dr Suhartono M.Kom @ Magister Informatika UIN Maulana Malik Ibrahim@ (http://suhartono.lecturer.uin-malang.ac.id/riwayat-hidup/)
Pengantar
vignette menjelaskan cara menggunakan paket tidybayes untuk mengekstraksi tidy data frame dari model brms, dan penarikan dari distribusi posterior variabel model, kecocokan, dan prediksi dari brms :: brm. Untuk pengenalan yang lebih umum tentang tidybayes, lihat vignette(“tidybayes”).
Untuk menjalankan vignette membutuhkan package sebagai berikut :
library(magrittr)
library(dplyr)##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(forcats)
library(modelr)
library(ggdist)
library(tidybayes)
library(ggplot2)
library(cowplot)
library(emmeans)
library(broom)##
## Attaching package: 'broom'## The following object is masked from 'package:modelr':
##
## bootstraplibrary(rstan)## Loading required package: StanHeaders## rstan (Version 2.21.2, GitRev: 2e1f913d3ca3)## For execution on a local, multicore CPU with excess RAM we recommend calling
## options(mc.cores = parallel::detectCores()).
## To avoid recompilation of unchanged Stan programs, we recommend calling
## rstan_options(auto_write = TRUE)##
## Attaching package: 'rstan'## The following object is masked from 'package:magrittr':
##
## extractlibrary(rstanarm)## Loading required package: Rcpp## This is rstanarm version 2.21.1## - See https://mc-stan.org/rstanarm/articles/priors for changes to default priors!## - Default priors may change, so it's safest to specify priors, even if equivalent to the defaults.## - For execution on a local, multicore CPU with excess RAM we recommend calling## options(mc.cores = parallel::detectCores())##
## Attaching package: 'rstanarm'## The following object is masked from 'package:rstan':
##
## loolibrary(brms)## Loading 'brms' package (version 2.14.4). Useful instructions
## can be found by typing help('brms'). A more detailed introduction
## to the package is available through vignette('brms_overview').##
## Attaching package: 'brms'## The following objects are masked from 'package:rstanarm':
##
## dirichlet, exponential, get_y, lasso, ngrps## The following object is masked from 'package:rstan':
##
## loo## The following objects are masked from 'package:tidybayes':
##
## dstudent_t, pstudent_t, qstudent_t, rstudent_t## The following objects are masked from 'package:ggdist':
##
## dstudent_t, pstudent_t, qstudent_t, rstudent_t## The following object is masked from 'package:stats':
##
## arlibrary(bayesplot)## This is bayesplot version 1.8.0## - Online documentation and vignettes at mc-stan.org/bayesplot## - bayesplot theme set to bayesplot::theme_default()## * Does _not_ affect other ggplot2 plots## * See ?bayesplot_theme_set for details on theme settinglibrary(MCMCglmm)## Loading required package: Matrix## Loading required package: coda##
## Attaching package: 'coda'## The following object is masked from 'package:rstan':
##
## traceplot## Loading required package: ape##
## Attaching package: 'MCMCglmm'## The following object is masked from 'package:brms':
##
## melibrary(RColorBrewer)
theme_set(theme_tidybayes() + panel_border())Opsi ini membantu Stan bekerja lebih cepat:
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())Kita mengambil dataset :
set.seed(5)
n = 10
n_condition = 5
ABC =
tibble(
condition = rep(c("A","B","C","D","E"), n),
response = rnorm(n * 5, c(0,1,2,1,-1), 0.5)
)head(ABC, 10)## # A tibble: 10 x 2
## condition response
## <chr> <dbl>
## 1 A -0.420
## 2 B 1.69
## 3 C 1.37
## 4 D 1.04
## 5 E -0.144
## 6 A -0.301
## 7 B 0.764
## 8 C 1.68
## 9 D 0.857
## 10 E -0.931ABC %>%
ggplot(aes(y = condition, x = response)) +
geom_point()
Untukmengkompile rstan model digunakan sebagai berikut :
stancode <- 'data {
int<lower=1> n;
int<lower=1> n_condition;
int<lower=1, upper=n_condition> condition[n];
real response[n];
}
parameters {
real overall_mean;
vector[n_condition] condition_zoffset;
real<lower=0> response_sd;
real<lower=0> condition_mean_sd;
}
transformed parameters {
vector[n_condition] condition_mean;
condition_mean = overall_mean + condition_zoffset * condition_mean_sd;
}
model {
response_sd ~ cauchy(0, 1); // => half-cauchy(0, 1)
condition_mean_sd ~ cauchy(0, 1); // => half-cauchy(0, 1)
overall_mean ~ normal(0, 5);
condition_zoffset ~ normal(0, 1); // => condition_mean ~ normal(overall_mean, condition_mean_sd)
for (i in 1:n) {
response[i] ~ normal(condition_mean[condition[i]], response_sd);
}
}'
ABC_stan <- stan_model(model_code = stancode, verbose = TRUE)##
## TRANSLATING MODEL '150269bbfe815db848542eb814742386' FROM Stan CODE TO C++ CODE NOW.
## successful in parsing the Stan model '150269bbfe815db848542eb814742386'.
## OS: x86_64, linux-gnu; rstan: 2.21.2; Rcpp: 1.0.6; inline: 0.3.17
## >> setting environment variables:
## PKG_LIBS = '/home/suhartono/R/x86_64-pc-linux-gnu-library/3.6/rstan/lib//libStanServices.a' -L'/home/suhartono/R/x86_64-pc-linux-gnu-library/3.6/StanHeaders/lib/' -lStanHeaders -L'/home/suhartono/R/x86_64-pc-linux-gnu-library/3.6/RcppParallel/lib/' -ltbb
## PKG_CPPFLAGS = -I"/home/suhartono/R/x86_64-pc-linux-gnu-library/3.6/Rcpp/include/" -I"/home/suhartono/R/x86_64-pc-linux-gnu-library/3.6/RcppEigen/include/" -I"/home/suhartono/R/x86_64-pc-linux-gnu-library/3.6/RcppEigen/include/unsupported" -I"/home/suhartono/R/x86_64-pc-linux-gnu-library/3.6/BH/include" -I"/home/suhartono/R/x86_64-pc-linux-gnu-library/3.6/StanHeaders/include/src/" -I"/home/suhartono/R/x86_64-pc-linux-gnu-library/3.6/StanHeaders/include/" -I"/home/suhartono/R/x86_64-pc-linux-gnu-library/3.6/RcppParallel/include/" -I"/home/suhartono/R/x86_64-pc-linux-gnu-library/3.6/rstan/include" -DEIGEN_NO_DEBUG -DBOOST_DISABLE_ASSERTS -DBOOST_PENDING_INTEGER_LOG2_HPP -DSTAN_THREADS -DBOOST_NO_AUTO_PTR -include '/home/suhartono/R/x86_64-pc-linux-gnu-library/3.6/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp' -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1
## >> Program source :
##
## 1 :
## 2 : // includes from the plugin
## 3 : // [[Rcpp::plugins(cpp14)]]
## 4 :
## 5 :
## 6 : // user includes
## 7 : #include <Rcpp.h>
## 8 : #include <rstan/io/rlist_ref_var_context.hpp>
## 9 : #include <rstan/io/r_ostream.hpp>
## 10 : #include <rstan/stan_args.hpp>
## 11 : #include <boost/integer/integer_log2.hpp>
## 12 : // Code generated by Stan version 2.21.0
## 13 :
## 14 : #include <stan/model/model_header.hpp>
## 15 :
## 16 : namespace model10837057f64f_150269bbfe815db848542eb814742386_namespace {
## 17 :
## 18 : using std::istream;
## 19 : using std::string;
## 20 : using std::stringstream;
## 21 : using std::vector;
## 22 : using stan::io::dump;
## 23 : using stan::math::lgamma;
## 24 : using stan::model::prob_grad;
## 25 : using namespace stan::math;
## 26 :
## 27 : static int current_statement_begin__;
## 28 :
## 29 : stan::io::program_reader prog_reader__() {
## 30 : stan::io::program_reader reader;
## 31 : reader.add_event(0, 0, "start", "model10837057f64f_150269bbfe815db848542eb814742386");
## 32 : reader.add_event(27, 25, "end", "model10837057f64f_150269bbfe815db848542eb814742386");
## 33 : return reader;
## 34 : }
## 35 :
## 36 : class model10837057f64f_150269bbfe815db848542eb814742386
## 37 : : public stan::model::model_base_crtp<model10837057f64f_150269bbfe815db848542eb814742386> {
## 38 : private:
## 39 : int n;
## 40 : int n_condition;
## 41 : std::vector<int> condition;
## 42 : std::vector<double> response;
## 43 : public:
## 44 : model10837057f64f_150269bbfe815db848542eb814742386(rstan::io::rlist_ref_var_context& context__,
## 45 : std::ostream* pstream__ = 0)
## 46 : : model_base_crtp(0) {
## 47 : ctor_body(context__, 0, pstream__);
## 48 : }
## 49 :
## 50 : model10837057f64f_150269bbfe815db848542eb814742386(stan::io::var_context& context__,
## 51 : unsigned int random_seed__,
## 52 : std::ostream* pstream__ = 0)
## 53 : : model_base_crtp(0) {
## 54 : ctor_body(context__, random_seed__, pstream__);
## 55 : }
## 56 :
## 57 : void ctor_body(stan::io::var_context& context__,
## 58 : unsigned int random_seed__,
## 59 : std::ostream* pstream__) {
## 60 : typedef double local_scalar_t__;
## 61 :
## 62 : boost::ecuyer1988 base_rng__ =
## 63 : stan::services::util::create_rng(random_seed__, 0);
## 64 : (void) base_rng__; // suppress unused var warning
## 65 :
## 66 : current_statement_begin__ = -1;
## 67 :
## 68 : static const char* function__ = "model10837057f64f_150269bbfe815db848542eb814742386_namespace::model10837057f64f_150269bbfe815db848542eb814742386";
## 69 : (void) function__; // dummy to suppress unused var warning
## 70 : size_t pos__;
## 71 : (void) pos__; // dummy to suppress unused var warning
## 72 : std::vector<int> vals_i__;
## 73 : std::vector<double> vals_r__;
## 74 : local_scalar_t__ DUMMY_VAR__(std::numeric_limits<double>::quiet_NaN());
## 75 : (void) DUMMY_VAR__; // suppress unused var warning
## 76 :
## 77 : try {
## 78 : // initialize data block variables from context__
## 79 : current_statement_begin__ = 2;
## 80 : context__.validate_dims("data initialization", "n", "int", context__.to_vec());
## 81 : n = int(0);
## 82 : vals_i__ = context__.vals_i("n");
## 83 : pos__ = 0;
## 84 : n = vals_i__[pos__++];
## 85 : check_greater_or_equal(function__, "n", n, 1);
## 86 :
## 87 : current_statement_begin__ = 3;
## 88 : context__.validate_dims("data initialization", "n_condition", "int", context__.to_vec());
## 89 : n_condition = int(0);
## 90 : vals_i__ = context__.vals_i("n_condition");
## 91 : pos__ = 0;
## 92 : n_condition = vals_i__[pos__++];
## 93 : check_greater_or_equal(function__, "n_condition", n_condition, 1);
## 94 :
## 95 : current_statement_begin__ = 4;
## 96 : validate_non_negative_index("condition", "n", n);
## 97 : context__.validate_dims("data initialization", "condition", "int", context__.to_vec(n));
## 98 : condition = std::vector<int>(n, int(0));
## 99 : vals_i__ = context__.vals_i("condition");
## 100 : pos__ = 0;
## 101 : size_t condition_k_0_max__ = n;
## 102 : for (size_t k_0__ = 0; k_0__ < condition_k_0_max__; ++k_0__) {
## 103 : condition[k_0__] = vals_i__[pos__++];
## 104 : }
## 105 : size_t condition_i_0_max__ = n;
## 106 : for (size_t i_0__ = 0; i_0__ < condition_i_0_max__; ++i_0__) {
## 107 : check_greater_or_equal(function__, "condition[i_0__]", condition[i_0__], 1);
## 108 : check_less_or_equal(function__, "condition[i_0__]", condition[i_0__], n_condition);
## 109 : }
## 110 :
## 111 : current_statement_begin__ = 5;
## 112 : validate_non_negative_index("response", "n", n);
## 113 : context__.validate_dims("data initialization", "response", "double", context__.to_vec(n));
## 114 : response = std::vector<double>(n, double(0));
## 115 : vals_r__ = context__.vals_r("response");
## 116 : pos__ = 0;
## 117 : size_t response_k_0_max__ = n;
## 118 : for (size_t k_0__ = 0; k_0__ < response_k_0_max__; ++k_0__) {
## 119 : response[k_0__] = vals_r__[pos__++];
## 120 : }
## 121 :
## 122 :
## 123 : // initialize transformed data variables
## 124 : // execute transformed data statements
## 125 :
## 126 : // validate transformed data
## 127 :
## 128 : // validate, set parameter ranges
## 129 : num_params_r__ = 0U;
## 130 : param_ranges_i__.clear();
## 131 : current_statement_begin__ = 8;
## 132 : num_params_r__ += 1;
## 133 : current_statement_begin__ = 9;
## 134 : validate_non_negative_index("condition_zoffset", "n_condition", n_condition);
## 135 : num_params_r__ += n_condition;
## 136 : current_statement_begin__ = 10;
## 137 : num_params_r__ += 1;
## 138 : current_statement_begin__ = 11;
## 139 : num_params_r__ += 1;
## 140 : } catch (const std::exception& e) {
## 141 : stan::lang::rethrow_located(e, current_statement_begin__, prog_reader__());
## 142 : // Next line prevents compiler griping about no return
## 143 : throw std::runtime_error("*** IF YOU SEE THIS, PLEASE REPORT A BUG ***");
## 144 : }
## 145 : }
## 146 :
## 147 : ~model10837057f64f_150269bbfe815db848542eb814742386() { }
## 148 :
## 149 :
## 150 : void transform_inits(const stan::io::var_context& context__,
## 151 : std::vector<int>& params_i__,
## 152 : std::vector<double>& params_r__,
## 153 : std::ostream* pstream__) const {
## 154 : typedef double local_scalar_t__;
## 155 : stan::io::writer<double> writer__(params_r__, params_i__);
## 156 : size_t pos__;
## 157 : (void) pos__; // dummy call to supress warning
## 158 : std::vector<double> vals_r__;
## 159 : std::vector<int> vals_i__;
## 160 :
## 161 : current_statement_begin__ = 8;
## 162 : if (!(context__.contains_r("overall_mean")))
## 163 : stan::lang::rethrow_located(std::runtime_error(std::string("Variable overall_mean missing")), current_statement_begin__, prog_reader__());
## 164 : vals_r__ = context__.vals_r("overall_mean");
## 165 : pos__ = 0U;
## 166 : context__.validate_dims("parameter initialization", "overall_mean", "double", context__.to_vec());
## 167 : double overall_mean(0);
## 168 : overall_mean = vals_r__[pos__++];
## 169 : try {
## 170 : writer__.scalar_unconstrain(overall_mean);
## 171 : } catch (const std::exception& e) {
## 172 : stan::lang::rethrow_located(std::runtime_error(std::string("Error transforming variable overall_mean: ") + e.what()), current_statement_begin__, prog_reader__());
## 173 : }
## 174 :
## 175 : current_statement_begin__ = 9;
## 176 : if (!(context__.contains_r("condition_zoffset")))
## 177 : stan::lang::rethrow_located(std::runtime_error(std::string("Variable condition_zoffset missing")), current_statement_begin__, prog_reader__());
## 178 : vals_r__ = context__.vals_r("condition_zoffset");
## 179 : pos__ = 0U;
## 180 : validate_non_negative_index("condition_zoffset", "n_condition", n_condition);
## 181 : context__.validate_dims("parameter initialization", "condition_zoffset", "vector_d", context__.to_vec(n_condition));
## 182 : Eigen::Matrix<double, Eigen::Dynamic, 1> condition_zoffset(n_condition);
## 183 : size_t condition_zoffset_j_1_max__ = n_condition;
## 184 : for (size_t j_1__ = 0; j_1__ < condition_zoffset_j_1_max__; ++j_1__) {
## 185 : condition_zoffset(j_1__) = vals_r__[pos__++];
## 186 : }
## 187 : try {
## 188 : writer__.vector_unconstrain(condition_zoffset);
## 189 : } catch (const std::exception& e) {
## 190 : stan::lang::rethrow_located(std::runtime_error(std::string("Error transforming variable condition_zoffset: ") + e.what()), current_statement_begin__, prog_reader__());
## 191 : }
## 192 :
## 193 : current_statement_begin__ = 10;
## 194 : if (!(context__.contains_r("response_sd")))
## 195 : stan::lang::rethrow_located(std::runtime_error(std::string("Variable response_sd missing")), current_statement_begin__, prog_reader__());
## 196 : vals_r__ = context__.vals_r("response_sd");
## 197 : pos__ = 0U;
## 198 : context__.validate_dims("parameter initialization", "response_sd", "double", context__.to_vec());
## 199 : double response_sd(0);
## 200 : response_sd = vals_r__[pos__++];
## 201 : try {
## 202 : writer__.scalar_lb_unconstrain(0, response_sd);
## 203 : } catch (const std::exception& e) {
## 204 : stan::lang::rethrow_located(std::runtime_error(std::string("Error transforming variable response_sd: ") + e.what()), current_statement_begin__, prog_reader__());
## 205 : }
## 206 :
## 207 : current_statement_begin__ = 11;
## 208 : if (!(context__.contains_r("condition_mean_sd")))
## 209 : stan::lang::rethrow_located(std::runtime_error(std::string("Variable condition_mean_sd missing")), current_statement_begin__, prog_reader__());
## 210 : vals_r__ = context__.vals_r("condition_mean_sd");
## 211 : pos__ = 0U;
## 212 : context__.validate_dims("parameter initialization", "condition_mean_sd", "double", context__.to_vec());
## 213 : double condition_mean_sd(0);
## 214 : condition_mean_sd = vals_r__[pos__++];
## 215 : try {
## 216 : writer__.scalar_lb_unconstrain(0, condition_mean_sd);
## 217 : } catch (const std::exception& e) {
## 218 : stan::lang::rethrow_located(std::runtime_error(std::string("Error transforming variable condition_mean_sd: ") + e.what()), current_statement_begin__, prog_reader__());
## 219 : }
## 220 :
## 221 : params_r__ = writer__.data_r();
## 222 : params_i__ = writer__.data_i();
## 223 : }
## 224 :
## 225 : void transform_inits(const stan::io::var_context& context,
## 226 : Eigen::Matrix<double, Eigen::Dynamic, 1>& params_r,
## 227 : std::ostream* pstream__) const {
## 228 : std::vector<double> params_r_vec;
## 229 : std::vector<int> params_i_vec;
## 230 : transform_inits(context, params_i_vec, params_r_vec, pstream__);
## 231 : params_r.resize(params_r_vec.size());
## 232 : for (int i = 0; i < params_r.size(); ++i)
## 233 : params_r(i) = params_r_vec[i];
## 234 : }
## 235 :
## 236 :
## 237 : template <bool propto__, bool jacobian__, typename T__>
## 238 : T__ log_prob(std::vector<T__>& params_r__,
## 239 : std::vector<int>& params_i__,
## 240 : std::ostream* pstream__ = 0) const {
## 241 :
## 242 : typedef T__ local_scalar_t__;
## 243 :
## 244 : local_scalar_t__ DUMMY_VAR__(std::numeric_limits<double>::quiet_NaN());
## 245 : (void) DUMMY_VAR__; // dummy to suppress unused var warning
## 246 :
## 247 : T__ lp__(0.0);
## 248 : stan::math::accumulator<T__> lp_accum__;
## 249 : try {
## 250 : stan::io::reader<local_scalar_t__> in__(params_r__, params_i__);
## 251 :
## 252 : // model parameters
## 253 : current_statement_begin__ = 8;
## 254 : local_scalar_t__ overall_mean;
## 255 : (void) overall_mean; // dummy to suppress unused var warning
## 256 : if (jacobian__)
## 257 : overall_mean = in__.scalar_constrain(lp__);
## 258 : else
## 259 : overall_mean = in__.scalar_constrain();
## 260 :
## 261 : current_statement_begin__ = 9;
## 262 : Eigen::Matrix<local_scalar_t__, Eigen::Dynamic, 1> condition_zoffset;
## 263 : (void) condition_zoffset; // dummy to suppress unused var warning
## 264 : if (jacobian__)
## 265 : condition_zoffset = in__.vector_constrain(n_condition, lp__);
## 266 : else
## 267 : condition_zoffset = in__.vector_constrain(n_condition);
## 268 :
## 269 : current_statement_begin__ = 10;
## 270 : local_scalar_t__ response_sd;
## 271 : (void) response_sd; // dummy to suppress unused var warning
## 272 : if (jacobian__)
## 273 : response_sd = in__.scalar_lb_constrain(0, lp__);
## 274 : else
## 275 : response_sd = in__.scalar_lb_constrain(0);
## 276 :
## 277 : current_statement_begin__ = 11;
## 278 : local_scalar_t__ condition_mean_sd;
## 279 : (void) condition_mean_sd; // dummy to suppress unused var warning
## 280 : if (jacobian__)
## 281 : condition_mean_sd = in__.scalar_lb_constrain(0, lp__);
## 282 : else
## 283 : condition_mean_sd = in__.scalar_lb_constrain(0);
## 284 :
## 285 : // transformed parameters
## 286 : current_statement_begin__ = 14;
## 287 : validate_non_negative_index("condition_mean", "n_condition", n_condition);
## 288 : Eigen::Matrix<local_scalar_t__, Eigen::Dynamic, 1> condition_mean(n_condition);
## 289 : stan::math::initialize(condition_mean, DUMMY_VAR__);
## 290 : stan::math::fill(condition_mean, DUMMY_VAR__);
## 291 :
## 292 : // transformed parameters block statements
## 293 : current_statement_begin__ = 15;
## 294 : stan::math::assign(condition_mean, add(overall_mean, multiply(condition_zoffset, condition_mean_sd)));
## 295 :
## 296 : // validate transformed parameters
## 297 : const char* function__ = "validate transformed params";
## 298 : (void) function__; // dummy to suppress unused var warning
## 299 :
## 300 : current_statement_begin__ = 14;
## 301 : size_t condition_mean_j_1_max__ = n_condition;
## 302 : for (size_t j_1__ = 0; j_1__ < condition_mean_j_1_max__; ++j_1__) {
## 303 : if (stan::math::is_uninitialized(condition_mean(j_1__))) {
## 304 : std::stringstream msg__;
## 305 : msg__ << "Undefined transformed parameter: condition_mean" << "(" << j_1__ << ")";
## 306 : stan::lang::rethrow_located(std::runtime_error(std::string("Error initializing variable condition_mean: ") + msg__.str()), current_statement_begin__, prog_reader__());
## 307 : }
## 308 : }
## 309 :
## 310 : // model body
## 311 :
## 312 : current_statement_begin__ = 18;
## 313 : lp_accum__.add(cauchy_log<propto__>(response_sd, 0, 1));
## 314 : current_statement_begin__ = 19;
## 315 : lp_accum__.add(cauchy_log<propto__>(condition_mean_sd, 0, 1));
## 316 : current_statement_begin__ = 20;
## 317 : lp_accum__.add(normal_log<propto__>(overall_mean, 0, 5));
## 318 : current_statement_begin__ = 21;
## 319 : lp_accum__.add(normal_log<propto__>(condition_zoffset, 0, 1));
## 320 : current_statement_begin__ = 22;
## 321 : for (int i = 1; i <= n; ++i) {
## 322 :
## 323 : current_statement_begin__ = 23;
## 324 : lp_accum__.add(normal_log<propto__>(get_base1(response, i, "response", 1), get_base1(condition_mean, get_base1(condition, i, "condition", 1), "condition_mean", 1), response_sd));
## 325 : }
## 326 :
## 327 : } catch (const std::exception& e) {
## 328 : stan::lang::rethrow_located(e, current_statement_begin__, prog_reader__());
## 329 : // Next line prevents compiler griping about no return
## 330 : throw std::runtime_error("*** IF YOU SEE THIS, PLEASE REPORT A BUG ***");
## 331 : }
## 332 :
## 333 : lp_accum__.add(lp__);
## 334 : return lp_accum__.sum();
## 335 :
## 336 : } // log_prob()
## 337 :
## 338 : template <bool propto, bool jacobian, typename T_>
## 339 : T_ log_prob(Eigen::Matrix<T_,Eigen::Dynamic,1>& params_r,
## 340 : std::ostream* pstream = 0) const {
## 341 : std::vector<T_> vec_params_r;
## 342 : vec_params_r.reserve(params_r.size());
## 343 : for (int i = 0; i < params_r.size(); ++i)
## 344 : vec_params_r.push_back(params_r(i));
## 345 : std::vector<int> vec_params_i;
## 346 : return log_prob<propto,jacobian,T_>(vec_params_r, vec_params_i, pstream);
## 347 : }
## 348 :
## 349 :
## 350 : void get_param_names(std::vector<std::string>& names__) const {
## 351 : names__.resize(0);
## 352 : names__.push_back("overall_mean");
## 353 : names__.push_back("condition_zoffset");
## 354 : names__.push_back("response_sd");
## 355 : names__.push_back("condition_mean_sd");
## 356 : names__.push_back("condition_mean");
## 357 : }
## 358 :
## 359 :
## 360 : void get_dims(std::vector<std::vector<size_t> >& dimss__) const {
## 361 : dimss__.resize(0);
## 362 : std::vector<size_t> dims__;
## 363 : dims__.resize(0);
## 364 : dimss__.push_back(dims__);
## 365 : dims__.resize(0);
## 366 : dims__.push_back(n_condition);
## 367 : dimss__.push_back(dims__);
## 368 : dims__.resize(0);
## 369 : dimss__.push_back(dims__);
## 370 : dims__.resize(0);
## 371 : dimss__.push_back(dims__);
## 372 : dims__.resize(0);
## 373 : dims__.push_back(n_condition);
## 374 : dimss__.push_back(dims__);
## 375 : }
## 376 :
## 377 : template <typename RNG>
## 378 : void write_array(RNG& base_rng__,
## 379 : std::vector<double>& params_r__,
## 380 : std::vector<int>& params_i__,
## 381 : std::vector<double>& vars__,
## 382 : bool include_tparams__ = true,
## 383 : bool include_gqs__ = true,
## 384 : std::ostream* pstream__ = 0) const {
## 385 : typedef double local_scalar_t__;
## 386 :
## 387 : vars__.resize(0);
## 388 : stan::io::reader<local_scalar_t__> in__(params_r__, params_i__);
## 389 : static const char* function__ = "model10837057f64f_150269bbfe815db848542eb814742386_namespace::write_array";
## 390 : (void) function__; // dummy to suppress unused var warning
## 391 :
## 392 : // read-transform, write parameters
## 393 : double overall_mean = in__.scalar_constrain();
## 394 : vars__.push_back(overall_mean);
## 395 :
## 396 : Eigen::Matrix<double, Eigen::Dynamic, 1> condition_zoffset = in__.vector_constrain(n_condition);
## 397 : size_t condition_zoffset_j_1_max__ = n_condition;
## 398 : for (size_t j_1__ = 0; j_1__ < condition_zoffset_j_1_max__; ++j_1__) {
## 399 : vars__.push_back(condition_zoffset(j_1__));
## 400 : }
## 401 :
## 402 : double response_sd = in__.scalar_lb_constrain(0);
## 403 : vars__.push_back(response_sd);
## 404 :
## 405 : double condition_mean_sd = in__.scalar_lb_constrain(0);
## 406 : vars__.push_back(condition_mean_sd);
## 407 :
## 408 : double lp__ = 0.0;
## 409 : (void) lp__; // dummy to suppress unused var warning
## 410 : stan::math::accumulator<double> lp_accum__;
## 411 :
## 412 : local_scalar_t__ DUMMY_VAR__(std::numeric_limits<double>::quiet_NaN());
## 413 : (void) DUMMY_VAR__; // suppress unused var warning
## 414 :
## 415 : if (!include_tparams__ && !include_gqs__) return;
## 416 :
## 417 : try {
## 418 : // declare and define transformed parameters
## 419 : current_statement_begin__ = 14;
## 420 : validate_non_negative_index("condition_mean", "n_condition", n_condition);
## 421 : Eigen::Matrix<double, Eigen::Dynamic, 1> condition_mean(n_condition);
## 422 : stan::math::initialize(condition_mean, DUMMY_VAR__);
## 423 : stan::math::fill(condition_mean, DUMMY_VAR__);
## 424 :
## 425 : // do transformed parameters statements
## 426 : current_statement_begin__ = 15;
## 427 : stan::math::assign(condition_mean, add(overall_mean, multiply(condition_zoffset, condition_mean_sd)));
## 428 :
## 429 : if (!include_gqs__ && !include_tparams__) return;
## 430 : // validate transformed parameters
## 431 : const char* function__ = "validate transformed params";
## 432 : (void) function__; // dummy to suppress unused var warning
## 433 :
## 434 : // write transformed parameters
## 435 : if (include_tparams__) {
## 436 : size_t condition_mean_j_1_max__ = n_condition;
## 437 : for (size_t j_1__ = 0; j_1__ < condition_mean_j_1_max__; ++j_1__) {
## 438 : vars__.push_back(condition_mean(j_1__));
## 439 : }
## 440 : }
## 441 : if (!include_gqs__) return;
## 442 : } catch (const std::exception& e) {
## 443 : stan::lang::rethrow_located(e, current_statement_begin__, prog_reader__());
## 444 : // Next line prevents compiler griping about no return
## 445 : throw std::runtime_error("*** IF YOU SEE THIS, PLEASE REPORT A BUG ***");
## 446 : }
## 447 : }
## 448 :
## 449 : template <typename RNG>
## 450 : void write_array(RNG& base_rng,
## 451 : Eigen::Matrix<double,Eigen::Dynamic,1>& params_r,
## 452 : Eigen::Matrix<double,Eigen::Dynamic,1>& vars,
## 453 : bool include_tparams = true,
## 454 : bool include_gqs = true,
## 455 : std::ostream* pstream = 0) const {
## 456 : std::vector<double> params_r_vec(params_r.size());
## 457 : for (int i = 0; i < params_r.size(); ++i)
## 458 : params_r_vec[i] = params_r(i);
## 459 : std::vector<double> vars_vec;
## 460 : std::vector<int> params_i_vec;
## 461 : write_array(base_rng, params_r_vec, params_i_vec, vars_vec, include_tparams, include_gqs, pstream);
## 462 : vars.resize(vars_vec.size());
## 463 : for (int i = 0; i < vars.size(); ++i)
## 464 : vars(i) = vars_vec[i];
## 465 : }
## 466 :
## 467 : std::string model_name() const {
## 468 : return "model10837057f64f_150269bbfe815db848542eb814742386";
## 469 : }
## 470 :
## 471 :
## 472 : void constrained_param_names(std::vector<std::string>& param_names__,
## 473 : bool include_tparams__ = true,
## 474 : bool include_gqs__ = true) const {
## 475 : std::stringstream param_name_stream__;
## 476 : param_name_stream__.str(std::string());
## 477 : param_name_stream__ << "overall_mean";
## 478 : param_names__.push_back(param_name_stream__.str());
## 479 : size_t condition_zoffset_j_1_max__ = n_condition;
## 480 : for (size_t j_1__ = 0; j_1__ < condition_zoffset_j_1_max__; ++j_1__) {
## 481 : param_name_stream__.str(std::string());
## 482 : param_name_stream__ << "condition_zoffset" << '.' << j_1__ + 1;
## 483 : param_names__.push_back(param_name_stream__.str());
## 484 : }
## 485 : param_name_stream__.str(std::string());
## 486 : param_name_stream__ << "response_sd";
## 487 : param_names__.push_back(param_name_stream__.str());
## 488 : param_name_stream__.str(std::string());
## 489 : param_name_stream__ << "condition_mean_sd";
## 490 : param_names__.push_back(param_name_stream__.str());
## 491 :
## 492 : if (!include_gqs__ && !include_tparams__) return;
## 493 :
## 494 : if (include_tparams__) {
## 495 : size_t condition_mean_j_1_max__ = n_condition;
## 496 : for (size_t j_1__ = 0; j_1__ < condition_mean_j_1_max__; ++j_1__) {
## 497 : param_name_stream__.str(std::string());
## 498 : param_name_stream__ << "condition_mean" << '.' << j_1__ + 1;
## 499 : param_names__.push_back(param_name_stream__.str());
## 500 : }
## 501 : }
## 502 :
## 503 : if (!include_gqs__) return;
## 504 : }
## 505 :
## 506 :
## 507 : void unconstrained_param_names(std::vector<std::string>& param_names__,
## 508 : bool include_tparams__ = true,
## 509 : bool include_gqs__ = true) const {
## 510 : std::stringstream param_name_stream__;
## 511 : param_name_stream__.str(std::string());
## 512 : param_name_stream__ << "overall_mean";
## 513 : param_names__.push_back(param_name_stream__.str());
## 514 : size_t condition_zoffset_j_1_max__ = n_condition;
## 515 : for (size_t j_1__ = 0; j_1__ < condition_zoffset_j_1_max__; ++j_1__) {
## 516 : param_name_stream__.str(std::string());
## 517 : param_name_stream__ << "condition_zoffset" << '.' << j_1__ + 1;
## 518 : param_names__.push_back(param_name_stream__.str());
## 519 : }
## 520 : param_name_stream__.str(std::string());
## 521 : param_name_stream__ << "response_sd";
## 522 : param_names__.push_back(param_name_stream__.str());
## 523 : param_name_stream__.str(std::string());
## 524 : param_name_stream__ << "condition_mean_sd";
## 525 : param_names__.push_back(param_name_stream__.str());
## 526 :
## 527 : if (!include_gqs__ && !include_tparams__) return;
## 528 :
## 529 : if (include_tparams__) {
## 530 : size_t condition_mean_j_1_max__ = n_condition;
## 531 : for (size_t j_1__ = 0; j_1__ < condition_mean_j_1_max__; ++j_1__) {
## 532 : param_name_stream__.str(std::string());
## 533 : param_name_stream__ << "condition_mean" << '.' << j_1__ + 1;
## 534 : param_names__.push_back(param_name_stream__.str());
## 535 : }
## 536 : }
## 537 :
## 538 : if (!include_gqs__) return;
## 539 : }
## 540 :
## 541 : }; // model
## 542 :
## 543 : } // namespace
## 544 :
## 545 : typedef model10837057f64f_150269bbfe815db848542eb814742386_namespace::model10837057f64f_150269bbfe815db848542eb814742386 stan_model;
## 546 :
## 547 : #ifndef USING_R
## 548 :
## 549 : stan::model::model_base& new_model(
## 550 : stan::io::var_context& data_context,
## 551 : unsigned int seed,
## 552 : std::ostream* msg_stream) {
## 553 : stan_model* m = new stan_model(data_context, seed, msg_stream);
## 554 : return *m;
## 555 : }
## 556 :
## 557 : #endif
## 558 :
## 559 :
## 560 :
## 561 : #include <rstan_next/stan_fit.hpp>
## 562 :
## 563 : struct stan_model_holder {
## 564 : stan_model_holder(rstan::io::rlist_ref_var_context rcontext,
## 565 : unsigned int random_seed)
## 566 : : rcontext_(rcontext), random_seed_(random_seed)
## 567 : {
## 568 : }
## 569 :
## 570 : //stan::math::ChainableStack ad_stack;
## 571 : rstan::io::rlist_ref_var_context rcontext_;
## 572 : unsigned int random_seed_;
## 573 : };
## 574 :
## 575 : Rcpp::XPtr<stan::model::model_base> model_ptr(stan_model_holder* smh) {
## 576 : Rcpp::XPtr<stan::model::model_base> model_instance(new stan_model(smh->rcontext_, smh->random_seed_), true);
## 577 : return model_instance;
## 578 : }
## 579 :
## 580 : Rcpp::XPtr<rstan::stan_fit_base> fit_ptr(stan_model_holder* smh) {
## 581 : return Rcpp::XPtr<rstan::stan_fit_base>(new rstan::stan_fit(model_ptr(smh), smh->random_seed_), true);
## 582 : }
## 583 :
## 584 : std::string model_name(stan_model_holder* smh) {
## 585 : return model_ptr(smh).get()->model_name();
## 586 : }
## 587 :
## 588 : RCPP_MODULE(stan_fit4model10837057f64f_150269bbfe815db848542eb814742386_mod){
## 589 : Rcpp::class_<stan_model_holder>("stan_fit4model10837057f64f_150269bbfe815db848542eb814742386")
## 590 : .constructor<rstan::io::rlist_ref_var_context, unsigned int>()
## 591 : .method("model_ptr", &model_ptr)
## 592 : .method("fit_ptr", &fit_ptr)
## 593 : .method("model_name", &model_name)
## 594 : ;
## 595 : }
## 596 :
## 597 :
## 598 : // declarations
## 599 : extern "C" {
## 600 : SEXP file10837406d0ec( ) ;
## 601 : }
## 602 :
## 603 : // definition
## 604 : SEXP file10837406d0ec() {
## 605 : return Rcpp::wrap("150269bbfe815db848542eb814742386");
## 606 : }compose_data(ABC)## $condition
## [1] 1 2 3 4 5 1 2 3 4 5 1 2 3 4 5 1 2 3 4 5 1 2 3 4 5 1 2 3 4 5 1 2 3 4 5 1 2 3
## [39] 4 5 1 2 3 4 5 1 2 3 4 5
##
## $n_condition
## [1] 5
##
## $response
## [1] -0.42042774 1.69217967 1.37225407 1.03507138 -0.14427956 -0.30145399
## [7] 0.76391681 1.68231434 0.85711318 -0.93094589 0.61381517 0.59911027
## [13] 1.45980370 0.92123282 -1.53588002 -0.06949307 0.70134345 0.90801662
## [19] 1.12040863 -1.12967770 0.45025597 1.47093470 2.73398095 1.35338054
## [25] -0.59049553 -0.14674092 1.70929454 2.74938691 0.67145895 -1.42639772
## [31] 0.15795752 1.55484708 3.10773029 1.60855182 -0.26038911 0.47578692
## [37] 0.49523368 0.99976363 0.11890706 -1.07130406 0.77503018 0.59878841
## [43] 1.96271054 1.94783398 -1.22828447 0.28111168 0.55649574 1.76987771
## [49] 0.63783576 -1.03460558
##
## $n
## [1] 50m = sampling(ABC_stan, data = compose_data(ABC), control = list(adapt_delta=0.99))print(m, pars = c("overall_mean", "condition_mean_sd", "condition_mean", "response_sd"))## Inference for Stan model: 150269bbfe815db848542eb814742386.
## 4 chains, each with iter=2000; warmup=1000; thin=1;
## post-warmup draws per chain=1000, total post-warmup draws=4000.
##
## mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
## overall_mean 0.62 0.03 0.70 -0.65 0.29 0.62 0.98 1.73 525 1.01
## condition_mean_sd 1.24 0.02 0.55 0.62 0.89 1.11 1.42 2.62 525 1.01
## condition_mean[1] 0.20 0.00 0.18 -0.15 0.08 0.20 0.31 0.54 4543 1.00
## condition_mean[2] 1.00 0.00 0.18 0.65 0.89 1.00 1.12 1.34 4735 1.00
## condition_mean[3] 1.84 0.00 0.18 1.48 1.72 1.84 1.96 2.19 4117 1.00
## condition_mean[4] 1.02 0.00 0.18 0.66 0.90 1.02 1.14 1.36 4359 1.00
## condition_mean[5] -0.88 0.00 0.18 -1.22 -1.00 -0.89 -0.77 -0.54 3927 1.00
## response_sd 0.56 0.00 0.06 0.46 0.52 0.56 0.60 0.70 1602 1.00
##
## Samples were drawn using NUTS(diag_e) at Wed Jan 20 07:14:11 2021.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at
## convergence, Rhat=1).str(rstan::extract(m))## List of 6
## $ overall_mean : num [1:4000(1d)] 0.948 0.765 1.313 0.762 0.784 ...
## ..- attr(*, "dimnames")=List of 1
## .. ..$ iterations: NULL
## $ condition_zoffset: num [1:4000, 1:5] -0.771 -0.767 -0.631 -0.688 -0.167 ...
## ..- attr(*, "dimnames")=List of 2
## .. ..$ iterations: NULL
## .. ..$ : NULL
## $ response_sd : num [1:4000(1d)] 0.63 0.498 0.572 0.556 0.572 ...
## ..- attr(*, "dimnames")=List of 1
## .. ..$ iterations: NULL
## $ condition_mean_sd: num [1:4000(1d)] 0.888 0.958 1.686 1.017 1.06 ...
## ..- attr(*, "dimnames")=List of 1
## .. ..$ iterations: NULL
## $ condition_mean : num [1:4000, 1:5] 0.2632 0.0303 0.25 0.0621 0.6069 ...
## ..- attr(*, "dimnames")=List of 2
## .. ..$ iterations: NULL
## .. ..$ : NULL
## $ lp__ : num [1:4000(1d)] 1.04082 1.62556 2.20505 0.65647 -0.00369 ...
## ..- attr(*, "dimnames")=List of 1
## .. ..$ iterations: NULLm %<>% recover_types(ABC)m %>%
spread_draws(condition_mean[condition]) %>%
head(10)## # A tibble: 10 x 5
## # Groups: condition [1]
## condition condition_mean .chain .iteration .draw
## <chr> <dbl> <int> <int> <int>
## 1 A -0.0199 1 1 1
## 2 A 0.647 1 2 2
## 3 A 0.145 1 3 3
## 4 A 0.384 1 4 4
## 5 A 0.298 1 5 5
## 6 A 0.494 1 6 6
## 7 A 0.138 1 7 7
## 8 A 0.252 1 8 8
## 9 A 0.272 1 9 9
## 10 A -0.0143 1 10 10m %>%
recover_types(ABC) %>%
spread_draws(condition_mean[condition]) %>%
head(10)## # A tibble: 10 x 5
## # Groups: condition [1]
## condition condition_mean .chain .iteration .draw
## <chr> <dbl> <int> <int> <int>
## 1 A -0.0199 1 1 1
## 2 A 0.647 1 2 2
## 3 A 0.145 1 3 3
## 4 A 0.384 1 4 4
## 5 A 0.298 1 5 5
## 6 A 0.494 1 6 6
## 7 A 0.138 1 7 7
## 8 A 0.252 1 8 8
## 9 A 0.272 1 9 9
## 10 A -0.0143 1 10 10Perbandingan dengan model lain melalui kompatibilitas dengan broom
Fungsi terjemahan seperti ggdist :: to_broom_names (), ggdist :: from_broom_names (), ggdist :: to_ggmcmc_names (), dll. Dapat digunakan untuk menerjemahkan antara bingkai data format rapi umum dengan skema penamaan yang berbeda. Ini memudahkan, misalnya untuk membandingkan ringkasan poin dan interval antara keluaran tidybayes dan model yang didukung oleh sapu :: tidy.
Misalnya, bandingkan dengan biasa disebut ordinary least squares (OLS) regression:
m_linear = lm(response ~ condition, data = ABC)linear_results = m_linear %>%
emmeans(~ condition) %>%
tidy(conf.int = TRUE) %>%
mutate(model = "OLS")
linear_results## # A tibble: 5 x 8
## condition estimate df conf.low conf.high statistic p.value model
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <chr>
## 1 A 0.182 45 -0.167 0.530 1.05 3.00e- 1 OLS
## 2 B 1.01 45 0.665 1.36 5.85 5.13e- 7 OLS
## 3 C 1.87 45 1.53 2.22 10.8 4.15e-14 OLS
## 4 D 1.03 45 0.678 1.38 5.93 3.97e- 7 OLS
## 5 E -0.935 45 -1.28 -0.586 -5.40 2.41e- 6 OLSKita juga dapat memperoleh kecocokan yang sesuai dari model bayes :
bayes_results = m %>%
spread_draws(condition_mean[condition]) %>%
median_qi(estimate = condition_mean) %>%
to_broom_names() %>%
mutate(model = "Bayes")
bayes_results## # A tibble: 5 x 8
## condition estimate conf.low conf.high .width .point .interval model
## <chr> <dbl> <dbl> <dbl> <dbl> <chr> <chr> <chr>
## 1 A 0.195 -0.149 0.540 0.95 median qi Bayes
## 2 B 1.00 0.647 1.34 0.95 median qi Bayes
## 3 C 1.84 1.48 2.19 0.95 median qi Bayes
## 4 D 1.02 0.655 1.36 0.95 median qi Bayes
## 5 E -0.885 -1.22 -0.543 0.95 median qi BayesIni membuatnya mudah untuk menggabungkan dua bingkai data yang rapi bersama-sama menggunakan bind_rows, dan memplotnya:
bind_rows(list(linear_results, bayes_results)) %>%
mutate(condition = fct_rev(condition)) %>%
ggplot(aes(y = condition, x = estimate, xmin = conf.low, xmax = conf.high, color = model)) +
geom_pointinterval(position = position_dodge(width = .3))
Daftar pustaka :