### Set up environment

```{r preamble}
library(dplyr)
library(ape)
library(geiger)
library(circlize)
library(corHMM)

# XXX set working directory to script source first
REPO_HOME <- paste(getwd(), "/../", sep = "")

# generic function to load single SCCS tabls
read.sccs.variable <- function(sccs.code) {
  sccs.file <- sprintf('%s/data/d-place/%s.csv',REPO_HOME,sccs.code)
  sccs.matrix <- select(
    read.table( sccs.file, header = T, sep = ',', quote = "\"" ),
    language_glottocode,
    code
  )
  sccs.df <- as.data.frame(sccs.matrix)
  row.names(sccs.df) <- sccs.df$language_glottocode
  names(sccs.df) <- c('language_glottocode',sccs.code)
  return(sccs.df)
}
```

### Domestic animals (SCCS)

The Standard Cross-Cultural Sample (SCCS) is a database that records many
variables across societies. Thanks to [D-PLACE](https://d-place.org), these
data are easily available and integrated with phylogenies of cultural evolution.

| *Species*                                                         | *Cluster* | *SCCS code*                                                            |
|-------------------------------------------------------------------|-----------|------------------------------------------------------------------------|
| _Capra hircus aegagrus_                                           | 7         | [1486. GOATS PRESENT](https://d-place.org/parameters/SCCS1486)         |
| _Ovis aries orientalis_                                           | 5         | [1487. SHEEP PRESENT](https://d-place.org/parameters/SCCS1487)         |
| _Sus scrofa_                                                      | 14        | [1488. PIGS PRESENT](https://d-place.org/parameters/SCCS1488)          |
| _Lama glama guanicoe_, _Vicagna vicugna_                          | 5, 5      | [1521. LLAMA PRESENT](https://d-place.org/parameters/SCCS1521)         |
| _Bubalus bubalis arnee_                                           | 2         | [1520. BUFFALO PRESENT](https://d-place.org/parameters/SCCS1520)       |
| _Bos grunniens mutus_                                             | 7         | [1519. YAKS PRESENT](https://d-place.org/parameters/SCCS1519)          |
| _Rangifer tarandus_                                               | 26        | [1518. REINDEER PRESENT](https://d-place.org/parameters/SCCS1518)      |
| _Camelus bactrianus_, _Camelus dromedarius_                       | 7, 7      | [1517. CAMELS PRESENT](https://d-place.org/parameters/SCCS1517)        |
| _Equus africanus_                                                 | NA        | [1516. DONKEYS/MULES PRESENT](https://d-place.org/parameters/SCCS1516) |
| _Equus przewalskii_                                               | 7         | [1515. HORSES PRESENT](https://d-place.org/parameters/SCCS1515)        |
| _Bos frontalis gaurus_, _Bos javanicus_, _Bos taurus primigenius_ | 2, 5, 5   | [1514. CATTLE PRESENT](https://d-place.org/parameters/SCCS1514)        |

```{r domestic_interactions}
# mapping from data set to taxa
sccs.code <- list(
  SCCS1486 = c('Capra_hircus_aegagrus'),
  SCCS1487 = c('Ovis_aries_orientalis'),
  SCCS1488 = c('Sus_scrofa'),
  SCCS1521 = c('Lama_glama_guanicoe', 'Vicagna_vicugna'),
  SCCS1520 = c('Bubalus_bubalis_arnee'),
  SCCS1519 = c('Bos_grunniens_mutus'),
  SCCS1518 = c('Rangifer_tarandus'),
  SCCS1517 = c('Camelus_bactrianus', 'Camelus_dromedarius'),
  SCCS1516 = c('Equus_africanus'),
  SCCS1515 = c('Equus_przewalskii'),
  SCCS1514 = c('Bos_frontalis_gaurus', 'Bos_javanicus', 'Bos_taurus_primigenius')
)

animal.present.df <- data.frame( language_glottocode = character())
for ( sc in names(sccs.code) ) {
  temp.df <- read.sccs.variable(sc)
  
  # copy values over to all taxa
  values <- temp.df[,2]
  for ( taxon in sccs.code[[sc]] ) {
    temp.df[[taxon]] <- values
  }
  rm(values, taxon) # clean
  
  # subset to remove opaque identifier column header
  temp.df <- temp.df[,c('language_glottocode',sccs.code[[sc]])]
  
  # do an outer join on the glottocodes, avoid coercion warning
  temp.df$language_glottocode <- as.character(temp.df$language_glottocode)
  animal.present.df <- full_join(animal.present.df, temp.df, by = c('language_glottocode'))
  rm(temp.df) # clean
}
rm(sc) # clean
```

### Subsistence (SCCS)

We are going to load data on the main mode of subsistence (pastoralism versus
agriculture). We will include the following data sets:

- Principal Subsistence Category [SCCS820](https://d-place.org/parameters/SCCS820)
	1. Gathering (n=10)
	2. Hunting (n=16)
	3. Fishing (n=23)
	4. Incipient Agriculture (n=18)
	5. Domestic Animals (n=16)
	6. Extensive Agriculture (n=46)
	7. Intensive Agriculture (n=55)
	8. Trade (n=2)
- Subsistence Economy: Dominant Mode [SCCS833](https://d-place.org/parameters/SCCS833)
	1. Advanced Agriculture (n=56)
	2. Horticulture (n=18)
	3. Simple or Shifting Cultivation (n=51)
	4. Domestic Animals (n=15)
	5. Exchange (n=1)
	6. Fishing (n=17)
	7. Gathering (n=13)
	8. Hunting (n=15)
- Subsistence economy: dominant activity [EA042](https://d-place.org/parameters/EA042)
	1. Gathering (n=103)
	2. Fishing (n=121)
	3. Hunting (n=77)
	4. Pastoralism (n=77)
	5. Casual agriculture (n=0)
	6. Extensive agriculture (n=477)
	7. Intensive agriculture (n=282)
	8. Two or more sources (n=67)
	9. Agriculture, type unknown (n=86)
- Primary Source of Subsistence [SCCS1716](https://d-place.org/parameters/SCCS1716)
	1. intensive agriculture (n=26)
	2. extensive agriculture (n=37)
	3. animal husbandry (n=9)
	4. fishing (n=11)
	5. hunting (n=9)
	6. gathering (n=8)
	7. trade (n=0)
	8. wage labor (n=1)
- Subsistence economy: dominant activity [Note, identical to EA042] [SCCS246](https://d-place.org/parameters/SCCS246)
	1. Gathering (n=15)
	2. Fishing (n=19)
	3. Hunting (n=12)
	4. Pastoralism (n=16)
	5. Casual agriculture (n=0)
	6. Extensive agriculture (n=60)
	7. Intensive agriculture (n=57)
	8. Two or more sources (n=7)
	9. Agriculture, type unknown (n=0)
- Scale 3- Agriculture [SCCS151](https://d-place.org/parameters/SCCS151)
	1. None (n=38)
	2. 10% food supply (n=17)
	3. 10 %; secondary (n=11)
	4. Primary; not intensive (n=63)
	5. Primary; intensive (n=57)

```{r subsistence}
frameize <- function(frame1,frame2,predicate) {
  result <- data.frame(
    language_glottocode = c(
      as.character(frame1$language_glottocode),
      as.character(frame2$language_glottocode)
    )
  )
  result[[predicate]] <- c(
    frame1[[predicate]],
    frame2[[predicate]]
  )
  return(result)
}

# read csv files
subsistence.SCCS820  <- read.sccs.variable('SCCS820')
subsistence.SCCS833  <- read.sccs.variable('SCCS833')
#subsistence.EA042    <- read.sccs.variable('EA042')
subsistence.SCCS1716 <- read.sccs.variable('SCCS1716')
subsistence.SCCS246  <- read.sccs.variable('SCCS246')
subsistence.SCCS151  <- read.sccs.variable('SCCS151')

# extract agriculturalism
subsistence.SCCS820.agriculturalism  <- filter(subsistence.SCCS820, SCCS820 == 4 | SCCS820 == 6 | SCCS820 == 7 )
subsistence.SCCS833.agriculturalism  <- filter(subsistence.SCCS833, SCCS833 < 4)
#subsistence.EA042.agriculturalism    <- filter(subsistence.EA042, EA042 > 4 & EA042 < 8)
subsistence.SCCS1716.agriculturalism <- filter(subsistence.SCCS1716, SCCS1716 < 3)
subsistence.SCCS246.agriculturalism  <- filter(subsistence.SCCS246, SCCS246 > 4 & SCCS246 < 8 )
subsistence.SCCS151.agriculturalism  <- filter(subsistence.SCCS151, SCCS151 > 3)

subsistence.SCCS820.agriculturalism$subsistence.SCCS820   <- 'agriculture'
subsistence.SCCS833.agriculturalism$subsistence.SCCS833   <- 'agriculture'
#subsistence.EA042.agriculturalism$subsistence.EA042       <- 'agriculture'
subsistence.SCCS1716.agriculturalism$subsistence.SCCS1716 <- 'agriculture'
subsistence.SCCS246.agriculturalism$subsistence.SCCS246   <- 'agriculture'
subsistence.SCCS151.agriculturalism$subsistence.SCCS151   <- 'agriculture'

# extract pastoralism
subsistence.SCCS820.pastoralism  <- filter(subsistence.SCCS820, SCCS820 != 4 & SCCS820 != 6 & SCCS820 != 7 )
subsistence.SCCS833.pastoralism  <- filter(subsistence.SCCS833, SCCS833 > 3 )
#subsistence.EA042.pastoralism    <- filter(subsistence.EA042, EA042 == 4)
subsistence.SCCS1716.pastoralism <- filter(subsistence.SCCS1716, SCCS1716 > 2)
subsistence.SCCS246.pastoralism  <- filter(subsistence.SCCS246, SCCS246 < 5)
subsistence.SCCS151.pastoralism  <- filter(subsistence.SCCS151, SCCS151 < 4)

subsistence.SCCS820.pastoralism$subsistence.SCCS820   <- 'other'
subsistence.SCCS833.pastoralism$subsistence.SCCS833   <- 'other'
#subsistence.EA042.pastoralism$subsistence.EA042       <- 'other'
subsistence.SCCS1716.pastoralism$subsistence.SCCS1716 <- 'other'
subsistence.SCCS246.pastoralism$subsistence.SCCS246   <- 'other'
subsistence.SCCS151.pastoralism$subsistence.SCCS151   <- 'other'


# merge
subsistence.SCCS820  <- frameize(subsistence.SCCS820.pastoralism,  subsistence.SCCS820.agriculturalism,  'subsistence.SCCS820')
subsistence.SCCS833  <- frameize(subsistence.SCCS833.pastoralism,  subsistence.SCCS833.agriculturalism,  'subsistence.SCCS833')
#subsistence.EA042    <- frameize(subsistence.EA042.pastoralism,    subsistence.EA042.agriculturalism,    'subsistence.EA042')
subsistence.SCCS1716 <- frameize(subsistence.SCCS1716.pastoralism, subsistence.SCCS1716.agriculturalism, 'subsistence.SCCS1716')
subsistence.SCCS246  <- frameize(subsistence.SCCS246.pastoralism,  subsistence.SCCS246.agriculturalism,  'subsistence.SCCS246')
subsistence.SCCS151  <- frameize(subsistence.SCCS151.pastoralism,  subsistence.SCCS151.agriculturalism,  'subsistence.SCCS151')
rm(
  subsistence.SCCS820.pastoralism,
  subsistence.SCCS820.agriculturalism,
  subsistence.SCCS833.pastoralism,
  subsistence.SCCS833.agriculturalism,
#  subsistence.EA042.pastoralism,
#  subsistence.EA042.agriculturalism,
  subsistence.SCCS1716.pastoralism,
  subsistence.SCCS1716.agriculturalism,
  subsistence.SCCS246.pastoralism,
  subsistence.SCCS246.agriculturalism,
  subsistence.SCCS151.pastoralism,  
  subsistence.SCCS151.agriculturalism
)

```


### Demographics (SCCS)

We interpret niche cluster 7 as a less productive, drier environment. Hence, 
these have relatively low population densities. We read these data:

- Population Density (from Pryor Data) [SCCS1130](https://d-place.org/parameters/SCCS1130)
	2. less than 1 per square mile (n=51), `low`
	3. 1-4.9 per square mile (n=28), `high`
	4. 5-24.9 per square mile (n=35), `high`
	5. 25-99.9 per square mile (n=37), `high`
	6. 99-499.9 per square mile (n=24), `high`
	7. 500 or more per square mile (n=11), `high`


```{r pop_density}
# update column names after loading
pop.density.SCCS1130 <- read.sccs.variable('SCCS1130')
names(pop.density.SCCS1130) <- c('language_glottocode', 'pop.density.SCCS1130')

# split in low and high density subsets
smallpop <- filter(pop.density.SCCS1130, pop.density.SCCS1130 == 2)
smallpop$pop.density <- 'low'
bigpop <- filter(pop.density.SCCS1130, pop.density.SCCS1130 > 2)
bigpop$pop.density <- 'high'

# concatenate into new data frame
pop.density.SCCS1130 <- data.frame(
  language_glottocode = c(
    as.character(smallpop$language_glottocode),
    as.character(bigpop$language_glottocode)
  ),
  pop.density.SCCS1130 = c(
    smallpop$pop.density,
    bigpop$pop.density
  )
)
rm(smallpop,bigpop)
plot(pop.density.SCCS1130$pop.density.SCCS1130, main = "Population density")
```

### Sedentism (SCCS)

We interpret niche cluster 5 as more productive, and therefore more conducive
to sedentism. We download the following data set:

- Scale 2 - Fixity of Residence [SCCS150](https://d-place.org/parameters/SCCS150)

The data is ordinal, with the following values:

1. Nomadic (n=28)
2. Seminomadic (n=21)
3. Semisedentary (n=20)
4. Sedentary; impermanent (n=15)
5. Sedentary (n=102)

We will lump the first four states together and contrast those with the only
fully sedentary state.

```{r sedentism}
sedentism.SCCS150 <- read.sccs.variable('SCCS150')

# make a subframe with all records that are not truly sedentary (i.e. anything
# below 3), and assign these the string value 'impermanent'
sedentism.SCCS150.nomadic <- filter(sedentism.SCCS150, SCCS150 < 3)
sedentism.SCCS150.nomadic$SCCS150 <- "impermanent"

# make a subframe with only those records that are truly sedentary (i.e. value
# is > 2), and assign these the string value 'sedentary'
sedentism.SCCS150.sedentary <- filter(sedentism.SCCS150, SCCS150 > 2)
sedentism.SCCS150.sedentary$SCCS150 <- "sedentary"

# concatenate the subframes into a new frame. Don't turn strings into factors
# just yet so that the glottocodes continue to be strings, then coerce the
# sedentary/impermanent states to factors.
sedentism.SCCS150 <- data.frame( 
  language_glottocode = c( 
    as.character(sedentism.SCCS150.nomadic$language_glottocode), 
    as.character(sedentism.SCCS150.sedentary$language_glottocode) 
  ), 
  sedentism.SCCS150 = c(
    sedentism.SCCS150.nomadic$SCCS150,
    sedentism.SCCS150.sedentary$SCCS150
  )
)

# clean up namespace, free memory
rm(sedentism.SCCS150.nomadic, sedentism.SCCS150.sedentary)
plot(sedentism.SCCS150$sedentism.SCCS150, main = "Fixity of residence")
```

### Stratification (SCCS)

Sedentary societies with high population densities should be conducive to local
stratification. We download the following data:

- Scale 10 - Social Stratification [SCCS158](https://d-place.org/parameters/SCCS158),
  which is ordinal and has the following levels:
  1. Egalitarian (n=65), `egalitarian`
  2. Hereditary slavery	(n=52), `egalitarian`
  3. 2 social classes, no castes/slavery (n=19), `stratified`
  4. 2 social classes, castes/slavery (n=20), `stratified`
  5. 3 social classes or castes, w/ or w/out slavery (n=30), `stratified`

```{r stratification}
stratification.SCCS158 <- read.sccs.variable('SCCS158')

# make a subframe only with records that are truly egalitarian
# and assign these the string (i.e. factor) value 'egalitarian'
stratification.SCCS158.egalitarian <- filter( stratification.SCCS158, SCCS158 < 3 )
stratification.SCCS158.egalitarian$SCCS158 <- "egalitarian"

# make a subframe with all records that are not truly egalitarian (values 2..5),
# and assign these the string (i.e. factor) value 'stratified'
stratification.SCCS158.stratified <- filter( stratification.SCCS158, SCCS158 > 2 )
stratification.SCCS158.stratified$SCCS158 <- "stratified"

# concatenate the subframes into a new frame. Don't make factors here so that
# we keep the glottocodes as strings, but coerce the egalitarian/stratification
# string values afterwards.
stratification.SCCS158 <- data.frame(
  language_glottocode = c(
    as.character(stratification.SCCS158.egalitarian$language_glottocode),
    as.character(stratification.SCCS158.stratified$language_glottocode)
  ),
  stratification.SCCS158 = c(
    stratification.SCCS158.egalitarian$SCCS158,
    stratification.SCCS158.stratified$SCCS158
  )
)

# clean up namespace, free memory
rm(
  stratification.SCCS158.egalitarian,
  stratification.SCCS158.stratified
)

# same for SCCS1751
stratification.SCCS1751 <- read.sccs.variable('SCCS1751')
stratification.SCCS1751.egalitarian <- filter( stratification.SCCS1751, SCCS1751 < 3 )
stratification.SCCS1751.egalitarian$SCCS1751 <- 'egalitarian'
stratification.SCCS1751.stratified <- filter( stratification.SCCS1751, SCCS1751 > 2 )
stratification.SCCS1751.stratified$SCCS1751 <- 'stratified'
stratification.SCCS1751 <- data.frame(
  language_glottocode = c(
    as.character(stratification.SCCS1751.egalitarian$language_glottocode),
    as.character(stratification.SCCS1751.stratified$language_glottocode)
  ),
  stratification.SCCS1751 = c(
    stratification.SCCS1751.egalitarian$SCCS1751,
    stratification.SCCS1751.stratified$SCCS1751
  )
)
rm(stratification.SCCS1751.egalitarian,stratification.SCCS1751.stratified)

# same for Social Stratification in the Local Community [SCCS710]
stratification.SCCS710 <- read.sccs.variable('SCCS710')
stratification.SCCS710.egalitarian <- filter( stratification.SCCS710, SCCS710 < 2 )
stratification.SCCS710.egalitarian$SCCS710 <- 'egalitarian'
stratification.SCCS710.stratified <- filter( stratification.SCCS710, SCCS710 > 1 )
stratification.SCCS710.stratified$SCCS710 <- 'stratified'
stratification.SCCS710 <- data.frame(
  language_glottocode = c(
    as.character(stratification.SCCS710.egalitarian$language_glottocode),
    as.character(stratification.SCCS710.stratified$language_glottocode)
  ),
  stratification.SCCS710 = c(
    stratification.SCCS710.egalitarian$SCCS710,
    stratification.SCCS710.stratified$SCCS710
  )
)
rm(stratification.SCCS710.egalitarian,stratification.SCCS710.stratified)

# Class differentiation: primary [EA066]
stratification.EA066 <- read_csv(sprintf('%s/data/d-place/EA066.csv',REPO_HOME))
stratification.EA066 <- unique(stratification.EA066)
stratification.EA066 <- select(left_join(animal.present.df,stratification.EA066),language_glottocode,code)
stratification.EA066.egalitarian <- filter(stratification.EA066, code < 3)
stratification.EA066.egalitarian$code <- 'egalitarian'
stratification.EA066.egalitarian <- unique(stratification.EA066.egalitarian)
stratification.EA066.stratified <- filter(stratification.EA066, code > 2)
stratification.EA066.stratified$code <- 'stratified'
stratification.EA066.stratified <- unique(stratification.EA066.stratified)
stratification.EA066 <- data.frame(
  language_glottocode = c(
    as.character(stratification.EA066.egalitarian$language_glottocode),
    as.character(stratification.EA066.stratified$language_glottocode)
  ),
  stratification.EA066 = c(
    stratification.EA066.egalitarian$code,
    stratification.EA066.stratified$code
  )
)
stratification.EA066 <- filter(stratification.EA066, language_glottocode != 'boro1274' & language_glottocode != 'kiku1240')
rm(stratification.EA066.egalitarian,stratification.EA066.stratified)

# plot results
plot(stratification.SCCS158$stratification.SCCS158, main = "Social stratification [SCCS158]")

```

### Join data

We left-join the different data frames on the `language_glottocode` columns:

```{r leftjoin}
# join data sets
joined.df <- dplyr::left_join(pop.density.SCCS1130, sedentism.SCCS150, by=c('language_glottocode'))
joined.df <- dplyr::left_join(joined.df, stratification.SCCS158, by=c('language_glottocode'))
joined.df <- dplyr::left_join(joined.df, stratification.SCCS1751, by=c('language_glottocode'))
joined.df <- dplyr::left_join(joined.df, stratification.SCCS710, by=c('language_glottocode'))
joined.df <- dplyr::left_join(joined.df, stratification.EA066, by=c('language_glottocode'))
joined.df <- dplyr::left_join(joined.df, subsistence.SCCS1716, by=c('language_glottocode'))
joined.df <- dplyr::left_join(joined.df, subsistence.SCCS246, by=c('language_glottocode'))
joined.df <- dplyr::left_join(joined.df, subsistence.SCCS820, by=c('language_glottocode'))
joined.df <- dplyr::left_join(joined.df, subsistence.SCCS833, by=c('language_glottocode'))
joined.df <- dplyr::left_join(joined.df, subsistence.SCCS151, by=c('language_glottocode'))

# assign row names
societies <- row.names(joined.df) <- joined.df$language_glottocode

# grouping of domesticated taxa by niche cluster
clusters <- list(
  cluster7 = c(
    'Capra_hircus_aegagrus',
    'Camelus_bactrianus', 
    'Camelus_dromedarius',
    'Equus_przewalskii',
    'Bos_grunniens_mutus'
  ),
  cluster5 = c(
    'Ovis_aries_orientalis',
    'Lama_glama_guanicoe', 
    'Vicagna_vicugna',
    'Bos_javanicus',
    'Bos_taurus_primigenius'
  ),
  cluster2 = c(
    'Bubalus_bubalis_arnee', # 5
    'Bos_frontalis_gaurus'   # 5
  ),
  cluster14 = c('Sus_scrofa'), # 5
  cluster26 = c('Rangifer_tarandus'), # 7
  cluster0 = c('Equus_africanus') # 7
)

# create and populate matrix
summarized.matrix <- matrix(
  nrow = length(animal.present.df$language_glottocode), 
  ncol = length(names(clusters))+1,
  dimnames = list(animal.present.df$language_glottocode,c(names(clusters),'clusters'))
)
row.names(animal.present.df) <- animal.present.df$language_glottocode
for ( soc in animal.present.df$language_glottocode ) {
  for ( cluster in names(clusters) ) {
    taxa <- clusters[[cluster]]
    total <- sum(animal.present.df[soc,taxa],na.rm=T)
    if ( total == 0 ) {
      summarized.matrix[soc,cluster] <- 0
    } else {
      summarized.matrix[soc,cluster] <- 1
    }
  }
}
summarized.df <- as.data.frame(summarized.matrix)
summarized.df$language_glottocode <- row.names(summarized.df)

# filter on steppe/pasture
steppe <- filter(filter(summarized.df, cluster7 == 1 | cluster26 == 1),
                 cluster5 == 0 & cluster2 == 0 & cluster14 == 0)
steppe$clusters <- 'steppe'
pasture <- filter(summarized.df, cluster5 == 1 | cluster2 == 1 | cluster14 == 1)
pasture$clusters <- 'pasture'
filtered <- data.frame(
  language_glottocode = c(
    steppe$language_glottocode,
    pasture$language_glottocode
  ),
  habitat = c(
    steppe$clusters,
    pasture$clusters
  )
)

joined.df <- dplyr::left_join(filtered,joined.df,by=c('language_glottocode'))

# clean up columns
joined.df <- select(
  joined.df,
  language_glottocode,
  pop.density.SCCS1130,
  sedentism.SCCS150,
  stratification.SCCS158,
  stratification.SCCS1751,
  stratification.SCCS710,
  stratification.EA066,
  subsistence.SCCS1716,
  subsistence.SCCS246,
  subsistence.SCCS820,
  subsistence.SCCS833,
  subsistence.SCCS151,
  habitat
)
rm(filtered,steppe,pasture,summarized.df,summarized.matrix,total,taxa,soc,cluster)
plot(joined.df$habitat)
```

### Load language tree

To correct for Galton's Problem in the comparative analysis of societies we must
correct for their relatedness. For this we use the following tree:

- [Phylogeny Global Classification (Glottolog 3.3)](https://d-place.org/phylogenys/glottolog_global)

```{r tree}
tree.file <- sprintf('%s/data/d-place/glottolog_global.trees',REPO_HOME)
tree <- read.tree( file = tree.file )

# intersection between tree tips and societies in data matrix
societies <- tree$tip.label[
  match(
    joined.df$language_glottocode, 
    tree$tip.label, 
    nomatch = 0
  )
]

# prune tree
tree <- keep.tip(tree, societies)

# prune matrix
row.names(joined.df) <- joined.df$language_glottocode
joined.df <- joined.df[societies,]

plot( tree, show.tip.label = F, use.edge.length = F, type = "fan" )
```
### Hypothesis testing: abiotic niche clusters as predictors 

Given the data and the tree, we can now perform a succession of hypothesis 
tests. Firstly, we will test the relationship between the niche clusters
(the provisional names we've chosen for this are `steppe` and `pasture`), which
we expect to predict the population density: the steppe supports much lower
population densities than the more fertile lands. The relationship doesn't work
the other way, the population density does not predict the niche cluster; hence,
the former is the dependent variable:

```{r niche_pop_density}
# the variables need to be named integer vectors whose values can only be 0 or 1
habitat <- as.integer(joined.df$habitat)-1
names(habitat) <- joined.df$language_glottocode
pop.density <- as.integer(joined.df$pop.density.SCCS1130)-1
names(pop.density) <- joined.df$language_glottocode

# note the dependent variable:
pagel <- fitPagel(tree, habitat, pop.density, dep.var = "y")
pagel
plot(pagel)
```

The relationship is highly significant: p=0.000008296388**

Next, we check to see how the niche clusters correlate with sedentism. The 
`steppe` species are roaming animals with large ranges (note: it would be worth
checking this using species traits), compatible with a nomadic lifestyle.
The `pasture` species are less migratory and live in habitats where agriculture
can be productive enough to enable sedentism. Hence, the lifestyle depends on
the niche cluster of domesticated animals:

```{r niche_sedentism}
# the variables need to be named integer vectors whose values can only be 0 or 1
sedentism <- as.integer(joined.df$sedentism.SCCS150)-1
names(sedentism) <- joined.df$language_glottocode

# note the dependent variable:
pagel <- fitPagel(tree, habitat, sedentism, dep.var = "y")
pagel
plot(pagel)
```

The relationship is significant: p=0.01743421*

The relationship between the niche clusters and sedentism is mediated by the
subsistence mode(s). A sedentary lifestyle is possible almost exclusively in the
context of agriculture. Conversely, impermanent settlement is the rule when the
food source needs to be followed, e.g. as in pastoral nomadism, but also when 
`steppe` animals are complementary to hunting or fishing. Hence, we expect the
same relationship between habitat and subsistence as between the former and 
sedentism:

```{r niche_subsistence}
# the variables need to be named integer vectors whose values can only be 0 or 1
subsistence <- as.integer(as.factor(joined.df$subsistence.SCCS833))-1
names(subsistence) <- joined.df$language_glottocode

# note the dependent variable:
pagel <- fitPagel(tree, habitat, subsistence, dep.var = "y")
pagel
plot(pagel)
```

The relationship is indeed significant (the patterns in the data are such that
the p value is actually identical): p=0.01743421*

Another, more direct way to go about this is to directly test the relationship
where sedentism is made possible by (i.e. depends on) subsistence, which is 
highly significant:

```{r sedentism_subsistence}
pagel <- fitPagel(tree, subsistence, sedentism, dep.var = "y")
pagel
plot(pagel)
```

p=0.0000000001493953**

Finally, there is evidence to suggest that social stratification correlates with
animal husbandry (in general). Here we test whether this differs between the 
different niche clusters such that there is a direct correlation between 
clusters and social complexity is something we can test:

```{r niche_stratification}
# the variables need to be named integer vectors whose values can only be 0 or 1
stratification <- as.integer(joined.df$stratification.SCCS158)-1 
names(stratification) <- joined.df$language_glottocode

pagel <- fitPagel(tree,habitat,stratification, dep.var = "y")
pagel
plot(pagel)
```

It turns out that this is not a significant correlation (p=0.05345692), though
it's not far off.

```{r fitpagel}


pagel <- fitPagel(tree,subsistence,stratification, dep.var = "y")
pagel
plot(pagel)
```

### Hypothesis testing: population density classes as predictors

Having established that the abiotic niche clusters strongly correlate with
population density, we now assess how the latter in turn correlates with 
the other variables: sedentism, subsistence, and stratification. To begin with,
population density and sedentism are linked: the societies with the highest
population densities are sedentary. Or, conversely, a society that can afford
to be sedentary due to its subsistence strategy can probably achieve high
densities. Let's be agnostic about what the dependent variable is:

```{r pop_sedentism}
pagel <- fitPagel(tree,pop.density,sedentism)
pagel
plot(pagel)
```

The two are strongly correlated: p=0.005100353**

How does population density interact with subsistence? The causal link would
appear to be more straightforward: _thanks to_ agriculture, the population
density can go up:

```{r pop_subsistence}
pagel <- fitPagel(tree,pop.density,subsistence,dep.var = "x")
pagel
plot(pagel)
```

Tested in this way, the two are highly correlated: p=0.0008374419**, which is
an order of magnitude higher than the link between population density and
sedentism (which are linked by way of the subsistence strategy, probably).

What about the relationship between population density and stratification? Let's
model stratification as something that depends on population density:

```{r pop_stratification}
pagel <- fitPagel(tree,pop.density,stratification, dep.var = "x")
pagel
plot(pagel)
```

Modeled in this way (or, indeed, in any way), the two are not correlated:
p=0.1084162.

### Hypothesis testing: sedentism as a predictor variable

Sedentism and subsistence should be tightly linked: agriculture is in most forms
a sedentary practice. Conversely, a sedentary lifestyle can mostly only be 
afforded thanks to agriculture, at least in comparison with pastoralism, hunting
and gathering, or fishing.