README.Rmd

---
output: github_document
---

<!-- README.md is generated from README.Rmd. Please edit that file -->

```{r setup, include = FALSE}
knitr::opts_chunk$set(
  warning = FALSE,
  message = FALSE,
  collapse = TRUE,
  comment = "#>",
  fig.path = "man/figures/README-",
  fig.retina = 2,
  out.width = "100%"
)
```
# signatselect: Identifying signatures and strengths of selection  

[![Travis build status](https://travis-ci.org/benmarwick/signatselect.svg?branch=master)](https://travis-ci.org/benmarwick/signatselect) [![Binder](http://mybinder.org/badge_logo.svg)](http://beta.mybinder.org/v2/gh/benmarwick/signatselect/master?urlpath=rstudio) 

The goal of signatselect is to provide two functions useful for investigating change over time in artefact assemblages (and genetic time-series data):

- `fit()` the frequency increment test as simple statistical test to aid in the detection and quantification of selective processes in the archaeological record. This is adapted directly from the algorithm published in Feder, A. F., Kryazhimskiy, S., & Plotkin, J. B. (2014). Identifying signatures of selection in genetic time series. _Genetics_, 196(2), 509-522. <https://doi.org/10.1534/genetics.113.158220> and Newberry's implementation of Feder's FIT for R at <https://github.com/mnewberry/ldrift> for their paper Newberry, et al. (2017). Detecting evolutionary forces in language change. _Nature_, 551(7679), 223–226. <https://doi.org/10.1038/nature24455> 

- `tsinfer()` to estimate the population size and the selection coefficient favoring one variant over another from time-series variant-frequency data. This is adapted from Feder et al. (2014) and Hezekiah Akiva Bacovcin's <https://github.com/bacovcin/tsinfer-R>, who adapted it from Sergey Kryazhimskiy's OCaml language version of `tsinfer` at <https://github.com/skryazhi/tsinfer>. 

## Installation

You can install the development version of signatselect from GitHub with:

``` r
# install.packages("pak")
pak::pkg_install("benmarwick/signatselect")
```

## Examples

### The Frequency Increment Test

Here is an example of the `fit()`, the frequency increment test:

```{r fit-example}
# install.packages("tidyverse")
suppressPackageStartupMessages(library(tidyverse))

  # data slightly modified from Feder et al. Table S2 
  feder_table_2 <- 
  tibble(time = c(415 , 505 , 585 , 665 , 745 , 825 , 910),
         freq = c(0.06956522, 0.23125000, 0.62352941, 0.78494624, 0.93333333, 0.97979798, 0.98979592))
```

Let's take a look:

```{r fit-plot, out.width="50%", fig.align='center'}
  ggplot(feder_table_2,
         aes(time, freq)) +
  geom_line() +
  geom_point(size = 5) +
  theme_minimal()
```

There's a trend of increasing frequencies, but is it a result of selection? Let's see:

```{r fit-test}
library(signatselect)

    fit(
      time = feder_table_2$time,
      v = feder_table_2$freq
    )
```

The result of the FIT, with the low p-value, indicates that selection is occuring in this time series.

How about the null situation, can we credibly detect a situation with no selection? Here's a random uniform distribution of a variant: 

```{r h0-fit-example, out.width="50%", fig.align='center'}
  no_selection <- 
  tibble(time = c(415 , 505 , 585 , 665 , 745 , 825 , 910),
         freq = runif(7))

  ggplot(no_selection,
         aes(time, freq)) +
  geom_line() +
  geom_point(size = 5) +
  theme_minimal()
```

```{r no-fit-test}
    fit(
      time = no_selection$time,
      v = no_selection$freq
    )
```

And we see a high p-value, indicating no selection in this time series.

### Infer population size and selection coefficient from time-series variant-frequency data

Here is an example of `tsinfer()` to estimate the population size and the selection coefficient favoring one variant over another from time-series variant-frequency data. Here's some sample data: 

- `tvec` contains sample times   
- `bvec` contains the number of samples of the focal variant (must be integers) and,   
- `nvec` containes total number of samples at each time point (must be integers)   

```{r}
  # data adapted from https://github.com/skryazhi/tsinfer
  tvec = c(0, 10, 20)
  bvec = c(2000, 4000, 6000)
  nvec = c(10000, 10000, 10000)
```

Now we compute the test result:
  
```{r}
  tsinfer_output <- 
  tsinfer(
    tvec = tvec,
    bvec = bvec,
    nvec = nvec,
    verbose = FALSE
  )

# and take a look at the result 
tsinfer_output
```

The selection coefficient for non-neutral model is in `s`, and so the value here is `r round(tsinfer_output$s, 3)`. The population size for non-neutral model is in `alpha`, and here is `r round(tsinfer_output$alpha, 4)`

## An archaeological application

### The frequency increment test

Here is an example of using the FIT to identify pottery types that indicate selection. We are using frequencies of different decorative motifs in the Merzbach assemblage, Neolithic Germany (Crema et al. 2016; many other papers). We can load the data from the `evoarchdata` package on GitHub:

```{r}
# pak::pkg_install("benmarwick/evoarchdata")
library(evoarchdata)
data("ceramics_lbk_merzbach")

# take a look
ceramics_lbk_merzbach
```

Here's a sample of some of the pottery decorations, from [Shennan and Wilkinson (2001)](https://www.jstor.org/stable/2694174):

```{r out.width="50%", fig.align='center', echo = FALSE}
knitr::include_graphics("man/figures/README-decorative-motifs.png")
```

Here's an overview of how each motif changes over time in this assemblage:

```{r}
# get ordered factor of decoration types so we can order the plots nicely 

decoration_types <- 
names(ceramics_lbk_merzbach)[-1] %>%
  enframe() %>% 
  separate(value, into = c('a', 'b'), 2) %>% 
  mutate(b = parse_number(b)) %>% 
  arrange(b) %>% 
  unite(decorations, c(a,b), sep = "") %>% 
  pull(decorations)

# see how the freqs of each change over time
ceramics_lbk_merzbach_long <-
  ceramics_lbk_merzbach %>%
  gather(variable, value, -Phase) %>% 
  mutate(Phase = fct_relevel(Phase, ceramics_lbk_merzbach$Phase)) %>% 
  mutate(variable = fct_relevel(variable, decoration_types))

# plot
ggplot(ceramics_lbk_merzbach_long,
       aes(Phase,
           value)) +
  geom_line(aes(group = 1)) +
  geom_point() +
  facet_wrap(~variable,
             scales = "free_y") +
  theme_minimal(base_size = 8) +
  ggtitle(str_glue('Ceramic decoration frequency data from Merzbach, Germany'))
```

There are many decoration types, so let's narrow it down to ones that have a maximum frequency of at least `r max_n <- 50; max_n`. 

```{r}
ceramics_lbk_merzbach_long_subset <-
  ceramics_lbk_merzbach_long %>% 
  group_by(variable) %>% 
  filter(max(value) > max_n)

# keep these decorations
decorations_to_keep <- unique(as.character(ceramics_lbk_merzbach_long_subset$variable))

# plot
ggplot(ceramics_lbk_merzbach_long_subset,
       aes(Phase,
           value)) +
  geom_line(aes(group = 1)) +
  geom_point() +
  facet_wrap(~variable,
             scales = "free_y") +
  theme_minimal(base_size = 8)
```

We can also view these in the classic archaeological visualiation, the battleship plot:

```{r out.width="50%", fig.align='center'}
phases <-  unique(ceramics_lbk_merzbach_long_subset$Phase)

library(plotrix)

ceramics_lbk_merzbach_matrix <- 
ceramics_lbk_merzbach_long_subset %>% 
  spread(variable, value) %>% 
  select(-Phase) %>% 
  as.matrix() 

row.names(ceramics_lbk_merzbach_matrix) <- phases

battleship.plot(ceramics_lbk_merzbach_matrix,
                yaxlab = phases,
                col = "grey",
                main=str_glue("Merzbach ceramic decoration types (max N > {max_n})"))
```

To prepare the data for the FIT we need to reshape it to a long form, compute frequency as ratio of count of type of a interest to all other types, and drop decoration types with less than threee time points (we need a min of three time points to compute the FIT).

```{r}
# reshape data for each decoration type to go into test:
ceramics_lbk_merzbach_prop <- 
  ceramics_lbk_merzbach %>% 
  select(Phase, decorations_to_keep)

df <- ceramics_lbk_merzbach_prop[ , 2:ncol(ceramics_lbk_merzbach_prop)]
time <- utils:::.roman2numeric(ceramics_lbk_merzbach_prop$Phase)

# compute frequency as ratio of count of type of interest to all other types
list_of_dfs <- vector("list", ncol(df))
names(list_of_dfs) <- names(df)

for(i in 1:ncol(df)){
  tmp <-
    data.frame(time = time,
               count_this_one = df[[i]],
               count_others = rowSums(df[, (seq_len(ncol(df)))[-i]   ]))
  
  # compute frequency
  tmp$frequency = with(tmp, count_this_one / count_others)
  
  # collect results and exclude rows with zero counts for this type i
  list_of_dfs[[i]] <- tmp[which(tmp$count_this_one != 0 ), ]
}

# we need a min of three time points to compute the FIT, so drop decoration types with less than 3
list_of_dfs_three_or_more <- 
  keep(list_of_dfs, ~nrow(.x) >= 3)
```
  
We can prepare safe version of the FIT so we can use it in loops without breaking out of the loop when there is an error for one iteration

```{r}
fit_safely <- 
  safely(fit, 
         otherwise = data.frame(fit_stat = NA,
                                fit_p = NA))
```

Now we can compute the FIT for each pottery decoration type:

```{r, fig.align='center'}
# apply test to each decoration type
df_fit_test_results <-
  list_of_dfs_three_or_more %>%
  bind_rows(.id = "type") %>%
  nest(-type) %>%
  mutate(fit_test = map(data,
                        ~fit_safely(time = .x$time,
                                    v =    .x$frequency))) %>%
  mutate(fit_p = map(fit_test, ~.x$result %>% bind_rows)) %>%
  unnest(fit_p) %>%
  mutate(sig = ifelse(fit_p <= 0.05, "selection", "neutral"))

ceramics_lbk_merzbach_long_sig <-
  ceramics_lbk_merzbach_long_subset %>%
  ungroup %>% 
  left_join(df_fit_test_results %>% 
              select(type, sig), by = c("variable" = "type")) %>%
  mutate(Phase_num = utils:::.roman2numeric(as.character(Phase))) %>% 
  mutate(variable = fct_relevel(variable, decoration_types)) %>% 
  arrange(variable, Phase_num)

# plot to indicate which styles show selection and which do not. 
ggplot(ceramics_lbk_merzbach_long_sig,
       aes(Phase_num,
           value,
           colour = sig,
           shape = sig,
           group = variable)) +
  geom_point(size = 3) +
  scale_color_viridis_d(name = "", 
                        begin = 0.25, 
                        end = 0.75) +
  geom_line() +
  facet_wrap(~variable,
             scales = "free_y") +
  guides(shape = FALSE) +
  theme_minimal(base_size = 8) +
  ggtitle(str_glue('Application of the FIT to decoration frequency data from Merzbach.\nShowing only decoration types that have a maximum frequency of at least {max_n}'))
```

We can subset each time series to see if we can identify episodes of selection among decoration types that might not show overall selection. First we need a function to isolate a data frame of rolling groups of three time points:

```{r}
# Function to get df with rolling groups of n
# from https://stackoverflow.com/a/5543632/1036500
n <- 5
df_with_rolling_idx <- function(df, window = n) {
  nr <- nrow(df)
  w <- window       # window size
  i <- 1:nr         # indices of the rows
  iw <-
    embed(i, w)[, w:1]   # matrix of rolling-window indices of length w
  wnum <- rep(1:nrow(iw), each = w)   # window number
  inds <-
    i[c(t(iw))]         # the indices flattened, to use below
  zw <- sapply(df, '[', inds)
  zw <- transform(data.frame(zw), w = wnum)
  return(zw)
}
```

Now we can compute the FIT on sections of the time series for each decoration type to identify time coordinates where selection has occurred, even when the overall series does not indicate selection:

```{r}
merzbach_long_sig_mid_time_point <- 
list_of_dfs %>% 
  bind_rows(.id = "type") %>% # to get rolling window of n
  nest(-type) %>%
  mutate(rolled = map(data, df_with_rolling_idx)) %>% 
  unnest(rolled) %>% 
  mutate(unid = str_glue('{type}_{w}')) %>% 
  nest(-unid) %>% 
  mutate(fit_test = map(data,
                        ~fit_safely(time = .x$time,
                                    v =    .x$frequency))) %>%
  mutate(fit_p = map(fit_test, ~.x$result %>% bind_rows)) %>%
  unnest(fit_p) %>%
  mutate(sig = ifelse(fit_p <= 0.05, "selection", "neutral")) %>% 
  unnest(data) 

# make type a factor so we can order the plots nicely 
merzbach_long_sig_mid_time_point$type <- 
  fct_relevel(merzbach_long_sig_mid_time_point$type,
              decoration_types[decoration_types %in% merzbach_long_sig_mid_time_point$type])

# plot with overall time-series results also
# harmonize some variable names first

ceramics_lbk_merzbach_long_sig_to_plot_with_others <- 
  ceramics_lbk_merzbach_long_sig %>% 
  rename( time = Phase_num,
          count_this_one = value,
          type = variable) %>% 
  filter(count_this_one != 0) %>% 
  arrange(type, time) %>% 
  mutate(type = fct_relevel(type, decoration_types))

sig_decorations <- 
  ceramics_lbk_merzbach_long_sig_to_plot_with_others %>% 
  filter(sig == "selection") %>% 
  pull(type) %>% 
  as.character() %>% 
  unique()

# here we have the plot showing overall selection, and point-wise selection
ggplot()  +
  geom_line(data = merzbach_long_sig_mid_time_point %>% 
              filter(sig == "selection"),
            aes(time,
                count_this_one,
                group = type),
            size = 5,
            colour = "grey80",
            lineend = "round") +
  geom_point(data = merzbach_long_sig_mid_time_point %>% 
               filter(sig == "selection"),
             aes(time,
                 count_this_one,
                 group = type),
             size = 5,
             colour = "grey80") +
  geom_line(data = ceramics_lbk_merzbach_long_sig_to_plot_with_others,
            aes(time,
                count_this_one, 
                group = type,
                colour = sig)) +
  geom_point(data = ceramics_lbk_merzbach_long_sig_to_plot_with_others,
            aes(time,
                count_this_one, 
                group = type,
                colour = sig,
                shape = sig))  +
  scale_color_viridis_d(name = "", 
                        begin = 0.25, 
                        end = 0.75) + 
  guides(shape = FALSE) +
  facet_wrap( ~ type, scales = "free_y") +
  theme_minimal(base_size = 8) +
  ggtitle(str_glue('Application of the FIT to decoration frequency data from Merzbach.\nShading highlights the data points where FIT identifies selection'))

```

The output is consistent with previous conclusions of the existence of selective forces acting on decorated LBK pottery from these sites [(Kandler and Shennan 2013)](https://doi.org/10.1016/j.jtbi.2013.03.006). We further note that selection is most widespread on the decorations in the later phases of site occupation, an observation also made by [Shennan and Wilkinson (2001)](https://www.jstor.org/stable/2694174) who describe a pro-novelty bias in the later phases, possibly reflecting a concern to establish distinct local identities once the area had filled with people.  

### Estimating the population size and the selection coefficient

We can compute the selection coefficient and related statistics for the decoration types that we previously identified as having a signal of selection. 

```{r}
decoration_types_tsinfer_output <- 
list_of_dfs_three_or_more %>% 
  keep(., names(.) %in% sig_decorations) %>% 
  map(~tsinfer(
    tvec = .x$time,
    bvec = .x$count_this_one,
    nvec = .x$count_others,
    verbose = FALSE
  )) %>% 
  bind_rows(.id = "type")

# inspect the output
decoration_types_tsinfer_output
```

And we can visualise it:

```{r}
# join with frequency data
ceramics_lbk_merzbach_long_sig_to_plot_with_others_select <- 
ceramics_lbk_merzbach_long_sig_to_plot_with_others %>% 
  filter(type %in% sig_decorations) %>% 
  left_join(decoration_types_tsinfer_output) %>% 
  filter(s > 0) %>% 
  mutate(sel_coef = s) %>% 
  mutate(plot_label = str_glue('{type}, s = {round(sel_coef, 3)}'))

# here we have the plot showing overall selection, and point-wise selection
library(scales)
ggplot(ceramics_lbk_merzbach_long_sig_to_plot_with_others_select)  +
  geom_line(aes(time,
                count_this_one, 
                group = plot_label,
                size = (sel_coef),
                alpha = rescale(sel_coef))) +
  geom_point(aes(time,
                count_this_one, 
                group = plot_label,
                size =  rescale(sel_coef),
                alpha = rescale(sel_coef)))  +
  facet_wrap( ~ plot_label, scales = "free_y") +
  scale_alpha_continuous(guide = FALSE) +
  scale_size(NULL, range = c(1, 5)) + 
  theme_minimal(base_size = 8) +
  theme(legend.position="none") + 
  ggtitle(str_glue('Computation of the selection coefficient on decoration frequency data from Merzbach'))
```

Decoration type B20 returns the highest selection coefficient. Looking back at the illustration of the decoration types above, we note that this is one of the more visually complex decoration types. 

Please note that the `signatselect` project is released with a [Contributor Code of Conduct](CODE_OF_CONDUCT.md). By contributing to this project, you agree to abide by its terms.


```{r echo = FALSE, eval=FALSE}

# schematic for FIT

  library(ggplot2)
  library(viridis)
  n <- 30
  mean1 <- 0
  mean2 <- 0.5
  sds <- 0.1
  # create dataset with n brownian distributed points
  # so rescaled frequencies are independent and approximately 
  # normally distributed with mean 0
  df1 <- data.frame(x = 1:n, 
                    y = rnorm(n, mean1, sds))
  # create a ggplot2 scatterplot
  p1 <- ggplot(df1, aes(x, y)) + 
    geom_point(size = 3,
               colour = viridis(1)) +
    geom_line( colour = viridis(1)) +
    ylim(-1,1) +
    geom_hline(yintercept = mean1, 
               colour = viridis(1)) +
    xlab("") +
    ylab(expression(Y[i])) +
    theme_minimal() +
    ggtitle(expression(paste("Rescaled variant frequencies under ", bold("drift")))) 
    
  # add marginal density along the y axis
  m1 <- 
  ggExtra::ggMarginal(p1,
                      type = "density", 
                      margins = "y", 
                      size = 3, 
                      colour = viridis(1),
                      fill = viridis(1))
  
  # create dataset with n brownian distributed points
  df2 <- data.frame(x = 1:n, 
                   y = rnorm(n, mean2, sds))
  # create a ggplot2 scatterplot
  p2 <- ggplot(df2, aes(x, y)) + 
    geom_point(size = 3,
               colour = viridis(5)[4]) +
    geom_line(colour = viridis(5)[4]) +
    geom_hline(yintercept = mean2, 
               colour = viridis(5)[4]) +
    ylim(-1,1) +
    theme_minimal() +
    xlab("Sampled time points") +
    ylab(expression(Y[i])) +
    ggtitle(expression(paste("Rescaled variant frequencies under ", bold("selection")))) + 
    geom_segment(aes(
      x =      0,
      xend =   0,
      y =      0,
      yend =   0.5,
    ),
    size = 1,
    colour = viridis(5)[1],
    arrow = arrow(length = unit(0.25, "cm"))
    ) 
  # add marginal density along the y axis
  m2 <- 
  ggExtra::ggMarginal(p2,
                      type = "density", 
                      margins = "y", 
                      size = 3, 
                      colour = viridis(5)[4],
                      fill = viridis(5)[4])


library(cowplot)
plot_grid(m1 , m2, ncol = 1)
ggsave("man/figures/README-fit-schematic.png",
       w = 7, h = 5)


```