Sample size simulation. See Chapter 16 in Regression and Other Stories.
library("rprojroot")
root<-has_file(".ROS-Examples-root")$make_fix_file()
library("rstanarm")
N <- 1000
sigma <- 10
y <- rnorm(N, 0, sigma)
x1 <- sample(c(-0.5,0.5), N, replace=TRUE)
x2 <- sample(c(-0.5,0.5), N, replace=TRUE)
fake <- data.frame(c(y,x1,x2))
fit_1a <- stan_glm(y ~ x1, data = fake, refresh = 0)
fit_1b <- stan_glm(y ~ x1 + x2 + x1:x2, data = fake, refresh = 0)
print(fit_1a)
stan_glm
family: gaussian [identity]
formula: y ~ x1
observations: 1000
predictors: 2
------
Median MAD_SD
(Intercept) 0.5 0.3
x1 -0.1 0.6
Auxiliary parameter(s):
Median MAD_SD
sigma 9.8 0.2
------
* For help interpreting the printed output see ?print.stanreg
* For info on the priors used see ?prior_summary.stanreg
print(fit_1b)
stan_glm
family: gaussian [identity]
formula: y ~ x1 + x2 + x1:x2
observations: 1000
predictors: 4
------
Median MAD_SD
(Intercept) 0.5 0.3
x1 -0.1 0.6
x2 -0.2 0.6
x1:x2 1.0 1.2
Auxiliary parameter(s):
Median MAD_SD
sigma 9.8 0.2
------
* For help interpreting the printed output see ?print.stanreg
* For info on the priors used see ?prior_summary.stanreg
x1 <- sample(c(0,1), N, replace=TRUE)
x2 <- sample(c(0,1), N, replace=TRUE)
fit_2a <- stan_glm(y ~ x1, data = fake, refresh = 0)
fit_2b <- stan_glm(y ~ x1 + x2 + x1:x2, data = fake, refresh = 0)
print(fit_2a)
stan_glm
family: gaussian [identity]
formula: y ~ x1
observations: 1000
predictors: 2
------
Median MAD_SD
(Intercept) 0.3 0.4
x1 0.4 0.6
Auxiliary parameter(s):
Median MAD_SD
sigma 9.8 0.2
------
* For help interpreting the printed output see ?print.stanreg
* For info on the priors used see ?prior_summary.stanreg
print(fit_2b)
stan_glm
family: gaussian [identity]
formula: y ~ x1 + x2 + x1:x2
observations: 1000
predictors: 4
------
Median MAD_SD
(Intercept) 0.4 0.6
x1 0.9 0.8
x2 -0.1 0.9
x1:x2 -1.0 1.2
Auxiliary parameter(s):
Median MAD_SD
sigma 9.8 0.2
------
* For help interpreting the printed output see ?print.stanreg
* For info on the priors used see ?prior_summary.stanreg
x1 <- sample(c(-1,1), N, replace=TRUE)
x2 <- sample(c(-1,1), N, replace=TRUE)
fit_3a <- stan_glm(y ~ x1, data = fake, refresh = 0)
fit_3b <- stan_glm(y ~ x1 + x2 + x1:x2, data = fake, refresh = 0)
print(fit_3a)
stan_glm
family: gaussian [identity]
formula: y ~ x1
observations: 1000
predictors: 2
------
Median MAD_SD
(Intercept) 0.5 0.3
x1 0.0 0.3
Auxiliary parameter(s):
Median MAD_SD
sigma 9.8 0.2
------
* For help interpreting the printed output see ?print.stanreg
* For info on the priors used see ?prior_summary.stanreg
print(fit_3b)
stan_glm
family: gaussian [identity]
formula: y ~ x1 + x2 + x1:x2
observations: 1000
predictors: 4
------
Median MAD_SD
(Intercept) 0.5 0.3
x1 0.0 0.3
x2 0.5 0.3
x1:x2 -0.2 0.3
Auxiliary parameter(s):
Median MAD_SD
sigma 9.8 0.2
------
* For help interpreting the printed output see ?print.stanreg
* For info on the priors used see ?prior_summary.stanreg