Intro to STAN Homework Set #1

After our Intro to Stan lecture I think it would be valuable to have you go through a similar exercise. Let’s test a second research question.

Research question: Is sea ice extent declining in the Southern Hemisphere over time? Is the same pattern happening in the Antarctic as in the Arctic? Fit a Stan model to find out!

Make sure you follow the steps we used in class.

What do your Stan model results indicate so far?

Load and Inspect Data

library(rstan)
## Warning: package 'rstan' was built under R version 3.5.3
## Loading required package: ggplot2
## Warning: package 'ggplot2' was built under R version 3.5.3
## Loading required package: StanHeaders
## Warning: package 'StanHeaders' was built under R version 3.5.3
## rstan (Version 2.18.2, GitRev: 2e1f913d3ca3)
## For execution on a local, multicore CPU with excess RAM we recommend calling
## options(mc.cores = parallel::detectCores()).
## To avoid recompilation of unchanged Stan programs, we recommend calling
## rstan_options(auto_write = TRUE)
## For improved execution time, we recommend calling
## Sys.setenv(LOCAL_CPPFLAGS = '-march=native')
## although this causes Stan to throw an error on a few processors.
library(gdata)
## Warning: package 'gdata' was built under R version 3.5.2
## gdata: Unable to locate valid perl interpreter
## gdata: 
## gdata: read.xls() will be unable to read Excel XLS and XLSX files
## gdata: unless the 'perl=' argument is used to specify the location
## gdata: of a valid perl intrpreter.
## gdata: 
## gdata: (To avoid display of this message in the future, please
## gdata: ensure perl is installed and available on the executable
## gdata: search path.)
## gdata: Unable to load perl libaries needed by read.xls()
## gdata: to support 'XLX' (Excel 97-2004) files.
## 
## gdata: Unable to load perl libaries needed by read.xls()
## gdata: to support 'XLSX' (Excel 2007+) files.
## 
## gdata: Run the function 'installXLSXsupport()'
## gdata: to automatically download and install the perl
## gdata: libaries needed to support Excel XLS and XLSX formats.
## 
## Attaching package: 'gdata'
## The following object is masked from 'package:stats':
## 
##     nobs
## The following object is masked from 'package:utils':
## 
##     object.size
## The following object is masked from 'package:base':
## 
##     startsWith
library(bayesplot)
## Warning: package 'bayesplot' was built under R version 3.5.3
## This is bayesplot version 1.7.0
## - Online documentation and vignettes at mc-stan.org/bayesplot
## - bayesplot theme set to bayesplot::theme_default()
##    * Does _not_ affect other ggplot2 plots
##    * See ?bayesplot_theme_set for details on theme setting
library(data.table)
## Warning: package 'data.table' was built under R version 3.5.2
## 
## Attaching package: 'data.table'
## The following objects are masked from 'package:gdata':
## 
##     first, last

Let’s load the data and take a look at the data:

seaice <- fread("C:/Users/jings/Desktop/2019spring/ANLY505/seaice.csv")
head(seaice)
##    year extent_north extent_south
## 1: 1979       12.328       11.700
## 2: 1980       12.337       11.230
## 3: 1981       12.127       11.435
## 4: 1982       12.447       11.640
## 5: 1983       12.332       11.389
## 6: 1984       11.910       11.454

Plot the data

To explore the answer to the question “Is sea ice extent declining in the Southern Hemisphere over time?”, first we can make a figure.

plot(extent_north ~ year, pch = 20, data = seaice)

Run a general linear model using lm()

Now, let’s run a general linear model using lm().

lm1 <- lm(extent_north ~ year, data = seaice)
lm1
## 
## Call:
## lm(formula = extent_north ~ year, data = seaice)
## 
## Coefficients:
## (Intercept)         year  
##   120.50304     -0.05457

Prepare the data, re-run the lm() and extract summary statistics

We can add that model fit to our plot:

summary(lm1)
## 
## Call:
## lm(formula = extent_north ~ year, data = seaice)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.49925 -0.17713  0.04898  0.16923  0.32829 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 120.503036   6.267203   19.23   <2e-16 ***
## year         -0.054574   0.003137  -17.40   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2205 on 37 degrees of freedom
## Multiple R-squared:  0.8911, Adjusted R-squared:  0.8881 
## F-statistic: 302.7 on 1 and 37 DF,  p-value: < 2.2e-16
plot(extent_north ~ year, pch = 20, data = seaice)
abline(lm1, col = 2, lty = 2, lw = 3)

Recall linear model: y = ?? + ?????x + error Let’s rename the variables and index the years from 1 to 39, the reason for doing that is because we don’t specifically interested in knowing if the sea ice changing between the years 1979 and 2017. What we really wants to know is that if the sea ice changing from the start of our dataset to the end of our dataset. In this case we want to index the year to 0. We don’t need our model to estimate what sea ice was like in the year 500, or 600, just over the duration of our dataset. So we set up our year data to index from 1 to 30 years.

x <- I(seaice$year - 1978)
y <- seaice$extent_north
N <- length(seaice$year)

Then, we rerun this linear model with our new data.

lm1 <- lm(y ~ x)
summary(lm1)
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.49925 -0.17713  0.04898  0.16923  0.32829 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 12.555888   0.071985   174.4   <2e-16 ***
## x           -0.054574   0.003137   -17.4   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2205 on 37 degrees of freedom
## Multiple R-squared:  0.8911, Adjusted R-squared:  0.8881 
## F-statistic: 302.7 on 1 and 37 DF,  p-value: < 2.2e-16

We can also extract some of the key summary statistics from our simple model, so that we can compare them with the outputs of the Stan models later.

lm_alpha <- summary(lm1)$coeff[1]#intersect
lm_beta <- summary(lm1)$coeff[2]#slope
lm_sigma <- sigma(lm1)#residual error

Now let’s turn that into a dataframe for inputting into a Stan model. Data passed to Stan needs to be a list of named objects. The names given here need to match the variable names used in the models.

stan_data <- list(N = N, x = x, y = y)

Write the Stan model

Now, let’s write stan models with the typical three blocks, which are:

“data” block “parameters” block “model” block

write("// Stan model for simple linear regression

data {
 int < lower = 1 > N; // Sample size
 vector[N] x; // Predictor
 vector[N] y; // Outcome
}

parameters {
 real alpha; // Intercept
 real beta; // Slope (regression coefficients)
 real < lower = 0 > sigma; // Error SD
}

model {
 y ~ normal(alpha + x * beta , sigma);
}

generated quantities {
} // The posterior predictive distribution",

"stan_model1.stan")

First thing after write the model is to check if the written Stan model is a file.

stanc("stan_model1.stan")
## $status
## [1] TRUE
## 
## $model_cppname
## [1] "model37ec544b2dac_stan_model1"
## 
## $cppcode
## [1] "// Code generated by Stan version 2.18.1\n\n#include <stan/model/model_header.hpp>\n\nnamespace model37ec544b2dac_stan_model1_namespace {\n\nusing std::istream;\nusing std::string;\nusing std::stringstream;\nusing std::vector;\nusing stan::io::dump;\nusing stan::math::lgamma;\nusing stan::model::prob_grad;\nusing namespace stan::math;\n\nstatic int current_statement_begin__;\n\nstan::io::program_reader prog_reader__() {\n    stan::io::program_reader reader;\n    reader.add_event(0, 0, \"start\", \"model37ec544b2dac_stan_model1\");\n    reader.add_event(22, 20, \"end\", \"model37ec544b2dac_stan_model1\");\n    return reader;\n}\n\nclass model37ec544b2dac_stan_model1 : public prob_grad {\nprivate:\n    int N;\n    vector_d x;\n    vector_d y;\npublic:\n    model37ec544b2dac_stan_model1(stan::io::var_context& context__,\n        std::ostream* pstream__ = 0)\n        : prob_grad(0) {\n        ctor_body(context__, 0, pstream__);\n    }\n\n    model37ec544b2dac_stan_model1(stan::io::var_context& context__,\n        unsigned int random_seed__,\n        std::ostream* pstream__ = 0)\n        : prob_grad(0) {\n        ctor_body(context__, random_seed__, pstream__);\n    }\n\n    void ctor_body(stan::io::var_context& context__,\n                   unsigned int random_seed__,\n                   std::ostream* pstream__) {\n        typedef double local_scalar_t__;\n\n        boost::ecuyer1988 base_rng__ =\n          stan::services::util::create_rng(random_seed__, 0);\n        (void) base_rng__;  // suppress unused var warning\n\n        current_statement_begin__ = -1;\n\n        static const char* function__ = \"model37ec544b2dac_stan_model1_namespace::model37ec544b2dac_stan_model1\";\n        (void) function__;  // dummy to suppress unused var warning\n        size_t pos__;\n        (void) pos__;  // dummy to suppress unused var warning\n        std::vector<int> vals_i__;\n        std::vector<double> vals_r__;\n        local_scalar_t__ DUMMY_VAR__(std::numeric_limits<double>::quiet_NaN());\n        (void) DUMMY_VAR__;  // suppress unused var warning\n\n        // initialize member variables\n        try {\n            current_statement_begin__ = 4;\n            context__.validate_dims(\"data initialization\", \"N\", \"int\", context__.to_vec());\n            N = int(0);\n            vals_i__ = context__.vals_i(\"N\");\n            pos__ = 0;\n            N = vals_i__[pos__++];\n            current_statement_begin__ = 5;\n            validate_non_negative_index(\"x\", \"N\", N);\n            context__.validate_dims(\"data initialization\", \"x\", \"vector_d\", context__.to_vec(N));\n            validate_non_negative_index(\"x\", \"N\", N);\n            x = vector_d(static_cast<Eigen::VectorXd::Index>(N));\n            vals_r__ = context__.vals_r(\"x\");\n            pos__ = 0;\n            size_t x_i_vec_lim__ = N;\n            for (size_t i_vec__ = 0; i_vec__ < x_i_vec_lim__; ++i_vec__) {\n                x[i_vec__] = vals_r__[pos__++];\n            }\n            current_statement_begin__ = 6;\n            validate_non_negative_index(\"y\", \"N\", N);\n            context__.validate_dims(\"data initialization\", \"y\", \"vector_d\", context__.to_vec(N));\n            validate_non_negative_index(\"y\", \"N\", N);\n            y = vector_d(static_cast<Eigen::VectorXd::Index>(N));\n            vals_r__ = context__.vals_r(\"y\");\n            pos__ = 0;\n            size_t y_i_vec_lim__ = N;\n            for (size_t i_vec__ = 0; i_vec__ < y_i_vec_lim__; ++i_vec__) {\n                y[i_vec__] = vals_r__[pos__++];\n            }\n\n            // validate, data variables\n            current_statement_begin__ = 4;\n            check_greater_or_equal(function__,\"N\",N,1);\n            current_statement_begin__ = 5;\n            current_statement_begin__ = 6;\n            // initialize data variables\n\n\n            // validate transformed data\n\n            // validate, set parameter ranges\n            num_params_r__ = 0U;\n            param_ranges_i__.clear();\n            current_statement_begin__ = 10;\n            ++num_params_r__;\n            current_statement_begin__ = 11;\n            ++num_params_r__;\n            current_statement_begin__ = 12;\n            ++num_params_r__;\n        } catch (const std::exception& e) {\n            stan::lang::rethrow_located(e, current_statement_begin__, prog_reader__());\n            // Next line prevents compiler griping about no return\n            throw std::runtime_error(\"*** IF YOU SEE THIS, PLEASE REPORT A BUG ***\");\n        }\n    }\n\n    ~model37ec544b2dac_stan_model1() { }\n\n\n    void transform_inits(const stan::io::var_context& context__,\n                         std::vector<int>& params_i__,\n                         std::vector<double>& params_r__,\n                         std::ostream* pstream__) const {\n        stan::io::writer<double> writer__(params_r__,params_i__);\n        size_t pos__;\n        (void) pos__; // dummy call to supress warning\n        std::vector<double> vals_r__;\n        std::vector<int> vals_i__;\n\n        if (!(context__.contains_r(\"alpha\")))\n            throw std::runtime_error(\"variable alpha missing\");\n        vals_r__ = context__.vals_r(\"alpha\");\n        pos__ = 0U;\n        context__.validate_dims(\"initialization\", \"alpha\", \"double\", context__.to_vec());\n        double alpha(0);\n        alpha = vals_r__[pos__++];\n        try {\n            writer__.scalar_unconstrain(alpha);\n        } catch (const std::exception& e) { \n            throw std::runtime_error(std::string(\"Error transforming variable alpha: \") + e.what());\n        }\n\n        if (!(context__.contains_r(\"beta\")))\n            throw std::runtime_error(\"variable beta missing\");\n        vals_r__ = context__.vals_r(\"beta\");\n        pos__ = 0U;\n        context__.validate_dims(\"initialization\", \"beta\", \"double\", context__.to_vec());\n        double beta(0);\n        beta = vals_r__[pos__++];\n        try {\n            writer__.scalar_unconstrain(beta);\n        } catch (const std::exception& e) { \n            throw std::runtime_error(std::string(\"Error transforming variable beta: \") + e.what());\n        }\n\n        if (!(context__.contains_r(\"sigma\")))\n            throw std::runtime_error(\"variable sigma missing\");\n        vals_r__ = context__.vals_r(\"sigma\");\n        pos__ = 0U;\n        context__.validate_dims(\"initialization\", \"sigma\", \"double\", context__.to_vec());\n        double sigma(0);\n        sigma = vals_r__[pos__++];\n        try {\n            writer__.scalar_lb_unconstrain(0,sigma);\n        } catch (const std::exception& e) { \n            throw std::runtime_error(std::string(\"Error transforming variable sigma: \") + e.what());\n        }\n\n        params_r__ = writer__.data_r();\n        params_i__ = writer__.data_i();\n    }\n\n    void transform_inits(const stan::io::var_context& context,\n                         Eigen::Matrix<double,Eigen::Dynamic,1>& params_r,\n                         std::ostream* pstream__) const {\n      std::vector<double> params_r_vec;\n      std::vector<int> params_i_vec;\n      transform_inits(context, params_i_vec, params_r_vec, pstream__);\n      params_r.resize(params_r_vec.size());\n      for (int i = 0; i < params_r.size(); ++i)\n        params_r(i) = params_r_vec[i];\n    }\n\n\n    template <bool propto__, bool jacobian__, typename T__>\n    T__ log_prob(vector<T__>& params_r__,\n                 vector<int>& params_i__,\n                 std::ostream* pstream__ = 0) const {\n\n        typedef T__ local_scalar_t__;\n\n        local_scalar_t__ DUMMY_VAR__(std::numeric_limits<double>::quiet_NaN());\n        (void) DUMMY_VAR__;  // suppress unused var warning\n\n        T__ lp__(0.0);\n        stan::math::accumulator<T__> lp_accum__;\n\n        try {\n            // model parameters\n            stan::io::reader<local_scalar_t__> in__(params_r__,params_i__);\n\n            local_scalar_t__ alpha;\n            (void) alpha;  // dummy to suppress unused var warning\n            if (jacobian__)\n                alpha = in__.scalar_constrain(lp__);\n            else\n                alpha = in__.scalar_constrain();\n\n            local_scalar_t__ beta;\n            (void) beta;  // dummy to suppress unused var warning\n            if (jacobian__)\n                beta = in__.scalar_constrain(lp__);\n            else\n                beta = in__.scalar_constrain();\n\n            local_scalar_t__ sigma;\n            (void) sigma;  // dummy to suppress unused var warning\n            if (jacobian__)\n                sigma = in__.scalar_lb_constrain(0,lp__);\n            else\n                sigma = in__.scalar_lb_constrain(0);\n\n\n            // transformed parameters\n\n\n\n            // validate transformed parameters\n\n            const char* function__ = \"validate transformed params\";\n            (void) function__;  // dummy to suppress unused var warning\n\n            // model body\n\n            current_statement_begin__ = 16;\n            lp_accum__.add(normal_log<propto__>(y, add(alpha,multiply(x,beta)), sigma));\n\n        } catch (const std::exception& e) {\n            stan::lang::rethrow_located(e, current_statement_begin__, prog_reader__());\n            // Next line prevents compiler griping about no return\n            throw std::runtime_error(\"*** IF YOU SEE THIS, PLEASE REPORT A BUG ***\");\n        }\n\n        lp_accum__.add(lp__);\n        return lp_accum__.sum();\n\n    } // log_prob()\n\n    template <bool propto, bool jacobian, typename T_>\n    T_ log_prob(Eigen::Matrix<T_,Eigen::Dynamic,1>& params_r,\n               std::ostream* pstream = 0) const {\n      std::vector<T_> vec_params_r;\n      vec_params_r.reserve(params_r.size());\n      for (int i = 0; i < params_r.size(); ++i)\n        vec_params_r.push_back(params_r(i));\n      std::vector<int> vec_params_i;\n      return log_prob<propto,jacobian,T_>(vec_params_r, vec_params_i, pstream);\n    }\n\n\n    void get_param_names(std::vector<std::string>& names__) const {\n        names__.resize(0);\n        names__.push_back(\"alpha\");\n        names__.push_back(\"beta\");\n        names__.push_back(\"sigma\");\n    }\n\n\n    void get_dims(std::vector<std::vector<size_t> >& dimss__) const {\n        dimss__.resize(0);\n        std::vector<size_t> dims__;\n        dims__.resize(0);\n        dimss__.push_back(dims__);\n        dims__.resize(0);\n        dimss__.push_back(dims__);\n        dims__.resize(0);\n        dimss__.push_back(dims__);\n    }\n\n    template <typename RNG>\n    void write_array(RNG& base_rng__,\n                     std::vector<double>& params_r__,\n                     std::vector<int>& params_i__,\n                     std::vector<double>& vars__,\n                     bool include_tparams__ = true,\n                     bool include_gqs__ = true,\n                     std::ostream* pstream__ = 0) const {\n        typedef double local_scalar_t__;\n\n        vars__.resize(0);\n        stan::io::reader<local_scalar_t__> in__(params_r__,params_i__);\n        static const char* function__ = \"model37ec544b2dac_stan_model1_namespace::write_array\";\n        (void) function__;  // dummy to suppress unused var warning\n        // read-transform, write parameters\n        double alpha = in__.scalar_constrain();\n        double beta = in__.scalar_constrain();\n        double sigma = in__.scalar_lb_constrain(0);\n        vars__.push_back(alpha);\n        vars__.push_back(beta);\n        vars__.push_back(sigma);\n\n        // declare and define transformed parameters\n        double lp__ = 0.0;\n        (void) lp__;  // dummy to suppress unused var warning\n        stan::math::accumulator<double> lp_accum__;\n\n        local_scalar_t__ DUMMY_VAR__(std::numeric_limits<double>::quiet_NaN());\n        (void) DUMMY_VAR__;  // suppress unused var warning\n\n        try {\n\n\n\n            // validate transformed parameters\n\n            // write transformed parameters\n            if (include_tparams__) {\n            }\n            if (!include_gqs__) return;\n            // declare and define generated quantities\n\n\n\n            // validate generated quantities\n\n            // write generated quantities\n        } catch (const std::exception& e) {\n            stan::lang::rethrow_located(e, current_statement_begin__, prog_reader__());\n            // Next line prevents compiler griping about no return\n            throw std::runtime_error(\"*** IF YOU SEE THIS, PLEASE REPORT A BUG ***\");\n        }\n    }\n\n    template <typename RNG>\n    void write_array(RNG& base_rng,\n                     Eigen::Matrix<double,Eigen::Dynamic,1>& params_r,\n                     Eigen::Matrix<double,Eigen::Dynamic,1>& vars,\n                     bool include_tparams = true,\n                     bool include_gqs = true,\n                     std::ostream* pstream = 0) const {\n      std::vector<double> params_r_vec(params_r.size());\n      for (int i = 0; i < params_r.size(); ++i)\n        params_r_vec[i] = params_r(i);\n      std::vector<double> vars_vec;\n      std::vector<int> params_i_vec;\n      write_array(base_rng,params_r_vec,params_i_vec,vars_vec,include_tparams,include_gqs,pstream);\n      vars.resize(vars_vec.size());\n      for (int i = 0; i < vars.size(); ++i)\n        vars(i) = vars_vec[i];\n    }\n\n    static std::string model_name() {\n        return \"model37ec544b2dac_stan_model1\";\n    }\n\n\n    void constrained_param_names(std::vector<std::string>& param_names__,\n                                 bool include_tparams__ = true,\n                                 bool include_gqs__ = true) const {\n        std::stringstream param_name_stream__;\n        param_name_stream__.str(std::string());\n        param_name_stream__ << \"alpha\";\n        param_names__.push_back(param_name_stream__.str());\n        param_name_stream__.str(std::string());\n        param_name_stream__ << \"beta\";\n        param_names__.push_back(param_name_stream__.str());\n        param_name_stream__.str(std::string());\n        param_name_stream__ << \"sigma\";\n        param_names__.push_back(param_name_stream__.str());\n\n        if (!include_gqs__ && !include_tparams__) return;\n\n        if (include_tparams__) {\n        }\n\n\n        if (!include_gqs__) return;\n    }\n\n\n    void unconstrained_param_names(std::vector<std::string>& param_names__,\n                                   bool include_tparams__ = true,\n                                   bool include_gqs__ = true) const {\n        std::stringstream param_name_stream__;\n        param_name_stream__.str(std::string());\n        param_name_stream__ << \"alpha\";\n        param_names__.push_back(param_name_stream__.str());\n        param_name_stream__.str(std::string());\n        param_name_stream__ << \"beta\";\n        param_names__.push_back(param_name_stream__.str());\n        param_name_stream__.str(std::string());\n        param_name_stream__ << \"sigma\";\n        param_names__.push_back(param_name_stream__.str());\n\n        if (!include_gqs__ && !include_tparams__) return;\n\n        if (include_tparams__) {\n        }\n\n\n        if (!include_gqs__) return;\n    }\n\n}; // model\n\n}\n\ntypedef model37ec544b2dac_stan_model1_namespace::model37ec544b2dac_stan_model1 stan_model;\n\n"
## 
## $model_name
## [1] "stan_model1"
## 
## $model_code
## [1] "// Stan model for simple linear regression\n\ndata {\n int < lower = 1 > N; // Sample size\n vector[N] x; // Predictor\n vector[N] y; // Outcome\n}\n\nparameters {\n real alpha; // Intercept\n real beta; // Slope (regression coefficients)\n real < lower = 0 > sigma; // Error SD\n}\n\nmodel {\n y ~ normal(alpha + x * beta , sigma);\n}\n\ngenerated quantities {\n} // The posterior predictive distribution"
## attr(,"model_name2")
## [1] "stan_model1"

Then, we can save that file path.

stan_model1 <- "stan_model1.stan"

Run the Stan model

fit <- stan(file = stan_model1, data = stan_data, warmup = 500, iter = 1000, chains = 4, cores = 2, thin = 1)

fit
## Inference for Stan model: stan_model1.
## 4 chains, each with iter=1000; warmup=500; thin=1; 
## post-warmup draws per chain=500, total post-warmup draws=2000.
## 
##        mean se_mean   sd  2.5%   25%   50%   75% 97.5% n_eff Rhat
## alpha 12.56    0.00 0.08 12.41 12.51 12.56 12.61 12.71   852    1
## beta  -0.05    0.00 0.00 -0.06 -0.06 -0.05 -0.05 -0.05   926    1
## sigma  0.23    0.00 0.03  0.18  0.21  0.22  0.24  0.29  1052    1
## lp__  37.41    0.05 1.29 34.14 36.84 37.76 38.36 38.87   670    1
## 
## Samples were drawn using NUTS(diag_e) at Mon Jun 10 23:52:47 2019.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at 
## convergence, Rhat=1).

From this output we can quickly assess model convergence by looking at the Rhat values for each parameter. When these are at or near 1, the chains have converged.

We can also look at the full posterior of our parameters by extracting them from the model object.

posterior <- extract(fit)
str(posterior)
## List of 4
##  $ alpha: num [1:2000(1d)] 12.6 12.6 12.5 12.6 12.5 ...
##   ..- attr(*, "dimnames")=List of 1
##   .. ..$ iterations: NULL
##  $ beta : num [1:2000(1d)] -0.056 -0.0559 -0.0537 -0.0566 -0.049 ...
##   ..- attr(*, "dimnames")=List of 1
##   .. ..$ iterations: NULL
##  $ sigma: num [1:2000(1d)] 0.216 0.21 0.25 0.215 0.232 ...
##   ..- attr(*, "dimnames")=List of 1
##   .. ..$ iterations: NULL
##  $ lp__ : num [1:2000(1d)] 38.8 38.8 38 37.6 36.3 ...
##   ..- attr(*, "dimnames")=List of 1
##   .. ..$ iterations: NULL

Compare your results to our results to “lm”

Let’s compare to our previous estimate with “lm”:

plot(y ~ x, pch = 20)

abline(lm1, col = 2, lty = 2, lw = 3)
abline( mean(posterior$alpha), mean(posterior$beta), col = 6, lw = 2)

The result is identical to the lm output. This is because we are using a simple model, and have put non-informative priors on our parameters.

One way to visualize the variability in our estimation of the regression line is to plot multiple estimates from the posterior.

plot(y ~ x, pch = 20)

for (i in 1:500) {
 abline(posterior$alpha[i], posterior$beta[i], col = "gray", lty = 1)
}

abline(mean(posterior$alpha), mean(posterior$beta), col = 6, lw = 2)

## Changing Our Prior

write("// Stan model for simple linear regression

data {
 int < lower = 1 > N; // Sample size
 vector[N] x; // Predictor
 vector[N] y; // Outcome
}

parameters {
 real alpha; // Intercept
 real beta; // Slope (regression coefficients)
 real < lower = 0 > sigma; // Error SD
}

model {
 alpha ~ normal(10, 0.1);
 beta ~ normal(1, 0.1);
 y ~ normal(alpha + x * beta , sigma);
}

generated quantities {}",

"stan_model2.stan")
stan_model2 <- "stan_model2.stan"
for (i in 1:500) {
  abline(posterior$alpha[i], posterior$beta[i], col="gray", lty=1)
} 
fit2 <- stan(stan_model2, data = stan_data, warmup = 500, iter = 1000, chains = 4, cores = 2, thin = 1)

posterior2 <- extract(fit2)


plot(y ~ x, pch = 20)
abline(lm_alpha, lm_beta, col = 4, lty = 2, lw = 2)
abline(mean(posterior2$alpha), mean(posterior2$beta), col = 3, lw = 2)
abline(mean(posterior$alpha), mean(posterior$beta), col = 36, lw = 3)

Now Let’s do a convergence diagnostic

For traceplots, we can view them directly from the posterior:

plot(posterior$alpha, type = "l")

plot(posterior$beta, type = "l")

plot(posterior$sigma, type = "l")

For poor convergence Try running a model for only 50 iterations and check the traceplots.

fit_bad <- stan(stan_model1, data = stan_data, warmup = 25, iter = 50, chains = 4, cores = 2, thin = 1)
## Warning: There were 11 divergent transitions after warmup. Increasing adapt_delta above 0.8 may help. See
## http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
## Warning: There were 4 chains where the estimated Bayesian Fraction of Missing Information was low. See
## http://mc-stan.org/misc/warnings.html#bfmi-low
## Warning: Examine the pairs() plot to diagnose sampling problems
posterior_bad <- extract(fit_bad)
plot(posterior_bad$alpha, type = "l")

plot(posterior_bad$beta, type = "l")

plot(posterior_bad$sigma, type = "l")

We can get summaries of the parameters through the posterior directly. Let’s also plot the non-Bayesian linear model values to make sure our model is doing what we think it is.

par(mfrow = c(1,3))

plot(density(posterior$alpha), main = "Alpha")
abline(v = lm_alpha, col = 4, lty = 2)

plot(density(posterior$beta), main = "Beta")
abline(v = lm_beta, col = 4, lty = 2)

plot(density(posterior$sigma), main = "Sigma")
abline(v = lm_sigma, col = 4, lty = 2)

From the posterior we can directly calculate the probability of any parameter being over or under a certain value of interest.

sum(posterior$beta>0)/length(posterior$beta)
## [1] 0
sum(posterior$beta>0.2)/length(posterior$beta)
## [1] 0

While we can work with the posterior directly, rstan has a lot of useful functions built-in.

traceplot(fit)

We can also look at the posterior densities & histograms.

stan_dens(fit)

stan_hist(fit)
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

And we can generate plots which indicate the mean parameter estimates and any credible intervals we may be interested in. Note that the 95% credible intervals for the beta and sigma parameters are very small, thus you only see the dots.

plot(fit, show_density = FALSE, ci_level = 0.5, outer_level = 0.95, fill_color = "salmon")
## ci_level: 0.5 (50% intervals)
## outer_level: 0.95 (95% intervals)

We’ll fit this model and compare it to the mean estimate using the uniform priors.

Posterior Predictive Checks

write("// Stan model for simple linear regression

data {
 int < lower = 1 > N; // Sample size
 vector[N] x; // Predictor
 vector[N] y; // Outcome
}

parameters {
 real alpha; // Intercept
 real beta; // Slope (regression coefficients)
 real < lower = 0 > sigma; // Error SD
}

model {
 y ~ normal(x * beta + alpha, sigma);
}

generated quantities {
 real y_rep[N];

 for (n in 1:N) {
 y_rep[n] = normal_rng(x[n] * beta + alpha, sigma);
 }

}",

"stan_model2_GQ.stan")

stan_model2_GQ <- "stan_model2_GQ.stan"
fit3 <- stan(stan_model2_GQ, data = stan_data, iter = 1000, chains = 4, cores = 2, thin = 1)

y_rep <- as.matrix(fit3, pars = "y_rep")
dim(y_rep)
## [1] 2000   39

Comparing density of y with densities of y over 200 posterior draws.

ppc_dens_overlay(y, y_rep[1:200, ])

Here we see data (dark blue) fit well with our posterior predictions.

We can also use this to compare estimates of summary statistics.

ppc_stat(y = y, yrep = y_rep, stat = "mean")
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

We can investigate mean posterior prediction per datapoint vs the observed value for each datapoint

ppc_scatter_avg(y = y, yrep = y_rep)

Here is a list of currently available plots

available_ppc()
## bayesplot PPC module:
##   ppc_bars
##   ppc_bars_grouped
##   ppc_boxplot
##   ppc_data
##   ppc_dens
##   ppc_dens_overlay
##   ppc_ecdf_overlay
##   ppc_error_binned
##   ppc_error_hist
##   ppc_error_hist_grouped
##   ppc_error_scatter
##   ppc_error_scatter_avg
##   ppc_error_scatter_avg_vs_x
##   ppc_freqpoly
##   ppc_freqpoly_grouped
##   ppc_hist
##   ppc_intervals
##   ppc_intervals_data
##   ppc_intervals_grouped
##   ppc_loo_intervals
##   ppc_loo_pit
##   ppc_loo_pit_overlay
##   ppc_loo_pit_qq
##   ppc_loo_ribbon
##   ppc_ribbon
##   ppc_ribbon_data
##   ppc_ribbon_grouped
##   ppc_rootogram
##   ppc_scatter
##   ppc_scatter_avg
##   ppc_scatter_avg_grouped
##   ppc_stat
##   ppc_stat_2d
##   ppc_stat_freqpoly_grouped
##   ppc_stat_grouped
##   ppc_violin_grouped
color_scheme_view(c("blue", "gray", "green", "pink", "purple",
 "red","teal","yellow"))

color_scheme_view("mix-blue-red")