Relatives and Neighbours
Muhammad Rehan

Relatives and Neighbours

Eight follow-along R tutorials on modelling phylogenetic and spatial autocorrelation in cross-linguistic data. Every model is fitted with an existing package: no custom Stan programs, no long data wrangling.

Start with Plate 1
World map of the 500 languages in the tutorial sampleLanguages with ejective consonants form regional clusters: the Caucasus, the Ethiopian highlands, southern Africa, the North American Cordillera, Mesoamerica and the Andes.

The sample: 500 languages

Each symbol is one language in the tutorial sample. Inventories from PHOIBLE 2.0, locations from Glottolog, Equal Earth projection.

The plates

The map above is the whole problem in one picture. Languages with ejectives are not scattered at random: they sit together in the Caucasus, the Ethiopian highlands, southern Africa, the North American Cordillera and the Andes. Neighbours share features through contact and relatives share them through inheritance, so 500 languages carry far less than 500 languages’ worth of independent evidence. A regression that ignores this is too sure of itself.

Each plate fixes that in one way, with one tool. They share a single running example, so the same slope can be followed from the first model to the last. Read them in order the first time; after that each one stands alone.

  1. Diagnose
  2. 1Seeing the problemA plain logistic regression, a simulation of how often it cries wolf, Moran's I and phylogenetic signal.base R, spdep, phytools · under a minute
  3. Model it in brms
  4. 2Phylogeny as a covariance matrixA family random intercept, then a phylogenetic random effect with gr(cov = A), and how to read phylogenetic signal off the fit.brms · about 3 minutes, mostly compiling
  5. 3Space as a smooth surfaceAn approximate Gaussian process over the coordinates with gp(), how to choose its settings, and a map of the surface it finds.brms · about 5 minutes, plus 5 for an optional check
  6. 4Space as a network of neighboursA neighbour graph from spdep, then ICAR and BYM2 models with car(), and a check on how much the graph matters.spdep and brms · about 6 minutes
  7. 5Both at onceA joint model with a phylogenetic effect and a Gaussian process, how the two share the credit, and one figure comparing every slope.brms · about 8 minutes
  8. Other engines
  9. 6The same models in INLAiid, generic0 and bym2 effects, an SPDE surface on the sphere, a joint model, and a check that catches the approximation failing.INLA · about 3 minutes
  10. 7The same models as a GAMRandom-effect, Markov random field and spline-on-the-sphere smooths in gam(), alone and together.mgcv · under a minute in total
  11. 8The classical toolkit: PGLS and SARphylolm() for continuous and binary outcomes on a tree, errorsarlm() and lagsarlm() on a neighbour graph, and when these are enough.phylolm and spatialreg · a few seconds

Before you start

You need R 4.3 or later. Install the packages once. The first block covers Plates 1, 7 and 8; the second and third add the two Bayesian engines.

install.packages(c("ape", "phytools", "phylolm", "nlme", "sf", "spdep",
                   "spatialreg", "rnaturalearth", "rnaturalearthdata",
                   "ggplot2", "mgcv"))

# the tree figures on Plates 1 and 2 use ggtree, from Bioconductor
install.packages("BiocManager")
BiocManager::install("ggtree")

# brms (Plates 2 to 5) compiles its models with Stan, through cmdstanr
install.packages("brms")
install.packages("cmdstanr",
                 repos = c("https://stan-dev.r-universe.dev", getOption("repos")))
cmdstanr::install_cmdstan()

# INLA (Plate 6) lives in its own repository
install.packages("INLA", dependencies = TRUE,
                 repos = c(getOption("repos"),
                           INLA = "https://inla.r-inla-download.org/R/stable"))

On Windows, cmdstanr needs Rtools; cmdstanr::check_cmdstan_toolchain() tells you whether yours is ready.

Every plate starts by reading the same two files, straight from this site:

library(ape)

base <- "https://rehan-muh.github.io/tutorials/data/"
d    <- read.csv(paste0(base, "languages.csv"))
tree <- read.tree(paste0(base, "glottolog_tree.nwk"))

The figures on these pages use a small ggplot2 theme. Your plots will show the same thing in ggplot’s default look; if you want them to match, download atlas_theme.R and source() it.

The running example

In 2013 Caleb Everett reported that languages spoken at high altitude are more likely to have ejective consonants, sounds such as [kʼ] and [tʼ] made by compressing air above a closed glottis. The proposed reason is articulatory: thinner air makes the compression cheaper. The claim is a good teaching case because the raw association is strong, the mechanism is plausible, and the data are exactly the kind where relatives and neighbours can manufacture a correlation by themselves.

The sample is 500 languages drawn at random from the 2,026 in PHOIBLE 2.0 that Glottolog locates on the map. For each language the file has:

Column Meaning
glottocode, name Glottolog identifier and name
family, macroarea Top-level Glottolog family and macroarea
lat, lon Glottolog’s point location
elevation Metres above sea level at that point, from ETOPO1
ejectives 1 if the PHOIBLE inventory has at least one ejective consonant
n_consonants, n_vowels Sizes of the consonant and vowel inventories
tone 1 if the inventory lists tones

The tree is Glottolog’s classification turned into a phylogeny. Each family is a clade shaped by Glottolog’s subgroups. Every level of the classification counts as one step, the family itself included, and all languages end at the same height, so two languages share as much of the tree as they share levels of classification. Families join only at the root: Glottolog makes no claim about relationships between families, so languages in different families share no history in this tree.

Three limits are worth knowing before you read any estimate off these pages. A language is one point, though its speakers occupy an area. The elevation is the height at that point on a grid of about two kilometres. And a classification with conventional branch lengths is a rough stand-in for a dated phylogeny. The tutorials teach the models; they are not a verdict on Everett’s hypothesis.

The script that builds both files from the sources is make_data.R. The files themselves are languages.csv and glottolog_tree.nwk.

Sources

  • Moran, Steven and Daniel McCloy (eds.) 2019. PHOIBLE 2.0. Jena: Max Planck Institute for the Science of Human History. CC BY-SA 3.0.
  • Hammarström, Harald, Robert Forkel, Martin Haspelmath and Sebastian Bank. Glottolog. Leipzig: Max Planck Institute for Evolutionary Anthropology. CC BY 4.0.
  • Amante, Christopher and Barry Eakins 2009. ETOPO1 1 arc-minute global relief model. NOAA. Public domain.
  • Everett, Caleb 2013. Evidence for direct geographic influences on linguistic sounds: the case of ejectives. PLoS ONE 8(6): e65275.