---
title: "Code for Aguirre, Blasi & Sanchez 2020"
author: "Damian E. Blasi"
date: "2/5/2020"
output: pdf_document
---

Load libraries and set options
```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
require(tidyverse)
require(party)
require(ggplot2)
require(GoodmanKruskal)
require(plyr)
require(Rtsne)
require(pals)
theme_set(theme_bw())
```

Loading and visualizing data
---
```{r load,include=T}
pan<-read.csv("/Users/damian/Desktop/Projects/Pan flutes/Final version/SuppFile2.csv",sep=",",
              stringsAsFactors = F,header=T,skip=1)
str(pan)
table(pan$Area_fine)
table(pan$Area_coarse)
```

Normalize NAs
```{r,cache=T}
pan<-as.data.frame(mapply(function(x) 
  ifelse(x %in% c("?","-"),NA,x),pan))
```

We recode number of pipes to coincide with its four quartiles
```{r ordinalization,cache=T}
pan[["No.pipes.S1"]]<-sapply(as.numeric(pan[["No.pipes.S1"]]),function(x)
  ifelse(is.na(x),NA,
          ifelse(x<=12,0,
           ifelse(x<=26,1,
            ifelse(x<=28,2,
             3)))))
```

Check overall missingness
```{r nas,cache=T}
nas<-apply(pan[,c(7:ncol(pan))],2,function(x) sum(is.na(x))/length(x))
fmax<-apply(pan[,c(7:ncol(pan))],2,function(x) max(table(x))/length(x[!is.na(x)]))
```

Plot fraction of missing data per feature
```{r,cache=T}
ggplot(data.frame(missing=nas,
                  id=sapply(names(nas),function(x) x),
                  threshold=sapply(nas,function(x) x>0.25)),
       aes(x=id,y=missing,color=threshold))+
  geom_point()+
  xlab("")+
  ylab("Fraction of missing data")+
  coord_flip()+
  theme(legend.position = "null")
```


Plot frequency of most frequent level per feature
```{r,cache=T}
ggplot(data.frame(freqmax=fmax,
                  id=sapply(names(fmax),function(x) x),
                  threshold=sapply(fmax,function(x) x>0.9)),
       aes(x=id,y=freqmax,color=threshold))+
  geom_point()+
  xlab("")+
  ylab("Frequency of the most frequent value")+
  coord_flip()+
  theme(legend.position = "null")
```

There are some features with a strong tendency of having the same coded value.

Plot nominal associations to assess the colinearity between traits
```{r goodmankruskal,cache=T}
pan_gk<-GKtauDataframe(pan[,c(7:ncol(pan))])
assoc_feat<-(adply(pan_gk,1:2,function(x) ifelse(x>0.5,x,0)) %>% subset((V1>0)&(X1!=X2)))
print(assoc_feat)
```

There are `r nrow(assoc_feat)` dependencies where one variable predicts the other with a precision of over 50%, covering `r length(unique(c(assoc_feat$X1,assoc_feat$X2)))` unique variables. This entails that colinearity will play an important role in our results.

Overall structure in panpipe data
---

We attempt at determining the provenance group where a flute comes from based on its morphological features. We consider how unique trait combinations are
```{r uniquetraits, cache=T}
pan_traits<-colnames(pan[,c(7:ncol(pan))])
unique_pan<-ddply(pan[,c("Ethnological.group",pan_traits)],pan_traits,
                  function(x) data.frame(l=length(unique(x$Ethnological.group)),
                                         e=paste0(unique(x$Ethnological.group),collapse = "|")))
```

There are `r sum(unique_pan$l==1)` instruments that are found only within one provenance. The rest have the following distribvution of provenances:
```{r where_unique,cache=T}
unique_pan$e[unique_pan$l>1]
```

We compare this against 500 randomized versions of the pan flutes
```{r,cache=T}
nsim=500

baseline_unique<-function(DATA) {
  DATA<-as.data.frame(lapply(DATA,function(x) sample(x)))
  DATA$Ethnological.group<-DATA$Ethnological.group
  unique_rand<-ddply(DATA[,c("Ethnological.group",pan_traits)],pan_traits,
                 function(x) data.frame(l=length(unique(x$Ethnological.group)),
                                         e=paste0(unique(x$Ethnological.group),collapse = "|")))
  return(sum(unique_rand$l==1))}

comparison_unique<-rdply(nsim,function(x) baseline_unique(pan[,c("Ethnological.group",pan_traits)]))
```

Most randomizations (`r sum(comparison_unique$V1==401)/nrow(comparison_unique)`) yield 401 distinct instruments. Now we move to compare the distribution of Gower dissimilarities between instruments - both in the actual dataset and the randomized control.

```{r,cache=T}
nsim=500

baseline_comp<-function(DATA) {
  D<-as.data.frame(lapply(DATA,function(x) sample(x)))
  D$Ethnological.group<-DATA$Ethnological.group
  unique_rand<-ddply(D[,c("Ethnological.group",pan_traits)],pan_traits,
                 function(x) data.frame(l=length(unique(x$Ethnological.group)),
                                         e=paste0(unique(x$Ethnological.group),collapse = "|")))
  return(as.vector(cluster::daisy(D))) }

comparison_random<-rdply(nsim,function(x) baseline_comp(pan[,c("Ethnological.group",pan_traits)]))
comparison_random<-data.frame(d=as.vector(as.matrix(comparison_random[,c(2:ncol(comparison_random))])),
                              type="random",stringsAsFactors = F)

distances_empirical<-cluster::daisy(pan[,pan_traits])
comparison_empirical<-data.frame(d=as.vector(distances_empirical),type="empirical",stringsAsFactors = F)

distances<-rbind(comparison_random,comparison_empirical)
```

We plot the distribution of dissimilarities
```{r randomization_plot,cache=F}
ggplot(distances,aes(x=d,y=..density..,group=type,fill=type,color=type))+
  geom_histogram(bins=15)+
  scale_x_continuous(limits=c(0,0.6))+
  facet_wrap(~type,scales = "free",nrow=2,
             labeller = as_labeller(c("empirical"="Panpipes","random"="Randomized control")))+
  labs(x="Gower dissimilarity",y="Density")+
  theme(legend.position = "none")+
    scale_color_manual(values=c("cornflowerblue","coral2"))+
  scale_fill_manual(values=c("cornflowerblue","coral2"))

ggsave("panpipes_distances.tiff",dpi=300)
```

Now we evaluate how do panpipes cluster in the trait space bty appealing to a t-SNE algorithm applied to the Gower dissimilarities between panpipes. Notice that we consider a range of perplexity values (a free parameter in the t-SNE algorithm that relates to the number of effective neighbors expected for each observation.)
```{r tsne,cache=T}
set.seed(26111985)

perplexities<-c(2,5,10,25,50,100)
plots_tsne<-vector("list",length=length(perplexities))

for(ppx in 1:length(perplexities)) {
  tsne_p <- Rtsne(distances_empirical, dims = 2,
              perplexity=perplexities[ppx], verbose=FALSE, max_iter = 500,
              check_duplicates=F)
  if(ppx==1) tsne<-data.frame(Dim1=tsne_p$Y[,1],
                              Dim2=tsne_p$Y[,2],
                              Area=pan$Area_fine,
                              Perplexity=perplexities[ppx])
  else{tsne<-rbind(tsne,data.frame(Dim1=tsne_p$Y[,1],
                              Dim2=tsne_p$Y[,2],
                              Area=pan$Area_fine,
                              Perplexity=perplexities[ppx]))}}

  tsne$Macro<-sapply(as.character(tsne$Area),
                   function(x) ifelse(x %in% c("Africa","Europe","China","SE_Asia"),"Old World",
                                      ifelse(x %in% c("SA_notAndes","Andes"),"South America","Pacific")))
  
   tsne$Area<-sapply(as.character(tsne$Area),
                    function(x) ifelse(x=="NewGuinea","New Guinea",
                                       ifelse(x=="Papuans","West Papua",
                                              ifelse(x=="SE_Asia","SE Asia",
                                                     ifelse(x=="SA_notAndes",
                                                     "SA not Andes",x)))))
  cvxhull <- tsne %>%
  group_by(Macro) %>%
  slice(chull(Dim1,Dim2))
  
  ## Plotting
  ggplot(tsne,aes(x=Dim1,y=Dim2,color=Area,label=Area))+
  scale_color_manual(values=as.vector(kelly(13)[c(3:13)]))+
  geom_point(alpha=0.7,size=2.5)+
    facet_wrap(~Perplexity,scales = "free")

ggsave("tsne_panpipes.tiff",width=7.5,dpi=300)
```

Traceability of provenance given traits
---
Now we focus on potential indicators of cultural transmission. First we put together the Old World (for data size considerations)
```{r,cache=T}
pan$Region<-sapply(as.character(pan$Area_coarse),function(x)
  ifelse(x %in% c("Africa","Asia","Europe"),"Old World",x))
table(pan$Region)
```

We remove trait levels that are present in a single region. We code NA for such levels.
```{r prune_vars,cache=T}
pan_redux<-pan[,c(7:ncol(pan))]
for(i in 1:(ncol(pan_redux)-1)){
tab<-table(pan_redux[,c("Region",colnames(pan_redux)[i])])  
check<-apply(tab,2,function(x) ifelse(sum(x!=0)==1,"fix","ok"))  
pan_redux[,colnames(pan_redux)[i]]<-sapply(pan_redux[,colnames(pan_redux)[i]],
                                           function(x) ifelse(x %in% names(check)[check=="fix"],NA,as.character(x)))
pan_redux[,colnames(pan_redux)[i]]<-as.factor(pan_redux[,colnames(pan_redux)[i]])}
pan_redux<-pan_redux[,sapply(pan_redux, function(x)!all(is.na(x)))]
pan_redux$Region<-as.factor(pan_redux$Region)
```

We focus only on pan pipes produced out of vegetal sources, since those are the only ones present around the globe
```{r material,cache=T}
print(summary(pan_redux$Region))
pan_redux<-subset(pan_redux,Material==2)
print(summary(pan_redux$Region))
pan_redux<-pan_redux[,colnames(pan_redux)!="Material"]
print(summary(pan_redux$Region))
```

We train an efficient weighted random forest (or more precisely, a random conditional inference forest) to predict the region of origin of a pipe given its features
```{r rforest,cache=T}
weights_pan<-sapply(pan_redux$Region,
                     function(x) nrow(pan_redux)/nrow(pan_redux[pan_redux$Region==x,]))
rf_redux<-cforest(Region~.,
                  controls=
                    cforest_control(ntree = 500,
                                    mtry=5,
                                    mincriterion=qnorm(0.9),
                                    fraction = 0.632,
                                    testtype ="Teststatistic",
                                    teststat="max",
                                    replace=TRUE,
                                    trace=F,
                                    savesplitstats=F,
                                    minsplit=20,
                                    minbucket=8),
                  data=pan_redux,
                  weights = weights_pan)
```

Produce the confusion matrix
```{r confusion_redux,cache=T}
cmat_redux<-table(pan_redux$Region, predict(rf_redux))
print(cmat_redux)
```

Statistics per region
```{r,cache=T}
stats_region<-ldply(c(1:nrow(cmat_redux)),
                  function(x) data.frame(id=rownames(cmat_redux)[x],
                                         N=sum(cmat_redux[,x]),
                                         d=cmat_redux[x,x],
                                         A=sum(cmat_redux[x,]))) %>% transform(p=d/N,
                                                                          r=d/A)

stats_region
```

Obtain all the panpipes that are misclassified
```{r misclas,cache=T}
pan_redux$predReg<-predict(rf_redux)
missclass<-cbind(pan[rownames(pan_redux[pan_redux$Region!=pan_redux$predReg,]),],
                 data.frame(predReg=pan_redux$predReg[pan_redux$Region!=pan_redux$predReg]))
missclass
```
