<

Stan: Bayesian Modeling Examples

statistics
Sample scripts and examples for Bayesian modeling using Stan.
Author

Sir Panda

Published

January 1, 2018

Keywords

Stan, Bayesian statistics, statistical modeling, MCMC, probabilistic programming, sensitivity analysis, uncertainty analysis, robust analysis

Please don’t mind this post. I use this to try out various highlighting styles for my code and formats.

{.algorithm} Algorithm 1 Bootstrap Confidence Interval
Input: Data X={x1,x2,...,xn}, statistic θ(X), confidence level α
Output: Confidence interval [θL,θU]
1. for b=1 to B do - Draw bootstrap sample Xb∗ by sampling n observations from X with replacement - Calculate θb∗=θ(Xb∗) 2. end for 3. Sort {θ1∗,θ2∗,...,θB∗} in ascending order 4. Set θL=quantile(θ∗,α/2) 5. Set θU=quantile(θ∗,1−α/2) 6. return [θL,θU] ``` 
\begin{algorithm}
    \caption{Quicksort}
    \begin{algorithmic}
    \PROCEDURE{Quicksort}{$A, p, r$}
        \IF{$p < r$}
            \STATE $q = $ \CALL{Partition}{$A, p, r$}
            \STATE \CALL{Quicksort}{$A, p, q - 1$}
            \STATE \CALL{Quicksort}{$A, q + 1, r$}
        \ENDIF
    \ENDPROCEDURE
    \PROCEDURE{Partition}{$A, p, r$}
        \STATE $x = A[r]$
        \STATE $i = p - 1$
        \FOR{$j = p$ \TO $r - 1$}
            \IF{$A[j] < x$}
                \STATE $i = i + 1$
                \STATE exchange $A[i]$ with $A[j]$
            \ENDIF
        \ENDFOR
        \STATE exchange $A[i + 1]$ with $A[r]$
        \RETURN $i + 1$
    \ENDPROCEDURE
    \end{algorithmic}
    \end{algorithm}
{js, echo = FALSE}
window.addEventListener('load', (event) => {
  elem = document.querySelectorAll("algorithm")
  elem.forEach(e => {
    pseudocode.renderElement(e, {
    indentSize: '1.3em',
    commentDelimiter: '//',
    lineNumber: true,
    lineNumberPunc: ':',
    noEnd: false,
    captionCount: undefined});
  })
});
window.addEventListener('load', (event) => {
  elem = document.querySelectorAll(".algorithm")
  elem.forEach(e => {
    pseudocode.renderElement(e, {
    indentSize: '1.3em',
    commentDelimiter: '//',
    lineNumber: true,
    lineNumberPunc: ':',
    noEnd: false,
    captionCount: undefined});
  })
});
YPost⁡=β0+β1Group⁡+β2Base⁡+β3Age⁡+β4Z⁡+β5R1⁡+β6R2⁡+ϵY^{\operatorname{Post}} = \beta_{0} + \beta_{1}^{\operatorname{Group}} + \beta_{2}^{\operatorname{Base}} + \beta_{3}^{\operatorname{Age}} + \beta_{4}^{\operatorname{Z}} + \beta_{5}^{\operatorname{R1}} + \beta_{6}^{\operatorname{R2}} + \epsilon

YPost=β0+β1Group+β2Base+β3Age+β4Z+β5R1+β6R2+ϵ

R Markdown

This is an R Markdown document. Markdown is a simple formatting syntax for authoring HTML, PDF, and MS Word documents. For more details on using R Markdown see http://rmarkdown.rstudio.com.

When you click the Knit button a document will be generated that includes both content as well as the output of any embedded R code chunks within the document. You can embed an R code chunk like this:

Stan

data {
  int<lower=0> J;         // number of schools
  real y[J];              // estimated treatment effects
  real<lower=0> sigma[J]; // standard error of effect estimates
}
parameters {
  real mu;                // population treatment effect
  real<lower=0> tau;      // standard deviation in treatment effects
  vector[J] eta;          // unscaled deviation from mu by school
}
transformed parameters {
  vector[J] theta = mu + tau * eta;        // school treatment effects
}
model {
  target += normal_lpdf(eta | 0, 1);       // prior log-density
  target += normal_lpdf(y | theta, sigma); // log-likelihood
}
df1 <- read.csv("../../../ts.csv")
df1$y <- ts(df1$Sales)
df1$ds <- as.Date(df1$Time.Increment)

splits <- initial_time_split(df1, prop = 0.5)
train <- training(splits)
test <- testing(splits)
interactive <- TRUE
# Forecasting with auto.arima
library("forecast")
md <- auto.arima(train$y)
fc <- forecast(md, h = 12)

model_fit_arima_no_boost <- arima_reg() %>%
  set_engine(engine = "auto_arima") %>%
  fit(y ~ ds, data = training(splits))

# Model 2: arima_boost ----
model_fit_arima_boosted <- arima_boost(
  min_n = 2,
  learn_rate = 0.015
) %>%
  set_engine(engine = "auto_arima_xgboost") %>%
  fit(y ~ ds + as.numeric(ds) + factor(month(ds, label = TRUE),
    ordered = F
  ),
  data = training(splits)
  )

# Model 3: ets ----
model_fit_ets <- exp_smoothing() %>%
  set_engine(engine = "ets") %>%
  fit(y ~ ds, data = training(splits))

model_fit_lm <- linear_reg() %>%
  set_engine("lm") %>%
  fit(y ~ as.numeric(ds) + factor(month(ds, label = TRUE),
    ordered = FALSE
  ),
  data = training(splits)
  )

# Model 4: prophet ----
model_fit_prophet <- prophet_reg() %>%
  set_engine(engine = "prophet") %>%
  fit(y ~ ds, data = training(splits))
# Model 6: earth ----
model_spec_mars <- mars(mode = "regression") %>%
  set_engine("earth")


recipe_spec <- recipe(y ~ ds, data = training(splits)) %>%
  step_date(ds, features = "month", ordinal = FALSE) %>%
  step_mutate(date_num = as.numeric(ds)) %>%
  step_normalize(date_num) %>%
  step_rm(ds)

wflw_fit_mars <- workflow() %>%
  add_recipe(recipe_spec) %>%
  add_model(model_spec_mars) %>%
  fit(training(splits))

models_tbl <- modeltime_table(
  model_fit_arima_no_boost,
  model_fit_arima_boosted,
  model_fit_ets,
  model_fit_prophet,
  model_fit_lm,
  wflw_fit_mars
)

calibration_tbl <- models_tbl %>%
  modeltime_calibrate(new_data = testing(splits))

calibration_tbl %>%
  modeltime_forecast(
    new_data    = testing(splits),
    actual_data = df1
  ) %>%
  plot_modeltime_forecast(
    .legend_max_width = 25, # For mobile screens
    .interactive      = interactive
  )
Oct 2020Jan 2021Apr 2021Jul 2021Oct 2021Jan 2022Apr 2022−20k020k40k60k80k100k
LegendACTUAL1_ARIMA(0,1,0)3_ETS(A,N,N)4_PROPHET6_EARTHForecast Plot
calibration_tbl %>%
  modeltime_accuracy() %>%
  table_modeltime_accuracy(
    .interactive = FALSE
  )
Accuracy Table
.model_id .model_desc .type mae mape mase smape rmse rsq
1 ARIMA(0,1,0) Test 3200 14.0 0.74 13.6 4252 NA
2 ARIMA(0,1,0) W/ XGBOOST ERRORS NA NA NA NA NA NA NA
3 ETS(A,N,N) Test 3200 14.0 0.74 13.6 4252 NA
4 PROPHET Test 20688 96.2 4.77 62.4 22234 0.03
5 LM NA NA NA NA NA NA NA
6 EARTH Test 3074 12.7 0.71 12.8 4743 0.01 
refit_tbl <- calibration_tbl %>%
  modeltime_refit(data = df1)

refit_tbl %>%
  modeltime_forecast(h = "4 years", actual_data = df1) %>%
  plot_modeltime_forecast(
    .legend_max_width = 10, # For mobile screens
    .interactive      = TRUE
  )
202120222023202420252026−20k020k40k60k80k100k
LegendACTUAL1_ARIMA...2_ARIMA...3_ETS(A...4_PROPHET5_LM6_EARTHForecast Plot
plot(greybox::forecast(smooth::adam(df1$y, h = 12, holdout = TRUE)))


plot(greybox::forecast(smooth::es(df1$y, h = 12, holdout = TRUE,
    silent = FALSE)))
#> Forming the pool of models based on... ANN , AAN , Estimation progress:    100 %... Done!


s1 <- bayesforecast::stan_naive(ts = df1$y, chains = 4, iter = 4000,
    cores = 8)

plot(s1) + theme_bw()

check_residuals(s1) + theme_light()

#> NULL
#> LING FOR MODEL 'Sarima' NOW (CHAIN 4).
#> Chain 4: 
#> Chain 4: Gradient evaluation took 3.3e-05 seconds
#> Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.33 seconds.
#> Chain 4: Adjust your expectations accordingly!
#> Chain 4: 
#> Chain 4: 
#> Chain 4: Iteration:    1 / 4000 [  0%]  (Warmup)
#> Chain 2: Iteration: 1200 / 4000 [ 30%]  (Warmup)
#> Chain 1: Iteration: 2400 / 4000 [ 60%]  (Sampling)
#> Chain 3: Iteration:  800 / 4000 [ 20%]  (Warmup)
#> Chain 2: Iteration: 1600 / 4000 [ 40%]  (Warmup)
#> Chain 1: Iteration: 2800 / 4000 [ 70%]  (Sampling)
#> Chain 4: Iteration:  400 / 4000 [ 10%]  (Warmup)
#> Chain 3: Iteration: 1200 / 4000 [ 30%]  (Warmup)
#> Chain 2: Iteration: 2000 / 4000 [ 50%]  (Warmup)
#> Chain 2: Iteration: 2001 / 4000 [ 50%]  (Sampling)
#> Chain 1: Iteration: 3200 / 4000 [ 80%]  (Sampling)
#> Chain 3: Iteration: 1600 / 4000 [ 40%]  (Warmup)
#> Chain 4: Iteration:  800 / 4000 [ 20%]  (Warmup)
#> Chain 2: Iteration: 2400 / 4000 [ 60%]  (Sampling)
#> Chain 1: Iteration: 3600 / 4000 [ 90%]  (Sampling)
#> Chain 3: Iteration: 2000 / 4000 [ 50%]  (Warmup)
#> Chain 3: Iteration: 2001 / 4000 [ 50%]  (Sampling)
#> Chain 4: Iteration: 1200 / 4000 [ 30%]  (Warmup)
#> Chain 2: Iteration: 2800 / 4000 [ 70%]  (Sampling)
#> Chain 3: Iteration: 2400 / 4000 [ 60%]  (Sampling)
#> Chain 1: Iteration: 4000 / 4000 [100%]  (Sampling)
#> Chain 1: 
#> Chain 1:  Elapsed Time: 0.022 seconds (Warm-up)
#> Chain 1:                0.023 seconds (Sampling)
#> Chain 1:                0.045 seconds (Total)
#> Chain 1: 
#> Chain 4: Iteration: 1600 / 4000 [ 40%]  (Warmup)
#> Chain 2: Iteration: 3200 / 4000 [ 80%]  (Sampling)
#> Chain 3: Iteration: 2800 / 4000 [ 70%]  (Sampling)
#> Chain 4: Iteration: 2000 / 4000 [ 50%]  (Warmup)
#> Chain 4: Iteration: 2001 / 4000 [ 50%]  (Sampling)
#> Chain 3: Iteration: 3200 / 4000 [ 80%]  (Sampling)
#> Chain 2: Iteration: 3600 / 4000 [ 90%]  (Sampling)
#> Chain 4: Iteration: 2400 / 4000 [ 60%]  (Sampling)
#> Chain 3: Iteration: 3600 / 4000 [ 90%]  (Sampling)
#> Chain 2: Iteration: 4000 / 4000 [100%]  (Sampling)
#> Chain 2: 
#> Chain 2:  Elapsed Time: 0.024 seconds (Warm-up)
#> Chain 2:                0.024 seconds (Sampling)
#> Chain 2:                0.048 seconds (Total)
#> Chain 2: 
#> Chain 3: Iteration: 4000 / 4000 [100%]  (Sampling)
#> Chain 3: 
#> Chain 3:  Elapsed Time: 0.023 seconds (Warm-up)
#> Chain 3:                0.018 seconds (Sampling)
#> Chain 3:                0.041 seconds (Total)
#> Chain 3: 
#> Chain 4: Iteration: 2800 / 4000 [ 70%]  (Sampling)
#> Chain 4: Iteration: 3200 / 4000 [ 80%]  (Sampling)
#> Chain 4: Iteration: 3600 / 4000 [ 90%]  (Sampling)
#> Chain 4: Iteration: 4000 / 4000 [100%]  (Sampling)
#> Chain 4: 
#> Chain 4:  Elapsed Time: 0.025 seconds (Warm-up)
#> Chain 4:                0.022 seconds (Sampling)
#> Chain 4:                0.047 seconds (Total)
#> Chain 4:
autoplot(object = forecast(s1, h = 12, biasadj = TRUE, PI = TRUE),
    include = 100) + theme_bw()

autoplot(object = forecast(s1, h = 52, biasadj = TRUE, PI = TRUE),
    include = 100) + theme_bw()

ztable(summary(s1))
mean se 5% 95% ess Rhat
mu0 0.02 0.03 -4.86 4.94 7744 1
sigma0 4664.66 7.50 3716.95 5872.17 7872 1
loglik -200.57 0.01 -202.77 -199.71 7663 1 
(meanf(df1$y))
(forecast(s1, h = 12))

autoplot(object = forecast(s1, h = 6, biasadj = TRUE, PI = TRUE),
    include = 100) + theme_bw()


df <- as.data.frame(df1)
m <- prophet(df1, growth = "linear", yearly.seasonality = "auto")
future <- make_future_dataframe(m, periods = 60, freq = "months",
    include_history = TRUE)
forecast <- predict(m, future)
plot(m, forecast) + theme_bw()


plot(greybox::forecast(smooth::auto.adam(df1$y, h = 12, holdout = TRUE,
    ic = "AICc", regressors = "select")))

adamAutoARIMAAir <- auto.adam(df1$y, h = 50)

plot(greybox::forecast(adam(df1$y, main = "Parametric prediction interval",
    h = 5, sim = 1000, lev = 0.99, interval = "prediction")))

plot(greybox::forecast(auto.adam(df1$y, h = 5, interval = "complete",
    nsim = 100, main = "Complete prediction interval")))


plot(greybox::forecast(adam(df1$y, c("CCN", "ANN", "AAN", "AAdN"),
    h = 10, holdout = TRUE, ic = "AICc")))


plot(greybox::forecast(adam(df1$y, model = "NNN", lags = 7, orders = c(0,
    1, 1), constant = TRUE, h = 5, interval = "complete", nsim = 100)))


plot(greybox::forecast((adam(df1$y, model = "NNN", lags = c(24,
    24 * 7, 24 * 365), orders = list(ar = c(3, 2, 2, 2), i = c(2,
    1, 1, 1), ma = c(3, 2, 2, 2), select = TRUE), initial = "backcasting"))))


adamARIMA
#> Error:
#> ! object 'adamARIMA' not found
# Apply models
adamPoolBJ <- vector("list", 3)
adamPoolBJ[[1]] <- adam(df1$y, "ZZN", h = 10, holdout = TRUE,
    ic = "BICc")
adamPoolBJ[[2]] <- adam(df1$y, "NNN", orders = list(ar = 3, i = 2,
    ma = 3, select = TRUE), h = 10, holdout = TRUE, ic = "BICc")
adamPoolBJ[[3]] <- adam(df1$y, "MMN", h = 10, holdout = TRUE,
    ic = "BICc", regressors = "select")

# Extract BICc values
adamsICs <- sapply(adamPoolBJ, BICc)

# Calculate weights
adamsICWeights <- adamsICs - min(adamsICs)
adamsICWeights[] <- exp(-0.5 * adamsICWeights)/sum(exp(-0.5 *
    adamsICWeights))
names(adamsICWeights) <- c("ETS", "ARIMA", "ETSX")
round(adamsICWeights, 3)

adamPoolBJForecasts <- vector("list", 3)
# Produce forecasts from the three models
for (i in 1:3) {
    adamPoolBJForecasts[[i]] <- forecast(adamPoolBJ[[i]], h = 10,
        interval = "pred")
}
# Produce combined conditional means and prediction
# intervals
finalForecast <- cbind(sapply(adamPoolBJForecasts, "[[", "mean") %*%
    adamsICWeights, sapply(adamPoolBJForecasts, "[[", "lower") %*%
    adamsICWeights, sapply(adamPoolBJForecasts, "[[", "upper") %*%
    adamsICWeights)
# Give the appropriate names
colnames(finalForecast) <- c("Mean", "Lower bound (2.5%)", "Upper bound (97.5%)")
# Transform the table in the ts format (for convenience)
finalForecast <- ts(finalForecast, start = start(adamPoolBJForecasts[[i]]$mean))
finalForecast

graphmaker(df1$y, finalForecast[, 1], lower = finalForecast[,
    2], upper = finalForecast[, 3], level = 0.95)


plot(forecast(forecastHybrid::hybridModel(df1$y, models = "aen",
    weights = "equal", cvHorizon = 8, num.cores = 4), h = 5))


plot(forecast(forecastHybrid::hybridModel(df1$y, models = "fnst",
    weights = "equal", errorMethod = "RMSE")))


oesModel <- oes(df1$y, model = "YYY", occurrence = "auto")

y0 <- stan_sarima(df1$y, refresh = 0, verbose = FALSE, open_progress = FALSE)
autoplot(forecast(y0, 5, 0.99), "red") + ggplot2::theme_bw()
#> Error in `tail.data.frame()`:
#> ! invalid 'n' - must be numeric, possibly NA.

df1$month <- lubridate::month(df1$ds)
gam1 <- mgcv::gam(Sales ~ s(month, bs = "cr", k = 12), data = df1,
    family = gaussian, correlation = SARIMA(form = ~month, p = 1),
    method = "REML") |>
    mgcv::plot.gam(lwd = 3, lty = 1, col = "#d46c5b")

ts_plot(df1, title = "US Monthly Natural Gas Consumption", Ytitle = "Billion Cubic Feet")
Oct 2020Jan 2021Apr 2021Jul 2021Oct 2021Jan 2022Apr 202205k10k15k20k25k30k35k
SalesOrdersTicket.SizemonthUS Monthly Natural Gas ConsumptionBillion Cubic Feet
months <- c("2022-01-01", "2022-03-01", "2022-06-01")
ubereats <- c(7327.55, 4653.53, 4833.21)
doordash <- c(1304.54, 2000.35, 1643.58)
grubhub <- c(1199.85, 941.68, 623.27)
total <- c(7222.85, 7464.68, 7100.06)

df <- data.frame(months, ubereats, doordash, grubhub, total)
df$months <- as.Date(df$months)


ts_plot(df, title = "US Monthly Natural Gas Consumption", Ytitle = "Billion Cubic Feet")
Jan 2022Feb 2022Mar 2022Apr 2022May 2022Jun 20221000200030004000500060007000
ubereatsdoordashgrubhubtotalUS Monthly Natural Gas ConsumptionBillion Cubic Feet
df$total <- ts(df$total)
plot(greybox::forecast(adam(df$total, h = 5, holdout = TRUE,
    ic = "AICc", regressors = "select")))
#> Error in `h()`:
#> ! error in evaluating the argument 'x' in selecting a method for function 'plot': The number of in-sample observations is not positive. Cannot do anything.

m <- prophet(df, growth = "linear", yearly.seasonality = "auto")
#> Error in `fit.prophet()`:
#> ! Dataframe must have columns 'ds' and 'y' with the dates and values respectively.
future <- make_future_dataframe(m, periods = 60, freq = "months",
    include_history = TRUE)
forecast <- predict(m, future)
plot(m, forecast) + theme_bw()

# Plotting actual vs. fitted and forecasted
test_forecast(actual = df1$y, forecast.obj = fc, test = test$y)
#> Error in `recycle_columns()`:
#> ! Tibble columns must have compatible sizes.
#> • Size 21: Column `x`.
#> • Size 22: Columns `y` and `text`.
#> ℹ Only values of size one are recycled.
plot(fc)

import delimited "../../../ts.csv", clear

Trace plot of imputed datasets.

Python

from pandas import DataFrame
import statsmodels.api as sm
import matplotlib as plot

Stock_Market = {'Year': [2017,2017,2017,2017,2017,
                         2017,2017,2017,2017,2017,
                         2017,2017,2016,2016,2016,
                         2016,2016,2016,2016,2016,
                         2016,2016,2016,2016],
                'Month': [12, 11,10,9,8,7,6,5,4,
                          3,2,1,12,11,10,9,8,7,6,
                          5,4,3,2,1],
                'Interest_Rate':[2.75,2.5,2.5,2.5,2.5,2.5,
                                 2.5,2.25,2.25,2.25,2,2,2,1.75,1.75,
                                 1.75,1.75,1.75,1.75,1.75,1.75,1.75,1.75,1.75],
                'Unemployment_Rate':[5.3,5.3,5.3,5.3,5.4,5.6,
                                     5.5,5.5,5.5,5.6,5.7,5.9,6,5.9,5.8,6.1,
                                     6.2,6.1,6.1,6.1,5.9,6.2,6.2,6.1],
                'Stock_Index_Price': [1464,1394,1357,1293,1256,1254,1234,1195,1159,1167,1130,
                                      1075,1047,965,943,958,971,949,884,866,876,822,704,719]
                }

df = DataFrame(Stock_Market,columns=['Year','Month','Interest_Rate',
                                     'Unemployment_Rate','Stock_Index_Price'])

X = df[['Interest_Rate','Unemployment_Rate']]

# here we have 2 variables for the multiple linear regression. If you just want to use one variable for simple linear regression, then use X = df['Interest_Rate'] for example

Y = df['Stock_Index_Price']

X = sm.add_constant(X) # adding a constant

model = sm.OLS(Y, X).fit()
predictions = model.predict(X)

print_model = model.summary()
print(print_model)
#>                             OLS Regression Results                            
#> ==============================================================================
#> Dep. Variable:      Stock_Index_Price   R-squared:                       0.898
#> Model:                            OLS   Adj. R-squared:                  0.888
#> Method:                 Least Squares   F-statistic:                     92.07
#> Date:                Mon, 02 Feb 2026   Prob (F-statistic):           4.04e-11
#> Time:                        23:46:44   Log-Likelihood:                -134.61
#> No. Observations:                  24   AIC:                             275.2
#> Df Residuals:                      21   BIC:                             278.8
#> Df Model:                           2                                         
#> Covariance Type:            nonrobust                                         
#> =====================================================================================
#>                         coef    std err          t      P>|t|      [0.025      0.975]
#> -------------------------------------------------------------------------------------
#> const              1798.4040    899.248      2.000      0.059     -71.685    3668.493
#> Interest_Rate       345.5401    111.367      3.103      0.005     113.940     577.140
#> Unemployment_Rate  -250.1466    117.950     -2.121      0.046    -495.437      -4.856
#> ==============================================================================
#> Omnibus:                        2.691   Durbin-Watson:                   0.530
#> Prob(Omnibus):                  0.260   Jarque-Bera (JB):                1.551
#> Skew:                          -0.612   Prob(JB):                        0.461
#> Kurtosis:                       3.226   Cond. No.                         394.
#> ==============================================================================
#> 
#> Notes:
#> [1] Standard Errors assume that the covariance matrix of the errors is correctly specified.

R

fit <- glm(mpg ~ cyl + disp, mtcars, family = gaussian())
# show the theoretical model
equatiomatic::extract_eq(fit)

E(mpg)=α+β1(cyl)+β2(disp)

Stata

sysuse auto2, clear
parallel initialize 8, f
mfp: glm price mpg
twoway (fpfitci price mpg, estcmd(glm) fcolor(dkorange%20) alcolor(%40))  || scatter price mpg, mcolor(dkorange) scale(0.75)
graph export "mfp.png", replace

Trace plot of imputed datasets. 
clear
set obs 100
/**
Generate variables from a normal distribution like before
*/
generate x = rnormal(0, 1)
generate y = rnormal(0, 1)
/**
We set up our model here

*/
parallel initialize 8, f
bayesmh y x, likelihood(normal({var})) prior({var}, normal(0, 10)) ///
prior({y:}, normal(0, 10)) rseed(1031) saving(coutput_pred, replace) mcmcsize(1000)
/**
We use the bayespredict command to make predictions from the model
*/
bayespredict (mean:@mean({_resid})) (var:@variance({_resid})), ///
rseed(1031) saving(coutput_pred, replace)
/**
Then we calculate the posterior predictive P-values
*/
bayesstats ppvalues {mean} using coutput_pred
library(lme4)
library(simr)

# Toy model
fm = lmer(y ~ x + (x | g), data = simdata)

# Extend sample size of `g`
fm_extended_g = extend(fm, along = 'g', n = 4)
# 4 levels of g
pwcurve_4g = powerCurve(fm_extended_g, fixed('x'), along = 'g', breaks = 4,
                        nsim = 50, seed = 123,
                        # No progress bar
                        progress = FALSE)
# 6 levels of g
# Create a destination object using any of the power curves above.
all_pwcurve = pwcurve_4g
# Combine results
all_pwcurve$ps = c(pwcurve_4g$ps[1])

# Combine the different numbers of levels.
all_pwcurve$xval = c(pwcurve_4g$nlevels)

print(all_pwcurve)
#> Power for predictor 'x', (95% confidence interval),
#> by number of levels in g:
#>       4: 44.00% (29.99, 58.75) - 40 rows
#> 
#> Time elapsed: 0 h 0 m 3 s

plot(all_pwcurve, xlab = 'Levels of g')

DataCamp notebooks:

Citation

BibTeX citation:
@misc{panda2018,
  author = {Panda, Sir},
  title = {Stan: {Bayesian} {Modeling} {Examples}},
  date = {2018-01-01},
  url = {https://lesslikely.com/posts/statistics/stan},
  langid = {en},
  abstract = {Sample scripts and examples for Bayesian modeling using
    Stan.}
}
For attribution, please cite this work as:
1. Panda S. (2018). ‘Stan: Bayesian Modeling Examples’. Less Likely. https://lesslikely.com/posts/statistics/stan.