Update lectuers for 3 dec

Author: Lauren Talluto

index.md

@@ -9,6 +9,7 @@ output:
     css: rmd_style.css
     self_contained: false
     lib_dir: lib
+    mathjax: https://cdn.jsdelivr.net/npm/mathjax@3/es5/tex-mml-chtml.js
   beamer_presentation: default
 ---
 
@@ -243,10 +244,9 @@ grpl
 
 ## 2. Sample from the joint posterior
 
-Let's try two models. First, we will see how NDVI affects diversity.
+First we compile the model and prepare the data.
+
 
-::: {.columns}
-:::: {.column}
 
 ```{r eval = FALSE}
 library(rstan)
@@ -270,8 +270,10 @@ X_scaled = scale(X)
 
 ```
 
-::::
-:::: {.column}
+
+## 2. Sample from the joint posterior
+
+Next we fit an initial model using only NVDI as a predictor.
 
 ```{r, cache = TRUE}
 
@@ -283,7 +285,8 @@ standat1 = list(
 	X = X_scaled[which_forest, "NDVI", drop=FALSE],
 	mu_a = 0,
 	mu_b = 0,
-	# these prior scales are still SUPER vague (exp(20) is a possible intercept under this prior!)
+	# these prior scales are still SUPER vague
+	# exp(20) is a possible intercept under this prior!
 	sigma_a = 10,
 	sigma_b = 5
 )
@@ -294,8 +297,7 @@ fit1 = sampling(bird_glm, data = standat1, iter=5000,
 
 ```
 
-::::
-:::
+
 
 ## 2. Sample from the joint posterior
 
@@ -315,9 +317,6 @@ print(fit1, par = c("a", "B"))
 
 Next, we choose two variables, and try them using birds of different habitat types.
 
-::: {.columns}
-:::: {.column}
-
 
 ```{r, message=FALSE, cache = TRUE}
 # Second, looking at how two variables influence birds of different types
@@ -339,8 +338,9 @@ X_2 = cbind(X_2,
 head(X_2)
 ```
 
-::::
-:::: {.column}
+## 2. Sample from the joint posterior
+
+Next, we choose two variables, and try them using birds of different habitat types.
 
 
 ```{r, message=FALSE, cache = TRUE}
@@ -360,8 +360,7 @@ standat2 = list(
 fit2 = sampling(bird_glm, data = standat2, iter=5000, 
    chains = 4, refresh = 0)
 ```
-::::
-:::
+
 
 ## 2. Sample from the joint posterior
 
@@ -608,15 +607,16 @@ quantile(phi, c(0.05, 0.95))
 
 * How does richness respond to the individual variables, holding other variables constant?
 
-```{r, message=FALSE, echo = FALSE, fig.width=10}
+```{r, message=FALSE, echo = FALSE, fig.width=14}
 ## generate a dataset to predict a line for all combinations
 xx = seq(-1.8,1.8, length.out = 200)
-predict_dat = rbind(data.table(forcover = xx, ndvi = 0, open=0, gen = 0, op_fc = 0, op_nd = 0, ge_fc = 0, ge_ndvi = 0), 
-			data.table(forcover = xx, ndvi = 0, open=1, gen = 0, op_fc = xx, op_nd = 0, ge_fc = 0, ge_ndvi = 0),
-			data.table(forcover = xx, ndvi = 0, open=0, gen = 1, op_fc = 0, op_nd = 0, ge_fc = xx, ge_ndvi = 0),
-			data.table(forcover = 0, ndvi = xx, open=0, gen = 0, op_fc = 0, op_nd = 0, ge_fc = 0, ge_ndvi = 0),
-			data.table(forcover = 0, ndvi = xx, open=1, gen = 0, op_fc = 0, op_nd = 1, ge_fc = 0, ge_ndvi = 0),
-			data.table(forcover = 0, ndvi = xx, open=0, gen = 1, op_fc = 0, op_nd = 0, ge_fc = 0, ge_ndvi = 1))
+predict_dat = rbind(
+	data.table(forcover = xx, ndvi = 0, open=0, gen = 0, op_fc = 0, op_nd = 0, ge_fc = 0, ge_ndvi = 0), 
+	data.table(forcover = xx, ndvi = 0, open=1, gen = 0, op_fc = xx, op_nd = 0, ge_fc = 0, ge_ndvi = 0),
+	data.table(forcover = xx, ndvi = 0, open=0, gen = 1, op_fc = 0, op_nd = 0, ge_fc = xx, ge_ndvi = 0),
+	data.table(forcover = 0, ndvi = xx, open=0, gen = 0, op_fc = 0, op_nd = 0, ge_fc = 0, ge_ndvi = 0),
+	data.table(forcover = 0, ndvi = xx, open=1, gen = 0, op_fc = 0, op_nd = 1, ge_fc = 0, ge_ndvi = 0),
+	data.table(forcover = 0, ndvi = xx, open=0, gen = 1, op_fc = 0, op_nd = 0, ge_fc = 0, ge_ndvi = 1))
 predict_dat = cbind(intercept=1, predict_dat)
 
 ## this computes a posterior distribution for E(y) | predict_dat
@@ -631,21 +631,25 @@ predict_dat$y_lower = apply(y, 2, quantile, 0.05)
 
 ## Create some nice labels for ggplot
 predict_dat$bird = with(predict_dat, ifelse(open == 1, "open", ifelse(gen == 1, "generalist", "forest")))
-predict_dat$panel = "Forest Cover | NDVI=0"
-predict_dat$panel[601:1200] = "NDVI | Forest Cover=0"
-
-## Create a combined x-variable for ggplot; this works because forcover
-## and ndvi are transformed to mean = 0, and the predictions are conditional; 
-## if one variable is nonzero, the other must be zero
-predict_dat$x = predict_dat$forcover + predict_dat$ndvi
-
-ggplot(predict_dat, aes(x=x, y=y_med, col=bird)) + 
-	geom_ribbon(aes(x=x, ymin=y_lower, ymax=y_upper, fill=bird), alpha=0.5) + 
-	geom_line() + facet_grid(.~panel) + 
-	theme_minimal() + ylab("Species Richness") + xlab(expression(sigma)) + 
+
+
+p1 = ggplot(predict_dat[1:600,], aes(x=forcover, y=y_med, col=bird)) + 
+	geom_ribbon(aes(x=forcover, ymin=y_lower, ymax=y_upper, fill=bird), alpha=0.5) + 
+	geom_line() + xlab("Forest Cover | NDVI=0 (sd)") + 
+	theme_minimal() + ylab("Species Richness") + 
 	labs(fill="Type of bird", colour="Type of bird") + 
 	xlim(-1.8, 1.8)
 
+p2 = ggplot(predict_dat[601:1200,], aes(x=ndvi, y=y_med, col=bird)) + 
+	geom_ribbon(aes(x=ndvi, ymin=y_lower, ymax=y_upper, fill=bird), alpha=0.5) + 
+	geom_line() + xlab("NDVI | Forest Cover=0 (sd)") + 
+	theme_minimal() + ylab("Species Richness") + 
+	labs(fill="Type of bird", colour="Type of bird") + 
+	xlim(-1.8, 1.8)
+
+grid.arrange(p1, p2, nrow = 1)
+
+
 ```