index.qmd

---
title: "Secret weapon workshop"
editor_options: 
  chunk_output_type: console
---

## What is the secret weapon? 

In many ways, ecology is the discipline of replication. 
Most ecological studies have multiple species, multiple sites, or multiple time points in their study -- almost by definition. 
As a result, our analyses are almost always designed to include replication. The most common tool that people use is the multi-level model (AKA the hierarchical model). 
These models are popular, useful, and very flexible; but they have the weakness of being very complex, difficult to interpret, and often difficult to write correctly.

Aside from the technical complexity, there's also conceptual complexity when ecological data has many dimensions of space, time, species, etc. It can be quite hard to visualize and understand what's happening. Also, we may have hypotheses that operate at multiple levels. We may have one hypothesis for what happens within any group, for example, across individuals or species, and another hypothesis for what happens across sites or over time.

Enter the secret weapon. The secret weapon is the name [Andrew Gelman ](https://statmodeling.stat.columbia.edu/2005/03/07/the_secret_weap/)gives to this simple idea. Instead of conducting one large analysis, divide your data into multiple small groups, fit a model to each small group, and compare one or more summary statistics across groups.

Then we can either analyze the summary statistics in a second study, which is a technique sometimes called "statistics on statistics", or go on to build a multilevel model. We may even decide that this is a sufficient endpoint for our study.

This is a technique that everyone in this room can use right now and which I often recommend to anybody who's beginning a quantitative analysis project.

This workshop will be divided into 3 parts. First, we'll do a worked example using the `mite` data from the package `vegan`. In the second, we're going to split into an exercise section where we will look at different examples from some long-term datasets that I've found. And third, we will have time to work on our own datasets. 


## Worked example

### Data preparation -- Center and scale your predictors!

This example uses the famous `mite` dataset, collected by Daniel Borcard in 1989 and included in the famous `vegan` package.

[Slides introducing the dataset](https://bios2.github.io/hiermod/slides/01_Data/#/oribatid-mite)

First, we'll load the data and do some reorganization and preparation.

```{r}
data(mite, package = "vegan")
data("mite.env", package = "vegan")
library(tidyverse)
library(cmdstanr)

# combine data and environment
mite_data_long <- bind_cols(mite.env, mite) |> 
  tibble::rownames_to_column(var = "site_id") |> 
  pivot_longer(Brachy:Trimalc2, names_to = "spp", values_to = "abd")
```

Because we're ecologists, we'll do one of the most classic ecological models: a Binomial GLM, modelling the presence or absence of a species based on the water content of the soil.

$$
\begin{align}
\text{Presence} &\sim \text{Binomial}\left(\frac{e^a}{1 + e^a}, 1\right) \\
a &= \beta_0 + \beta_xx
\end{align}
$$

Where $x$ is some predictor variable and $\frac{e^a}{1 + e^a}$ is the logit link function.
To keep things simple and univariate, let's consider only water for $x$. 

First, a quick word about centering and scaling a predictor variable:

1.  I center the predictor by subtracting the mean. This changes the *intercept* of my linear predictor. It becomes the mean log-odds of occurrance when the water content is average
2.  I divide water content by 100. The dataset has units of **grams per Litre** of water (see `?vegan::mite.env` for more details). This is fine, but I don't think mites are able to sense differences as precise as a millimeter of water either way. by dividing by 10 I transform this into centilitres, which is more informative.

```{r}
mite_data_long_transformed <- mite_data_long |> 
  mutate(presabs = as.numeric(abd>0),
         # center predictors
         water = (WatrCont - mean(WatrCont)) / 100
         )
```


This makes the model:

$$
\begin{align}
\text{Presence} &\sim \text{Binomial}\left(\frac{e^a}{1 + e^a}, 1\right) \\
a &= \beta_0 + \beta_{\text{water}}\frac{\text{water} - \overline{\text{water}}}{100}
\end{align}
$$


### The secret weapon: visually

We can already get a very good visual expression of the Secret Weapon approach using the power of ggplot2 and `stat_smooth()`

```{r}
mite_data_long_transformed |> 
  ggplot(aes(x = water, y = presabs)) + 
  geom_point() + 
  stat_smooth(method = "glm", method.args = list(family = "binomial")) + 
  facet_wrap(~spp)
```

some things to notice about this figure:

-   the x-axis scale has been transformed from "grams per litre" to "centilitres away from average"
-   there is a ton of variation in how different species respond to water!

### The secret weapon: in tidyverse code

We're going to perform the following steps on these data:

1. divide the dataset up into different species
1. fit _precisely the same_ model to each one
1. extract the coefficients
1. visualize these

To do this, we'll be using some handy tricks from the tidyverse, which let us implement a strategy called "split-apply-combine"

#### Split-Apply-Combine

["split-apply-combine"](https://vita.had.co.nz/papers/plyr.pdf) describes a really general workflow in statistics / data science, which was described my Rstudio's Hadley Wickham in the paper linked above. The idea is to divide data into groups, transform each, and then put the results together. The Secret Weapon is one example of this approach.

There are many possible ways to do this in practice. We are using a technique here from the tidyverse, which you can [read more about here](https://dplyr.tidyverse.org/articles/rowwise.html).  A for-loop would also work for this kind of problem!


```{r}
mite_many_glms <- mite_data_long_transformed |> 
  nest_by(spp) |> 
  mutate(logistic_regressions = list(
    glm(presabs ~ water,
        family = "binomial",
        data = data))) |> 
  mutate(coefs = list(broom::tidy(logistic_regressions)))
```

```{r}
mite_many_glm_coefs <- mite_many_glms |> 
  select(-data, -logistic_regressions) |> 
  unnest(coefs)

mite_many_glm_coefs |> 
  ggplot(aes(x = estimate, y = spp,
             xmin = estimate - std.error,
             xmax = estimate + std.error)) + 
  geom_pointrange() + 
  facet_wrap(~term, scales = "free")
```

As you can see, some of these estimates are high, others low. We could also plot these as histograms to see this distribution.

```{r}
mite_many_glm_coefs |> 
  ggplot(aes(x = estimate)) + 
  geom_histogram(binwidth = .5) + 
  facet_wrap(~term, scales = "free")

```

Once again, the two parameters of this model represent:

-   *Intercept* The probability (in log-odds) of a species being present at the average water concentration. some species are common, others are rare.
-   *water* this is the change in probability (in log-odds) as water increases by one centilitre per litre of substrate.



#### For the Bayesians

The above is adapted from a tutorial I wrote for a course in Bayesian statistics. If you would like to see how this is done in the Bayesian software Stan, [the course notes are here](https://bios2.github.io/hiermod/topics/correlated_effects/)


## Practice examples

Here are some exercises that might make good practice for using the Secret Weapon:

* The `mite.env` data has another continuous predictor, `SubsDens`. Adapt the example above for this variable. Which has stronger effects on species presence-absence: water, or substrate density?
* We used the tidyverse to implement the Secret Weapon. This would also work very well in a for-loop! Write a for-loop version of the same workflow.
* We did a binomial GLM for presence absence. We could also use a linear model to estimate evenness in each plot, based on a rank-abundance plot, like this:

```{r}
rank_df <- mite_data_long_transformed |> 
  filter(abd != 0) |> 
  group_by(site_id) |> 
  mutate(spp_rank = dense_rank(desc(abd)))

rank_df |> 
  ungroup() |> 
  ## drop some sites just for appearance
  filter(site_id %in% c("1", "15", "47")) |> 
  ggplot(aes(x = spp_rank, y = log(abd))) + 
  facet_wrap(~site_id) + 
  geom_point() + 
  stat_smooth(method = "lm")

```

The slope of these lines is one way to think about evenness -- more negative slopes mean a less even community.  Use this approach to ask the question: Does evenness change with water concentration?

## Take it to a brand-new dataset (or to your own science!)

Install the LTER data sampler package

```{r eval=FALSE}
# install.packages("remotes")
remotes::install_github("lter/lterdatasampler")
```

```{r}
library(lterdatasampler)
```

Check out the [documentation here](https://lter.github.io/lterdatasampler/)

There are LOTS of useful datasets in here, which we can use to practice the Secret Weapon. Here are some suggestions to get us started:

#### Seedling stem length vs Stem dry mass

```{r}
knitr::kable(head(hbr_maples))
```

#### Growth rates of bison over time

Alternative idea: average differences between males and females across years.

```{r}
knitr::kable(head(knz_bison))
```