ANLY505 - Intro to Stan

Week 4

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("D:\\seaice.csv")
head(seaice)

##   锘縴ear 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 ~ 锘縴ear, pch = 20, data = seaice)

3. Run a general linear model using lm()

#write the code

lm1 <- lm(extent_north~锘縴ear, data=seaice)
lm1

## 
## Call:
## lm(formula = extent_north ~ 锘縴ear, data = seaice)
## 
## Coefficients:
## (Intercept)      锘縴ear  
##   120.50304     -0.05457

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$锘縴ear - 1978)
y <- seaice$extent_south
N <- length(seaice$year)

lm1 <- lm(y~x)
summary(lm1)

## 
## 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.02318

plot(extent_north ~ 锘縴ear, pch = 20, data = seaice)
abline(lm1, col = 2, lty = 2, lw = 3)

5. Write the Stan model

#write the code

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

 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

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

7. Extract the posterior estimates into a list so we can plot them

#code here

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

#code here

9. Plot multiple estimates from the posterior

#code here