ANLY505 - Intro to Stan
Week 4
Kirti Tambe
2020-06-18
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
#place the code here
seaice <- read.csv("~/Desktop/ANLY 505 - Modeling Simulation & Game Theory/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.4542. Plot the data
#plot data
plot(extent_south ~ year, pch = 20, col = 'blue',main = "Sea ice extent in the southern hemisphere over time ", data = seaice)
3. Run a general linear model using lm()
model1 <- lm(extent_south ~ year, data = seaice)
summary(model1)##
## 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)
model2 <- lm(y ~ x)
summary(model2)##
## 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(model2)$coeff[1]
lm_beta <- summary(model2)$coeff[2]
lm_sigma <- sigma(model2)
stan_input_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
#code here
stan_model1 <- "stan_model1.stan"
fitmod <- stan(file = stan_model1, data = stan_input_data, warmup = 400, iter = 1500, 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
## clang -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG -I"/Library/Frameworks/R.framework/Versions/3.6/Resources/library/Rcpp/include/" -I"/Library/Frameworks/R.framework/Versions/3.6/Resources/library/RcppEigen/include/" -I"/Library/Frameworks/R.framework/Versions/3.6/Resources/library/RcppEigen/include/unsupported" -I"/Library/Frameworks/R.framework/Versions/3.6/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/3.6/Resources/library/StanHeaders/include/src/" -I"/Library/Frameworks/R.framework/Versions/3.6/Resources/library/StanHeaders/include/" -I"/Library/Frameworks/R.framework/Versions/3.6/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 -isysroot /Library/Developer/CommandLineTools/SDKs/MacOSX.sdk -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.6/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/3.6/Resources/library/RcppEigen/include/Eigen/Dense:1:
## In file included from /Library/Frameworks/R.framework/Versions/3.6/Resources/library/RcppEigen/include/Eigen/Core:88:
## /Library/Frameworks/R.framework/Versions/3.6/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.6/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.6/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/3.6/Resources/library/RcppEigen/include/Eigen/Dense:1:
## /Library/Frameworks/R.framework/Versions/3.6/Resources/library/RcppEigen/include/Eigen/Core:96:10: fatal error: 'complex' file not found
## #include <complex>
## ^~~~~~~~~
## 3 errors generated.
## make: *** [foo.o] Error 1fitmod## Inference for Stan model: stan_model1.
## 4 chains, each with iter=1500; warmup=400; thin=1;
## post-warmup draws per chain=1100, total post-warmup draws=4400.
##
## mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
## alpha 11.42 0.00 0.13 11.16 11.33 11.42 11.50 11.69 1224 1
## beta 0.01 0.00 0.01 0.00 0.01 0.01 0.02 0.02 1603 1
## sigma 0.40 0.00 0.05 0.31 0.36 0.39 0.43 0.50 2062 1
## lp__ 16.25 0.03 1.30 12.86 15.71 16.57 17.20 17.73 1498 1
##
## Samples were drawn using NUTS(diag_e) at Thu Jun 18 17:14:41 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(fitmod)
str(posterior)## List of 4
## $ alpha: num [1:4400(1d)] 11.3 11.3 11.5 11.4 11.3 ...
## ..- attr(*, "dimnames")=List of 1
## .. ..$ iterations: NULL
## $ beta : num [1:4400(1d)] 0.01596 0.01309 0.00702 0.01262 0.01906 ...
## ..- attr(*, "dimnames")=List of 1
## .. ..$ iterations: NULL
## $ sigma: num [1:4400(1d)] 0.392 0.423 0.41 0.319 0.439 ...
## ..- attr(*, "dimnames")=List of 1
## .. ..$ iterations: NULL
## $ lp__ : num [1:4400(1d)] 17.6 16.2 17.1 16.5 16.5 ...
## ..- attr(*, "dimnames")=List of 1
## .. ..$ iterations: NULL8. Compare your results to our results to “lm”
plot(y ~ x, pch = 20, main = "Comparison: lm and STAN")
abline(model2, col = "navy", pch = 22, lty = 2, lw = 3)
abline(mean(posterior$alpha), mean(posterior$beta), col = "red", lw = 1)
legend("topleft", c("Linear Model","Stan Model"), fill = c("blue","red"))
9. Plot multiple estimates from the posterior
plot(y ~ x, pch = 20, main = "lm and All STAN: Comparison")
for (i in 1:500) {
abline(posterior$alpha[i], posterior$beta[i], col = "grey", lty = 1)
}
abline(model2, col = "navy", pch=22, lty = 2, lw = 3)
abline( mean(posterior$alpha), mean(posterior$beta), col = "red", lw = 1)
legend("topleft",c("Linear Model","Mean Stan Model", "All stan Models"), fill = c("blue","red","grey"))