Add files via upload

aammd · Jun 21, 2020 · 37d010d · 37d010d
1 parent dac8827
commit 37d010d
Show file tree

Hide file tree

Showing 9 changed files with 455 additions and 54 deletions.
diff --git a/StanCourseVLB/index.Rmd b/StanCourseVLB/index.Rmd
@@ -345,7 +345,7 @@ model{
                     x[t]*beta));
   } else {
     y[t] ~ bernoulli_logit(alpha[1] +
-                    alpha[2]*sex[t]
+                    alpha[2]*sex[t] + 
                     alpha[3]*bsize[t] +
                     x[t]*beta);
   }
@@ -385,13 +385,13 @@ Changes: have an index for individual shark and $\beta_1$ is now a hierarchical
 ```{stan, output.var = "ex4", eval=F}
 data{
    ...
-  int no.sharks;
+  int nsharks;
   vector[TT] sharkid;
 }
 
 parameters{
   ...
-  vector[no.sharks] beta;
+  vector[nsharks] beta;
   real betamu;
   real<lower=0> betasig;
 }
@@ -466,7 +466,7 @@ https://mc-stan.org/projpred/articles/quickstart.html
 Back to the fixed effects for a minute: 
 
 ```{r, eval=F}
-n <- 4000
+n <- 7300
 D <- 5
 p0 <- 2 # prior guess for the number of relevant variables
 tau0 <- p0/(D-p0) * 1/sqrt(n) 
@@ -477,7 +477,7 @@ prior_coeff <- hs(global_scale = tau0, slab_scale = 1)
 logregcov.horseshoe <- stan_glm(formula = Presence~., 
                       family = binomial(link="logit"), 
                       data = sharkfun, 
-                      prior = hs(global_scale = 0.03049858, 
+                      prior = hs(global_scale = 0.007802743, 
                                  slab_scale = 1))
 
 logregcov.defaultpriors <- stan_glm(formula = Presence~., 

diff --git a/StanCourseVLB/index.html b/StanCourseVLB/index.html
@@ -330,7 +330,7 @@
                     x[t]*beta));
   } else {
     y[t] ~ bernoulli_logit(alpha[1] +
-                    alpha[2]*sex[t]
+                    alpha[2]*sex[t] + 
                     alpha[3]*bsize[t] +
                     x[t]*beta);
   }
@@ -371,13 +371,13 @@
 ```stan
 data{
    ...
-  int no.sharks;
+  int nsharks;
   vector[TT] sharkid;
 }
 
 parameters{
   ...
-  vector[no.sharks] beta;
+  vector[nsharks] beta;
   real betamu;
   real&lt;lower=0&gt; betasig;
 }
@@ -455,7 +455,7 @@
 
 
 ```r
-n &lt;- 4000
+n &lt;- 7300
 D &lt;- 5
 p0 &lt;- 2 # prior guess for the number of relevant variables
 tau0 &lt;- p0/(D-p0) * 1/sqrt(n) 
@@ -466,7 +466,7 @@
 logregcov.horseshoe &lt;- stan_glm(formula = Presence~., 
                       family = binomial(link="logit"), 
                       data = sharkfun, 
-                      prior = hs(global_scale = 0.03049858, 
+                      prior = hs(global_scale = 0.007802743, 
                                  slab_scale = 1))
 
 logregcov.defaultpriors &lt;- stan_glm(formula = Presence~., 

diff --git a/StanCourseVLB/photos/realworld_abstract_chart.png b/StanCourseVLB/photos/realworld_abstract_chart.png
diff --git a/StanCourseVLB/photos/stanforums.png b/StanCourseVLB/photos/stanforums.png
diff --git a/StanCourseVLB/runningstancode_vlb1.Rmd b/StanCourseVLB/runningstancode_vlb1.Rmd
@@ -34,7 +34,7 @@ alpha[3] <- 0.01
 ## number of `sharks` = 20
 
 sex <- rbinom(20, size=1, prob=0.5)
-size <- rpois(20, lambda=200) ## in cm
+size <- rpois(20, lambda=220) ## in cm
 
 y <- list()
 for(n in 1:20){
@@ -140,19 +140,21 @@ prior_coeff <- hs(global_scale = tau0, slab_scale = 1)
 #----------------------------------------------------------------------
 ## FITTING THE MODELS
 
-logregcov.horseshoe <- stan_glm(formula = Presence~., 
+logregcov.horseshoe <- stan_glm(formula = Presence ~ Sex + BSize + timecos + timesin, 
                       family = binomial(link="logit"), 
-                      data = sharkfun%>%dplyr::select(-Shark), 
+                      data = sharkfun%>%dplyr::select(-Shark)%>%mutate(BSize = Size - 220), 
                       prior = hs(global_scale = 0.007802743, 
                                  slab_scale = 1))
 
-logregcov.defaultpriors <- stan_glm(formula = Presence~., 
+#launch_shinystan(logregcov.horseshoe)
+
+logregcov.defaultpriors <- stan_glm(formula = Presence ~ Sex + BSize + timecos + timesin, 
                       family = binomial(link="logit"), 
-                      data = sharkfun%>%dplyr::select(-Shark))
+                      data = sharkfun%>%dplyr::select(-Shark)%>%mutate(BSize = Size - 220))
 
-logregcov.n01 <- stan_glm(formula = Presence~., 
+logregcov.n01 <- stan_glm(formula = Presence ~ Sex + BSize + timecos + timesin, 
                       family = binomial(link="logit"), 
-                      data = sharkfun%>%dplyr::select(-Shark), 
+                      data = sharkfun%>%dplyr::select(-Shark)%>%mutate(BSize = Size - 220), 
                       prior = normal(0,1))
 
 
@@ -164,14 +166,14 @@ cvs.horseshoe <- cv_varsel(logregcov.horseshoe,
                            cv_method='kfold', K=5)
 
 cvplot.horseshoe <- varsel_plot(cvs.horseshoe, 
-                                stats=c('elpd, acc'))  + 
+                                stats=c('elpd', 'acc'))  + 
   theme_minimal() + theme(text=element_text(size=15), 
                           legend.position = "none") + 
   ggtitle("Classification Accuracy (5-fold CV)")
 
 cvplot.horseshoe
 
-suggest_size(cvs)
+suggest_size(cvs.horseshoe)
 
 vs <- varsel(logregcov.horseshoe, method='forward')
 vs$vind
@@ -184,7 +186,38 @@ cvs.horseshoe$vind
 ```
 
 
+## Trying to have the covariates on approximately the same scale...let's see how that works out
+```{r}
+
+logregscaledcov.horseshoe <- stan_glm(formula = Presence ~ Sex + BSize + timecos + timesin, 
+                      family = binomial(link="logit"), 
+                      data = sharkfun%>%dplyr::select(-Shark)%>%mutate(BSize = (Size - 220)/30), 
+                      prior = hs(global_scale = 0.007802743, 
+                                 slab_scale = 1))
+
+#----------------------------------------------------------------------
+## K-FOLD CROSS-VALIDATION
+
+cvs.horseshoe2 <- cv_varsel(logregscaledcov.horseshoe, 
+                           method='forward', 
+                           cv_method='kfold', K=5)
+
+cvplot.horseshoe2 <- varsel_plot(cvs.horseshoe2, 
+                                stats=c('elpd', 'acc'))  + 
+  theme_minimal() + theme(text=element_text(size=15), 
+                          legend.position = "none") + 
+  ggtitle("Classification Accuracy (5-fold CV)")
+
+cvplot.horseshoe2
+
+suggest_size(cvs.horseshoe2)
 
+vs2 <- varsel(logregscaledcov.horseshoe, method='forward')
+vs2$vind
+varsel_plot(vs2, stats=c('elpd', 'acc'), deltas=F)
 
+#----------
 
+cvs.horseshoe2$vind
 
+```
diff --git a/StanCourseVLB/stanhour3.Rmd b/StanCourseVLB/stanhour3.Rmd
@@ -111,9 +111,9 @@ class: center
 
 Combinations: 
 
-( $\nu_1 = 1, \nu_2=1$ ), ( $\nu_1=2, \nu_2=1$ ) , ( $\alpha=1, \beta=2$ )
+( $\nu_1 = 1, \nu_2=1$ ), ( $\nu_1=2, \nu_2=1$ ) , ( $\nu_1=1, \nu_2=2$ )
 
-( $\alpha=10, \beta=1$ ), ( $\alpha=1, \beta=10$ ), ( $\alpha=10, \beta=10$ )
+( $\nu_1=10, \nu_2=1$ ), ( $\nu_1=1, \nu_2=10$ ), ( $\nu_1=10, \nu_2=10$ )
 
 ```{r, echo=FALSE, fig.height=5}
 
@@ -415,8 +415,10 @@ model{
  
  beta ~ normal(betamu, betasig);
 
- y[t] ~ bernoulli_logit(alpha[1] + alpha[2]*sex[t] +
-                       alpha[3]*bsize[t] + x[t]*beta[sharkid]);
+ for(t in 1:TT){
+   y[t] ~ bernoulli_logit(alpha[1] + alpha[2]*sex[t] +
+                       alpha[3]*bsize[t] + x[t]*beta[sharkid[t]]);
+ }
 
 }
 
@@ -470,24 +472,23 @@ n <- 7300
 D <- 5
 p0 <- 2 # prior guess for the number of relevant variables
 tau0 <- p0/(D-p0) * 1/sqrt(n) 
-# regularized horseshoe prior
-prior_coeff <- hs(global_scale = tau0, slab_scale = 1) 
-
+#prior_coeff <- hs(global_scale = tau0, slab_scale = 1) 
 
 logregcov.horseshoe <- stan_glm(formula = Presence~., 
                       family = binomial(link="logit"), 
-                      data = sharkfun, 
-                      prior = hs(global_scale = 0.007802743, 
-                                 slab_scale = 1))
+                      data = sharkfun%>%dplyr::select(-Shark)%>%
+                        mutate(BSize = Size - 220), 
+                      prior = hs(global_scale = 0.007802743, slab_scale = 1))
 
 logregcov.defaultpriors <- stan_glm(formula = Presence~., 
                       family = binomial(link="logit"), 
-                      data = sharkfun)
+                      data = sharkfun%>%dplyr::select(-Shark)%>%
+                        mutate(BSize = Size - 220))
 
 logregcov.n01 <- stan_glm(formula = Presence~., 
                       family = binomial(link="logit"), 
-                      data = sharkfun, 
-                      prior = normal(0,1))
+                      data = sharkfun%>%%>%dplyr::select(-Shark)
+                       mutate(BSize = Size - 220), prior = normal(0,1))
 
 
 ```

diff --git a/StanCourseVLB/stanhour3.html b/StanCourseVLB/stanhour3.html
@@ -1,5 +1,5 @@
 <!DOCTYPE html>
-<html xmlns="http://www.w3.org/1999/xhtml" lang="" xml:lang="">
+<html lang="" xml:lang="">
   <head>
     <title>Stan @ ISEC 2020</title>
     <meta charset="utf-8" />
@@ -107,9 +107,9 @@
 
 Combinations: 
 
-( `\(\nu_1 = 1, \nu_2=1\)` ), ( `\(\nu_1=2, \nu_2=1\)` ) , ( `\(\alpha=1, \beta=2\)` )
+( `\(\nu_1 = 1, \nu_2=1\)` ), ( `\(\nu_1=2, \nu_2=1\)` ) , ( `\(\nu_1=1, \nu_2=2\)` )
 
-( `\(\alpha=10, \beta=1\)` ), ( `\(\alpha=1, \beta=10\)` ), ( `\(\alpha=10, \beta=10\)` )
+( `\(\nu_1=10, \nu_2=1\)` ), ( `\(\nu_1=1, \nu_2=10\)` ), ( `\(\nu_1=10, \nu_2=10\)` )
 
 ![](stanhour3_files/figure-html/unnamed-chunk-3-1.png)&lt;!-- --&gt;
 
@@ -330,7 +330,7 @@
                     x[t]*beta));
   } else {
     y[t] ~ bernoulli_logit(alpha[1] +
-                    alpha[2]*sex[t]
+                    alpha[2]*sex[t] + 
                     alpha[3]*bsize[t] +
                     x[t]*beta);
   }
@@ -371,13 +371,13 @@
 ```stan
 data{
    ...
-  int no.sharks;
+  int nsharks;
   vector[TT] sharkid;
 }
 
 parameters{
   ...
-  vector[no.sharks] beta;
+  vector[nsharks] beta;
   real betamu;
   real&lt;lower=0&gt; betasig;
 }
@@ -402,8 +402,10 @@
 
  beta ~ normal(betamu, betasig);
 
- y[t] ~ bernoulli_logit(alpha[1] + alpha[2]*sex[t] +
-                       alpha[3]*bsize[t] + x[t]*beta[sharkid]);
+ for(t in 1:TT){
+   y[t] ~ bernoulli_logit(alpha[1] + alpha[2]*sex[t] +
+                       alpha[3]*bsize[t] + x[t]*beta[sharkid[t]]);
+ }
 
 }
 
@@ -455,28 +457,27 @@
 
 
 ```r
-n &lt;- 4000
+n &lt;- 7300
 D &lt;- 5
 p0 &lt;- 2 # prior guess for the number of relevant variables
 tau0 &lt;- p0/(D-p0) * 1/sqrt(n) 
-# regularized horseshoe prior
-prior_coeff &lt;- hs(global_scale = tau0, slab_scale = 1) 
-
+#prior_coeff &lt;- hs(global_scale = tau0, slab_scale = 1) 
 
 logregcov.horseshoe &lt;- stan_glm(formula = Presence~., 
                       family = binomial(link="logit"), 
-                      data = sharkfun, 
-                      prior = hs(global_scale = 0.03049858, 
-                                 slab_scale = 1))
+                      data = sharkfun%&gt;%dplyr::select(-Shark)%&gt;%
+                        mutate(BSize = Size - 220), 
+                      prior = hs(global_scale = 0.007802743, slab_scale = 1))
 
 logregcov.defaultpriors &lt;- stan_glm(formula = Presence~., 
                       family = binomial(link="logit"), 
-                      data = sharkfun)
+                      data = sharkfun%&gt;%dplyr::select(-Shark)%&gt;%
+                        mutate(BSize = Size - 220))
 
 logregcov.n01 &lt;- stan_glm(formula = Presence~., 
                       family = binomial(link="logit"), 
-                      data = sharkfun, 
-                      prior = normal(0,1))
+                      data = sharkfun%&gt;%%&gt;%dplyr::select(-Shark)
+                       mutate(BSize = Size - 220), prior = normal(0,1))
 ```
 
 ---
@@ -602,9 +603,25 @@
     }
     deleted = true;
   });
-})();</script>
-
-<script>
+})();
+(function() {
+  "use strict"
+  // Replace <script> tags in slides area to make them executable
+  var scripts = document.querySelectorAll(
+    '.remark-slides-area .remark-slide-container script'
+  );
+  if (!scripts.length) return;
+  for (var i = 0; i < scripts.length; i++) {
+    var s = document.createElement('script');
+    var code = document.createTextNode(scripts[i].textContent);
+    s.appendChild(code);
+    var scriptAttrs = scripts[i].attributes;
+    for (var j = 0; j < scriptAttrs.length; j++) {
+      s.setAttribute(scriptAttrs[j].name, scriptAttrs[j].value);
+    }
+    scripts[i].parentElement.replaceChild(s, scripts[i]);
+  }
+})();
 (function() {
   var links = document.getElementsByTagName('a');
   for (var i = 0; i < links.length; i++) {
@@ -613,7 +630,23 @@
     }
   }
 })();
-</script>
+// adds .remark-code-has-line-highlighted class to <pre> parent elements
+// of code chunks containing highlighted lines with class .remark-code-line-highlighted
+(function(d) {
+  const hlines = d.querySelectorAll('.remark-code-line-highlighted');
+  const preParents = [];
+  const findPreParent = function(line, p = 0) {
+    if (p > 1) return null; // traverse up no further than grandparent
+    const el = line.parentElement;
+    return el.tagName === "PRE" ? el : findPreParent(el, ++p);
+  };
+
+  for (let line of hlines) {
+    let pre = findPreParent(line);
+    if (pre && !preParents.includes(pre)) preParents.push(pre);
+  }
+  preParents.forEach(p => p.classList.add("remark-code-has-line-highlighted"));
+})(document);</script>
 
 <script>
 slideshow._releaseMath = function(el) {