Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
UseR 2019 tutorial - Get up to speed with Bayes test script
# Prior to the tutorial make sure that the script below runs without error on your R installation.
# You first need to install the follwoing packages:
# install.packages(c("rstanarm", "prophet", "CausalImpact"))
library(rstanarm)
library(prophet)
library(CausalImpact)
# This will test that rstanarm works
# Don't be alarmed if you get a warning about "divergent transitions "
fit <- stan_lm(mpg ~ wt + qsec + am, data = mtcars, prior = R2(0.75))
plot(fit, prob = 0.8)
# This tests that prophet is working
history <- data.frame(ds = seq(as.Date('2015-01-01'), as.Date('2016-01-01'), by = 'd'),
y = sin(1:366/200) + rnorm(366)/10)
m <- prophet(history)
future <- make_future_dataframe(m, periods = 365)
forecast <- predict(m, future)
plot(m, forecast)
# This tests that CausalImpact is working
# First simulating some data
x1 <- 100 + arima.sim(model = list(ar = 0.999), n = 52)
y <- 1.2 * x1 + rnorm(52)
y[41:52] <- y[41:52] + 10
data <- cbind(y, x1)
pre.period <- c(1, 40)
post.period <- c(41, 52)
# Then running CausalImpact
impact <- CausalImpact(data, pre.period, post.period)
plot(impact)
@SamuelWiqvist
Copy link

SamuelWiqvist commented Jul 3, 2019

@rasmusab Thanks for your comments!

I finally got rstanarm to work by before installing it running sudo apt install libssl-dev (see the comment related to Ubuntu here https://github.com/rocker-org/rocker/issues/124).

@johanrex
Copy link

johanrex commented Jul 4, 2019

@rasmusab Is this the expected output? There is some output regarding "Disabling yearly seasonality." and warnings regarding "Removed 52 rows containing missing values (geom_path)."

> library(rstanarm)
> library(prophet)
> library(CausalImpact)
> 
> # This will test that rstanarm works
> # Don't be alarmed if you get a warning about "divergent transitions "
> fit <- stan_lm(mpg ~ wt + qsec + am, data = mtcars, prior = R2(0.75))

SAMPLING FOR MODEL 'lm' NOW (CHAIN 1).
Chain 1: 
Chain 1: Gradient evaluation took 0 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1: 
Chain 1: 
Chain 1: Iteration:    1 / 2000 [  0%]  (Warmup)
Chain 1: Iteration:  200 / 2000 [ 10%]  (Warmup)
Chain 1: Iteration:  400 / 2000 [ 20%]  (Warmup)
Chain 1: Iteration:  600 / 2000 [ 30%]  (Warmup)
Chain 1: Iteration:  800 / 2000 [ 40%]  (Warmup)
Chain 1: Iteration: 1000 / 2000 [ 50%]  (Warmup)
Chain 1: Iteration: 1001 / 2000 [ 50%]  (Sampling)
Chain 1: Iteration: 1200 / 2000 [ 60%]  (Sampling)
Chain 1: Iteration: 1400 / 2000 [ 70%]  (Sampling)
Chain 1: Iteration: 1600 / 2000 [ 80%]  (Sampling)
Chain 1: Iteration: 1800 / 2000 [ 90%]  (Sampling)
Chain 1: Iteration: 2000 / 2000 [100%]  (Sampling)
Chain 1: 
Chain 1:  Elapsed Time: 3.998 seconds (Warm-up)
Chain 1:                2.404 seconds (Sampling)
Chain 1:                6.402 seconds (Total)
Chain 1: 

SAMPLING FOR MODEL 'lm' NOW (CHAIN 2).
Chain 2: 
Chain 2: Gradient evaluation took 0 seconds
Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
Chain 2: Adjust your expectations accordingly!
Chain 2: 
Chain 2: 
Chain 2: Iteration:    1 / 2000 [  0%]  (Warmup)
Chain 2: Iteration:  200 / 2000 [ 10%]  (Warmup)
Chain 2: Iteration:  400 / 2000 [ 20%]  (Warmup)
Chain 2: Iteration:  600 / 2000 [ 30%]  (Warmup)
Chain 2: Iteration:  800 / 2000 [ 40%]  (Warmup)
Chain 2: Iteration: 1000 / 2000 [ 50%]  (Warmup)
Chain 2: Iteration: 1001 / 2000 [ 50%]  (Sampling)
Chain 2: Iteration: 1200 / 2000 [ 60%]  (Sampling)
Chain 2: Iteration: 1400 / 2000 [ 70%]  (Sampling)
Chain 2: Iteration: 1600 / 2000 [ 80%]  (Sampling)
Chain 2: Iteration: 1800 / 2000 [ 90%]  (Sampling)
Chain 2: Iteration: 2000 / 2000 [100%]  (Sampling)
Chain 2: 
Chain 2:  Elapsed Time: 4.17 seconds (Warm-up)
Chain 2:                5.554 seconds (Sampling)
Chain 2:                9.724 seconds (Total)
Chain 2: 

SAMPLING FOR MODEL 'lm' NOW (CHAIN 3).
Chain 3: 
Chain 3: Gradient evaluation took 0 seconds
Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
Chain 3: Adjust your expectations accordingly!
Chain 3: 
Chain 3: 
Chain 3: Iteration:    1 / 2000 [  0%]  (Warmup)
Chain 3: Iteration:  200 / 2000 [ 10%]  (Warmup)
Chain 3: Iteration:  400 / 2000 [ 20%]  (Warmup)
Chain 3: Iteration:  600 / 2000 [ 30%]  (Warmup)
Chain 3: Iteration:  800 / 2000 [ 40%]  (Warmup)
Chain 3: Iteration: 1000 / 2000 [ 50%]  (Warmup)
Chain 3: Iteration: 1001 / 2000 [ 50%]  (Sampling)
Chain 3: Iteration: 1200 / 2000 [ 60%]  (Sampling)
Chain 3: Iteration: 1400 / 2000 [ 70%]  (Sampling)
Chain 3: Iteration: 1600 / 2000 [ 80%]  (Sampling)
Chain 3: Iteration: 1800 / 2000 [ 90%]  (Sampling)
Chain 3: Iteration: 2000 / 2000 [100%]  (Sampling)
Chain 3: 
Chain 3:  Elapsed Time: 3.83 seconds (Warm-up)
Chain 3:                2.556 seconds (Sampling)
Chain 3:                6.386 seconds (Total)
Chain 3: 

SAMPLING FOR MODEL 'lm' NOW (CHAIN 4).
Chain 4: 
Chain 4: Gradient evaluation took 0 seconds
Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
Chain 4: Adjust your expectations accordingly!
Chain 4: 
Chain 4: 
Chain 4: Iteration:    1 / 2000 [  0%]  (Warmup)
Chain 4: Iteration:  200 / 2000 [ 10%]  (Warmup)
Chain 4: Iteration:  400 / 2000 [ 20%]  (Warmup)
Chain 4: Iteration:  600 / 2000 [ 30%]  (Warmup)
Chain 4: Iteration:  800 / 2000 [ 40%]  (Warmup)
Chain 4: Iteration: 1000 / 2000 [ 50%]  (Warmup)
Chain 4: Iteration: 1001 / 2000 [ 50%]  (Sampling)
Chain 4: Iteration: 1200 / 2000 [ 60%]  (Sampling)
Chain 4: Iteration: 1400 / 2000 [ 70%]  (Sampling)
Chain 4: Iteration: 1600 / 2000 [ 80%]  (Sampling)
Chain 4: Iteration: 1800 / 2000 [ 90%]  (Sampling)
Chain 4: Iteration: 2000 / 2000 [100%]  (Sampling)
Chain 4: 
Chain 4:  Elapsed Time: 4 seconds (Warm-up)
Chain 4:                5.935 seconds (Sampling)
Chain 4:                9.935 seconds (Total)
Chain 4: 
Warning messages:
1: There were 13 divergent transitions after warmup. Increasing adapt_delta above 0.99 may help. See
http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup 
2: Examine the pairs() plot to diagnose sampling problems
 
> plot(fit, prob = 0.8)
> 
> # This tests that prophet is working
> history <- data.frame(ds = seq(as.Date('2015-01-01'), as.Date('2016-01-01'), by = 'd'),
+                       y = sin(1:366/200) + rnorm(366)/10)
> m <- prophet(history)
Disabling yearly seasonality. Run prophet with yearly.seasonality=TRUE to override this.
Disabling daily seasonality. Run prophet with daily.seasonality=TRUE to override this.
> future <- make_future_dataframe(m, periods = 365)
> forecast <- predict(m, future)
> plot(m, forecast)
> 
> # This tests that CausalImpact is working
> # First simulating some data
> x1 <- 100 + arima.sim(model = list(ar = 0.999), n = 52)
> y <- 1.2 * x1 + rnorm(52)
> y[41:52] <- y[41:52] + 10
> data <- cbind(y, x1)
> pre.period <- c(1, 40)
> post.period <- c(41, 52)
> # Then running CausalImpact
> impact <- CausalImpact(data, pre.period, post.period)
> plot(impact)
Warning messages:
1: Removed 52 rows containing missing values (geom_path). 
2: Removed 104 rows containing missing values (geom_path). 
> 

@rasmusab
Copy link
Author

rasmusab commented Jul 4, 2019

@serifatf
Copy link

serifatf commented Jul 5, 2019

Thank you Rasmus, my script is working with graphs plotting fine but i understand what am running can you please explain what we are trying to achieve even if there is a document to study.

@rasmusab
Copy link
Author

rasmusab commented Jul 6, 2019

@serifatf The only purpose of this script is to see if you have correctly installed the three packages, nothing else :)

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment