---
title: "Supplementary Material 3: Random Forests"
subtitle: "Measuring and assessing indeterminacy and variation in the morphology-syntax distinction"
author: 
  - Sandra Auderset
  - Adam J.R. Tallman
# date: "6/27/2021"
output: 
  pdf_document: 
    fig_caption: yes
    latex_engine: xelatex
    keep_tex: yes
---

```{r setup, include=FALSE}
library(knitr)
opts_chunk$set(tidy.opts=list(width.cutoff=60),tidy=TRUE, echo=TRUE)

hook_output = knit_hooks$get('output')
knit_hooks$set(output = function(x, options) {
  # this hook is used only when the linewidth option is not NULL
  if (!is.null(n <- options$linewidth)) {
    x = knitr:::split_lines(x)
    # any lines wider than n should be wrapped
    if (any(nchar(x) > n)) x = strwrap(x, width = n)
    x = paste(x, collapse = '\n')
  }
  hook_output(x, options)
})
```

In this file we provide the code as well as some additional explanation and plots of the Random Forest models (cf. Section 6 of the paper).


## Data Preparation

First, we load the necessary packages and set the seed to make the analysis fully reproducible (see the Rmd document for details).
```{r preamble, include=FALSE}
# load packages
library(groundhog)
packages <- c('here', 'tidyverse', 'randomForest', 'gdata', 'ggpubr')
groundhog.library(packages, "2022-02-04")
```

```{r set seed}
# set seed
set.seed(200)
```

Next we load data (available at DOI), excluding columns we do not need for the analysis, i.e. we exclude the variables that are summarized in Exponence Complexity and certain meta information about the languages such as affiliation and glottocode.
```{r data}
db <- read_csv(here("SM1_Database.csv")) %>% select(Language, ClassAuthorSimple:EXPcomplex)
glimpse(db)
```

Since there are too few data points for `incorporand', we exclude this category from the analysis. We convert all variables to factors and subset the columns relevant for the analysis. Finally, we recode exponence complexity as a binary variable with values 'low' (1) and 'high' (over 1).

```{r prep}
db_rf <- db %>% filter(ClassAuthorSimple!="incorporand") %>% 
  droplevels() %>% 
  mutate(EXPcomplex = if_else(EXPcomplex==1, "low", "high")) %>%
  mutate_if(is.character, as.factor)
glimpse(db_rf)
```

We split the data set by language.
```{r langsets}
db_names <- db_rf %>% group_by(Language) %>% group_keys() %>% pull(1)
db_lang <- db_rf %>% group_split(Language) %>% 
  setNames(db_names) %>% 
  map(as_tibble) %>% 
  map(., ~ (.x %>% select(-Language))) %>%
  drop.levels()
glimpse(db_lang)
```

We create a custom function for calculating baseline, accuracy and difference.

```{r rffunc}
rf_acc <- function(x){
  cm <- as.matrix(x$confusion[ , -3])
  bl <- round(max(rowSums(cm))/sum(cm), 4)
  acc <- round(sum(diag(cm))/sum(cm), 4)
  diff <- round(acc-bl, 4)
  returnlist <- c("baseline" = bl, "accuracy" = acc, "difference" = diff)
}
``` 


# Author Classification

We apply a Random Forest model for each language separately with author classification as the dependent variable and all the wordhood criterial variables (exponence complexity, prominence projection, coding elaboration, fixedness, and interruptability) as predictors. More formally, the model is specified as:

```{=latex}
\begin{quote}
	Author Classification ~ Exponence Complexity + Fixedness + Boundedness + Interruptability + Prominence + Coding Elaboration
\end{quote}
```


This one and all the following RF models are implemented aggregating over 10'000 trees with two variables tried at each split. Each model outputs error rates for each level of the variable in a confusion matrix and an overall error rate for the model (called out-of-bag error or estimate, OOB for short). The former provides the classification accuracy for words vs. clitics vs. affixes and the latter how accurate the classification is overall.
The model also provides a plot displaying Mean Decrease in Accuracy, which can be related to the relative importance of variables. Broadly speaking, certain variables are better and more consistent at classifying the dependent variable correctly than others across all the decision trees and this is reflected by MDA. 
Note that the actual value does not matter, since it is not comparable across models. What is important is the ranking of the independent variables and whether or not there are clear breaks between them. 

To this we add a calculation of the baseline, the accuracy, and the difference between these two. The idea is to account for the skewness of the data, since the more skewed the data are, the easier it is to get a correct classification by chance. For example, imagine two languages X and Y with 100 data points each. Language X has 34 affixes, 33 clitics, and 33 words, while language Y has 15 affixes, 5 clitics, and 85 words. If we classify all elements as words in language Y, we get 85\% right. But no matter what category we choose for language X, we do not achieve over 34\% correct classification. 
The relative skewness by language needs to be taken into account when assessing whether the RF model performed well or not. The calculation of each of these measures is very simple. The baseline is exactly what we just discussed. More simply put, it represents the maximal row-sum of the confusion matrix divided by the number of data points. The accuracy represents the sum of correct predictions (i.e. the diagonal in the confusion matrix) divided by the number of data points, i.e. it reflects the proportion of correct predictions. 
The difference is calculated by subtracting the baseline measure from the accuracy, i.e. it reflects how much better the RF model performed over chance.

```{r rf_auth, linewidth=90}
rf_aut <- lapply(1:length(db_lang), function(a) randomForest(ClassAuthorSimple~., data=db_lang[[a]], ntree=10000, mtry=2, importance=TRUE))
names(rf_aut) <- names(db_lang)
rf_aut
```

The OOB error rates range from 0 to over 23\%, basically reflecting skewness as discussed above. CAY and Movima, which both have OOB errors of 0\% also both have very skewed distributions while Chácobo, which has the most even distribution across the three categories also has a quite high error. 
It thus more instructive to look at the difference between baseline and accuracy. The highest performing models -- in Chácobo and Ashéninka Perené -- are performing over 30\% better than chance. In Wãnsöjöt and Cavineña, however, the model is performing close to chance level and in Hup it is underperforming the baseline. 
Overall, we see that even though many models achieve an impressionistically high accuracy, the performance is not better than the baseline (i.e. chance). 


```{r rf_aut_measures}
rf_aut_measures <- lapply(1:length(rf_aut), function(m) rf_acc(rf_aut[[m]]))
names(rf_aut_measures) <- names(rf_aut)
rf_aut_measures
```


```{r rf_aut_plotsprep, results='hide', message=FALSE, fig.show='hide', linewidth=60}
# create data frames for variable importance plots
rf_aut_imp <- lapply(1:length(rf_aut), function(i) 
  as_tibble(varImpPlot(rf_aut[[i]]), rownames = "Variable") %>% slice(-1) %>% mutate(MeanDecreaseAccuracy = round(MeanDecreaseAccuracy)))
```

```{r rf_aut_plots, fig.width=3, fig.height=3 , message=FALSE, linewidth=60}
# display plots
rf_aut_plots <- lapply(1:length(rf_aut_imp), function(p)
  ggplot(rf_aut_imp[[p]], aes(x = MeanDecreaseAccuracy, y = reorder(Variable, MeanDecreaseAccuracy))) +
    geom_point(size = 3) + 
    theme_bw() + 
    theme_update(panel.grid.minor.x = element_blank(), 
     panel.grid.major.y = element_line(color = "black", linetype = "dotted"), axis.title.y = element_blank(), axis.text.y = element_text(face = "bold")) + 
    ggtitle(names(rf_aut)[[p]]))
names(rf_aut_plots) <- names(rf_aut)
rf_aut_plots
```

The MDA plots show that authors generally rely on one to three variables to classify elements into affixes, clitics, and words. This can be seen by the break that is present in each plot between the upper one to three variables to the lower ones. The plots also show that different authors rely on different variables, although prominence projection appears to be the most important. 
Both Wãnsöjöt and Hup show the first break already after one variable. This means the authors rely predominantly on prosodic prominence and coding elaboration, respectively. Cavineña shows a break after three variables, fixedness, coding elaboration, and prominence. However, it appears that relying on either just one or all of three variables results in the model performing at about chance level.
For the languages studied here, it seems that relying on two variables results in the most consistent classifications by authors. In all of the remaining languages, prominence projection is one of these variables, but the others vary. 
To summarize, the results show that authors rely on different variables and different numbers of variables, but prominence often plays a role when classifying elements into words, affixes, and clitics. 


# Exponence Complexity

We now apply a Random Forest model for each language separately with exponence complexity as the dependent variable and the other wordhood criterial variables (prominence projection, coding elaboration, fixedness, and interruptability) as predictors. More formally, the model is specified as:

```{=latex}
\begin{quote}
	Exponence Complexity ~ Fixedness + Boundedness + Interruptability + Prominence + Coding Elaboration
\end{quote}
```

We split the data by language, this time excluding author classification.
```{r langsets2}
db_names <- db_rf %>% group_by(Language) %>% group_keys() %>% pull(1)
db_lang2 <- db_rf %>% group_split(Language) %>% 
  setNames(db_names) %>% 
  map(as_tibble) %>% 
  map(., ~ (.x %>% select(-Language, -ClassAuthorSimple))) %>%
  drop.levels()
glimpse(db_lang2)
```

We apply the random forest models to this new data frame as specified above. 
```{r rf_ex, linewidth=90}
rf_ex <- lapply(1:length(db_lang2), function(a) randomForest(EXPcomplex~., data=db_lang2[[a]], ntree=10000, mtry=2, importance=TRUE))
names(rf_ex) <- names(db_lang2)
rf_ex
```

We calculate baseline, accuraccy and difference.
```{r rf_ex_measures}
rf_ex_measures <- lapply(1:length(rf_ex), function(m) rf_acc(rf_ex[[m]]))
names(rf_ex_measures) <- names(rf_ex)
rf_ex_measures
```


```{r rf_ex_plotsprep, results='hide', message=FALSE, fig.show='hide'}
# create data frames for variable importance plots
rf_ex_imp <- lapply(1:length(rf_ex), function(i) 
  as_tibble(varImpPlot(rf_ex[[i]]), rownames = "Variable") %>% mutate(MeanDecreaseAccuracy = round(MeanDecreaseAccuracy)))
```

And we display the variable importance plots.

```{r rf_ex_plots, linewidth=40, fig.width=3, fig.height=3 , message=FALSE}
# display plots
rf_ex_plots <- lapply(1:length(rf_ex_imp), function(p)
  ggplot(rf_ex_imp[[p]], aes(x = MeanDecreaseAccuracy, y = reorder(Variable, MeanDecreaseAccuracy))) +
    geom_point(size = 3) + 
    theme_bw() + 
    theme_update(panel.grid.minor.x = element_blank(), 
      panel.grid.major.y = element_line(color = "black", linetype = "dotted"),       axis.title.y = element_blank(), 
      axis.text.y = element_text(face = "bold")) + 
    ggtitle(names(rf_ex)[[p]]))
names(rf_ex_plots) <- names(rf_ex)
rf_ex_plots
```

```{r export Movima and Chacobo plot, include=FALSE}
rf_paper_plots <- ggarrange(rf_ex_plots$Movima, rf_ex_plots$Chacobo, nrow = 1, ncol = 2)
rf_paper_plots
ggsave("rfexcompplots.png", plot = rf_paper_plots, device = "png", width = 13, height = 7, units = "cm", dpi = 600)
```



