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("/Users/XiaoZhou/Downloads/seaice.csv", stringsAsFactors=F)
colnames(seaice) <- c("year", "extent_north", "extent_south")
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

2. Plot the data

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

3. Run a general linear model using lm()

#write the code
lm1 <- lm(extent_north ~ year, data = seaice)
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

4. Index the data, re-run the lm(), extract summary statistics and turn the indexed data into a dataframe to pass into Stan

#write the code here
x <- I(seaice$year - 1978)
y <- seaice$extent_north
N <- length(seaice$year)

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
lm_alpha <- summary(lm1)$coeff[1]
lm_beta <- summary(lm1)$coeff[2]
lm_sigma <- summary(lm1)$sigma

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

#code here
detectCores(all.tests = FALSE, logical = TRUE)
## [1] 4
options(mc.cores = parallel::detectCores())

fit <- stan(file = stan_model1, data = stan_data, warmup = 500, iter = 1000, chains = 4, cores = 4, 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.07 12.41 12.51 12.56 12.61 12.70   941 1.00
## beta  -0.05    0.00 0.00 -0.06 -0.06 -0.05 -0.05 -0.05  1000 1.00
## sigma  0.23    0.00 0.03  0.18  0.21  0.23  0.25  0.29   990 1.01
## lp__  37.41    0.05 1.22 34.41 36.85 37.69 38.31 38.83   707 1.01
## 
## Samples were drawn using NUTS(diag_e) at Wed Jan 15 13:08:43 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

#code here
posterior <- rstan::extract(fit)
str(posterior)
## List of 4
##  $ alpha: num [1:2000(1d)] 12.5 12.6 12.6 12.6 12.5 ...
##   ..- attr(*, "dimnames")=List of 1
##   .. ..$ iterations: NULL
##  $ beta : num [1:2000(1d)] -0.0498 -0.0585 -0.0579 -0.0536 -0.0533 ...
##   ..- attr(*, "dimnames")=List of 1
##   .. ..$ iterations: NULL
##  $ sigma: num [1:2000(1d)] 0.193 0.199 0.256 0.179 0.313 ...
##   ..- attr(*, "dimnames")=List of 1
##   .. ..$ iterations: NULL
##  $ lp__ : num [1:2000(1d)] 36.2 37.7 37.6 36.3 34.6 ...
##   ..- attr(*, "dimnames")=List of 1
##   .. ..$ iterations: NULL

8. Compare your results to our results to “lm”

plot(y ~ x, pch = 20)

abline(lm1, col = 2, lty = 2, lw = 3)
abline( mean(posterior$alpha), mean(posterior$beta), col = 5, 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 = "gray", lty =1)
}
abline(mean(posterior$alpha), mean(posterior$beta), col = 6, lw = 2)