diff --git a/articles/performance.html b/articles/performance.html index 4ba4cdf..794108d 100644 --- a/articles/performance.html +++ b/articles/performance.html @@ -74,14 +74,12 @@

Using default slice function
 res <- compare_eta_comptime_across_nvars(
-  n_vars = c(2, seq(from = 20, to = 200, by = 20)),
+  n_vars = c(2, seq(from = 40, to = 400, by = 40)),
   n_samples = 1,
   burnin = 0,
   w = 0.5)
-plot_eta_comptime(res, facet_by = "w")
-#> Don't know how to automatically pick scale for object of type <difftime>.
-#> Defaulting to continuous.
+plot_eta_comptime(res, facet_by = "w")

Could run over different slice widths if interested

diff --git a/articles/performance_files/figure-html/show-compare-1.png b/articles/performance_files/figure-html/show-compare-1.png
index fe2364f..8b58f33 100644
Binary files a/articles/performance_files/figure-html/show-compare-1.png and b/articles/performance_files/figure-html/show-compare-1.png differ
diff --git a/articles/pospkg_files/figure-html/normal-normal-1.png b/articles/pospkg_files/figure-html/normal-normal-1.png
index c4b7b05..0fe6efd 100644
Binary files a/articles/pospkg_files/figure-html/normal-normal-1.png and b/articles/pospkg_files/figure-html/normal-normal-1.png differ
diff --git a/articles/pospkg_files/figure-html/slice-elliptical-1.png b/articles/pospkg_files/figure-html/slice-elliptical-1.png
index ef9dba1..e96fc66 100644
Binary files a/articles/pospkg_files/figure-html/slice-elliptical-1.png and b/articles/pospkg_files/figure-html/slice-elliptical-1.png differ
diff --git a/articles/pospkg_files/figure-html/vary-df-1.png b/articles/pospkg_files/figure-html/vary-df-1.png
index 7e2b30d..5b5cef1 100644
Binary files a/articles/pospkg_files/figure-html/vary-df-1.png and b/articles/pospkg_files/figure-html/vary-df-1.png differ
diff --git a/articles/pospkg_files/figure-html/vary-w-1.png b/articles/pospkg_files/figure-html/vary-w-1.png
index b1cf70d..1b8daa8 100644
Binary files a/articles/pospkg_files/figure-html/vary-w-1.png and b/articles/pospkg_files/figure-html/vary-w-1.png differ
diff --git a/pkgdown.yml b/pkgdown.yml
index aa2f0a1..ef6cfcc 100644
--- a/pkgdown.yml
+++ b/pkgdown.yml
@@ -5,7 +5,7 @@ articles:
   customising: customising.html
   performance: performance.html
   pospkg: pospkg.html
-last_built: 2024-12-08T22:25Z
+last_built: 2024-12-08T23:32Z
 urls:
   reference: https://mathiaslj.github.io/mcmcglm/reference
   article: https://mathiaslj.github.io/mcmcglm/articles
diff --git a/reference/compare_eta_comptime_across_nvars.html b/reference/compare_eta_comptime_across_nvars.html
index b9e607d..263d6bc 100644
--- a/reference/compare_eta_comptime_across_nvars.html
+++ b/reference/compare_eta_comptime_across_nvars.html
@@ -145,39 +145,40 @@ 
Examplescompare_eta_comptime_across_nvars(n_vars = c(2, seq(from = 20, to = 60, by = 20)),
                                   n_samples = 100,
                                   burnin = 20)
-#> Sampling from posterior ■■■■■■■■                          22% | ETA:  4s
-#> Sampling from posterior ■■■■■■■■■■■■■■■■■■■■■■■■■■■■      91% | ETA:  0s
+#> Sampling from posterior ■■■■■■■                           21% | ETA:  4s
+#> Sampling from posterior ■■■■■■■■■■■■■■■■■■■■■■■■■■■       87% | ETA:  0s
 #> Sampling from posterior ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■  100% | ETA:  0s
-#> Sampling from posterior ■■■■■■■■■■■■■■■■■■■■■■            69% | ETA:  1s
+#> Sampling from posterior ■■■■■■■■■■■■■■■■■■■■              63% | ETA:  1s
 #> Sampling from posterior ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■  100% | ETA:  0s
-#> Sampling from posterior ■■■■                               9% | ETA: 20s
-#> Sampling from posterior ■■■■■■■■■■■■■■■■                  49% | ETA:  5s
-#> Sampling from posterior ■■■■■■■■■■■■■■■■■■■■■■■■■■■■      89% | ETA:  1s
+#> Sampling from posterior ■■■                                6% | ETA: 27s
+#> Sampling from posterior ■■■■■■■■■■■■■■■                   45% | ETA:  6s
+#> Sampling from posterior ■■■■■■■■■■■■■■■■■■■■■■■■■■■       85% | ETA:  1s
 #> Sampling from posterior ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■  100% | ETA:  0s
-#> Sampling from posterior ■■■■■                             13% | ETA: 15s
-#> Sampling from posterior ■■■■■■■■■■■■■■■■                  49% | ETA:  5s
-#> Sampling from posterior ■■■■■■■■■■■■■■■■■■■■■■■■■■■■      89% | ETA:  1s
+#> Sampling from posterior ■■■■                               9% | ETA: 19s
+#> Sampling from posterior ■■■■■■■■■■■■■■■                   45% | ETA:  6s
+#> Sampling from posterior ■■■■■■■■■■■■■■■■■■■■■■■■■■        84% | ETA:  1s
 #> Sampling from posterior ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■  100% | ETA:  0s
-#> Sampling from posterior ■■                                 2% | ETA:  2m
-#> Sampling from posterior ■■■■■                             12% | ETA: 38s
-#> Sampling from posterior ■■■■■■■■■■■                       32% | ETA: 17s
-#> Sampling from posterior ■■■■■■■■■■■■■■■■■                 53% | ETA: 10s
-#> Sampling from posterior ■■■■■■■■■■■■■■■■■■■■■■■■          76% | ETA:  4s
+#> Sampling from posterior ■                                  1% | ETA:  3m
+#> Sampling from posterior ■■■■                              10% | ETA: 44s
+#> Sampling from posterior ■■■■■■■■■■                        29% | ETA: 19s
+#> Sampling from posterior ■■■■■■■■■■■■■■■■                  51% | ETA: 10s
+#> Sampling from posterior ■■■■■■■■■■■■■■■■■■■■■■■           73% | ETA:  5s
+#> Sampling from posterior ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■    97% | ETA:  1s
 #> Sampling from posterior ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■  100% | ETA:  0s
-#> Sampling from posterior ■■■■                               9% | ETA: 30s
-#> Sampling from posterior ■■■■■■■■■■■                       33% | ETA: 12s
-#> Sampling from posterior ■■■■■■■■■■■■■■■■■■                58% | ETA:  7s
-#> Sampling from posterior ■■■■■■■■■■■■■■■■■■■■■■■■■■        83% | ETA:  2s
+#> Sampling from posterior ■■■                                6% | ETA: 39s
+#> Sampling from posterior ■■■■■■■■■■                        29% | ETA: 14s
+#> Sampling from posterior ■■■■■■■■■■■■■■■■■                 53% | ETA:  8s
+#> Sampling from posterior ■■■■■■■■■■■■■■■■■■■■■■■■          78% | ETA:  3s
 #> Sampling from posterior ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■  100% | ETA:  0s
 #>              time linear_predictor_calc n_vars n_samples beta_mean
-#> 1  0.3299866 secs                update      2       100         0
-#> 2  0.2936881 secs                 naive      2       100         0
-#> 3  3.6798954 secs                update     20       100         0
-#> 4  3.7449553 secs                 naive     20       100         0
-#> 5  8.7484875 secs                update     40       100         0
-#> 6  9.0414798 secs                 naive     40       100         0
-#> 7 17.2183602 secs                update     60       100         0
-#> 8 14.0114274 secs                 naive     60       100         0
+#> 1  0.3153636 secs                update      2       100         0
+#> 2  0.2818551 secs                 naive      2       100         0
+#> 3  3.7287543 secs                update     20       100         0
+#> 4  3.7928767 secs                 naive     20       100         0
+#> 5  8.8522904 secs                update     40       100         0
+#> 6  9.1041164 secs                 naive     40       100         0
+#> 7 17.2623024 secs                update     60       100         0
+#> 8 14.2019286 secs                 naive     60       100         0
 #>   beta_variance   family sd                 qslice_fun   w parallelised
 #> 1             1 gaussian  1 qslice::slice_stepping_out 0.5        FALSE
 #> 2             1 gaussian  1 qslice::slice_stepping_out 0.5        FALSE
diff --git a/reference/plot_eta_comptime-1.png b/reference/plot_eta_comptime-1.png
index 90e6d92..a9355a4 100644
Binary files a/reference/plot_eta_comptime-1.png and b/reference/plot_eta_comptime-1.png differ
diff --git a/reference/plot_eta_comptime.html b/reference/plot_eta_comptime.html
index dd91650..52683c9 100644
--- a/reference/plot_eta_comptime.html
+++ b/reference/plot_eta_comptime.html
@@ -76,33 +76,32 @@ 
Examples                                         n_samples = 100,
                                          burnin = 20)
 #> Sampling from posterior ■■■■■■                            18% | ETA:  5s
-#> Sampling from posterior ■■■■■■■■■■■■■■■■■■■■■■■■          76% | ETA:  1s
+#> Sampling from posterior ■■■■■■■■■■■■■■■■■■■■■■            69% | ETA:  1s
 #> Sampling from posterior ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■  100% | ETA:  0s
-#> Sampling from posterior ■■■■■■■■■■■■■■■■■■■               60% | ETA:  1s
+#> Sampling from posterior ■■■■■■■■■■■■■■■■■                 55% | ETA:  2s
 #> Sampling from posterior ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■  100% | ETA:  0s
-#> Sampling from posterior ■■■                                8% | ETA: 19s
-#> Sampling from posterior ■■■■■■■■■■■■■■                    44% | ETA:  6s
-#> Sampling from posterior ■■■■■■■■■■■■■■■■■■■■■■■■■■        82% | ETA:  2s
+#> Sampling from posterior ■■■                                5% | ETA: 26s
+#> Sampling from posterior ■■■■■■■■■■■■■                     41% | ETA:  6s
+#> Sampling from posterior ■■■■■■■■■■■■■■■■■■■■■■■■          78% | ETA:  2s
 #> Sampling from posterior ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■  100% | ETA:  0s
-#> Sampling from posterior ■■                                 2% | ETA:  1m
-#> Sampling from posterior ■■■■■■■■■■                        31% | ETA: 10s
-#> Sampling from posterior ■■■■■■■■■■■■■■■■■■■■■■■           74% | ETA:  3s
+#> Sampling from posterior ■                                  1% | ETA:  2m
+#> Sampling from posterior ■■■■■■■■■                         26% | ETA: 12s
+#> Sampling from posterior ■■■■■■■■■■■■■■■■■■■■■■            70% | ETA:  3s
 #> Sampling from posterior ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■  100% | ETA:  0s
 #> Sampling from posterior ■                                  1% | ETA:  4m
-#> Sampling from posterior ■■■                                6% | ETA:  1m
-#> Sampling from posterior ■■■■■■■■                          23% | ETA: 24s
-#> Sampling from posterior ■■■■■■■■■■■■■■                    43% | ETA: 14s
-#> Sampling from posterior ■■■■■■■■■■■■■■■■■■■■■             65% | ETA:  7s
-#> Sampling from posterior ■■■■■■■■■■■■■■■■■■■■■■■■■■■       88% | ETA:  2s
+#> Sampling from posterior ■■■                                5% | ETA:  1m
+#> Sampling from posterior ■■■■■■■                           21% | ETA: 26s
+#> Sampling from posterior ■■■■■■■■■■■■■                     41% | ETA: 14s
+#> Sampling from posterior ■■■■■■■■■■■■■■■■■■■■              63% | ETA:  8s
+#> Sampling from posterior ■■■■■■■■■■■■■■■■■■■■■■■■■■■       85% | ETA:  3s
 #> Sampling from posterior ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■  100% | ETA:  0s
-#> Sampling from posterior ■■                                 2% | ETA:  1m
-#> Sampling from posterior ■■■■■■■■■                         25% | ETA: 13s
-#> Sampling from posterior ■■■■■■■■■■■■■■■■                  49% | ETA:  8s
-#> Sampling from posterior ■■■■■■■■■■■■■■■■■■■■■■■           74% | ETA:  4s
+#> Sampling from posterior ■                                  1% | ETA:  2m
+#> Sampling from posterior ■■■■■■■■                          22% | ETA: 14s
+#> Sampling from posterior ■■■■■■■■■■■■■■■                   46% | ETA:  8s
+#> Sampling from posterior ■■■■■■■■■■■■■■■■■■■■■■            71% | ETA:  4s
+#> Sampling from posterior ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■    96% | ETA:  1s
 #> Sampling from posterior ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■  100% | ETA:  0s
 plot_eta_comptime(res)
-#> Don't know how to automatically pick scale for object of type <difftime>.
-#> Defaulting to continuous.
diff --git a/reference/plot_mcmcglm_across_tuningparams.html b/reference/plot_mcmcglm_across_tuningparams.html index c99928f..1801321 100644 --- a/reference/plot_mcmcglm_across_tuningparams.html +++ b/reference/plot_mcmcglm_across_tuningparams.html @@ -85,11 +85,11 @@

Examples family = "gaussian", data = dat_norm ) -#> Sampling from posterior ■■■■■■■■■■■■■■■ 47% | ETA: 1s +#> Sampling from posterior ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ 98% | ETA: 0s #> Sampling from posterior ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ 100% | ETA: 0s -#> Sampling from posterior ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ 92% | ETA: 0s +#> Sampling from posterior ■■■■■■■■■■■■■■ 45% | ETA: 1s #> Sampling from posterior ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ 100% | ETA: 0s -#> Sampling from posterior ■■■■■■■■■■■■■■■■■■■■■■■■■■ 82% | ETA: 0s +#> Sampling from posterior ■■■■■■■■■■■■■■■■■■■■■■■■■ 79% | ETA: 1s #> Sampling from posterior ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ 100% | ETA: 0s plot_mcmcglm_across_tuningparams(w05_mcmcglms) diff --git a/search.json b/search.json index c04ef7e..be83c55 100644 --- a/search.json +++ b/search.json @@ -1 +1 @@ -[{"path":"https://mathiaslj.github.io/mcmcglm/LICENSE.html","id":null,"dir":"","previous_headings":"","what":"MIT License","title":"MIT License","text":"Copyright (c) 2024 mcmcglm authors Permission hereby granted, free charge, person obtaining copy software associated documentation files (“Software”), deal Software without restriction, including without limitation rights use, copy, modify, merge, publish, distribute, sublicense, /sell copies Software, permit persons Software furnished , subject following conditions: copyright notice permission notice shall included copies substantial portions Software. SOFTWARE PROVIDED “”, WITHOUT WARRANTY KIND, EXPRESS IMPLIED, INCLUDING LIMITED WARRANTIES MERCHANTABILITY, FITNESS PARTICULAR PURPOSE NONINFRINGEMENT. EVENT SHALL AUTHORS COPYRIGHT HOLDERS LIABLE CLAIM, DAMAGES LIABILITY, WHETHER ACTION CONTRACT, TORT OTHERWISE, ARISING , CONNECTION SOFTWARE USE DEALINGS SOFTWARE.","code":""},{"path":"https://mathiaslj.github.io/mcmcglm/articles/customising.html","id":"adding-methods-for-new-families","dir":"Articles","previous_headings":"","what":"Adding methods for new families","title":"Adding methods for families","text":", easy extend package include family existing package. missing piece families method calculating log-density distribution. Since calculation log-density family implemented using S3 methods, user needs define function called log_density.*family_name* family_name specified relevant family. example showcasing single line code user needs write enable use inverse gaussian family","code":""},{"path":"https://mathiaslj.github.io/mcmcglm/articles/customising.html","id":"example-of-implementing-method-for-inverse-gaussian-family","dir":"Articles","previous_headings":"Adding methods for new families","what":"Example of implementing method for inverse gaussian family","title":"Adding methods for families","text":"First, generate data showcase example. user needs specify method computing density inverse gaussian distribution like : mu modelled mean μ=g−1(Xβ)\\mu=g^{-1}(X\\beta) GLM model, Y response variable. , mcmcglm() can called like : specify inverse.gaussian() family argument, make sure add additional arguments density log_likelihood_extra_args argument.","code":"n <- 1000 x1 <- rexp(n, 2) x2 <- rbinom(n, 1, .5) b0 <- 1 b1 <- 1.5 b2 <- 2 lin_pred <- b0+b1*x1+b2*x2 invgauss_fam <- inverse.gaussian() y_invnorm <- statmod::rinvgauss(n, mean = invgauss_fam$linkinv(lin_pred), shape = 1, dispersion = 1) dat_invnorm <- data.frame(Y = y_invnorm, X1 = x1, X2 = x2) log_density.inverse.gaussian <- function(family, mu, Y, ...) { statmod::dinvgauss(Y, mean = mu, ..., log = T) } invnorm <- mcmcglm::mcmcglm(formula = Y ~ ., family = inverse.gaussian(), log_likelihood_extra_args = list(shape = 1, dispersion = 1), data = dat_invnorm, beta_prior = distributional::dist_gamma(1, 1), w = 0.5) trace_plot(invnorm)"},{"path":[]},{"path":"https://mathiaslj.github.io/mcmcglm/articles/performance.html","id":"using-default-slice-function","dir":"Articles","previous_headings":"Graph of compute time as function of dimension","what":"Using default slice function","title":"Performance of package","text":"run different slice widths interested","code":"res <- compare_eta_comptime_across_nvars( n_vars = c(2, seq(from = 20, to = 200, by = 20)), n_samples = 1, burnin = 0, w = 0.5) plot_eta_comptime(res, facet_by = \"w\") #> Don't know how to automatically pick scale for object of type . #> Defaulting to continuous. ws <- c(0.5, 10) res <- lapply(ws, function(w) { compare_eta_comptime_across_nvars( w = w, n_vars = c(2, seq(from = 10, to = 50, by = 10)), n_samples = 50, burnin = 1, parallelise = TRUE) }) %>% dplyr::bind_rows()"},{"path":[]},{"path":"https://mathiaslj.github.io/mcmcglm/articles/pospkg.html","id":"family-of-response","dir":"Articles","previous_headings":"","what":"Family of response","title":"Examples of family, prior and slice sampler specifications","text":"family response can family implemented family class, principle (see help(\"family\")). now, implemented families package Gaussian Binomial Poisson Negative binomial However, small amount work needed user enable use family. process expanding available families described vignette(\"customising\"). Due implementation using family object, also means compatible link function can used family, like fx. logit probit binomial family. Also, unexported check_family function package ensures accepted values family argument package’s functions character, family function family class object. different uses seen throughout examples vignette also available example section mcmcglm() documentation.","code":""},{"path":"https://mathiaslj.github.io/mcmcglm/articles/pospkg.html","id":"data-simulation-for-examples","dir":"Articles","previous_headings":"Family of response","what":"Data simulation for examples","title":"Examples of family, prior and slice sampler specifications","text":"order show package provides samples target posterior distribution, simulate data scenario known β\\beta vector order see empirical distribution coefficient seems match value simulated data . start generating linear predictor η=Xβ\\eta=X\\beta known values β\\beta sampled values making design matrix XX. , example can apply relevant inverse link function g−1g^{-1} obtain modelled mean μ=g−1(Xβ)\\mu=g^{-1}(X\\beta) GLM model. can simulate response sampling distribution corresponding family using modelled mean μ\\mu relevant parameter distribution. example gaussian family available introductory example README.","code":"n <- 1000 x1 <- rnorm (n) x2 <- rbinom (n, 1, .5) b0 <- 1 b1 <- 1.5 b2 <- 2 lin_pred <- b0+b1*x1+b2*x2"},{"path":"https://mathiaslj.github.io/mcmcglm/articles/pospkg.html","id":"example-for-binomial-family","dir":"Articles","previous_headings":"Family of response","what":"Example for binomial family","title":"Examples of family, prior and slice sampler specifications","text":"binomial family, commonly used link function logit given logit(p)=lnp1−p\\begin{align} logit(p)=\\ln\\frac{p}{1-p} \\end{align} link functions also available, including probit, cauchit, log cloglog (see help(family) information). ’ll show example first binomial family logit link function.","code":""},{"path":"https://mathiaslj.github.io/mcmcglm/articles/pospkg.html","id":"logit-link-function","dir":"Articles","previous_headings":"Family of response > Example for binomial family","what":"Logit link function","title":"Examples of family, prior and slice sampler specifications","text":"First, define data fit mcmcglm inspect trace plot.","code":"logit_binom <- binomial(link = \"logit\") mu_logit <- logit_binom$linkinv(lin_pred) y_logit <- rbinom(n, size = 1, prob = mu_logit) dat_logit <- data.frame(Y = y_logit, X1 = x1, X2 = x2) logit <- mcmcglm(formula = Y ~ ., data = dat_logit, family = binomial(link = \"logit\"), beta_prior = distributional::dist_normal(mean = 0, sd = 1), w = 0.5) logit #> Object of class 'mcmcglm' #> #> Call: mcmcglm(formula = Y ~ ., family = binomial(link = \"logit\"), data = dat_logit, #> beta_prior = distributional::dist_normal(mean = 0, sd = 1), #> w = 0.5) #> #> Average of parameter samples: #> (Intercept) X1 X2 #> 1 0.8802159 1.461209 2.061524 trace_plot(logit)"},{"path":"https://mathiaslj.github.io/mcmcglm/articles/pospkg.html","id":"probit-link-function","dir":"Articles","previous_headings":"Family of response > Example for binomial family","what":"Probit link function","title":"Examples of family, prior and slice sampler specifications","text":"First, define data fit mcmcglm inspect trace plot.","code":"probit_binom <- binomial(link = \"probit\") mu_probit <- probit_binom$linkinv(lin_pred) y_probit <- rbinom(n, size = 1, prob = mu_probit) dat_probit <- data.frame(Y = y_probit, X1 = x1, X2 = x2) probit <- mcmcglm(formula = Y ~ ., data = dat_probit, family = binomial(link = \"probit\"), beta_prior = distributional::dist_normal(mean = 0, sd = 1), w = 0.5) probit #> Object of class 'mcmcglm' #> #> Call: mcmcglm(formula = Y ~ ., family = binomial(link = \"probit\"), #> data = dat_probit, beta_prior = distributional::dist_normal(mean = 0, #> sd = 1), w = 0.5) #> #> Average of parameter samples: #> (Intercept) X1 X2 #> 1 0.8318177 1.374661 2.162339 trace_plot(probit)"},{"path":"https://mathiaslj.github.io/mcmcglm/articles/pospkg.html","id":"example-for-poisson-family","dir":"Articles","previous_headings":"Family of response","what":"Example for poisson family","title":"Examples of family, prior and slice sampler specifications","text":"poisson family, commonly used link function log, one showcase example binomial family shows ’s easily possible specify link function. Define data fit mcmcglm inspect trace plot.","code":"mu_log <- exp(lin_pred) y_pois <- rpois(n, lambda = mu_log) dat_pois <- data.frame(Y = y_pois, X1 = x1, X2 = x2) pois <- mcmcglm(formula = Y ~ ., data = dat_pois, family = \"poisson\", beta_prior = distributional::dist_normal(mean = 0, sd = 1), w = 0.5) pois #> Object of class 'mcmcglm' #> #> Call: mcmcglm(formula = Y ~ ., family = \"poisson\", data = dat_pois, #> beta_prior = distributional::dist_normal(mean = 0, sd = 1), #> w = 0.5) #> #> Average of parameter samples: #> (Intercept) X1 X2 #> 1 0.9969416 1.504678 1.997051 trace_plot(pois)"},{"path":"https://mathiaslj.github.io/mcmcglm/articles/pospkg.html","id":"example-for-negative-binomial-family","dir":"Articles","previous_headings":"Family of response","what":"Example for negative binomial family","title":"Examples of family, prior and slice sampler specifications","text":"Besides families available stats package (via help(\"family\")), families implemented various packages. Fx. MASS::negative.binomial() function allows use negative binomial family various inference procedures. distribution interesting choice modelling count data Poisson distribution appropriate due equality first second moment distribution. showcase use family using package. fit mcmcglm inspect trace plot.","code":"mu_log <- exp(lin_pred) y_nbinom <- rnbinom(n, size = 1, mu = mu_log) dat_nbinom <- data.frame(Y = y_nbinom, X1 = x1, X2 = x2) theta <- 3 nbinom <- mcmcglm(formula = Y ~ ., data = dat_nbinom, family = MASS::negative.binomial(theta), beta_prior = distributional::dist_normal(mean = 0, sd = 1), w = 0.5) nbinom #> Object of class 'mcmcglm' #> #> Call: mcmcglm(formula = Y ~ ., family = MASS::negative.binomial(theta), #> data = dat_nbinom, beta_prior = distributional::dist_normal(mean = 0, #> sd = 1), w = 0.5) #> #> Average of parameter samples: #> (Intercept) X1 X2 #> 1 0.9608508 1.431976 2.059662 trace_plot(nbinom)"},{"path":"https://mathiaslj.github.io/mcmcglm/articles/pospkg.html","id":"specifying-the-prior-distribution","dir":"Articles","previous_headings":"","what":"Specifying the prior distribution","title":"Examples of family, prior and slice sampler specifications","text":"package uses distributional package specify prior distribution β\\beta parameter mentioned documentation mcmcglm(). distributional package provides vectorised distribution objects specified parameters methods evaluating densities, CDFs, etc. generate random samples distribution. use package allows user great flexibility specifying varying independent univariate distributions multivariate distributions parameter vector illustrated examples . examples simply use gaussian family response introductory example README alternate prior specifications. Start defining data.","code":"y_norm <- rnorm(n, mean = lin_pred, sd = 1) dat_norm <- data.frame(Y = y_norm, X1 = x1, X2 = x2)"},{"path":"https://mathiaslj.github.io/mcmcglm/articles/pospkg.html","id":"example-of-iid-priors","dir":"Articles","previous_headings":"Specifying the prior distribution","what":"Example of iid priors","title":"Examples of family, prior and slice sampler specifications","text":"case iid priors, single univariate distribution specification needed: make example exact introductory example README, show posterior distribution affected prior mean different data(/likelihood) showing. course, distribution used using dist_ functions distributional package. , simply use normal distribution distributional::dist_normal(), might well distributional::dist_gamma(), distributional::dist_student_t(), etc. See Reference full list distributions.","code":"iid <- mcmcglm(formula = Y ~ ., family = \"gaussian\", data = dat_norm, beta_prior = distributional::dist_normal(1000, 1), w = 0.5) trace_plot(iid)"},{"path":"https://mathiaslj.github.io/mcmcglm/articles/pospkg.html","id":"example-of-independent-and-potentially-different-priors","dir":"Articles","previous_headings":"Specifying the prior distribution","what":"Example of independent and potentially different priors","title":"Examples of family, prior and slice sampler specifications","text":"Due S3 method implementation mcmcglm package, ’s possible specify beta_prior list distribution objects. illustrative purposes, use distributions mentioned example specify different prior β\\beta component.","code":"indep <- mcmcglm(formula = Y ~ ., family = \"gaussian\", data = dat_norm, beta_prior = list( distributional::dist_normal(mean = 0, sd = 1), distributional::dist_gamma(shape = 1, rate = 2), distributional::dist_student_t(df = 8, mu = 3, sigma = 1) ), w = 0.5) trace_plot(indep)"},{"path":"https://mathiaslj.github.io/mcmcglm/articles/pospkg.html","id":"example-of-multivariate-prior","dir":"Articles","previous_headings":"Specifying the prior distribution","what":"Example of multivariate prior","title":"Examples of family, prior and slice sampler specifications","text":"implemented multivariate distribution distributional::dist_multivariate_normal(), can use specify different moments potential correlation parameters model. example shown prior distribution specified β∼𝒩3((123),[10.50.20.510.50.20.51])\\begin{align} \\beta\\sim \\mathcal{N}_3\\left( \\begin{pmatrix} 1\\\\2\\\\3 \\end{pmatrix}, \\begin{bmatrix} 1 & 0.5 & 0.2\\\\ 0.5 & 1 & 0.5\\\\ 0.2 & 0.5 & 1 \\end{bmatrix} \\right) \\end{align}","code":"mult <- mcmcglm(formula = Y ~ ., family = \"gaussian\", data = dat_norm, beta_prior = distributional::dist_multivariate_normal( mu = list(1:3), sigma = list( matrix( c(1, 0.5, 0.2, 0.5, 1, 0.5, 0.2, 0.5, 1), ncol = 3, byrow = TRUE)) ), w = 0.5) trace_plot(mult)"},{"path":"https://mathiaslj.github.io/mcmcglm/articles/pospkg.html","id":"tuning-and-changing-the-slice-sampler-procedure","dir":"Articles","previous_headings":"","what":"Tuning (and changing) the slice sampler procedure","title":"Examples of family, prior and slice sampler specifications","text":"sampling procedure package uses compute graph gibbs (CGGibbs) sampler introduced Luu et al. 2024. done efficiently “update” linear predictor gibbs update single component parameter vector rather naively recalculating linear predictor using terms.","code":""},{"path":"https://mathiaslj.github.io/mcmcglm/articles/pospkg.html","id":"advantage-of-cggibbs","dir":"Articles","previous_headings":"Tuning (and changing) the slice sampler procedure","what":"Advantage of CGGibbs","title":"Examples of family, prior and slice sampler specifications","text":"See details article, ’ll quickly show simple elegant idea . Assuming kk’th iteration (KK total iterations) sampled new jj’th coordinate dd-dimensional parameter vector β\\beta, βj(k)\\beta_j^{(k)}, “current” value β(j)(k)=(β1(k),β2(k),...,βj(k),βj+1(k−1),...,βd(k−1))\\begin{align} \\beta_{(j)}^{(k)}=\\left(\\beta_1^{(k)}, \\beta_2^{(k)}, ..., \\beta_j^{(k)}, \\beta_{j+1}^{(k-1)}, ..., \\beta_{d}^{(k-1)}\\right) \\end{align} CGGibbs approach calculate linear predictor kk’th iteration component jj, η(j)(k)\\eta_{(j)}^{(k)}, η(j)(k)=η(j−1)(k)+∑=1nxij(βj(k)−βj(k−1))\\begin{align} \\eta_{(j)}^{(k)}=\\eta_{(j-1)}^{(k)}+\\sum_{=1}^n x_{ij} (\\beta_j^{(k)}-\\beta_j^{(k-1)}) \\end{align} naive approach ηj(k)=∑=1n∑l=1dxilβl(k−𝟏[j>l])\\begin{align} \\eta_j^{(k)}=\\sum_{=1}^n\\sum_{l=1}^d x_{il} \\beta_l^{(k-\\mathbf{1}[j>l])} \\end{align} means using CGGibbs, component update β\\beta, computation time O(n)O(n), compared O(dn)O(dn) naive approach.","code":""},{"path":"https://mathiaslj.github.io/mcmcglm/articles/pospkg.html","id":"advantages-of-slice-sampling-in-combination-with-cggibbs","dir":"Articles","previous_headings":"Tuning (and changing) the slice sampler procedure","what":"Advantages of slice sampling (in combination with CGGibbs)","title":"Examples of family, prior and slice sampler specifications","text":"linear predictor used calculate log-likelihood, summed log-density prior distribution obtain log-potential, log-posterior normalising constant. log-potential calculated algorithm order use slice sampling, enables us use combination family response prior distribution parameters rely conjugacy get closed-form distributions sample . , slice sampling avoids spending computational power numerical integration find normalising constant corresponding marginal distribution data sampling procedure needs know posterior proportionality; reason use log-potential. default method slice sampling “stepping-” “shrinkage” described Neal 2003.","code":""},{"path":"https://mathiaslj.github.io/mcmcglm/articles/pospkg.html","id":"using-the-default-slice-sampling-method","dir":"Articles","previous_headings":"Tuning (and changing) the slice sampler procedure","what":"Using the default slice sampling method","title":"Examples of family, prior and slice sampler specifications","text":"slice_ functions qslice package arguments x log_target. mcmcglm function finds entities new component β\\beta vector given iteration log-potential evaluated function new value parameter. Thus, task left user specify tuning parameters chosen algorithm, passed ... mcmcglm(). default slice sampler uses qslice::slice_stepping_out() function, mandatory non-default argument w slice width algorithm. Thus, readers might noticed examples within vignette README, specification w = x given calls mcmcglm(). Example completeness:","code":"norm <- mcmcglm(formula = Y ~ ., family = \"gaussian\", data = dat_norm, beta_prior = distributional::dist_normal(0, 1), w = 0.5)"},{"path":"https://mathiaslj.github.io/mcmcglm/articles/pospkg.html","id":"changing-the-default-slice-sampling-method","dir":"Articles","previous_headings":"Tuning (and changing) the slice sampler procedure","what":"Changing the default slice sampling method","title":"Examples of family, prior and slice sampler specifications","text":"argument qslice_fun mcmcglm() function allows user change default slice sampler qslice::slice_stepping_out() slice sampler qslice package. Thus, example using qslice::slice_elliptical() Murray et al. (2010), tuning parameters mu sigma:","code":"el <- mcmcglm(formula = Y ~ ., family = \"gaussian\", data = dat_norm, beta_prior = distributional::dist_normal(0, 1), qslice_fun = qslice::slice_elliptical, mu = 3, sigma = 2) trace_plot(el)"},{"path":"https://mathiaslj.github.io/mcmcglm/articles/pospkg.html","id":"investigating-effect-of-tuning-parameters-of-slice-sampler","dir":"Articles","previous_headings":"Tuning (and changing) the slice sampler procedure","what":"Investigating effect of tuning parameters of slice sampler","title":"Examples of family, prior and slice sampler specifications","text":"utility functions mcmcglm_across_tuningparams() plot_mcmcglm_across_tuningparams() can used create trace plots samples running algorithm different values tuning parameters.","code":""},{"path":"https://mathiaslj.github.io/mcmcglm/articles/pospkg.html","id":"varying-slice-width-for-the-default-slice-sampler","dir":"Articles","previous_headings":"Tuning (and changing) the slice sampler procedure > Investigating effect of tuning parameters of slice sampler","what":"Varying slice width for the default slice sampler","title":"Examples of family, prior and slice sampler specifications","text":"Fx. default choice slice sampler, qslice::slice_stepping_out(), tuning parameter slice width w. way create plot samples across different values slice width seen :","code":"w05_mcmcglms <- mcmcglm_across_tuningparams( seq(from = 0.5, by = 0.5, length.out = 4), tuning_parameter_name = \"w\", formula = Y ~ ., family = \"gaussian\", data = dat_norm, parallelise = FALSE) plot_mcmcglm_across_tuningparams(w05_mcmcglms)"},{"path":"https://mathiaslj.github.io/mcmcglm/articles/pospkg.html","id":"varying-degrees-of-freedom-in-a-generalized-elliptical-slice-sampler","dir":"Articles","previous_headings":"Tuning (and changing) the slice sampler procedure > Investigating effect of tuning parameters of slice sampler","what":"Varying degrees of freedom in a generalized elliptical slice sampler","title":"Examples of family, prior and slice sampler specifications","text":"Since utility function work generally single tuning parameter, example qslice::slice_elliptical sampling procedure described Nishihara et al. (2014), varying df parameter algorithm. Note single tuning parameter can “varied” time. example, keep sigma mu constant, vary df","code":"df1_mcmcglms <- mcmcglm_across_tuningparams( c(1, seq(from = 10, by = 10, length.out = 3)), tuning_parameter_name = \"df\", formula = Y ~ ., family = \"gaussian\", data = dat_norm, qslice_fun = qslice::slice_genelliptical, mu = 2, sigma = 1, parallelise = FALSE) plot_mcmcglm_across_tuningparams(df1_mcmcglms)"},{"path":"https://mathiaslj.github.io/mcmcglm/articles/pospkg.html","id":"normal-normal-sample-method","dir":"Articles","previous_headings":"Tuning (and changing) the slice sampler procedure","what":"Normal-normal sample method","title":"Examples of family, prior and slice sampler specifications","text":"Mostly testing purposes, package implementation sampling closed-form normal distribution case response prior normally distributed. Running requires specifying \"normal-normal\" argument sample_method mcmcglm() seen example :","code":"norm_norm <- mcmcglm(formula = Y ~ ., family = \"gaussian\", data = dat_norm, beta_prior = distributional::dist_normal(0, 1), sample_method = \"normal-normal\") trace_plot(norm_norm)"},{"path":"https://mathiaslj.github.io/mcmcglm/authors.html","id":null,"dir":"","previous_headings":"","what":"Authors","title":"Authors and Citation","text":"Mathias Lerbech Jeppesen. Author, maintainer. Julie Køhler Munkvad. Author.","code":""},{"path":"https://mathiaslj.github.io/mcmcglm/authors.html","id":"citation","dir":"","previous_headings":"","what":"Citation","title":"Authors and Citation","text":"Jeppesen M, Munkvad J (2024). mcmcglm: Gibbs Sampling Methods. R package version 0.0.0.9000, https://mathiaslj.github.io/mcmcglm/, https://github.com/mathiaslj/mcmcglm.","code":"@Manual{, title = {mcmcglm: Gibbs Sampling Methods}, author = {Mathias Lerbech Jeppesen and Julie Køhler Munkvad}, year = {2024}, note = {R package version 0.0.0.9000, https://mathiaslj.github.io/mcmcglm/}, url = {https://github.com/mathiaslj/mcmcglm}, }"},{"path":"https://mathiaslj.github.io/mcmcglm/index.html","id":"mcmcglm","dir":"","previous_headings":"","what":"Gibbs Sampling Methods","title":"Gibbs Sampling Methods","text":"mcmcglm package implements CGGibbs sampler article Gibbs sampling faster Hamiltonian Monte Carlo GLMs?, linear run time function number parameters GLM model due clever “update” linear predictor. See details section vignette package implemented way user can specify family response distribution prior β\\beta parameter, 𝔼[Y|Xβ]=g−1(Xβ)\\mathbb{E}[Y|X\\beta]=g^{-1}(X\\beta) link function gg specified family. See vignette(\"pospkg\").","code":""},{"path":"https://mathiaslj.github.io/mcmcglm/index.html","id":"installation","dir":"","previous_headings":"","what":"Installation","title":"Gibbs Sampling Methods","text":"Installation available GitHub :","code":"devtools::install_github(\"mathiaslj/mcmcglm\")"},{"path":"https://mathiaslj.github.io/mcmcglm/index.html","id":"example","dir":"","previous_headings":"","what":"Example","title":"Gibbs Sampling Methods","text":"first simulate data linear model use showcasing use mcmcglm gaussian family. use function mcmcglm similar interface glm function added mandatory specification prior parameter β\\beta default qslice::slice_stepping_out function, performs slice sampling described Neal 2003 slice width w needs specified. creates mcmcglm object prints summarising call function averages samples parameter GLM model.","code":"n <- 1000 x1 <- rnorm (n) x2 <- rbinom (n, 1, .5) b0 <- 1 b1 <- 1.5 b2 <- 2 lin_pred <- b0+b1*x1+b2*x2 y_norm <- rnorm(n, mean = lin_pred, sd = 1) dat_norm <- data.frame(Y = y_norm, X1 = x1, X2 = x2) norm <- mcmcglm(formula = Y ~ ., family = \"gaussian\", data = dat_norm, beta_prior = distributional::dist_normal(0, 1), w = 0.5) norm #> Object of class 'mcmcglm' #> #> Call: mcmcglm(formula = Y ~ ., family = \"gaussian\", data = dat_norm, #> beta_prior = distributional::dist_normal(0, 1), w = 0.5) #> #> Average of parameter samples: #> (Intercept) X1 X2 #> 1 1.011134 1.490459 2.026047"},{"path":"https://mathiaslj.github.io/mcmcglm/index.html","id":"investigating-results","dir":"","previous_headings":"Example","what":"Investigating results","title":"Gibbs Sampling Methods","text":"averages shown print method object can retrieved generic coef like : full data set samples can accessed samples function: trace plot can seen function trace_plot:","code":"coef(norm) #> (Intercept) X1 X2 #> 1 1.011134 1.490459 2.026047 head(samples(norm)) #> (Intercept) X1 X2 iteration burnin #> 1 0.6173367 -0.004541141 -0.09125636 0 TRUE #> 2 2.6508146 0.281295470 0.68343132 1 TRUE #> 3 0.8240996 0.324627659 2.30889073 2 TRUE #> 4 0.8170086 1.028326905 2.20351455 3 TRUE #> 5 0.8777350 1.592074284 2.16115289 4 TRUE #> 6 0.9092187 1.442872350 2.02913214 5 TRUE trace_plot(norm)"},{"path":"https://mathiaslj.github.io/mcmcglm/reference/coef.mcmcglm.html","id":null,"dir":"Reference","previous_headings":"","what":"S3 method for getting the average value of coefficients — coef.mcmcglm","title":"S3 method for getting the average value of coefficients — coef.mcmcglm","text":"S3 method getting average value coefficients","code":""},{"path":"https://mathiaslj.github.io/mcmcglm/reference/coef.mcmcglm.html","id":"ref-usage","dir":"Reference","previous_headings":"","what":"Usage","title":"S3 method for getting the average value of coefficients — coef.mcmcglm","text":"","code":"# S3 method for class 'mcmcglm' coef(x)"},{"path":"https://mathiaslj.github.io/mcmcglm/reference/coef.mcmcglm.html","id":"arguments","dir":"Reference","previous_headings":"","what":"Arguments","title":"S3 method for getting the average value of coefficients — coef.mcmcglm","text":"x mcmcglm object","code":""},{"path":"https://mathiaslj.github.io/mcmcglm/reference/compare_eta_comptime_across_nvars.html","id":null,"dir":"Reference","previous_headings":"","what":"Compare runtime using CGGibbs and naive approach to calculate linear predictor — compare_eta_comptime_across_nvars","title":"Compare runtime using CGGibbs and naive approach to calculate linear predictor — compare_eta_comptime_across_nvars","text":"comparison methods described vignette. user can specify different arguments alter specification comparison, possible specify values n_vars argument. possible change fact data simulated normal response normal explanatory variables.","code":""},{"path":"https://mathiaslj.github.io/mcmcglm/reference/compare_eta_comptime_across_nvars.html","id":"ref-usage","dir":"Reference","previous_headings":"","what":"Usage","title":"Compare runtime using CGGibbs and naive approach to calculate linear predictor — compare_eta_comptime_across_nvars","text":"","code":"compare_eta_comptime_across_nvars( n_vars, n = 100, beta_prior = distributional::dist_normal(0, 1), log_likelihood_extra_args = list(sd = 1), sample_method = c(\"slice_sampling\", \"normal-normal\"), qslice_fun = qslice::slice_stepping_out, ..., n_samples = 500, burnin = 100, parallelise = FALSE, n_cores = as.numeric(Sys.getenv(\"NUMBER_OF_PROCESSORS\")) - 1 )"},{"path":"https://mathiaslj.github.io/mcmcglm/reference/compare_eta_comptime_across_nvars.html","id":"arguments","dir":"Reference","previous_headings":"","what":"Arguments","title":"Compare runtime using CGGibbs and naive approach to calculate linear predictor — compare_eta_comptime_across_nvars","text":"n_vars numeric vector numbers normally distributed columns design matrix different runs n numeric number observations data analysed GLM beta_prior distribution object created function distributional package. fx. distributional::dist_normal(mean = 0, sd = 1). log_likelihood_extra_args named list arguments passed onto log_density function. Fx. specification log_likelihood_extra_args = list(sd = x) needed case family = \"gaussian\" sample_method character specifying method used sampling. default \"slice_sampling\" intended value cases, works specification family beta_prior. \"normal-normal\" uses conditional normal distribution sample case conjugate prior gaussian response beta_prior. Implemented testing purposes works niche case. qslice_fun function qslice package. Default qslice::slice_stepping_out uses slice sampler Neal 2003, functions available. ... arguments passed onto function specified qslice_fun. default qslice::slice_stepping_out w needs specified, fx. qslice::slice_elliptical, mu sigma need specified n_samples numeric number samples draw parameter(/variable) model burnin numeric number samples marked \"burnin\". Burnin samples included beta_mean calculation increase finite sample performance LLN estimate parallelise logical runs algorithm across values n_vars parallelised using future.apply::future_lapply n_cores numeric number cores use parallelisation. Default 1 less number available cores","code":""},{"path":"https://mathiaslj.github.io/mcmcglm/reference/compare_eta_comptime_across_nvars.html","id":"value","dir":"Reference","previous_headings":"","what":"Value","title":"Compare runtime using CGGibbs and naive approach to calculate linear predictor — compare_eta_comptime_across_nvars","text":"data.frame information computation time different values n_vars","code":""},{"path":"https://mathiaslj.github.io/mcmcglm/reference/compare_eta_comptime_across_nvars.html","id":"ref-examples","dir":"Reference","previous_headings":"","what":"Examples","title":"Compare runtime using CGGibbs and naive approach to calculate linear predictor — compare_eta_comptime_across_nvars","text":"","code":"# Compare the runtime for 2, 20, 40, 60 variables in the model compare_eta_comptime_across_nvars(n_vars = c(2, seq(from = 20, to = 60, by = 20)), n_samples = 100, burnin = 20) #> Sampling from posterior ■■■■■■■■ 22% | ETA: 4s #> Sampling from posterior ■■■■■■■■■■■■■■■■■■■■■■■■■■■■ 91% | ETA: 0s #> Sampling from posterior ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ 100% | ETA: 0s #> Sampling from posterior ■■■■■■■■■■■■■■■■■■■■■■ 69% | ETA: 1s #> Sampling from posterior ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ 100% | ETA: 0s #> Sampling from posterior ■■■■ 9% | ETA: 20s #> Sampling from posterior ■■■■■■■■■■■■■■■■ 49% | ETA: 5s #> Sampling from posterior ■■■■■■■■■■■■■■■■■■■■■■■■■■■■ 89% | ETA: 1s #> Sampling from posterior ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ 100% | ETA: 0s #> Sampling from posterior ■■■■■ 13% | ETA: 15s #> Sampling from posterior ■■■■■■■■■■■■■■■■ 49% | ETA: 5s #> Sampling from posterior ■■■■■■■■■■■■■■■■■■■■■■■■■■■■ 89% | ETA: 1s #> Sampling from posterior ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ 100% | ETA: 0s #> Sampling from posterior ■■ 2% | ETA: 2m #> Sampling from posterior ■■■■■ 12% | ETA: 38s #> Sampling from posterior ■■■■■■■■■■■ 32% | ETA: 17s #> Sampling from posterior ■■■■■■■■■■■■■■■■■ 53% | ETA: 10s #> Sampling from posterior ■■■■■■■■■■■■■■■■■■■■■■■■ 76% | ETA: 4s #> Sampling from posterior ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ 100% | ETA: 0s #> Sampling from posterior ■■■■ 9% | ETA: 30s #> Sampling from posterior ■■■■■■■■■■■ 33% | ETA: 12s #> Sampling from posterior ■■■■■■■■■■■■■■■■■■ 58% | ETA: 7s #> Sampling from posterior ■■■■■■■■■■■■■■■■■■■■■■■■■■ 83% | ETA: 2s #> Sampling from posterior ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ 100% | ETA: 0s #> time linear_predictor_calc n_vars n_samples beta_mean #> 1 0.3299866 secs update 2 100 0 #> 2 0.2936881 secs naive 2 100 0 #> 3 3.6798954 secs update 20 100 0 #> 4 3.7449553 secs naive 20 100 0 #> 5 8.7484875 secs update 40 100 0 #> 6 9.0414798 secs naive 40 100 0 #> 7 17.2183602 secs update 60 100 0 #> 8 14.0114274 secs naive 60 100 0 #> beta_variance family sd qslice_fun w parallelised #> 1 1 gaussian 1 qslice::slice_stepping_out 0.5 FALSE #> 2 1 gaussian 1 qslice::slice_stepping_out 0.5 FALSE #> 3 1 gaussian 1 qslice::slice_stepping_out 0.5 FALSE #> 4 1 gaussian 1 qslice::slice_stepping_out 0.5 FALSE #> 5 1 gaussian 1 qslice::slice_stepping_out 0.5 FALSE #> 6 1 gaussian 1 qslice::slice_stepping_out 0.5 FALSE #> 7 1 gaussian 1 qslice::slice_stepping_out 0.5 FALSE #> 8 1 gaussian 1 qslice::slice_stepping_out 0.5 FALSE"},{"path":"https://mathiaslj.github.io/mcmcglm/reference/log_density.html","id":null,"dir":"Reference","previous_headings":"","what":"S3 generic for calculating the log density of a distribution dispatched via a family — log_density","title":"S3 generic for calculating the log density of a distribution dispatched via a family — log_density","text":"methods parametrised \"mu\" takes additional arguments needed calculation ... argument","code":""},{"path":"https://mathiaslj.github.io/mcmcglm/reference/log_density.html","id":"ref-usage","dir":"Reference","previous_headings":"","what":"Usage","title":"S3 generic for calculating the log density of a distribution dispatched via a family — log_density","text":"","code":"log_density(family, mu, Y, ...)"},{"path":"https://mathiaslj.github.io/mcmcglm/reference/log_density.html","id":"arguments","dir":"Reference","previous_headings":"","what":"Arguments","title":"S3 generic for calculating the log density of a distribution dispatched via a family — log_density","text":"family description error distribution link function used model. can character string naming family function, family function result call family function. (See family details family functions.) mu numeric vector values \"main\" parameter distribution specified family argument Y numeric vector response variable evaluate density ... arguments passed relevant methods","code":""},{"path":"https://mathiaslj.github.io/mcmcglm/reference/log_density.html","id":"value","dir":"Reference","previous_headings":"","what":"Value","title":"S3 generic for calculating the log density of a distribution dispatched via a family — log_density","text":"numeric vector log_density values","code":""},{"path":"https://mathiaslj.github.io/mcmcglm/reference/log_density.html","id":"details","dir":"Reference","previous_headings":"","what":"Details","title":"S3 generic for calculating the log density of a distribution dispatched via a family — log_density","text":"Supported families gaussian, binomial, poisson negative.binomial. Implement S3 method add support new family. mcmcglm function called, mu modelled mean glm model (meaning inverse link linear predictor). Reference methods see parametrisation","code":""},{"path":"https://mathiaslj.github.io/mcmcglm/reference/log_likelihood.html","id":null,"dir":"Reference","previous_headings":"","what":"Calculate log likelihood parametrised by ","title":"Calculate log likelihood parametrised by ","text":"Calculate log likelihood parametrised \"mu\"","code":""},{"path":"https://mathiaslj.github.io/mcmcglm/reference/log_likelihood.html","id":"ref-usage","dir":"Reference","previous_headings":"","what":"Usage","title":"Calculate log likelihood parametrised by ","text":"","code":"log_likelihood(family, mu, Y, ...)"},{"path":"https://mathiaslj.github.io/mcmcglm/reference/log_likelihood.html","id":"arguments","dir":"Reference","previous_headings":"","what":"Arguments","title":"Calculate log likelihood parametrised by ","text":"family description error distribution link function used model. can character string naming family function, family function result call family function. (See family details family functions.) mu numeric vector values \"main\" parameter distribution specified family argument Y numeric vector response variable evaluate density ... arguments passed S3 generic log_density","code":""},{"path":"https://mathiaslj.github.io/mcmcglm/reference/log_likelihood.html","id":"value","dir":"Reference","previous_headings":"","what":"Value","title":"Calculate log likelihood parametrised by ","text":"numeric Value log-likelihood","code":""},{"path":"https://mathiaslj.github.io/mcmcglm/reference/log_likelihood.html","id":"details","dir":"Reference","previous_headings":"","what":"Details","title":"Calculate log likelihood parametrised by ","text":"See log_density details","code":""},{"path":"https://mathiaslj.github.io/mcmcglm/reference/log_likelihood.html","id":"ref-examples","dir":"Reference","previous_headings":"","what":"Examples","title":"Calculate log likelihood parametrised by ","text":"","code":"# Create a test data n <- 100 x1 <- rnorm (n) x2 <- rbinom (n, 1, .5) b0 <- 1 b1 <- 1.5 b2 <- 2 lin_pred <- b0+b1*x1+b2*x2 known_sigma <- 1 y_norm <- rnorm(n, mean = lin_pred, sd = known_sigma) model_matrix_norm <- as.matrix( data.frame(int = 1, X1 = x1, X2 = x2)) b_prior <- 1:3 mu <- model_matrix_norm %*% b_prior log_likelihood(family = gaussian, mu = mu, Y = y_norm, sd = known_sigma) #> [1] -186.0887"},{"path":"https://mathiaslj.github.io/mcmcglm/reference/log_potential_from_betaj.html","id":null,"dir":"Reference","previous_headings":"","what":"Calculate the log-potential (log-likelihood plus log-density of prior) — log_potential_from_betaj","title":"Calculate the log-potential (log-likelihood plus log-density of prior) — log_potential_from_betaj","text":"Calculates log-potential function new coordinate beta parameter vector. Done like use unexporte","code":""},{"path":"https://mathiaslj.github.io/mcmcglm/reference/log_potential_from_betaj.html","id":"ref-usage","dir":"Reference","previous_headings":"","what":"Usage","title":"Calculate the log-potential (log-likelihood plus log-density of prior) — log_potential_from_betaj","text":"","code":"log_potential_from_betaj( new_beta_j, j, current_beta, current_eta, Y, X, family, beta_prior, linear_predictor_calc = \"update\", ... )"},{"path":"https://mathiaslj.github.io/mcmcglm/reference/log_potential_from_betaj.html","id":"arguments","dir":"Reference","previous_headings":"","what":"Arguments","title":"Calculate the log-potential (log-likelihood plus log-density of prior) — log_potential_from_betaj","text":"new_beta_j numeric new value j'th component beta parameter vector j numeric index parameter vector current_beta current value beta parameter vector sampling procedure current_eta current value linear predictor corresponding current_beta value Y numeric vector response variable evaluate density X design matrix family description error distribution link function used model. can character string naming family function, family function result call family function. (See family details family functions.) beta_prior distribution object created function distributional package. fx. distributional::dist_normal(mean = 0, sd = 1). linear_predictor_calc character specifying method used calculate linear predictor step gibbs algorithm. Default \"update\", uses CGGibbs procedure described start section vignette. option \"naive\", usual ... arguments passed S3 generic log_density","code":""},{"path":"https://mathiaslj.github.io/mcmcglm/reference/log_potential_from_betaj.html","id":"value","dir":"Reference","previous_headings":"","what":"Value","title":"Calculate the log-potential (log-likelihood plus log-density of prior) — log_potential_from_betaj","text":"value log-potential changed j'th component current_beta new_beta_j","code":""},{"path":"https://mathiaslj.github.io/mcmcglm/reference/log_potential_from_betaj.html","id":"ref-examples","dir":"Reference","previous_headings":"","what":"Examples","title":"Calculate the log-potential (log-likelihood plus log-density of prior) — log_potential_from_betaj","text":"","code":"# Create a test data n <- 100 x1 <- rnorm (n) x2 <- rbinom (n, 1, .5) b0 <- 1 b1 <- 1.5 b2 <- 2 lin_pred <- b0+b1*x1+b2*x2 known_sigma <- 1 y_norm <- rnorm(n, mean = lin_pred, sd = known_sigma) model_matrix_norm <- as.matrix( data.frame(int = 1, X1 = x1, X2 = x2)) b_prior <- distributional::dist_normal(mean = 0, sd = 1) b_prior_init <- distributional::generate( b_prior, ncol(model_matrix_norm) )[[1]] eta_init <- model_matrix_norm %*% b_prior_init j <- 1 new_beta_j <- 4 log_potential_from_betaj(new_beta_j = new_beta_j, j = j, current_beta = b_prior_init, current_eta = eta_init, Y = y_norm, X = model_matrix_norm, family = gaussian, beta_prior = b_prior, sd = known_sigma) #> [1] -535.7085"},{"path":"https://mathiaslj.github.io/mcmcglm/reference/mcmcglm-package.html","id":null,"dir":"Reference","previous_headings":"","what":"mcmcglm: Gibbs Sampling Methods — mcmcglm-package","title":"mcmcglm: Gibbs Sampling Methods — mcmcglm-package","text":"Sample posterior distributions using variations slice within Gibbs samplers GLMs.","code":""},{"path":[]},{"path":"https://mathiaslj.github.io/mcmcglm/reference/mcmcglm-package.html","id":"author","dir":"Reference","previous_headings":"","what":"Author","title":"mcmcglm: Gibbs Sampling Methods — mcmcglm-package","text":"Maintainer: Mathias Lerbech Jeppesen mathiasljeppesen@outlook.com Authors: Julie Køhler Munkvad julie.m@hotmail.dk","code":""},{"path":"https://mathiaslj.github.io/mcmcglm/reference/mcmcglm.html","id":null,"dir":"Reference","previous_headings":"","what":"Efficient Gibbs sampling of posterior distribution of parameters in GLM — mcmcglm","title":"Efficient Gibbs sampling of posterior distribution of parameters in GLM — mcmcglm","text":"Obtain MCMC samples using slice sampling within Gibbs generalized linear models (GLMs) using compute graph Gibbs (CGGibbs) linear runtime number variables model matrix. Method described article Gibbs sampling faster Hamiltonian Monte Carlo GLMs?, see details .","code":""},{"path":"https://mathiaslj.github.io/mcmcglm/reference/mcmcglm.html","id":"ref-usage","dir":"Reference","previous_headings":"","what":"Usage","title":"Efficient Gibbs sampling of posterior distribution of parameters in GLM — mcmcglm","text":"","code":"mcmcglm( formula, family = gaussian, data, beta_prior = distributional::dist_normal(0, 1), log_likelihood_extra_args = list(sd = 1), linear_predictor_calc = c(\"update\", \"naive\"), sample_method = c(\"slice_sampling\", \"normal-normal\"), qslice_fun = qslice::slice_stepping_out, ..., n_samples = 500, burnin = 100 )"},{"path":"https://mathiaslj.github.io/mcmcglm/reference/mcmcglm.html","id":"arguments","dir":"Reference","previous_headings":"","what":"Arguments","title":"Efficient Gibbs sampling of posterior distribution of parameters in GLM — mcmcglm","text":"formula object class \"formula\" (one can coerced class): symbolic description model fitted. See details stats::glm family description error distribution link function used model. can character string naming family function, family function result call family function. (See family details family functions.) data optional data frame, list environment (object coercible .data.frame data frame) containing variables model. found data, variables taken environment(formula), typically environment function called. beta_prior distribution object created function distributional package. fx. distributional::dist_normal(mean = 0, sd = 1). log_likelihood_extra_args named list arguments passed onto log_density function. Fx. specification log_likelihood_extra_args = list(sd = x) needed case family = \"gaussian\" linear_predictor_calc character specifying method used calculate linear predictor step gibbs algorithm. Default \"update\", uses CGGibbs procedure described start section vignette. option \"naive\", usual sample_method character specifying method used sampling. default \"slice_sampling\" intended value cases, works specification family beta_prior. \"normal-normal\" uses conditional normal distribution sample case conjugate prior gaussian response beta_prior. Implemented testing purposes works niche case. qslice_fun function qslice package. Default qslice::slice_stepping_out uses slice sampler Neal 2003, functions available. ... arguments passed onto function specified qslice_fun. default qslice::slice_stepping_out w needs specified, fx. qslice::slice_elliptical, mu sigma need specified n_samples numeric number samples draw parameter(/variable) model burnin numeric number samples marked \"burnin\". Burnin samples included beta_mean calculation increase finite sample performance LLN estimate","code":""},{"path":"https://mathiaslj.github.io/mcmcglm/reference/mcmcglm.html","id":"value","dir":"Reference","previous_headings":"","what":"Value","title":"Efficient Gibbs sampling of posterior distribution of parameters in GLM — mcmcglm","text":"object class mcmcglm methods getting data.frame parameter samples, plotting, etc.","code":""},{"path":"https://mathiaslj.github.io/mcmcglm/reference/mcmcglm.html","id":"details","dir":"Reference","previous_headings":"","what":"Details","title":"Efficient Gibbs sampling of posterior distribution of parameters in GLM — mcmcglm","text":"uses updating scheme linear predictor draw Gibbs sampling coordinates parameter vector","code":""},{"path":"https://mathiaslj.github.io/mcmcglm/reference/mcmcglm.html","id":"references","dir":"Reference","previous_headings":"","what":"References","title":"Efficient Gibbs sampling of posterior distribution of parameters in GLM — mcmcglm","text":"Gibbs sampling faster Hamiltonian Monte Carlo GLMs?, Neal 2003","code":""},{"path":"https://mathiaslj.github.io/mcmcglm/reference/mcmcglm.html","id":"ref-examples","dir":"Reference","previous_headings":"","what":"Examples","title":"Efficient Gibbs sampling of posterior distribution of parameters in GLM — mcmcglm","text":"","code":"if (FALSE) { # \\dontrun{ # Create test data for different scenarios n <- 100 x1 <- rnorm (n) x2 <- rbinom (n, 1, .5) b0 <- 1 b1 <- 1.5 b2 <- 2 lin_pred <- b0+b1*x1+b2*x2 ############################################# # Different families and priors # For family \"gaussian\" and iid normal prior y_norm <- rnorm(n, mean = lin_pred, sd = 1) dat_norm <- data.frame(Y = y_norm, X1 = x1, X2 = x2) norm <- mcmcglm(formula = Y ~ ., data = dat_norm, beta_prior = distributional::dist_normal(0, 1), family = \"gaussian\", n_samples = 100, burnin = 10, sample_method = \"slice_sampling\", qslice_fun = qslice::slice_stepping_out, w = 0.5) norm # For family \"binomial\" with logit link and iid gamma distributed prior y_logit <- rbinom (n, 1, arm::invlogit(lin_pred)) dat_logit <- data.frame(Y = y_logit, X1 = x1, X2 = x2) logit <- mcmcglm(formula = Y ~ ., data = dat_logit, beta_prior = distributional::dist_gamma(shape = 1, rate = 0.5), family = binomial(link = \"logit\"), n_samples = 100, burnin = 10, sample_method = \"slice_sampling\", qslice_fun = qslice::slice_stepping_out, w = 0.8) logit # For family \"negative.binomial\" and multivariate normal specification of parameter priors y_log <- rnbinom(n, size = 1, mu = exp(lin_pred)) dat_log <- data.frame(Y = y_log, X1 = x1, X2 = x2) log <- mcmcglm(formula = Y ~ X1, data = dat_log, beta_prior = distributional::dist_multivariate_normal( mu = list(c(1, 2)), sigma = list(matrix(c(1, 0.5, 0.5, 1), ncol = 2)) ), family = MASS::negative.binomial(3), n_samples = 100, burnin = 10, sample_method = \"slice_sampling\", qslice_fun = qslice::slice_stepping_out, w = 0.8) log # For family \"negative.binomial\" and specification of different independent # priors for each parameter log2 <- mcmcglm(formula = Y ~ ., data = dat_log, beta_prior = list(distributional::dist_normal(0, 1), distributional::dist_gamma(1, 1), distributional::dist_exponential(2)), family = MASS::negative.binomial(3), n_samples = 100, burnin = 10, sample_method = \"slice_sampling\", qslice_fun = qslice::slice_stepping_out, w = 0.8) log2 ############################################# # Using a different slice function log3 <- mcmcglm(formula = Y ~ ., data = dat_log, beta_prior = list(distributional::dist_normal(0, 1), distributional::dist_gamma(1, 1), distributional::dist_exponential(2)), family = MASS::negative.binomial(3), n_samples = 100, burnin = 10, sample_method = \"slice_sampling\", qslice_fun = qslice::slice_elliptical, mu = 1.5, sigma = 2) log3 } # }"},{"path":"https://mathiaslj.github.io/mcmcglm/reference/mcmcglm_across_tuningparams.html","id":null,"dir":"Reference","previous_headings":"","what":"Get list of mcmcglms run across values of slice sampling tuning parameters — mcmcglm_across_tuningparams","title":"Get list of mcmcglms run across values of slice sampling tuning parameters — mcmcglm_across_tuningparams","text":"function simply performs lapply mcmcglm values provided tuning parameter values. See general usage vignette","code":""},{"path":"https://mathiaslj.github.io/mcmcglm/reference/mcmcglm_across_tuningparams.html","id":"ref-usage","dir":"Reference","previous_headings":"","what":"Usage","title":"Get list of mcmcglms run across values of slice sampling tuning parameters — mcmcglm_across_tuningparams","text":"","code":"mcmcglm_across_tuningparams( ..., tuning_parameter_name = \"w\", formula, family, data, beta_prior = distributional::dist_normal(0, 1), log_likelihood_extra_args = list(sd = 1), sample_method = c(\"slice_sampling\", \"normal-normal\"), qslice_fun = qslice::slice_stepping_out, n_samples = 500, burnin = 100, parallelise = FALSE, n_cores = as.numeric(Sys.getenv(\"NUMBER_OF_PROCESSORS\")) - 1 )"},{"path":"https://mathiaslj.github.io/mcmcglm/reference/mcmcglm_across_tuningparams.html","id":"arguments","dir":"Reference","previous_headings":"","what":"Arguments","title":"Get list of mcmcglms run across values of slice sampling tuning parameters — mcmcglm_across_tuningparams","text":"... list vector values tuning parameter tuning_parameter_name name tuning parameter. Fx. default qslice_fun qslice::slice_stepping_out, tuning parameter called w formula object class \"formula\" (one can coerced class): symbolic description model fitted. See details stats::glm family description error distribution link function used model. can character string naming family function, family function result call family function. (See family details family functions.) data optional data frame, list environment (object coercible .data.frame data frame) containing variables model. found data, variables taken environment(formula), typically environment function called. beta_prior distribution object created function distributional package. fx. distributional::dist_normal(mean = 0, sd = 1). log_likelihood_extra_args named list arguments passed onto log_density function. Fx. specification log_likelihood_extra_args = list(sd = x) needed case family = \"gaussian\" sample_method character specifying method used sampling. default \"slice_sampling\" intended value cases, works specification family beta_prior. \"normal-normal\" uses conditional normal distribution sample case conjugate prior gaussian response beta_prior. Implemented testing purposes works niche case. qslice_fun function qslice package. Default qslice::slice_stepping_out uses slice sampler Neal 2003, functions available. n_samples numeric number samples draw parameter(/variable) model burnin numeric number samples marked \"burnin\". Burnin samples included beta_mean calculation increase finite sample performance LLN estimate parallelise logical runs algorithm across values n_vars parallelised using future.apply::future_lapply n_cores numeric number cores use parallelisation. Default 1 less number available cores","code":""},{"path":"https://mathiaslj.github.io/mcmcglm/reference/mcmcglm_across_tuningparams.html","id":"ref-examples","dir":"Reference","previous_headings":"","what":"Examples","title":"Get list of mcmcglms run across values of slice sampling tuning parameters — mcmcglm_across_tuningparams","text":"","code":"# Create test data for different scenarios n <- 100 x1 <- rnorm (n) x2 <- rbinom (n, 1, .5) b0 <- 1 b1 <- 1.5 b2 <- 2 lin_pred <- b0+b1*x1+b2*x2 ############################################# # Different families and priors # For family \"gaussian\" and iid normal prior y_norm <- rnorm(n, mean = lin_pred, sd = 1) dat_norm <- data.frame(Y = y_norm, X1 = x1, X2 = x2) w05_mcmcglms <- mcmcglm_across_tuningparams( seq(from = 0.5, by = 0.5, length.out = 9), tuning_parameter_name = \"w\", formula = Y ~ ., family = \"gaussian\", data = dat_norm, n_samples = 10, burnin = 0 )"},{"path":"https://mathiaslj.github.io/mcmcglm/reference/plot_eta_comptime.html","id":null,"dir":"Reference","previous_headings":"","what":"Plot the results of compare_eta_comptime_across_nvars — plot_eta_comptime","title":"Plot the results of compare_eta_comptime_across_nvars — plot_eta_comptime","text":"Plot results compare_eta_comptime_across_nvars","code":""},{"path":"https://mathiaslj.github.io/mcmcglm/reference/plot_eta_comptime.html","id":"ref-usage","dir":"Reference","previous_headings":"","what":"Usage","title":"Plot the results of compare_eta_comptime_across_nvars — plot_eta_comptime","text":"","code":"plot_eta_comptime(eta_comptime_data, facet_by = \"qslice_fun\")"},{"path":"https://mathiaslj.github.io/mcmcglm/reference/plot_eta_comptime.html","id":"arguments","dir":"Reference","previous_headings":"","what":"Arguments","title":"Plot the results of compare_eta_comptime_across_nvars — plot_eta_comptime","text":"eta_comptime_data data.frame results computation times. Result call compare_eta_comptime_across_nvars facet_by character variable facet plots . Default \"qslice_fun\", enabling user combine results compare_eta_comptime_across_nvars across different slice functions plot easily. options fx. running compare_eta_comptime_across_nvars different values tuning parameter seeing affects runtime","code":""},{"path":"https://mathiaslj.github.io/mcmcglm/reference/plot_eta_comptime.html","id":"value","dir":"Reference","previous_headings":"","what":"Value","title":"Plot the results of compare_eta_comptime_across_nvars — plot_eta_comptime","text":"ggplot object","code":""},{"path":"https://mathiaslj.github.io/mcmcglm/reference/plot_eta_comptime.html","id":"ref-examples","dir":"Reference","previous_headings":"","what":"Examples","title":"Plot the results of compare_eta_comptime_across_nvars — plot_eta_comptime","text":"","code":"# Compare the runtime for 2, 20, 40, 60 variables in the model res <- compare_eta_comptime_across_nvars(n_vars = c(2, seq(from = 20, to = 60, by = 20)), n_samples = 100, burnin = 20) #> Sampling from posterior ■■■■■■ 18% | ETA: 5s #> Sampling from posterior ■■■■■■■■■■■■■■■■■■■■■■■■ 76% | ETA: 1s #> Sampling from posterior ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ 100% | ETA: 0s #> Sampling from posterior ■■■■■■■■■■■■■■■■■■■ 60% | ETA: 1s #> Sampling from posterior ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ 100% | ETA: 0s #> Sampling from posterior ■■■ 8% | ETA: 19s #> Sampling from posterior ■■■■■■■■■■■■■■ 44% | ETA: 6s #> Sampling from posterior ■■■■■■■■■■■■■■■■■■■■■■■■■■ 82% | ETA: 2s #> Sampling from posterior ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ 100% | ETA: 0s #> Sampling from posterior ■■ 2% | ETA: 1m #> Sampling from posterior ■■■■■■■■■■ 31% | ETA: 10s #> Sampling from posterior ■■■■■■■■■■■■■■■■■■■■■■■ 74% | ETA: 3s #> Sampling from posterior ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ 100% | ETA: 0s #> Sampling from posterior ■ 1% | ETA: 4m #> Sampling from posterior ■■■ 6% | ETA: 1m #> Sampling from posterior ■■■■■■■■ 23% | ETA: 24s #> Sampling from posterior ■■■■■■■■■■■■■■ 43% | ETA: 14s #> Sampling from posterior ■■■■■■■■■■■■■■■■■■■■■ 65% | ETA: 7s #> Sampling from posterior ■■■■■■■■■■■■■■■■■■■■■■■■■■■ 88% | ETA: 2s #> Sampling from posterior ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ 100% | ETA: 0s #> Sampling from posterior ■■ 2% | ETA: 1m #> Sampling from posterior ■■■■■■■■■ 25% | ETA: 13s #> Sampling from posterior ■■■■■■■■■■■■■■■■ 49% | ETA: 8s #> Sampling from posterior ■■■■■■■■■■■■■■■■■■■■■■■ 74% | ETA: 4s #> Sampling from posterior ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ 100% | ETA: 0s plot_eta_comptime(res) #> Don't know how to automatically pick scale for object of type . #> Defaulting to continuous."},{"path":"https://mathiaslj.github.io/mcmcglm/reference/plot_mcmcglm_across_tuningparams.html","id":null,"dir":"Reference","previous_headings":"","what":"Plot a list of mcmcglms showing varying tuning parameters in the title of the plots — plot_mcmcglm_across_tuningparams","title":"Plot a list of mcmcglms showing varying tuning parameters in the title of the plots — plot_mcmcglm_across_tuningparams","text":"See decription mcmcglm_across_tuningparams details functionality can used ","code":""},{"path":"https://mathiaslj.github.io/mcmcglm/reference/plot_mcmcglm_across_tuningparams.html","id":"ref-usage","dir":"Reference","previous_headings":"","what":"Usage","title":"Plot a list of mcmcglms showing varying tuning parameters in the title of the plots — plot_mcmcglm_across_tuningparams","text":"","code":"plot_mcmcglm_across_tuningparams(list_mcmcglms)"},{"path":"https://mathiaslj.github.io/mcmcglm/reference/plot_mcmcglm_across_tuningparams.html","id":"arguments","dir":"Reference","previous_headings":"","what":"Arguments","title":"Plot a list of mcmcglms showing varying tuning parameters in the title of the plots — plot_mcmcglm_across_tuningparams","text":"list_mcmcglms list mcmcglm objects. Intended result call mcmcglm_across_tuningparams","code":""},{"path":"https://mathiaslj.github.io/mcmcglm/reference/plot_mcmcglm_across_tuningparams.html","id":"ref-examples","dir":"Reference","previous_headings":"","what":"Examples","title":"Plot a list of mcmcglms showing varying tuning parameters in the title of the plots — plot_mcmcglm_across_tuningparams","text":"","code":"# Create test data for different scenarios n <- 100 x1 <- rnorm (n) x2 <- rbinom (n, 1, .5) b0 <- 1 b1 <- 1.5 b2 <- 2 lin_pred <- b0+b1*x1+b2*x2 ############################################# # Different families and priors # For family \"gaussian\" and iid normal prior y_norm <- rnorm(n, mean = lin_pred, sd = 1) dat_norm <- data.frame(Y = y_norm, X1 = x1, X2 = x2) w05_mcmcglms <- mcmcglm_across_tuningparams( seq(from = 0.5, by = 0.5, length.out = 4), tuning_parameter_name = \"w\", formula = Y ~ ., family = \"gaussian\", data = dat_norm ) #> Sampling from posterior ■■■■■■■■■■■■■■■ 47% | ETA: 1s #> Sampling from posterior ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ 100% | ETA: 0s #> Sampling from posterior ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ 92% | ETA: 0s #> Sampling from posterior ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ 100% | ETA: 0s #> Sampling from posterior ■■■■■■■■■■■■■■■■■■■■■■■■■■ 82% | ETA: 0s #> Sampling from posterior ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ 100% | ETA: 0s plot_mcmcglm_across_tuningparams(w05_mcmcglms)"},{"path":"https://mathiaslj.github.io/mcmcglm/reference/samples.html","id":null,"dir":"Reference","previous_headings":"","what":"Get the drawn samples from the object — samples","title":"Get the drawn samples from the object — samples","text":"Get drawn samples object","code":""},{"path":"https://mathiaslj.github.io/mcmcglm/reference/samples.html","id":"ref-usage","dir":"Reference","previous_headings":"","what":"Usage","title":"Get the drawn samples from the object — samples","text":"","code":"samples(x)"},{"path":"https://mathiaslj.github.io/mcmcglm/reference/samples.html","id":"arguments","dir":"Reference","previous_headings":"","what":"Arguments","title":"Get the drawn samples from the object — samples","text":"x mcmcglm object","code":""},{"path":"https://mathiaslj.github.io/mcmcglm/reference/trace_plot.html","id":null,"dir":"Reference","previous_headings":"","what":"Create a trace plot of the MCMC samples — trace_plot","title":"Create a trace plot of the MCMC samples — trace_plot","text":"Create trace plot MCMC samples","code":""},{"path":"https://mathiaslj.github.io/mcmcglm/reference/trace_plot.html","id":"ref-usage","dir":"Reference","previous_headings":"","what":"Usage","title":"Create a trace plot of the MCMC samples — trace_plot","text":"","code":"trace_plot(x, samples_drop = NULL)"},{"path":"https://mathiaslj.github.io/mcmcglm/reference/trace_plot.html","id":"arguments","dir":"Reference","previous_headings":"","what":"Arguments","title":"Create a trace plot of the MCMC samples — trace_plot","text":"x mcmcglm object samples_drop numeric specifying number initial samples exclude trace_plot improve axis zoom plot","code":""},{"path":"https://mathiaslj.github.io/mcmcglm/reference/update_linear_predictor.html","id":null,"dir":"Reference","previous_headings":"","what":"Update value of a linear predictor as function of a single coordinate change — update_linear_predictor","title":"Update value of a linear predictor as function of a single coordinate change — update_linear_predictor","text":"Function updating linear predictor n actions rather nn_vars actions naively matrix-vector product X %% beta","code":""},{"path":"https://mathiaslj.github.io/mcmcglm/reference/update_linear_predictor.html","id":"ref-usage","dir":"Reference","previous_headings":"","what":"Usage","title":"Update value of a linear predictor as function of a single coordinate change — update_linear_predictor","text":"","code":"update_linear_predictor(new_beta_j, current_beta_j, current_eta, X_j)"},{"path":"https://mathiaslj.github.io/mcmcglm/reference/update_linear_predictor.html","id":"arguments","dir":"Reference","previous_headings":"","what":"Arguments","title":"Update value of a linear predictor as function of a single coordinate change — update_linear_predictor","text":"new_beta_j numeric new value j'th component beta parameter vector current_beta_j numeric current value j'th component beta parameter vector current_eta current value linear predictor corresponding current_beta value X_j j'th column design matrix","code":""}] +[{"path":"https://mathiaslj.github.io/mcmcglm/LICENSE.html","id":null,"dir":"","previous_headings":"","what":"MIT License","title":"MIT License","text":"Copyright (c) 2024 mcmcglm authors Permission hereby granted, free charge, person obtaining copy software associated documentation files (“Software”), deal Software without restriction, including without limitation rights use, copy, modify, merge, publish, distribute, sublicense, /sell copies Software, permit persons Software furnished , subject following conditions: copyright notice permission notice shall included copies substantial portions Software. SOFTWARE PROVIDED “”, WITHOUT WARRANTY KIND, EXPRESS IMPLIED, INCLUDING LIMITED WARRANTIES MERCHANTABILITY, FITNESS PARTICULAR PURPOSE NONINFRINGEMENT. EVENT SHALL AUTHORS COPYRIGHT HOLDERS LIABLE CLAIM, DAMAGES LIABILITY, WHETHER ACTION CONTRACT, TORT OTHERWISE, ARISING , CONNECTION SOFTWARE USE DEALINGS SOFTWARE.","code":""},{"path":"https://mathiaslj.github.io/mcmcglm/articles/customising.html","id":"adding-methods-for-new-families","dir":"Articles","previous_headings":"","what":"Adding methods for new families","title":"Adding methods for families","text":", easy extend package include family existing package. missing piece families method calculating log-density distribution. Since calculation log-density family implemented using S3 methods, user needs define function called log_density.*family_name* family_name specified relevant family. example showcasing single line code user needs write enable use inverse gaussian family","code":""},{"path":"https://mathiaslj.github.io/mcmcglm/articles/customising.html","id":"example-of-implementing-method-for-inverse-gaussian-family","dir":"Articles","previous_headings":"Adding methods for new families","what":"Example of implementing method for inverse gaussian family","title":"Adding methods for families","text":"First, generate data showcase example. user needs specify method computing density inverse gaussian distribution like : mu modelled mean μ=g−1(Xβ)\\mu=g^{-1}(X\\beta) GLM model, Y response variable. , mcmcglm() can called like : specify inverse.gaussian() family argument, make sure add additional arguments density log_likelihood_extra_args argument.","code":"n <- 1000 x1 <- rexp(n, 2) x2 <- rbinom(n, 1, .5) b0 <- 1 b1 <- 1.5 b2 <- 2 lin_pred <- b0+b1*x1+b2*x2 invgauss_fam <- inverse.gaussian() y_invnorm <- statmod::rinvgauss(n, mean = invgauss_fam$linkinv(lin_pred), shape = 1, dispersion = 1) dat_invnorm <- data.frame(Y = y_invnorm, X1 = x1, X2 = x2) log_density.inverse.gaussian <- function(family, mu, Y, ...) { statmod::dinvgauss(Y, mean = mu, ..., log = T) } invnorm <- mcmcglm::mcmcglm(formula = Y ~ ., family = inverse.gaussian(), log_likelihood_extra_args = list(shape = 1, dispersion = 1), data = dat_invnorm, beta_prior = distributional::dist_gamma(1, 1), w = 0.5) trace_plot(invnorm)"},{"path":[]},{"path":"https://mathiaslj.github.io/mcmcglm/articles/performance.html","id":"using-default-slice-function","dir":"Articles","previous_headings":"Graph of compute time as function of dimension","what":"Using default slice function","title":"Performance of package","text":"run different slice widths interested","code":"res <- compare_eta_comptime_across_nvars( n_vars = c(2, seq(from = 40, to = 400, by = 40)), n_samples = 1, burnin = 0, w = 0.5) plot_eta_comptime(res, facet_by = \"w\") ws <- c(0.5, 10) res <- lapply(ws, function(w) { compare_eta_comptime_across_nvars( w = w, n_vars = c(2, seq(from = 10, to = 50, by = 10)), n_samples = 50, burnin = 1, parallelise = TRUE) }) %>% dplyr::bind_rows()"},{"path":[]},{"path":"https://mathiaslj.github.io/mcmcglm/articles/pospkg.html","id":"family-of-response","dir":"Articles","previous_headings":"","what":"Family of response","title":"Examples of family, prior and slice sampler specifications","text":"family response can family implemented family class, principle (see help(\"family\")). now, implemented families package Gaussian Binomial Poisson Negative binomial However, small amount work needed user enable use family. process expanding available families described vignette(\"customising\"). Due implementation using family object, also means compatible link function can used family, like fx. logit probit binomial family. Also, unexported check_family function package ensures accepted values family argument package’s functions character, family function family class object. different uses seen throughout examples vignette also available example section mcmcglm() documentation.","code":""},{"path":"https://mathiaslj.github.io/mcmcglm/articles/pospkg.html","id":"data-simulation-for-examples","dir":"Articles","previous_headings":"Family of response","what":"Data simulation for examples","title":"Examples of family, prior and slice sampler specifications","text":"order show package provides samples target posterior distribution, simulate data scenario known β\\beta vector order see empirical distribution coefficient seems match value simulated data . start generating linear predictor η=Xβ\\eta=X\\beta known values β\\beta sampled values making design matrix XX. , example can apply relevant inverse link function g−1g^{-1} obtain modelled mean μ=g−1(Xβ)\\mu=g^{-1}(X\\beta) GLM model. can simulate response sampling distribution corresponding family using modelled mean μ\\mu relevant parameter distribution. example gaussian family available introductory example README.","code":"n <- 1000 x1 <- rnorm (n) x2 <- rbinom (n, 1, .5) b0 <- 1 b1 <- 1.5 b2 <- 2 lin_pred <- b0+b1*x1+b2*x2"},{"path":"https://mathiaslj.github.io/mcmcglm/articles/pospkg.html","id":"example-for-binomial-family","dir":"Articles","previous_headings":"Family of response","what":"Example for binomial family","title":"Examples of family, prior and slice sampler specifications","text":"binomial family, commonly used link function logit given logit(p)=lnp1−p\\begin{align} logit(p)=\\ln\\frac{p}{1-p} \\end{align} link functions also available, including probit, cauchit, log cloglog (see help(family) information). ’ll show example first binomial family logit link function.","code":""},{"path":"https://mathiaslj.github.io/mcmcglm/articles/pospkg.html","id":"logit-link-function","dir":"Articles","previous_headings":"Family of response > Example for binomial family","what":"Logit link function","title":"Examples of family, prior and slice sampler specifications","text":"First, define data fit mcmcglm inspect trace plot.","code":"logit_binom <- binomial(link = \"logit\") mu_logit <- logit_binom$linkinv(lin_pred) y_logit <- rbinom(n, size = 1, prob = mu_logit) dat_logit <- data.frame(Y = y_logit, X1 = x1, X2 = x2) logit <- mcmcglm(formula = Y ~ ., data = dat_logit, family = binomial(link = \"logit\"), beta_prior = distributional::dist_normal(mean = 0, sd = 1), w = 0.5) logit #> Object of class 'mcmcglm' #> #> Call: mcmcglm(formula = Y ~ ., family = binomial(link = \"logit\"), data = dat_logit, #> beta_prior = distributional::dist_normal(mean = 0, sd = 1), #> w = 0.5) #> #> Average of parameter samples: #> (Intercept) X1 X2 #> 1 0.8802159 1.461209 2.061524 trace_plot(logit)"},{"path":"https://mathiaslj.github.io/mcmcglm/articles/pospkg.html","id":"probit-link-function","dir":"Articles","previous_headings":"Family of response > Example for binomial family","what":"Probit link function","title":"Examples of family, prior and slice sampler specifications","text":"First, define data fit mcmcglm inspect trace plot.","code":"probit_binom <- binomial(link = \"probit\") mu_probit <- probit_binom$linkinv(lin_pred) y_probit <- rbinom(n, size = 1, prob = mu_probit) dat_probit <- data.frame(Y = y_probit, X1 = x1, X2 = x2) probit <- mcmcglm(formula = Y ~ ., data = dat_probit, family = binomial(link = \"probit\"), beta_prior = distributional::dist_normal(mean = 0, sd = 1), w = 0.5) probit #> Object of class 'mcmcglm' #> #> Call: mcmcglm(formula = Y ~ ., family = binomial(link = \"probit\"), #> data = dat_probit, beta_prior = distributional::dist_normal(mean = 0, #> sd = 1), w = 0.5) #> #> Average of parameter samples: #> (Intercept) X1 X2 #> 1 0.8318177 1.374661 2.162339 trace_plot(probit)"},{"path":"https://mathiaslj.github.io/mcmcglm/articles/pospkg.html","id":"example-for-poisson-family","dir":"Articles","previous_headings":"Family of response","what":"Example for poisson family","title":"Examples of family, prior and slice sampler specifications","text":"poisson family, commonly used link function log, one showcase example binomial family shows ’s easily possible specify link function. Define data fit mcmcglm inspect trace plot.","code":"mu_log <- exp(lin_pred) y_pois <- rpois(n, lambda = mu_log) dat_pois <- data.frame(Y = y_pois, X1 = x1, X2 = x2) pois <- mcmcglm(formula = Y ~ ., data = dat_pois, family = \"poisson\", beta_prior = distributional::dist_normal(mean = 0, sd = 1), w = 0.5) pois #> Object of class 'mcmcglm' #> #> Call: mcmcglm(formula = Y ~ ., family = \"poisson\", data = dat_pois, #> beta_prior = distributional::dist_normal(mean = 0, sd = 1), #> w = 0.5) #> #> Average of parameter samples: #> (Intercept) X1 X2 #> 1 0.9969416 1.504678 1.997051 trace_plot(pois)"},{"path":"https://mathiaslj.github.io/mcmcglm/articles/pospkg.html","id":"example-for-negative-binomial-family","dir":"Articles","previous_headings":"Family of response","what":"Example for negative binomial family","title":"Examples of family, prior and slice sampler specifications","text":"Besides families available stats package (via help(\"family\")), families implemented various packages. Fx. MASS::negative.binomial() function allows use negative binomial family various inference procedures. distribution interesting choice modelling count data Poisson distribution appropriate due equality first second moment distribution. showcase use family using package. fit mcmcglm inspect trace plot.","code":"mu_log <- exp(lin_pred) y_nbinom <- rnbinom(n, size = 1, mu = mu_log) dat_nbinom <- data.frame(Y = y_nbinom, X1 = x1, X2 = x2) theta <- 3 nbinom <- mcmcglm(formula = Y ~ ., data = dat_nbinom, family = MASS::negative.binomial(theta), beta_prior = distributional::dist_normal(mean = 0, sd = 1), w = 0.5) nbinom #> Object of class 'mcmcglm' #> #> Call: mcmcglm(formula = Y ~ ., family = MASS::negative.binomial(theta), #> data = dat_nbinom, beta_prior = distributional::dist_normal(mean = 0, #> sd = 1), w = 0.5) #> #> Average of parameter samples: #> (Intercept) X1 X2 #> 1 0.9608508 1.431976 2.059662 trace_plot(nbinom)"},{"path":"https://mathiaslj.github.io/mcmcglm/articles/pospkg.html","id":"specifying-the-prior-distribution","dir":"Articles","previous_headings":"","what":"Specifying the prior distribution","title":"Examples of family, prior and slice sampler specifications","text":"package uses distributional package specify prior distribution β\\beta parameter mentioned documentation mcmcglm(). distributional package provides vectorised distribution objects specified parameters methods evaluating densities, CDFs, etc. generate random samples distribution. use package allows user great flexibility specifying varying independent univariate distributions multivariate distributions parameter vector illustrated examples . examples simply use gaussian family response introductory example README alternate prior specifications. Start defining data.","code":"y_norm <- rnorm(n, mean = lin_pred, sd = 1) dat_norm <- data.frame(Y = y_norm, X1 = x1, X2 = x2)"},{"path":"https://mathiaslj.github.io/mcmcglm/articles/pospkg.html","id":"example-of-iid-priors","dir":"Articles","previous_headings":"Specifying the prior distribution","what":"Example of iid priors","title":"Examples of family, prior and slice sampler specifications","text":"case iid priors, single univariate distribution specification needed: make example exact introductory example README, show posterior distribution affected prior mean different data(/likelihood) showing. course, distribution used using dist_ functions distributional package. , simply use normal distribution distributional::dist_normal(), might well distributional::dist_gamma(), distributional::dist_student_t(), etc. See Reference full list distributions.","code":"iid <- mcmcglm(formula = Y ~ ., family = \"gaussian\", data = dat_norm, beta_prior = distributional::dist_normal(1000, 1), w = 0.5) trace_plot(iid)"},{"path":"https://mathiaslj.github.io/mcmcglm/articles/pospkg.html","id":"example-of-independent-and-potentially-different-priors","dir":"Articles","previous_headings":"Specifying the prior distribution","what":"Example of independent and potentially different priors","title":"Examples of family, prior and slice sampler specifications","text":"Due S3 method implementation mcmcglm package, ’s possible specify beta_prior list distribution objects. illustrative purposes, use distributions mentioned example specify different prior β\\beta component.","code":"indep <- mcmcglm(formula = Y ~ ., family = \"gaussian\", data = dat_norm, beta_prior = list( distributional::dist_normal(mean = 0, sd = 1), distributional::dist_gamma(shape = 1, rate = 2), distributional::dist_student_t(df = 8, mu = 3, sigma = 1) ), w = 0.5) trace_plot(indep)"},{"path":"https://mathiaslj.github.io/mcmcglm/articles/pospkg.html","id":"example-of-multivariate-prior","dir":"Articles","previous_headings":"Specifying the prior distribution","what":"Example of multivariate prior","title":"Examples of family, prior and slice sampler specifications","text":"implemented multivariate distribution distributional::dist_multivariate_normal(), can use specify different moments potential correlation parameters model. example shown prior distribution specified β∼𝒩3((123),[10.50.20.510.50.20.51])\\begin{align} \\beta\\sim \\mathcal{N}_3\\left( \\begin{pmatrix} 1\\\\2\\\\3 \\end{pmatrix}, \\begin{bmatrix} 1 & 0.5 & 0.2\\\\ 0.5 & 1 & 0.5\\\\ 0.2 & 0.5 & 1 \\end{bmatrix} \\right) \\end{align}","code":"mult <- mcmcglm(formula = Y ~ ., family = \"gaussian\", data = dat_norm, beta_prior = distributional::dist_multivariate_normal( mu = list(1:3), sigma = list( matrix( c(1, 0.5, 0.2, 0.5, 1, 0.5, 0.2, 0.5, 1), ncol = 3, byrow = TRUE)) ), w = 0.5) trace_plot(mult)"},{"path":"https://mathiaslj.github.io/mcmcglm/articles/pospkg.html","id":"tuning-and-changing-the-slice-sampler-procedure","dir":"Articles","previous_headings":"","what":"Tuning (and changing) the slice sampler procedure","title":"Examples of family, prior and slice sampler specifications","text":"sampling procedure package uses compute graph gibbs (CGGibbs) sampler introduced Luu et al. 2024. done efficiently “update” linear predictor gibbs update single component parameter vector rather naively recalculating linear predictor using terms.","code":""},{"path":"https://mathiaslj.github.io/mcmcglm/articles/pospkg.html","id":"advantage-of-cggibbs","dir":"Articles","previous_headings":"Tuning (and changing) the slice sampler procedure","what":"Advantage of CGGibbs","title":"Examples of family, prior and slice sampler specifications","text":"See details article, ’ll quickly show simple elegant idea . Assuming kk’th iteration (KK total iterations) sampled new jj’th coordinate dd-dimensional parameter vector β\\beta, βj(k)\\beta_j^{(k)}, “current” value β(j)(k)=(β1(k),β2(k),...,βj(k),βj+1(k−1),...,βd(k−1))\\begin{align} \\beta_{(j)}^{(k)}=\\left(\\beta_1^{(k)}, \\beta_2^{(k)}, ..., \\beta_j^{(k)}, \\beta_{j+1}^{(k-1)}, ..., \\beta_{d}^{(k-1)}\\right) \\end{align} CGGibbs approach calculate linear predictor kk’th iteration component jj, η(j)(k)\\eta_{(j)}^{(k)}, η(j)(k)=η(j−1)(k)+∑=1nxij(βj(k)−βj(k−1))\\begin{align} \\eta_{(j)}^{(k)}=\\eta_{(j-1)}^{(k)}+\\sum_{=1}^n x_{ij} (\\beta_j^{(k)}-\\beta_j^{(k-1)}) \\end{align} naive approach ηj(k)=∑=1n∑l=1dxilβl(k−𝟏[j>l])\\begin{align} \\eta_j^{(k)}=\\sum_{=1}^n\\sum_{l=1}^d x_{il} \\beta_l^{(k-\\mathbf{1}[j>l])} \\end{align} means using CGGibbs, component update β\\beta, computation time O(n)O(n), compared O(dn)O(dn) naive approach.","code":""},{"path":"https://mathiaslj.github.io/mcmcglm/articles/pospkg.html","id":"advantages-of-slice-sampling-in-combination-with-cggibbs","dir":"Articles","previous_headings":"Tuning (and changing) the slice sampler procedure","what":"Advantages of slice sampling (in combination with CGGibbs)","title":"Examples of family, prior and slice sampler specifications","text":"linear predictor used calculate log-likelihood, summed log-density prior distribution obtain log-potential, log-posterior normalising constant. log-potential calculated algorithm order use slice sampling, enables us use combination family response prior distribution parameters rely conjugacy get closed-form distributions sample . , slice sampling avoids spending computational power numerical integration find normalising constant corresponding marginal distribution data sampling procedure needs know posterior proportionality; reason use log-potential. default method slice sampling “stepping-” “shrinkage” described Neal 2003.","code":""},{"path":"https://mathiaslj.github.io/mcmcglm/articles/pospkg.html","id":"using-the-default-slice-sampling-method","dir":"Articles","previous_headings":"Tuning (and changing) the slice sampler procedure","what":"Using the default slice sampling method","title":"Examples of family, prior and slice sampler specifications","text":"slice_ functions qslice package arguments x log_target. mcmcglm function finds entities new component β\\beta vector given iteration log-potential evaluated function new value parameter. Thus, task left user specify tuning parameters chosen algorithm, passed ... mcmcglm(). default slice sampler uses qslice::slice_stepping_out() function, mandatory non-default argument w slice width algorithm. Thus, readers might noticed examples within vignette README, specification w = x given calls mcmcglm(). Example completeness:","code":"norm <- mcmcglm(formula = Y ~ ., family = \"gaussian\", data = dat_norm, beta_prior = distributional::dist_normal(0, 1), w = 0.5)"},{"path":"https://mathiaslj.github.io/mcmcglm/articles/pospkg.html","id":"changing-the-default-slice-sampling-method","dir":"Articles","previous_headings":"Tuning (and changing) the slice sampler procedure","what":"Changing the default slice sampling method","title":"Examples of family, prior and slice sampler specifications","text":"argument qslice_fun mcmcglm() function allows user change default slice sampler qslice::slice_stepping_out() slice sampler qslice package. Thus, example using qslice::slice_elliptical() Murray et al. (2010), tuning parameters mu sigma:","code":"el <- mcmcglm(formula = Y ~ ., family = \"gaussian\", data = dat_norm, beta_prior = distributional::dist_normal(0, 1), qslice_fun = qslice::slice_elliptical, mu = 3, sigma = 2) trace_plot(el)"},{"path":"https://mathiaslj.github.io/mcmcglm/articles/pospkg.html","id":"investigating-effect-of-tuning-parameters-of-slice-sampler","dir":"Articles","previous_headings":"Tuning (and changing) the slice sampler procedure","what":"Investigating effect of tuning parameters of slice sampler","title":"Examples of family, prior and slice sampler specifications","text":"utility functions mcmcglm_across_tuningparams() plot_mcmcglm_across_tuningparams() can used create trace plots samples running algorithm different values tuning parameters.","code":""},{"path":"https://mathiaslj.github.io/mcmcglm/articles/pospkg.html","id":"varying-slice-width-for-the-default-slice-sampler","dir":"Articles","previous_headings":"Tuning (and changing) the slice sampler procedure > Investigating effect of tuning parameters of slice sampler","what":"Varying slice width for the default slice sampler","title":"Examples of family, prior and slice sampler specifications","text":"Fx. default choice slice sampler, qslice::slice_stepping_out(), tuning parameter slice width w. way create plot samples across different values slice width seen :","code":"w05_mcmcglms <- mcmcglm_across_tuningparams( seq(from = 0.5, by = 0.5, length.out = 4), tuning_parameter_name = \"w\", formula = Y ~ ., family = \"gaussian\", data = dat_norm, parallelise = FALSE) plot_mcmcglm_across_tuningparams(w05_mcmcglms)"},{"path":"https://mathiaslj.github.io/mcmcglm/articles/pospkg.html","id":"varying-degrees-of-freedom-in-a-generalized-elliptical-slice-sampler","dir":"Articles","previous_headings":"Tuning (and changing) the slice sampler procedure > Investigating effect of tuning parameters of slice sampler","what":"Varying degrees of freedom in a generalized elliptical slice sampler","title":"Examples of family, prior and slice sampler specifications","text":"Since utility function work generally single tuning parameter, example qslice::slice_elliptical sampling procedure described Nishihara et al. (2014), varying df parameter algorithm. Note single tuning parameter can “varied” time. example, keep sigma mu constant, vary df","code":"df1_mcmcglms <- mcmcglm_across_tuningparams( c(1, seq(from = 10, by = 10, length.out = 3)), tuning_parameter_name = \"df\", formula = Y ~ ., family = \"gaussian\", data = dat_norm, qslice_fun = qslice::slice_genelliptical, mu = 2, sigma = 1, parallelise = FALSE) plot_mcmcglm_across_tuningparams(df1_mcmcglms)"},{"path":"https://mathiaslj.github.io/mcmcglm/articles/pospkg.html","id":"normal-normal-sample-method","dir":"Articles","previous_headings":"Tuning (and changing) the slice sampler procedure","what":"Normal-normal sample method","title":"Examples of family, prior and slice sampler specifications","text":"Mostly testing purposes, package implementation sampling closed-form normal distribution case response prior normally distributed. Running requires specifying \"normal-normal\" argument sample_method mcmcglm() seen example :","code":"norm_norm <- mcmcglm(formula = Y ~ ., family = \"gaussian\", data = dat_norm, beta_prior = distributional::dist_normal(0, 1), sample_method = \"normal-normal\") trace_plot(norm_norm)"},{"path":"https://mathiaslj.github.io/mcmcglm/authors.html","id":null,"dir":"","previous_headings":"","what":"Authors","title":"Authors and Citation","text":"Mathias Lerbech Jeppesen. Author, maintainer. Julie Køhler Munkvad. Author.","code":""},{"path":"https://mathiaslj.github.io/mcmcglm/authors.html","id":"citation","dir":"","previous_headings":"","what":"Citation","title":"Authors and Citation","text":"Jeppesen M, Munkvad J (2024). mcmcglm: Gibbs Sampling Methods. R package version 0.0.0.9000, https://mathiaslj.github.io/mcmcglm/, https://github.com/mathiaslj/mcmcglm.","code":"@Manual{, title = {mcmcglm: Gibbs Sampling Methods}, author = {Mathias Lerbech Jeppesen and Julie Køhler Munkvad}, year = {2024}, note = {R package version 0.0.0.9000, https://mathiaslj.github.io/mcmcglm/}, url = {https://github.com/mathiaslj/mcmcglm}, }"},{"path":"https://mathiaslj.github.io/mcmcglm/index.html","id":"mcmcglm","dir":"","previous_headings":"","what":"Gibbs Sampling Methods","title":"Gibbs Sampling Methods","text":"mcmcglm package implements CGGibbs sampler article Gibbs sampling faster Hamiltonian Monte Carlo GLMs?, linear run time function number parameters GLM model due clever “update” linear predictor. See details section vignette package implemented way user can specify family response distribution prior β\\beta parameter, 𝔼[Y|Xβ]=g−1(Xβ)\\mathbb{E}[Y|X\\beta]=g^{-1}(X\\beta) link function gg specified family. See vignette(\"pospkg\").","code":""},{"path":"https://mathiaslj.github.io/mcmcglm/index.html","id":"installation","dir":"","previous_headings":"","what":"Installation","title":"Gibbs Sampling Methods","text":"Installation available GitHub :","code":"devtools::install_github(\"mathiaslj/mcmcglm\")"},{"path":"https://mathiaslj.github.io/mcmcglm/index.html","id":"example","dir":"","previous_headings":"","what":"Example","title":"Gibbs Sampling Methods","text":"first simulate data linear model use showcasing use mcmcglm gaussian family. use function mcmcglm similar interface glm function added mandatory specification prior parameter β\\beta default qslice::slice_stepping_out function, performs slice sampling described Neal 2003 slice width w needs specified. creates mcmcglm object prints summarising call function averages samples parameter GLM model.","code":"n <- 1000 x1 <- rnorm (n) x2 <- rbinom (n, 1, .5) b0 <- 1 b1 <- 1.5 b2 <- 2 lin_pred <- b0+b1*x1+b2*x2 y_norm <- rnorm(n, mean = lin_pred, sd = 1) dat_norm <- data.frame(Y = y_norm, X1 = x1, X2 = x2) norm <- mcmcglm(formula = Y ~ ., family = \"gaussian\", data = dat_norm, beta_prior = distributional::dist_normal(0, 1), w = 0.5) norm #> Object of class 'mcmcglm' #> #> Call: mcmcglm(formula = Y ~ ., family = \"gaussian\", data = dat_norm, #> beta_prior = distributional::dist_normal(0, 1), w = 0.5) #> #> Average of parameter samples: #> (Intercept) X1 X2 #> 1 1.011134 1.490459 2.026047"},{"path":"https://mathiaslj.github.io/mcmcglm/index.html","id":"investigating-results","dir":"","previous_headings":"Example","what":"Investigating results","title":"Gibbs Sampling Methods","text":"averages shown print method object can retrieved generic coef like : full data set samples can accessed samples function: trace plot can seen function trace_plot:","code":"coef(norm) #> (Intercept) X1 X2 #> 1 1.011134 1.490459 2.026047 head(samples(norm)) #> (Intercept) X1 X2 iteration burnin #> 1 0.6173367 -0.004541141 -0.09125636 0 TRUE #> 2 2.6508146 0.281295470 0.68343132 1 TRUE #> 3 0.8240996 0.324627659 2.30889073 2 TRUE #> 4 0.8170086 1.028326905 2.20351455 3 TRUE #> 5 0.8777350 1.592074284 2.16115289 4 TRUE #> 6 0.9092187 1.442872350 2.02913214 5 TRUE trace_plot(norm)"},{"path":"https://mathiaslj.github.io/mcmcglm/reference/coef.mcmcglm.html","id":null,"dir":"Reference","previous_headings":"","what":"S3 method for getting the average value of coefficients — coef.mcmcglm","title":"S3 method for getting the average value of coefficients — coef.mcmcglm","text":"S3 method getting average value coefficients","code":""},{"path":"https://mathiaslj.github.io/mcmcglm/reference/coef.mcmcglm.html","id":"ref-usage","dir":"Reference","previous_headings":"","what":"Usage","title":"S3 method for getting the average value of coefficients — coef.mcmcglm","text":"","code":"# S3 method for class 'mcmcglm' coef(x)"},{"path":"https://mathiaslj.github.io/mcmcglm/reference/coef.mcmcglm.html","id":"arguments","dir":"Reference","previous_headings":"","what":"Arguments","title":"S3 method for getting the average value of coefficients — coef.mcmcglm","text":"x mcmcglm object","code":""},{"path":"https://mathiaslj.github.io/mcmcglm/reference/compare_eta_comptime_across_nvars.html","id":null,"dir":"Reference","previous_headings":"","what":"Compare runtime using CGGibbs and naive approach to calculate linear predictor — compare_eta_comptime_across_nvars","title":"Compare runtime using CGGibbs and naive approach to calculate linear predictor — compare_eta_comptime_across_nvars","text":"comparison methods described vignette. user can specify different arguments alter specification comparison, possible specify values n_vars argument. possible change fact data simulated normal response normal explanatory variables.","code":""},{"path":"https://mathiaslj.github.io/mcmcglm/reference/compare_eta_comptime_across_nvars.html","id":"ref-usage","dir":"Reference","previous_headings":"","what":"Usage","title":"Compare runtime using CGGibbs and naive approach to calculate linear predictor — compare_eta_comptime_across_nvars","text":"","code":"compare_eta_comptime_across_nvars( n_vars, n = 100, beta_prior = distributional::dist_normal(0, 1), log_likelihood_extra_args = list(sd = 1), sample_method = c(\"slice_sampling\", \"normal-normal\"), qslice_fun = qslice::slice_stepping_out, ..., n_samples = 500, burnin = 100, parallelise = FALSE, n_cores = as.numeric(Sys.getenv(\"NUMBER_OF_PROCESSORS\")) - 1 )"},{"path":"https://mathiaslj.github.io/mcmcglm/reference/compare_eta_comptime_across_nvars.html","id":"arguments","dir":"Reference","previous_headings":"","what":"Arguments","title":"Compare runtime using CGGibbs and naive approach to calculate linear predictor — compare_eta_comptime_across_nvars","text":"n_vars numeric vector numbers normally distributed columns design matrix different runs n numeric number observations data analysed GLM beta_prior distribution object created function distributional package. fx. distributional::dist_normal(mean = 0, sd = 1). log_likelihood_extra_args named list arguments passed onto log_density function. Fx. specification log_likelihood_extra_args = list(sd = x) needed case family = \"gaussian\" sample_method character specifying method used sampling. default \"slice_sampling\" intended value cases, works specification family beta_prior. \"normal-normal\" uses conditional normal distribution sample case conjugate prior gaussian response beta_prior. Implemented testing purposes works niche case. qslice_fun function qslice package. Default qslice::slice_stepping_out uses slice sampler Neal 2003, functions available. ... arguments passed onto function specified qslice_fun. default qslice::slice_stepping_out w needs specified, fx. qslice::slice_elliptical, mu sigma need specified n_samples numeric number samples draw parameter(/variable) model burnin numeric number samples marked \"burnin\". Burnin samples included beta_mean calculation increase finite sample performance LLN estimate parallelise logical runs algorithm across values n_vars parallelised using future.apply::future_lapply n_cores numeric number cores use parallelisation. Default 1 less number available cores","code":""},{"path":"https://mathiaslj.github.io/mcmcglm/reference/compare_eta_comptime_across_nvars.html","id":"value","dir":"Reference","previous_headings":"","what":"Value","title":"Compare runtime using CGGibbs and naive approach to calculate linear predictor — compare_eta_comptime_across_nvars","text":"data.frame information computation time different values n_vars","code":""},{"path":"https://mathiaslj.github.io/mcmcglm/reference/compare_eta_comptime_across_nvars.html","id":"ref-examples","dir":"Reference","previous_headings":"","what":"Examples","title":"Compare runtime using CGGibbs and naive approach to calculate linear predictor — compare_eta_comptime_across_nvars","text":"","code":"# Compare the runtime for 2, 20, 40, 60 variables in the model compare_eta_comptime_across_nvars(n_vars = c(2, seq(from = 20, to = 60, by = 20)), n_samples = 100, burnin = 20) #> Sampling from posterior ■■■■■■■ 21% | ETA: 4s #> Sampling from posterior ■■■■■■■■■■■■■■■■■■■■■■■■■■■ 87% | ETA: 0s #> Sampling from posterior ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ 100% | ETA: 0s #> Sampling from posterior ■■■■■■■■■■■■■■■■■■■■ 63% | ETA: 1s #> Sampling from posterior ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ 100% | ETA: 0s #> Sampling from posterior ■■■ 6% | ETA: 27s #> Sampling from posterior ■■■■■■■■■■■■■■■ 45% | ETA: 6s #> Sampling from posterior ■■■■■■■■■■■■■■■■■■■■■■■■■■■ 85% | ETA: 1s #> Sampling from posterior ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ 100% | ETA: 0s #> Sampling from posterior ■■■■ 9% | ETA: 19s #> Sampling from posterior ■■■■■■■■■■■■■■■ 45% | ETA: 6s #> Sampling from posterior ■■■■■■■■■■■■■■■■■■■■■■■■■■ 84% | ETA: 1s #> Sampling from posterior ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ 100% | ETA: 0s #> Sampling from posterior ■ 1% | ETA: 3m #> Sampling from posterior ■■■■ 10% | ETA: 44s #> Sampling from posterior ■■■■■■■■■■ 29% | ETA: 19s #> Sampling from posterior ■■■■■■■■■■■■■■■■ 51% | ETA: 10s #> Sampling from posterior ■■■■■■■■■■■■■■■■■■■■■■■ 73% | ETA: 5s #> Sampling from posterior ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ 97% | ETA: 1s #> Sampling from posterior ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ 100% | ETA: 0s #> Sampling from posterior ■■■ 6% | ETA: 39s #> Sampling from posterior ■■■■■■■■■■ 29% | ETA: 14s #> Sampling from posterior ■■■■■■■■■■■■■■■■■ 53% | ETA: 8s #> Sampling from posterior ■■■■■■■■■■■■■■■■■■■■■■■■ 78% | ETA: 3s #> Sampling from posterior ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ 100% | ETA: 0s #> time linear_predictor_calc n_vars n_samples beta_mean #> 1 0.3153636 secs update 2 100 0 #> 2 0.2818551 secs naive 2 100 0 #> 3 3.7287543 secs update 20 100 0 #> 4 3.7928767 secs naive 20 100 0 #> 5 8.8522904 secs update 40 100 0 #> 6 9.1041164 secs naive 40 100 0 #> 7 17.2623024 secs update 60 100 0 #> 8 14.2019286 secs naive 60 100 0 #> beta_variance family sd qslice_fun w parallelised #> 1 1 gaussian 1 qslice::slice_stepping_out 0.5 FALSE #> 2 1 gaussian 1 qslice::slice_stepping_out 0.5 FALSE #> 3 1 gaussian 1 qslice::slice_stepping_out 0.5 FALSE #> 4 1 gaussian 1 qslice::slice_stepping_out 0.5 FALSE #> 5 1 gaussian 1 qslice::slice_stepping_out 0.5 FALSE #> 6 1 gaussian 1 qslice::slice_stepping_out 0.5 FALSE #> 7 1 gaussian 1 qslice::slice_stepping_out 0.5 FALSE #> 8 1 gaussian 1 qslice::slice_stepping_out 0.5 FALSE"},{"path":"https://mathiaslj.github.io/mcmcglm/reference/log_density.html","id":null,"dir":"Reference","previous_headings":"","what":"S3 generic for calculating the log density of a distribution dispatched via a family — log_density","title":"S3 generic for calculating the log density of a distribution dispatched via a family — log_density","text":"methods parametrised \"mu\" takes additional arguments needed calculation ... argument","code":""},{"path":"https://mathiaslj.github.io/mcmcglm/reference/log_density.html","id":"ref-usage","dir":"Reference","previous_headings":"","what":"Usage","title":"S3 generic for calculating the log density of a distribution dispatched via a family — log_density","text":"","code":"log_density(family, mu, Y, ...)"},{"path":"https://mathiaslj.github.io/mcmcglm/reference/log_density.html","id":"arguments","dir":"Reference","previous_headings":"","what":"Arguments","title":"S3 generic for calculating the log density of a distribution dispatched via a family — log_density","text":"family description error distribution link function used model. can character string naming family function, family function result call family function. (See family details family functions.) mu numeric vector values \"main\" parameter distribution specified family argument Y numeric vector response variable evaluate density ... arguments passed relevant methods","code":""},{"path":"https://mathiaslj.github.io/mcmcglm/reference/log_density.html","id":"value","dir":"Reference","previous_headings":"","what":"Value","title":"S3 generic for calculating the log density of a distribution dispatched via a family — log_density","text":"numeric vector log_density values","code":""},{"path":"https://mathiaslj.github.io/mcmcglm/reference/log_density.html","id":"details","dir":"Reference","previous_headings":"","what":"Details","title":"S3 generic for calculating the log density of a distribution dispatched via a family — log_density","text":"Supported families gaussian, binomial, poisson negative.binomial. Implement S3 method add support new family. mcmcglm function called, mu modelled mean glm model (meaning inverse link linear predictor). Reference methods see parametrisation","code":""},{"path":"https://mathiaslj.github.io/mcmcglm/reference/log_likelihood.html","id":null,"dir":"Reference","previous_headings":"","what":"Calculate log likelihood parametrised by ","title":"Calculate log likelihood parametrised by ","text":"Calculate log likelihood parametrised \"mu\"","code":""},{"path":"https://mathiaslj.github.io/mcmcglm/reference/log_likelihood.html","id":"ref-usage","dir":"Reference","previous_headings":"","what":"Usage","title":"Calculate log likelihood parametrised by ","text":"","code":"log_likelihood(family, mu, Y, ...)"},{"path":"https://mathiaslj.github.io/mcmcglm/reference/log_likelihood.html","id":"arguments","dir":"Reference","previous_headings":"","what":"Arguments","title":"Calculate log likelihood parametrised by ","text":"family description error distribution link function used model. can character string naming family function, family function result call family function. (See family details family functions.) mu numeric vector values \"main\" parameter distribution specified family argument Y numeric vector response variable evaluate density ... arguments passed S3 generic log_density","code":""},{"path":"https://mathiaslj.github.io/mcmcglm/reference/log_likelihood.html","id":"value","dir":"Reference","previous_headings":"","what":"Value","title":"Calculate log likelihood parametrised by ","text":"numeric Value log-likelihood","code":""},{"path":"https://mathiaslj.github.io/mcmcglm/reference/log_likelihood.html","id":"details","dir":"Reference","previous_headings":"","what":"Details","title":"Calculate log likelihood parametrised by ","text":"See log_density details","code":""},{"path":"https://mathiaslj.github.io/mcmcglm/reference/log_likelihood.html","id":"ref-examples","dir":"Reference","previous_headings":"","what":"Examples","title":"Calculate log likelihood parametrised by ","text":"","code":"# Create a test data n <- 100 x1 <- rnorm (n) x2 <- rbinom (n, 1, .5) b0 <- 1 b1 <- 1.5 b2 <- 2 lin_pred <- b0+b1*x1+b2*x2 known_sigma <- 1 y_norm <- rnorm(n, mean = lin_pred, sd = known_sigma) model_matrix_norm <- as.matrix( data.frame(int = 1, X1 = x1, X2 = x2)) b_prior <- 1:3 mu <- model_matrix_norm %*% b_prior log_likelihood(family = gaussian, mu = mu, Y = y_norm, sd = known_sigma) #> [1] -186.0887"},{"path":"https://mathiaslj.github.io/mcmcglm/reference/log_potential_from_betaj.html","id":null,"dir":"Reference","previous_headings":"","what":"Calculate the log-potential (log-likelihood plus log-density of prior) — log_potential_from_betaj","title":"Calculate the log-potential (log-likelihood plus log-density of prior) — log_potential_from_betaj","text":"Calculates log-potential function new coordinate beta parameter vector. Done like use unexporte","code":""},{"path":"https://mathiaslj.github.io/mcmcglm/reference/log_potential_from_betaj.html","id":"ref-usage","dir":"Reference","previous_headings":"","what":"Usage","title":"Calculate the log-potential (log-likelihood plus log-density of prior) — log_potential_from_betaj","text":"","code":"log_potential_from_betaj( new_beta_j, j, current_beta, current_eta, Y, X, family, beta_prior, linear_predictor_calc = \"update\", ... )"},{"path":"https://mathiaslj.github.io/mcmcglm/reference/log_potential_from_betaj.html","id":"arguments","dir":"Reference","previous_headings":"","what":"Arguments","title":"Calculate the log-potential (log-likelihood plus log-density of prior) — log_potential_from_betaj","text":"new_beta_j numeric new value j'th component beta parameter vector j numeric index parameter vector current_beta current value beta parameter vector sampling procedure current_eta current value linear predictor corresponding current_beta value Y numeric vector response variable evaluate density X design matrix family description error distribution link function used model. can character string naming family function, family function result call family function. (See family details family functions.) beta_prior distribution object created function distributional package. fx. distributional::dist_normal(mean = 0, sd = 1). linear_predictor_calc character specifying method used calculate linear predictor step gibbs algorithm. Default \"update\", uses CGGibbs procedure described start section vignette. option \"naive\", usual ... arguments passed S3 generic log_density","code":""},{"path":"https://mathiaslj.github.io/mcmcglm/reference/log_potential_from_betaj.html","id":"value","dir":"Reference","previous_headings":"","what":"Value","title":"Calculate the log-potential (log-likelihood plus log-density of prior) — log_potential_from_betaj","text":"value log-potential changed j'th component current_beta new_beta_j","code":""},{"path":"https://mathiaslj.github.io/mcmcglm/reference/log_potential_from_betaj.html","id":"ref-examples","dir":"Reference","previous_headings":"","what":"Examples","title":"Calculate the log-potential (log-likelihood plus log-density of prior) — log_potential_from_betaj","text":"","code":"# Create a test data n <- 100 x1 <- rnorm (n) x2 <- rbinom (n, 1, .5) b0 <- 1 b1 <- 1.5 b2 <- 2 lin_pred <- b0+b1*x1+b2*x2 known_sigma <- 1 y_norm <- rnorm(n, mean = lin_pred, sd = known_sigma) model_matrix_norm <- as.matrix( data.frame(int = 1, X1 = x1, X2 = x2)) b_prior <- distributional::dist_normal(mean = 0, sd = 1) b_prior_init <- distributional::generate( b_prior, ncol(model_matrix_norm) )[[1]] eta_init <- model_matrix_norm %*% b_prior_init j <- 1 new_beta_j <- 4 log_potential_from_betaj(new_beta_j = new_beta_j, j = j, current_beta = b_prior_init, current_eta = eta_init, Y = y_norm, X = model_matrix_norm, family = gaussian, beta_prior = b_prior, sd = known_sigma) #> [1] -535.7085"},{"path":"https://mathiaslj.github.io/mcmcglm/reference/mcmcglm-package.html","id":null,"dir":"Reference","previous_headings":"","what":"mcmcglm: Gibbs Sampling Methods — mcmcglm-package","title":"mcmcglm: Gibbs Sampling Methods — mcmcglm-package","text":"Sample posterior distributions using variations slice within Gibbs samplers GLMs.","code":""},{"path":[]},{"path":"https://mathiaslj.github.io/mcmcglm/reference/mcmcglm-package.html","id":"author","dir":"Reference","previous_headings":"","what":"Author","title":"mcmcglm: Gibbs Sampling Methods — mcmcglm-package","text":"Maintainer: Mathias Lerbech Jeppesen mathiasljeppesen@outlook.com Authors: Julie Køhler Munkvad julie.m@hotmail.dk","code":""},{"path":"https://mathiaslj.github.io/mcmcglm/reference/mcmcglm.html","id":null,"dir":"Reference","previous_headings":"","what":"Efficient Gibbs sampling of posterior distribution of parameters in GLM — mcmcglm","title":"Efficient Gibbs sampling of posterior distribution of parameters in GLM — mcmcglm","text":"Obtain MCMC samples using slice sampling within Gibbs generalized linear models (GLMs) using compute graph Gibbs (CGGibbs) linear runtime number variables model matrix. Method described article Gibbs sampling faster Hamiltonian Monte Carlo GLMs?, see details .","code":""},{"path":"https://mathiaslj.github.io/mcmcglm/reference/mcmcglm.html","id":"ref-usage","dir":"Reference","previous_headings":"","what":"Usage","title":"Efficient Gibbs sampling of posterior distribution of parameters in GLM — mcmcglm","text":"","code":"mcmcglm( formula, family = gaussian, data, beta_prior = distributional::dist_normal(0, 1), log_likelihood_extra_args = list(sd = 1), linear_predictor_calc = c(\"update\", \"naive\"), sample_method = c(\"slice_sampling\", \"normal-normal\"), qslice_fun = qslice::slice_stepping_out, ..., n_samples = 500, burnin = 100 )"},{"path":"https://mathiaslj.github.io/mcmcglm/reference/mcmcglm.html","id":"arguments","dir":"Reference","previous_headings":"","what":"Arguments","title":"Efficient Gibbs sampling of posterior distribution of parameters in GLM — mcmcglm","text":"formula object class \"formula\" (one can coerced class): symbolic description model fitted. See details stats::glm family description error distribution link function used model. can character string naming family function, family function result call family function. (See family details family functions.) data optional data frame, list environment (object coercible .data.frame data frame) containing variables model. found data, variables taken environment(formula), typically environment function called. beta_prior distribution object created function distributional package. fx. distributional::dist_normal(mean = 0, sd = 1). log_likelihood_extra_args named list arguments passed onto log_density function. Fx. specification log_likelihood_extra_args = list(sd = x) needed case family = \"gaussian\" linear_predictor_calc character specifying method used calculate linear predictor step gibbs algorithm. Default \"update\", uses CGGibbs procedure described start section vignette. option \"naive\", usual sample_method character specifying method used sampling. default \"slice_sampling\" intended value cases, works specification family beta_prior. \"normal-normal\" uses conditional normal distribution sample case conjugate prior gaussian response beta_prior. Implemented testing purposes works niche case. qslice_fun function qslice package. Default qslice::slice_stepping_out uses slice sampler Neal 2003, functions available. ... arguments passed onto function specified qslice_fun. default qslice::slice_stepping_out w needs specified, fx. qslice::slice_elliptical, mu sigma need specified n_samples numeric number samples draw parameter(/variable) model burnin numeric number samples marked \"burnin\". Burnin samples included beta_mean calculation increase finite sample performance LLN estimate","code":""},{"path":"https://mathiaslj.github.io/mcmcglm/reference/mcmcglm.html","id":"value","dir":"Reference","previous_headings":"","what":"Value","title":"Efficient Gibbs sampling of posterior distribution of parameters in GLM — mcmcglm","text":"object class mcmcglm methods getting data.frame parameter samples, plotting, etc.","code":""},{"path":"https://mathiaslj.github.io/mcmcglm/reference/mcmcglm.html","id":"details","dir":"Reference","previous_headings":"","what":"Details","title":"Efficient Gibbs sampling of posterior distribution of parameters in GLM — mcmcglm","text":"uses updating scheme linear predictor draw Gibbs sampling coordinates parameter vector","code":""},{"path":"https://mathiaslj.github.io/mcmcglm/reference/mcmcglm.html","id":"references","dir":"Reference","previous_headings":"","what":"References","title":"Efficient Gibbs sampling of posterior distribution of parameters in GLM — mcmcglm","text":"Gibbs sampling faster Hamiltonian Monte Carlo GLMs?, Neal 2003","code":""},{"path":"https://mathiaslj.github.io/mcmcglm/reference/mcmcglm.html","id":"ref-examples","dir":"Reference","previous_headings":"","what":"Examples","title":"Efficient Gibbs sampling of posterior distribution of parameters in GLM — mcmcglm","text":"","code":"if (FALSE) { # \\dontrun{ # Create test data for different scenarios n <- 100 x1 <- rnorm (n) x2 <- rbinom (n, 1, .5) b0 <- 1 b1 <- 1.5 b2 <- 2 lin_pred <- b0+b1*x1+b2*x2 ############################################# # Different families and priors # For family \"gaussian\" and iid normal prior y_norm <- rnorm(n, mean = lin_pred, sd = 1) dat_norm <- data.frame(Y = y_norm, X1 = x1, X2 = x2) norm <- mcmcglm(formula = Y ~ ., data = dat_norm, beta_prior = distributional::dist_normal(0, 1), family = \"gaussian\", n_samples = 100, burnin = 10, sample_method = \"slice_sampling\", qslice_fun = qslice::slice_stepping_out, w = 0.5) norm # For family \"binomial\" with logit link and iid gamma distributed prior y_logit <- rbinom (n, 1, arm::invlogit(lin_pred)) dat_logit <- data.frame(Y = y_logit, X1 = x1, X2 = x2) logit <- mcmcglm(formula = Y ~ ., data = dat_logit, beta_prior = distributional::dist_gamma(shape = 1, rate = 0.5), family = binomial(link = \"logit\"), n_samples = 100, burnin = 10, sample_method = \"slice_sampling\", qslice_fun = qslice::slice_stepping_out, w = 0.8) logit # For family \"negative.binomial\" and multivariate normal specification of parameter priors y_log <- rnbinom(n, size = 1, mu = exp(lin_pred)) dat_log <- data.frame(Y = y_log, X1 = x1, X2 = x2) log <- mcmcglm(formula = Y ~ X1, data = dat_log, beta_prior = distributional::dist_multivariate_normal( mu = list(c(1, 2)), sigma = list(matrix(c(1, 0.5, 0.5, 1), ncol = 2)) ), family = MASS::negative.binomial(3), n_samples = 100, burnin = 10, sample_method = \"slice_sampling\", qslice_fun = qslice::slice_stepping_out, w = 0.8) log # For family \"negative.binomial\" and specification of different independent # priors for each parameter log2 <- mcmcglm(formula = Y ~ ., data = dat_log, beta_prior = list(distributional::dist_normal(0, 1), distributional::dist_gamma(1, 1), distributional::dist_exponential(2)), family = MASS::negative.binomial(3), n_samples = 100, burnin = 10, sample_method = \"slice_sampling\", qslice_fun = qslice::slice_stepping_out, w = 0.8) log2 ############################################# # Using a different slice function log3 <- mcmcglm(formula = Y ~ ., data = dat_log, beta_prior = list(distributional::dist_normal(0, 1), distributional::dist_gamma(1, 1), distributional::dist_exponential(2)), family = MASS::negative.binomial(3), n_samples = 100, burnin = 10, sample_method = \"slice_sampling\", qslice_fun = qslice::slice_elliptical, mu = 1.5, sigma = 2) log3 } # }"},{"path":"https://mathiaslj.github.io/mcmcglm/reference/mcmcglm_across_tuningparams.html","id":null,"dir":"Reference","previous_headings":"","what":"Get list of mcmcglms run across values of slice sampling tuning parameters — mcmcglm_across_tuningparams","title":"Get list of mcmcglms run across values of slice sampling tuning parameters — mcmcglm_across_tuningparams","text":"function simply performs lapply mcmcglm values provided tuning parameter values. See general usage vignette","code":""},{"path":"https://mathiaslj.github.io/mcmcglm/reference/mcmcglm_across_tuningparams.html","id":"ref-usage","dir":"Reference","previous_headings":"","what":"Usage","title":"Get list of mcmcglms run across values of slice sampling tuning parameters — mcmcglm_across_tuningparams","text":"","code":"mcmcglm_across_tuningparams( ..., tuning_parameter_name = \"w\", formula, family, data, beta_prior = distributional::dist_normal(0, 1), log_likelihood_extra_args = list(sd = 1), sample_method = c(\"slice_sampling\", \"normal-normal\"), qslice_fun = qslice::slice_stepping_out, n_samples = 500, burnin = 100, parallelise = FALSE, n_cores = as.numeric(Sys.getenv(\"NUMBER_OF_PROCESSORS\")) - 1 )"},{"path":"https://mathiaslj.github.io/mcmcglm/reference/mcmcglm_across_tuningparams.html","id":"arguments","dir":"Reference","previous_headings":"","what":"Arguments","title":"Get list of mcmcglms run across values of slice sampling tuning parameters — mcmcglm_across_tuningparams","text":"... list vector values tuning parameter tuning_parameter_name name tuning parameter. Fx. default qslice_fun qslice::slice_stepping_out, tuning parameter called w formula object class \"formula\" (one can coerced class): symbolic description model fitted. See details stats::glm family description error distribution link function used model. can character string naming family function, family function result call family function. (See family details family functions.) data optional data frame, list environment (object coercible .data.frame data frame) containing variables model. found data, variables taken environment(formula), typically environment function called. beta_prior distribution object created function distributional package. fx. distributional::dist_normal(mean = 0, sd = 1). log_likelihood_extra_args named list arguments passed onto log_density function. Fx. specification log_likelihood_extra_args = list(sd = x) needed case family = \"gaussian\" sample_method character specifying method used sampling. default \"slice_sampling\" intended value cases, works specification family beta_prior. \"normal-normal\" uses conditional normal distribution sample case conjugate prior gaussian response beta_prior. Implemented testing purposes works niche case. qslice_fun function qslice package. Default qslice::slice_stepping_out uses slice sampler Neal 2003, functions available. n_samples numeric number samples draw parameter(/variable) model burnin numeric number samples marked \"burnin\". Burnin samples included beta_mean calculation increase finite sample performance LLN estimate parallelise logical runs algorithm across values n_vars parallelised using future.apply::future_lapply n_cores numeric number cores use parallelisation. Default 1 less number available cores","code":""},{"path":"https://mathiaslj.github.io/mcmcglm/reference/mcmcglm_across_tuningparams.html","id":"ref-examples","dir":"Reference","previous_headings":"","what":"Examples","title":"Get list of mcmcglms run across values of slice sampling tuning parameters — mcmcglm_across_tuningparams","text":"","code":"# Create test data for different scenarios n <- 100 x1 <- rnorm (n) x2 <- rbinom (n, 1, .5) b0 <- 1 b1 <- 1.5 b2 <- 2 lin_pred <- b0+b1*x1+b2*x2 ############################################# # Different families and priors # For family \"gaussian\" and iid normal prior y_norm <- rnorm(n, mean = lin_pred, sd = 1) dat_norm <- data.frame(Y = y_norm, X1 = x1, X2 = x2) w05_mcmcglms <- mcmcglm_across_tuningparams( seq(from = 0.5, by = 0.5, length.out = 9), tuning_parameter_name = \"w\", formula = Y ~ ., family = \"gaussian\", data = dat_norm, n_samples = 10, burnin = 0 )"},{"path":"https://mathiaslj.github.io/mcmcglm/reference/plot_eta_comptime.html","id":null,"dir":"Reference","previous_headings":"","what":"Plot the results of compare_eta_comptime_across_nvars — plot_eta_comptime","title":"Plot the results of compare_eta_comptime_across_nvars — plot_eta_comptime","text":"Plot results compare_eta_comptime_across_nvars","code":""},{"path":"https://mathiaslj.github.io/mcmcglm/reference/plot_eta_comptime.html","id":"ref-usage","dir":"Reference","previous_headings":"","what":"Usage","title":"Plot the results of compare_eta_comptime_across_nvars — plot_eta_comptime","text":"","code":"plot_eta_comptime(eta_comptime_data, facet_by = \"qslice_fun\")"},{"path":"https://mathiaslj.github.io/mcmcglm/reference/plot_eta_comptime.html","id":"arguments","dir":"Reference","previous_headings":"","what":"Arguments","title":"Plot the results of compare_eta_comptime_across_nvars — plot_eta_comptime","text":"eta_comptime_data data.frame results computation times. Result call compare_eta_comptime_across_nvars facet_by character variable facet plots . Default \"qslice_fun\", enabling user combine results compare_eta_comptime_across_nvars across different slice functions plot easily. options fx. running compare_eta_comptime_across_nvars different values tuning parameter seeing affects runtime","code":""},{"path":"https://mathiaslj.github.io/mcmcglm/reference/plot_eta_comptime.html","id":"value","dir":"Reference","previous_headings":"","what":"Value","title":"Plot the results of compare_eta_comptime_across_nvars — plot_eta_comptime","text":"ggplot object","code":""},{"path":"https://mathiaslj.github.io/mcmcglm/reference/plot_eta_comptime.html","id":"ref-examples","dir":"Reference","previous_headings":"","what":"Examples","title":"Plot the results of compare_eta_comptime_across_nvars — plot_eta_comptime","text":"","code":"# Compare the runtime for 2, 20, 40, 60 variables in the model res <- compare_eta_comptime_across_nvars(n_vars = c(2, seq(from = 20, to = 60, by = 20)), n_samples = 100, burnin = 20) #> Sampling from posterior ■■■■■■ 18% | ETA: 5s #> Sampling from posterior ■■■■■■■■■■■■■■■■■■■■■■ 69% | ETA: 1s #> Sampling from posterior ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ 100% | ETA: 0s #> Sampling from posterior ■■■■■■■■■■■■■■■■■ 55% | ETA: 2s #> Sampling from posterior ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ 100% | ETA: 0s #> Sampling from posterior ■■■ 5% | ETA: 26s #> Sampling from posterior ■■■■■■■■■■■■■ 41% | ETA: 6s #> Sampling from posterior ■■■■■■■■■■■■■■■■■■■■■■■■ 78% | ETA: 2s #> Sampling from posterior ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ 100% | ETA: 0s #> Sampling from posterior ■ 1% | ETA: 2m #> Sampling from posterior ■■■■■■■■■ 26% | ETA: 12s #> Sampling from posterior ■■■■■■■■■■■■■■■■■■■■■■ 70% | ETA: 3s #> Sampling from posterior ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ 100% | ETA: 0s #> Sampling from posterior ■ 1% | ETA: 4m #> Sampling from posterior ■■■ 5% | ETA: 1m #> Sampling from posterior ■■■■■■■ 21% | ETA: 26s #> Sampling from posterior ■■■■■■■■■■■■■ 41% | ETA: 14s #> Sampling from posterior ■■■■■■■■■■■■■■■■■■■■ 63% | ETA: 8s #> Sampling from posterior ■■■■■■■■■■■■■■■■■■■■■■■■■■■ 85% | ETA: 3s #> Sampling from posterior ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ 100% | ETA: 0s #> Sampling from posterior ■ 1% | ETA: 2m #> Sampling from posterior ■■■■■■■■ 22% | ETA: 14s #> Sampling from posterior ■■■■■■■■■■■■■■■ 46% | ETA: 8s #> Sampling from posterior ■■■■■■■■■■■■■■■■■■■■■■ 71% | ETA: 4s #> Sampling from posterior ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ 96% | ETA: 1s #> Sampling from posterior ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ 100% | ETA: 0s plot_eta_comptime(res)"},{"path":"https://mathiaslj.github.io/mcmcglm/reference/plot_mcmcglm_across_tuningparams.html","id":null,"dir":"Reference","previous_headings":"","what":"Plot a list of mcmcglms showing varying tuning parameters in the title of the plots — plot_mcmcglm_across_tuningparams","title":"Plot a list of mcmcglms showing varying tuning parameters in the title of the plots — plot_mcmcglm_across_tuningparams","text":"See decription mcmcglm_across_tuningparams details functionality can used ","code":""},{"path":"https://mathiaslj.github.io/mcmcglm/reference/plot_mcmcglm_across_tuningparams.html","id":"ref-usage","dir":"Reference","previous_headings":"","what":"Usage","title":"Plot a list of mcmcglms showing varying tuning parameters in the title of the plots — plot_mcmcglm_across_tuningparams","text":"","code":"plot_mcmcglm_across_tuningparams(list_mcmcglms)"},{"path":"https://mathiaslj.github.io/mcmcglm/reference/plot_mcmcglm_across_tuningparams.html","id":"arguments","dir":"Reference","previous_headings":"","what":"Arguments","title":"Plot a list of mcmcglms showing varying tuning parameters in the title of the plots — plot_mcmcglm_across_tuningparams","text":"list_mcmcglms list mcmcglm objects. Intended result call mcmcglm_across_tuningparams","code":""},{"path":"https://mathiaslj.github.io/mcmcglm/reference/plot_mcmcglm_across_tuningparams.html","id":"ref-examples","dir":"Reference","previous_headings":"","what":"Examples","title":"Plot a list of mcmcglms showing varying tuning parameters in the title of the plots — plot_mcmcglm_across_tuningparams","text":"","code":"# Create test data for different scenarios n <- 100 x1 <- rnorm (n) x2 <- rbinom (n, 1, .5) b0 <- 1 b1 <- 1.5 b2 <- 2 lin_pred <- b0+b1*x1+b2*x2 ############################################# # Different families and priors # For family \"gaussian\" and iid normal prior y_norm <- rnorm(n, mean = lin_pred, sd = 1) dat_norm <- data.frame(Y = y_norm, X1 = x1, X2 = x2) w05_mcmcglms <- mcmcglm_across_tuningparams( seq(from = 0.5, by = 0.5, length.out = 4), tuning_parameter_name = \"w\", formula = Y ~ ., family = \"gaussian\", data = dat_norm ) #> Sampling from posterior ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ 98% | ETA: 0s #> Sampling from posterior ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ 100% | ETA: 0s #> Sampling from posterior ■■■■■■■■■■■■■■ 45% | ETA: 1s #> Sampling from posterior ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ 100% | ETA: 0s #> Sampling from posterior ■■■■■■■■■■■■■■■■■■■■■■■■■ 79% | ETA: 1s #> Sampling from posterior ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ 100% | ETA: 0s plot_mcmcglm_across_tuningparams(w05_mcmcglms)"},{"path":"https://mathiaslj.github.io/mcmcglm/reference/samples.html","id":null,"dir":"Reference","previous_headings":"","what":"Get the drawn samples from the object — samples","title":"Get the drawn samples from the object — samples","text":"Get drawn samples object","code":""},{"path":"https://mathiaslj.github.io/mcmcglm/reference/samples.html","id":"ref-usage","dir":"Reference","previous_headings":"","what":"Usage","title":"Get the drawn samples from the object — samples","text":"","code":"samples(x)"},{"path":"https://mathiaslj.github.io/mcmcglm/reference/samples.html","id":"arguments","dir":"Reference","previous_headings":"","what":"Arguments","title":"Get the drawn samples from the object — samples","text":"x mcmcglm object","code":""},{"path":"https://mathiaslj.github.io/mcmcglm/reference/trace_plot.html","id":null,"dir":"Reference","previous_headings":"","what":"Create a trace plot of the MCMC samples — trace_plot","title":"Create a trace plot of the MCMC samples — trace_plot","text":"Create trace plot MCMC samples","code":""},{"path":"https://mathiaslj.github.io/mcmcglm/reference/trace_plot.html","id":"ref-usage","dir":"Reference","previous_headings":"","what":"Usage","title":"Create a trace plot of the MCMC samples — trace_plot","text":"","code":"trace_plot(x, samples_drop = NULL)"},{"path":"https://mathiaslj.github.io/mcmcglm/reference/trace_plot.html","id":"arguments","dir":"Reference","previous_headings":"","what":"Arguments","title":"Create a trace plot of the MCMC samples — trace_plot","text":"x mcmcglm object samples_drop numeric specifying number initial samples exclude trace_plot improve axis zoom plot","code":""},{"path":"https://mathiaslj.github.io/mcmcglm/reference/update_linear_predictor.html","id":null,"dir":"Reference","previous_headings":"","what":"Update value of a linear predictor as function of a single coordinate change — update_linear_predictor","title":"Update value of a linear predictor as function of a single coordinate change — update_linear_predictor","text":"Function updating linear predictor n actions rather nn_vars actions naively matrix-vector product X %% beta","code":""},{"path":"https://mathiaslj.github.io/mcmcglm/reference/update_linear_predictor.html","id":"ref-usage","dir":"Reference","previous_headings":"","what":"Usage","title":"Update value of a linear predictor as function of a single coordinate change — update_linear_predictor","text":"","code":"update_linear_predictor(new_beta_j, current_beta_j, current_eta, X_j)"},{"path":"https://mathiaslj.github.io/mcmcglm/reference/update_linear_predictor.html","id":"arguments","dir":"Reference","previous_headings":"","what":"Arguments","title":"Update value of a linear predictor as function of a single coordinate change — update_linear_predictor","text":"new_beta_j numeric new value j'th component beta parameter vector current_beta_j numeric current value j'th component beta parameter vector current_eta current value linear predictor corresponding current_beta value X_j j'th column design matrix","code":""}]