---
title: 'Supplementary Material 4: Cluster Validation'
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=80),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)
})
```


```{r preamble, include=FALSE}
# load packages
library(groundhog)
packages <- c('here', 'tidyverse', 'cluster', 'clValid', 'ggdendro', 'gdata', 'ggpubr')
groundhog.library(packages, "2022-02-04")
```

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


# Data Preparation

We first produce the simulated data set. We do this by developing four random list of 0s and 1s, a random list of 0s, 1s and 2s and a random list of integers from 0 to 10. The  0s and 1s for the binary variables were replaced 'y(es)'s and 'n(o)'s to correspond to interruptability/contiguity, fixedness/variability, simplex/complex and freedom/boundedness variables of our data set. The 0s, 1s and 2s were replaced with 'y(es)', 'both' and 'n(o)' respectively corresponding to prominence projection / projection deficiency. In our dataset exponence complexity varies from 1 to 10, and we thus create a variable that has a uniform distribution of numbers from 1 to 10 corresponding to exponence complexity. 

```{r simulate}
# simulate random data set
FREE <- as.character(sample(c("y", "n"), 100, replace=T))
INTERone <- as.character(sample(c("y", "n"), 100, replace=T))
CODelab <- as.character(sample(c("y", "n"), 100, replace=T))
PMfixed <- as.character(sample(c("y", "n"), 100, replace=T))
PRM <- as.character(sample(c("y", "both", "n"), 100, replace=T))
EXPcomplex <- sample(1:10, 100, replace=T, prob = 10:1)

sim_db <- as_tibble(cbind(FREE, INTERone, CODelab,PMfixed, PRM, EXPcomplex)) %>% mutate_if(is.character, as.factor)
glimpse(sim_db)
summary(sim_db)
```

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 load}
# load data
db <- read_csv(here("SM1_Database.csv")) %>% select(Language, PMfixed:EXPcomplex)
glimpse(db)
```


We split the data set by language and select the relevant columns. Then we convert all the character columns to factors to be able to compute the dissimilarity matrices.
```{r langsets}
db_names <- db %>% group_by(Language) %>% group_keys() %>% pull(1)
db_lang <- db %>% group_split(Language) %>% 
  setNames(db_names) %>% 
  map(as_tibble) %>% 
  map(., ~ (.x %>% select(-Language))) %>%
  map(., ~ (.x %>% mutate_if(is.character, as.factor))) %>%
  drop.levels()
glimpse(db_lang)
```


# Creating the (dis)similarity matrices

We use the Gower distance appropriate for mixed data types to create the dissimiliarity matrices that will serve as the input for the clustering.

```{r dissmatrix, message=FALSE, warning=FALSE}
db_dist <- lapply(1:length(db_lang), function(d) 
  daisy(db_lang[[d]] ,metric = "gower"))

# simulated data set
sim_dist <- daisy(sim_db, metric = "gower")
```


# Hierarchical clustering per language

We apply hierarchical clustering to the (dis)similarity matrices and plot the results as a dendrogram.
We do this furst with the simulated data set and then with the data of the sample.

```{r hclustsim}
# hierarchical clustering
sim_hclust <- hclust(sim_dist, method="ward.D2")
# plot results
sim_dendro <- ggdendrogram(sim_hclust, leaf_labels = FALSE, labels = FALSE)
sim_dendro
```

```{r simplot, include=FALSE}
ggsave("hclustrandom.png", sim_dendro, device = "png", width = 7, height = 6, units = "cm", dpi = 600)
```

```{r hclust, message=FALSE}
# hierarchical clustering
db_hclust <- lapply(db_dist, function(c)
  hclust(c, method = "ward.D2"))
names(db_hclust) <- names(db_lang)
# plot results
db_dendro <- lapply(names(db_hclust), function(p)
  ggdendrogram(db_hclust[[p]], leaf_labels = FALSE, labels = FALSE) + ylim(0, 5) + labs(title = p))
db_dendro
``` 

```{r export clusters, include=FALSE}
# arrange plots
dend_grid <- ggarrange(db_dendro[[1]], db_dendro[[2]], db_dendro[[3]], db_dendro[[4]], db_dendro[[5]], db_dendro[[6]], db_dendro[[7]], db_dendro[[8]],   ncol = 3, nrow = 3)
# export
ggsave("hclustplots.png", plot = dend_grid, device = "png", width = 20, height = 20, units = "cm", dpi = 600)
```


We also calculate the cophenetic coefficient and distance correlation, as well as the height difference between the first and second partition of the cluster. Do do this efficiently, we write a custom function.

```{r clustcalc}
coph_height <- function(x, y){
  coph <- cophenetic(x)
  cor <- round(cor(y, coph), 4)
  height <- x$height
  height_diff <- round(height[length(height)]-height[length(height)-1], 4)
  returnlist <- c("heightdiff" = height_diff, "coph" = cor)
}
```

We apply this to the simulated data set and to the language data. 

```{r coph}
# simulated data
sim_ch <- coph_height(sim_hclust, sim_dist)
sim_ch
# language data
db_ch <- lapply(1:length(db_hclust), function(d) 
  coph_height(db_hclust[[d]], db_dist[[d]]))
names(db_ch) <- names(db_lang)
db_ch
```


# Comparison with simulated languages

An anonymous reviewer suggests that the languages of our sample should better be compared to many simulated languages rather than one.
Thus we simulate 1000 languages, apply hierarchical clustering and calculate cophenetic distance and height difference of the first partition

```{r simulate 1000 languages, warning=FALSE}
# set up variables and loop over them
sim1000 <- lapply(1:1000, function(s) {
  FREE <- as.character(sample(c("y", "n"), 100, replace=T))
  INTERone <- as.character(sample(c("y", "n"), 100, replace=T))
  CODelab <- as.character(sample(c("y", "n"), 100, replace=T))
  PMfixed <- as.character(sample(c("y", "n"), 100, replace=T))
  PRM <- as.character(sample(c("y", "both", "n"), 100, replace=T))
  EXPcomplex <- sample(1:10, 100, replace=T, prob = 10:1)
  # combine to df
  sim_db <- as_tibble(cbind(FREE, INTERone, CODelab,PMfixed, PRM, EXPcomplex)) %>% mutate_if(is.character, as.factor)
  # compute distances
  sim_dist <- daisy(sim_db, metric = "gower")
  # apply hierarchical clustering
  sim_hclust <- hclust(sim_dist, method="ward.D2")
  # compute height and cophenetic distances
  sim_ch <- coph_height(sim_hclust, sim_dist)}
  )
# convert to df
sim1000_df <- as.data.frame(do.call(rbind, sim1000))
# means of simulated languages
summary(sim1000_df)
```


We then plot the cophenetic distances and height differences of the simulated data and the languages of our sample against each other.

```{r plot simulated vs. real, linewidth=60}
# convert language list into df and assign rownames as column
db_ch_df <- as.data.frame(do.call(rbind, db_ch))
db_ch_df <- rownames_to_column(db_ch_df, "Language")

# plot cophenetic distance and height difference
sim_comp <- ggplot(sim1000_df, aes(x = coph, y = heightdiff)) + geom_point(size = 2, shape = 1) +
  geom_point(data = db_ch_df, aes(x = coph, y = heightdiff), color = "red", size = 3) +
  geom_text(data = db_ch_df, aes(x = coph, y = heightdiff, label = Language), nudge_x = -0.02, nudge_y = -0.1) +
  xlim(0.4, 1.05) +
  labs(x = "Cophenetic Correlation", y = "Height Difference (morphology-syntax division)")
sim_comp
```

```{r export plot, include=FALSE}
# export
ggsave("simcomparison.png", plot = sim_comp, device = "png", width = 20, height = 15, units = "cm", dpi = 600)
```



# Reproducing Haspelmath's cluster graphic

Below we include the code we used to reproduce the graphics from Haspelmath 2011:64.

```{r Haspelmath graphic, warning=FALSE, linewidth=60}
# reproducing Haspelmath's graphic
# set up vectors
noclust <- runif(200, min = 0, max = 1)
clust1 <- rnorm(100, mean = 0.09, sd = 0.15)
clust2 <- rnorm(100, mean = 0.91, sd = 0.15)
# combine into df
noclust_db <- as_tibble(noclust)
clust_db <- as_tibble(clust1, clust2)
# make scatterplot
noclust_plot <- ggplot(noclust_db) + 
  geom_point(aes(x = noclust, y = 1), shape = 1, size = 3) + theme_linedraw() + theme(axis.title=element_blank(), axis.text.x=element_text(size = 12), axis.text.y=element_blank(), axis.ticks.y=element_blank(), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), plot.title = element_text(size = 16)) + 
  xlim(0, 1) +
  ylim(0.7, 1.3) +
  ggtitle("Random distribution")

clust_plot <- ggplot(clust_db) + 
  geom_point(aes(x = clust1, y = 1), shape = 1, size = 3, inherit.aes = FALSE) +
  geom_point(aes(x = clust2, y = 1), shape = 1, size = 3, inherit.aes = FALSE) +
  theme_linedraw() + 
  theme(axis.title=element_blank(), axis.text.x=element_text(size = 12), axis.text.y=element_blank(), axis.ticks.y=element_blank(), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), plot.title = element_text(size = 16)) + 
  xlim(0, 1) +
  ylim(0.7, 1.3) +
  ggtitle("Two 'clusters' distribution")

hm_plots <- ggarrange(noclust_plot, clust_plot, ncol = 1, nrow = 2)
hm_plots
```

```{r export cluster graphics, include=FALSE}
ggsave("clusterdistributions.png", hm_plots, dpi = 600)
```
