Sample size simulation. See Chapter 16 in Regression and Other Stories.


Load packages

library("rprojroot")
root<-has_file(".ROS-Examples-root")$make_fix_file()
library("rstanarm")

Simulated data 1: predictor range (-0.5, 0.5)

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 models

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

Simulated data 2: predictor range (0, 1)

x1 <- sample(c(0,1), N, replace=TRUE)
x2 <- sample(c(0,1), N, replace=TRUE)

Fit models

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

Simulated data 2: predictor range (-1, 1)

x1 <- sample(c(-1,1), N, replace=TRUE)
x2 <- sample(c(-1,1), N, replace=TRUE)

Fit models

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
LS0tCnRpdGxlOiAiUmVncmVzc2lvbiBhbmQgT3RoZXIgU3RvcmllczogU2FtcGxlIHNpemUgc2ltdWxhdGlvbiIKYXV0aG9yOiAiQW5kcmV3IEdlbG1hbiwgSmVubmlmZXIgSGlsbCwgQWtpIFZlaHRhcmkiCmRhdGU6ICJgciBmb3JtYXQoU3lzLkRhdGUoKSlgIgpvdXRwdXQ6CiAgaHRtbF9kb2N1bWVudDoKICAgIHRoZW1lOiByZWFkYWJsZQogICAgdG9jOiB0cnVlCiAgICB0b2NfZGVwdGg6IDIKICAgIHRvY19mbG9hdDogdHJ1ZQogICAgY29kZV9kb3dubG9hZDogdHJ1ZQotLS0KU2FtcGxlIHNpemUgc2ltdWxhdGlvbi4gU2VlIENoYXB0ZXIgMTYgaW4gUmVncmVzc2lvbiBhbmQgT3RoZXIgU3Rvcmllcy4KCi0tLS0tLS0tLS0tLS0KCgpgYGB7ciBzZXR1cCwgaW5jbHVkZT1GQUxTRX0Ka25pdHI6Om9wdHNfY2h1bmskc2V0KG1lc3NhZ2U9RkFMU0UsIGVycm9yPUZBTFNFLCB3YXJuaW5nPUZBTFNFLCBjb21tZW50PU5BKQpgYGAKCiMjIyMgTG9hZCBwYWNrYWdlcwoKYGBge3IgfQpsaWJyYXJ5KCJycHJvanJvb3QiKQpyb290PC1oYXNfZmlsZSgiLlJPUy1FeGFtcGxlcy1yb290IikkbWFrZV9maXhfZmlsZSgpCmxpYnJhcnkoInJzdGFuYXJtIikKYGBgCgojIyMjIFNpbXVsYXRlZCBkYXRhIDE6IHByZWRpY3RvciByYW5nZSAoLTAuNSwgMC41KQoKYGBge3IgfQpOIDwtIDEwMDAKc2lnbWEgPC0gMTAKeSA8LSBybm9ybShOLCAwLCBzaWdtYSkKeDEgPC0gc2FtcGxlKGMoLTAuNSwwLjUpLCBOLCByZXBsYWNlPVRSVUUpCngyIDwtIHNhbXBsZShjKC0wLjUsMC41KSwgTiwgcmVwbGFjZT1UUlVFKQpmYWtlIDwtIGRhdGEuZnJhbWUoYyh5LHgxLHgyKSkKYGBgCgojIyMjIEZpdCBtb2RlbHMKCmBgYHtyIH0KZml0XzFhIDwtIHN0YW5fZ2xtKHkgfiB4MSwgZGF0YSA9IGZha2UsIHJlZnJlc2ggPSAwKQpmaXRfMWIgPC0gc3Rhbl9nbG0oeSB+IHgxICsgeDIgKyB4MTp4MiwgZGF0YSA9IGZha2UsIHJlZnJlc2ggPSAwKQpwcmludChmaXRfMWEpCnByaW50KGZpdF8xYikKYGBgCgojIyMjIFNpbXVsYXRlZCBkYXRhIDI6IHByZWRpY3RvciByYW5nZSAoMCwgMSkKCmBgYHtyIH0KeDEgPC0gc2FtcGxlKGMoMCwxKSwgTiwgcmVwbGFjZT1UUlVFKQp4MiA8LSBzYW1wbGUoYygwLDEpLCBOLCByZXBsYWNlPVRSVUUpCmBgYAoKIyMjIyBGaXQgbW9kZWxzCgpgYGB7ciB9CmZpdF8yYSA8LSBzdGFuX2dsbSh5IH4geDEsIGRhdGEgPSBmYWtlLCByZWZyZXNoID0gMCkKZml0XzJiIDwtIHN0YW5fZ2xtKHkgfiB4MSArIHgyICsgeDE6eDIsIGRhdGEgPSBmYWtlLCByZWZyZXNoID0gMCkKcHJpbnQoZml0XzJhKQpwcmludChmaXRfMmIpCmBgYAoKIyMjIyBTaW11bGF0ZWQgZGF0YSAyOiBwcmVkaWN0b3IgcmFuZ2UgKC0xLCAxKQoKYGBge3IgfQp4MSA8LSBzYW1wbGUoYygtMSwxKSwgTiwgcmVwbGFjZT1UUlVFKQp4MiA8LSBzYW1wbGUoYygtMSwxKSwgTiwgcmVwbGFjZT1UUlVFKQpgYGAKCiMjIyMgRml0IG1vZGVscwoKYGBge3IgfQpmaXRfM2EgPC0gc3Rhbl9nbG0oeSB+IHgxLCBkYXRhID0gZmFrZSwgcmVmcmVzaCA9IDApCmZpdF8zYiA8LSBzdGFuX2dsbSh5IH4geDEgKyB4MiArIHgxOngyLCBkYXRhID0gZmFrZSwgcmVmcmVzaCA9IDApCnByaW50KGZpdF8zYSkKcHJpbnQoZml0XzNiKQpgYGAKCg==