ANLY505 - Intro to Stan
Week 4
Shuai Wang
2020-03-23
Intro to STAN Homework Part #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?
1. Load and Inspect Data
seaice <- read.csv("seaice.csv")
colnames(seaice) <- c("year", "extent_north", "extent_south")2. Plot the data
plot(extent_south ~ year, pch = 20, data = seaice)
3. Run a general linear model using lm()
lm1 <- lm(extent_south ~ year, data = seaice)
summary(lm1)##
## Call:
## lm(formula = extent_south ~ year, data = seaice)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.23372 -0.18142 0.01587 0.18465 0.88814
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -14.199551 10.925576 -1.300 0.2018
## year 0.012953 0.005468 2.369 0.0232 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3843 on 37 degrees of freedom
## Multiple R-squared: 0.1317, Adjusted R-squared: 0.1082
## F-statistic: 5.611 on 1 and 37 DF, p-value: 0.023184. Index the data, re-run the lm(), extract summary statistics and turn the indexed data into a dataframe to pass into Stan
x <- I(seaice$year - 1978)
y <- seaice$extent_south
N <- length(seaice$year)
lm2 <- lm(y ~ x)
summary(lm2)##
## Call:
## lm(formula = y ~ x)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.23372 -0.18142 0.01587 0.18465 0.88814
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 11.421555 0.125490 91.015 <2e-16 ***
## x 0.012953 0.005468 2.369 0.0232 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3843 on 37 degrees of freedom
## Multiple R-squared: 0.1317, Adjusted R-squared: 0.1082
## F-statistic: 5.611 on 1 and 37 DF, p-value: 0.02318lm_alpha <- summary(lm1)$coeff[1]
lm_beta <- summary(lm1)$coeff[2]
lm_sigma <- sigma(lm1)
stan_data <- list(N = N, x = x, y = y)5. Write the Stan model
#write the code
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")
stan_model1 <- "stan_model1.stan"6. Run the Stan model and inspect the results
fit <- stan(file = stan_model1, data = stan_data, warmup = 500, iter = 1000, chains = 4, cores = 2, thin = 1)## Trying to compile a simple C file## Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
## During startup - Warning messages:
## 1: Setting LC_CTYPE failed, using "C"
## 2: Setting LC_COLLATE failed, using "C"
## 3: Setting LC_TIME failed, using "C"
## 4: Setting LC_MESSAGES failed, using "C"
## 5: Setting LC_MONETARY failed, using "C"
## clang -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG -I"/Library/Frameworks/R.framework/Versions/3.5/Resources/library/Rcpp/include/" -I"/Library/Frameworks/R.framework/Versions/3.5/Resources/library/RcppEigen/include/" -I"/Library/Frameworks/R.framework/Versions/3.5/Resources/library/RcppEigen/include/unsupported" -I"/Library/Frameworks/R.framework/Versions/3.5/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/src/" -I"/Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/" -I"/Library/Frameworks/R.framework/Versions/3.5/Resources/library/rstan/include" -DEIGEN_NO_DEBUG -D_REENTRANT -DBOOST_DISABLE_ASSERTS -DBOOST_PENDING_INTEGER_LOG2_HPP -include stan/math/prim/mat/fun/Eigen.hpp -I/usr/local/include -fPIC -Wall -g -O2 -c foo.c -o foo.o
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/RcppEigen/include/Eigen/Dense:1:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/RcppEigen/include/Eigen/Core:88:
## /Library/Frameworks/R.framework/Versions/3.5/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:613:1: error: unknown type name 'namespace'
## namespace Eigen {
## ^
## /Library/Frameworks/R.framework/Versions/3.5/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:613:16: error: expected ';' after top level declarator
## namespace Eigen {
## ^
## ;
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/RcppEigen/include/Eigen/Dense:1:
## /Library/Frameworks/R.framework/Versions/3.5/Resources/library/RcppEigen/include/Eigen/Core:96:10: fatal error: 'complex' file not found
## #include <complex>
## ^~~~~~~~~
## 3 errors generated.
## make: *** [foo.o] Error 1fit## 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 11.42 0.00 0.13 11.16 11.34 11.43 11.51 11.68 833 1
## beta 0.01 0.00 0.01 0.00 0.01 0.01 0.02 0.02 918 1
## sigma 0.40 0.00 0.05 0.32 0.37 0.40 0.43 0.50 1000 1
## lp__ 16.33 0.05 1.26 13.13 15.76 16.67 17.23 17.72 696 1
##
## Samples were drawn using NUTS(diag_e) at Mon Mar 23 20:52:32 2020.
## 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).7. Extract the posterior estimates into a list so we can plot them
posterior <- extract(fit)
str(posterior)## List of 4
## $ alpha: num [1:2000(1d)] 11.5 11.3 11.4 11.3 11.6 ...
## ..- attr(*, "dimnames")=List of 1
## .. ..$ iterations: NULL
## $ beta : num [1:2000(1d)] 0.01228 0.01812 0.00854 0.01793 0.00561 ...
## ..- attr(*, "dimnames")=List of 1
## .. ..$ iterations: NULL
## $ sigma: num [1:2000(1d)] 0.384 0.466 0.355 0.344 0.407 ...
## ..- attr(*, "dimnames")=List of 1
## .. ..$ iterations: NULL
## $ lp__ : num [1:2000(1d)] 17.4 15.8 16.7 16.4 16.6 ...
## ..- attr(*, "dimnames")=List of 1
## .. ..$ iterations: NULL8. Compare your results to our results to “lm”
par(mfrow=c(1,2))
plot(extent_south ~ year, data= seaice, pch = 20)
abline(lm1 , col = 6, lw = 2)
plot(y~x , pch = 20)
abline(mean(posterior$alpha), mean(posterior$beta), col = 6, lw = 2)
9. Plot multiple estimates from the posterior
plot(y ~ x, pch = 20)
for (i in 1:500) {
abline(posterior$alpha[i], posterior$beta[i], col = "blue", lty =1)
}
abline(mean(posterior$alpha), mean(posterior$beta), col = 6, lw = 2)