This is a blog post that I’ve had in mind for awhile, and as I was working on it I realized that there is a bit of groundwork needed. Part of that groundwork was laid in a previous post where I showed a simple way to control for language family and area with categorical variables in a linear mixed model using Python.
Another part of the groundwork is understanding the effect of language area, which I touched on in a more recent post. One drawback of the macroarea approach is that it treats all languages within a given area as having the same relationship. Since we know that languages closer to each other in an area are more likely to share similarities than those farther apart, this is problematic. We can of course visualize various geographical relationships, but how do we know whether those relationships have anything to do with features in question? What we would prefer to model is how geographically close languages are to one another, and whether this proximity affects how typological features are distributed as well.
Modeling relationships using gaussian process
These kind of dependencies can be modeled using a Gaussian Process. I won’t go into detail in this post about how these work (if you are interested, here is a decent explanation for starters), but essentially for our purposes (geographic relatedness) they function similarly to a covariance matrix - a grid that identifies how each point is related (via distance) to the others.
For the remainder of this post, I will be
illustrating how you can control for language
area in typological studies by integrating a
gaussian process using GPS coordinates. This can
be considered a companion to the blog posts
linked above, as it uses the same data - a
single coded universal from Verkerk
et al (2025). In a future post I may try to
replicate their results using the
gpboost library, which I will be
describing below.
Getting the data
As with previous tutorials, I’m assuming that you’ve set up a Python 3.10 environment and have some familiarity with writing and running Python scripts (i.e. via a terminal). For a basic setup tutorial, you can follow these instructions.
Following my earlier post, we are again
selecting the purported universal #005KA
(“If a language has dominant SOV order and
the genitive follows the governing noun, then
the adjective likewise follows the noun.”).
This was found to be significant by Verkerk et
al in a Bayesian analysis without family/area
controls, but when language family and area were
controlled for, the effect was non-significant.
The datasheet for this feature can
be found here. This datasheet does not
contain GPS coordinates, family membership, or
macroarea identifiers, but since it does contain
glottocodes, we can link it to the Glottolog
dataset. The following code uses a
downloader function to get both these files and
write them to your local directory, then uses
pandas to read them into dataframes
and combine them.
import os, requests # for file handling
import pandas as pd # import pandas for tabular data
def get_file_from_url(url, localpath):
"""Download a url to a local path"""
# check whether the file exists already
if not os.path.isfile(localpath):
response = requests.get(url) # fetch the file content
if response.status_code == 200:
with open(localpath, "wb") as file:
file.write(response.content) # write content to the local file
print("Download complete!")
else:
print(f"Failed to download. Status code: {response.status_code}")
else:
print("File exists!")
# first we get the coded data for this feature
featfile = "0005KA-BT_data.txt" # local name for tab delimited file with relevant features
# use the raw GitHub URL for this file
url = "https://raw.githubusercontent.com/SimonGreenhill/TestingLinguisticUniversals/43da16f/results/0005KA/BT_data.txt"
get_file_from_url(url, featfile) # get the file from the url
# read the csv and name the columns
fdf = pd.read_csv(featfile, delimiter="\t", header=None, names=["glottocode", "GB133_GB065", "GB193"])
print(fdf.head()) # visualize the data in the terminal
# glottocode GB133_GB065 GB193
# 0 abun1252 0 1
# 1 area1240 0 1
# 2 axam1237 0 1
# 3 cash1251 1 0
# 4 chuk1273 0 0
# ...
feat = featfile.split("-")[0] # get the name of the feature from the filename
print(f"There are {len(fdf)} Glottocodes in the {feat} feature dataset.")
# There are 1787 Glottocodes in the 0005KA feature dataset.
# now we want to get the information on languages and their coordinates
glottofile = "Glottolog_Languages.csv" # local name for comma-separated file with glottocodes and other info
# the info for the languages is conveniently located in the same repo
url = "https://raw.githubusercontent.com/SimonGreenhill/TestingLinguisticUniversals/43da16f/results/Glottolog_Languages.csv"
get_file_from_url(url, glottofile) # get the file from the url
gldf = pd.read_csv(glottofile) # read in the dataset
print(gldf.head()) # check its structure
# ID Name
# 0 abkh1242 Abkhaz-Adyge
# 1 surm1244 Surmic
# 2 tama1329 Tamaic
# 3 yare1250 Yareban
# 4 monu1249 Bogia
# ...
print(f'There are {len(gldf)} languoids in the dataset')
# There are 21508 languoids in the datasetExamining and comparing the data
We can compare the two datasets
programmatically to ensure that the data we need
can be matched up. The following code prints out
some information about the columns, and then
does a transformation to ensure that columns
will match. We then get the intersection of the
glottocode column to combine the
data, and print out some information about the
results.
print(fdf.columns) # columns in the feature dataframe
print(gldf.columns) # columns in the glottolog dataframe
gldf = gldf.dropna(subset=['latitude']) # drop languoids without latitude
print(f"There are {len(gldf)} languoids in the dataset with latitude info")
# There are 8704 languoids in the dataset with latitude info
# combine the datasets, keeping only matching glottocodes
intersect = pd.merge(gldf, fdf, on='glottocode', how='inner') # merge
print(intersect.head()) # visualize
# glottocode name isocodes level ... Family_ID Family_name GB133_GB065 GB193
# 0 abad1241 Abadi kbt language ... aust1307 Austronesian 1 1
# 1 abee1242 Abé aba language ... atla1278 Atlantic-Congo 0 1
# 2 abkh1244 Abkhaz abk language ... abkh1242 Abkhaz-Adyge 0 0
# 3 abua1245 Abu' Arapesh aah language ... nucl1708 Nuclear Torricelli 0 1
# 4 abun1252 Abun kgr language ... NaN NaN 0 1
print(intersect.columns) # the combined column names
# Index(['glottocode', 'name', 'isocodes', 'level', 'macroarea', 'latitude',
# 'longitude', 'Family_ID', 'Family_name', 'GB133_GB065', 'GB193'],
# dtype='object')
print(f"There are {len(intersect)} languoids for assessment")
# There are 1786 languoids for assessment
# print out some information about the data
nfams = len(intersect['Family_ID'].value_counts()) # these are the unique language families
nisols = len(intersect[intersect['Family_ID'].isna()]) # rows in this dataset with "NaN" values are isolates
print(f"There are {len(intersect['macroarea'].value_counts())} macroareas")
print(f"There are {nfams} language families in the dataset")
print(f"There are {nisols} isolates in the dataset")
print(f"There are {nfams+nisols} total `families` including isolates")
# There are 6 macroareas
# There are 175 language families in the dataset
# There are 57 isolates in the dataset
# There are 232 total `families` including isolatesOne observation here is that there are a
number of isolates in the data,
i.e. languages that have no known relatives.
There are three basic ways to handle isolates -
remove them, make them a single group, or treat
each as unique. Each has various implications,
as we would expect a language family control to
have differing impact on various linguistic
features if there are more or less members in
each group. For our purposes, the simplest way
to deal with this is to remove these languages
from observation, as in the following code.
# remove the languoids without Family info from the dataset and create a new dataframe called `df`
df = intersect[~intersect['Family_ID'].isna()]
print(f"There are {len(df)} non-isolate languoids for assessment")
# There are 1729 non-isolate languoids for assessmentUsing the GPBoost library
Now that we have prepared the data, we can
finally start to run models to see actual
results. To illustrate, we are using the GPBoost
library, which supports statistical methods
such as tree-boosting, Gaussian processes, and
mixed-effects models. The real attraction of
this library is that it allows you either to use
methods separately or to combine multiple
methods at the same time, which means features
can be both categorical and continuous, and
gaussian processes can be integrated with linear
models.
Identifying the variables of interest
In the code below, I show how this works step
by step. First we import the necessary libraries
that we need. We then identify the independent
(GB133_GB065) and dependent
(GB193) variables in our data. We
can then instantiate various models to see
whether, if a language is SOV and the Genitive
follows the noun (GB133_GB065),
that language’s adjectives also follow
the noun (GB193). To make it a bit
easier, we write a function that takes a
gp_model and our dataframe along
with the dv/iv and prints the result.
import numpy as np
import gpboost as gpb
from scipy import stats
iv = 'GB133_GB065' # independent variable (GB133_GB065: SOV order + [noun > genitive])
dv = 'GB193' # dependent variable (GB193: noun > adj)
def get_results(gp_model, df=df, dv=dv, iv=iv):
y = df[dv].to_numpy() # here we set our dependent variable
X = df[iv] # here is our independent variable
# add an intercept as a baseline for modeling the independent variable
X_with_intercept = np.column_stack((np.ones(len(X)), X))
# fit the Linear Mixed Model
gp_model.fit(y=y, X=X_with_intercept)
print(gp_model.summary())Instantiating the base model
In the code below we instantiate a base linear model to observe the potential relationship between our iv and dv.
# instantiate the vanilla gp_model with arguments
gp_model = gpb.GPModel(
num_data=len(df), # use number of observations to get an i.i.d baseline
likelihood="bernoulli_logit", # since our dv/iv are binary (yes/no), we use logit likelihood
)
get_results(gp_model) # run the model and print the results
# =====================================================
# Model summary:
# Nb. observations: 1729
# Log-lik AIC BIC
# -1139.64 2283.28 2294.19
# -----------------------------------------------------
# -----------------------------------------------------
# Linear regression coefficients (fixed effects):
# Param. Std. err. z value P(>|z|)
# Covariate_1 0.6353 0.0605 10.5013 0.0
# Covariate_2 -0.7504 0.1065 -7.0438 0.0
# =====================================================The summary table that prints shows the
number of observations in the dataset along with
additional information about model fit. What we
are mainly interested in is the linear
regression info. Because our independent
variable requires a baseline intercept for
regression, this is calculated as
Covariate_1, suggesting that the
true baseline is confidently different from zero
(i.e. that for languages without the
SOV + noun>genitive profile,
there is a tendency for adjectives to follow the
noun). Our dependent variable is
Covariate_2, and here it indicates
that having an adjective following a noun is
negatively correlated with
SOV + noun>genitive. Both these
tendencies are significant.
The Macroarea model
However, this is just the baseline - we want
to see whether controlling for family or area
affects these results, so in the following code
we examine each separately and then together.
First we instantiate a model to determine fit
with the macroarea group, then with
the Family_ID group, and then with
both together.
# instantiate a model grouped by macroarea
gp_model = gpb.GPModel(
group_data=df['macroarea'].astype(str).to_numpy(),
likelihood="bernoulli_logit",
)
get_results(gp_model) # run the model and print the results
# Model summary:
# Nb. observations: 1729
# Nb. groups: 6 (Group_1)
# Log-lik AIC BIC
# -1081.68 2169.37 2185.73
# -----------------------------------------------------
# Covariance parameters (random effects):
# Param. Std. err.
# Group_1 0.3003 0.1831
# -----------------------------------------------------
# Linear regression coefficients (fixed effects):
# Param. Std. err. z value P(>|z|)
# Covariate_1 0.3794 0.2364 1.6049 0.1085
# Covariate_2 -0.3424 0.1189 -2.8801 0.0040
# =====================================================This table indicates that adding
macroarea as a grouping factor
greatly diminishes the likelihood for adjectives
to follow the noun in languages without the
SOV + noun>genitive profile. But
for languages with the
SOV + noun>genitive profile,
this is negatively correlated with adjectives
following a noun, and the tendency is still
significant.
The Language family model
Now we can try grouping by language family, as in the following code:
# instantiate a model grouped by language family
gp_model = gpb.GPModel(
group_data=df['Family_ID'].astype(str).to_numpy(),
likelihood="bernoulli_logit",
)
get_results(gp_model) # run the model and print the results
# =====================================================
# Model summary:
# Nb. observations: 1729
# Nb. groups: 175 (Group_1)
# Log-lik AIC BIC
# -957.75 1921.49 1937.86
# -----------------------------------------------------
# Covariance parameters (random effects):
# Param. Std. err.
# Group_1 6.3624 1.8441
# -----------------------------------------------------
# Linear regression coefficients (fixed effects):
# Param. Std. err. z value P(>|z|)
# Covariate_1 0.1541 0.2701 0.5704 0.5684
# Covariate_2 -0.2577 0.1920 -1.3426 0.1794
# =====================================================Here we see that the effect completely disappears when grouping by language family.
The Macroarea + Language family model
Similarly, the effect is also non-significant when adding language family as a second random effect:
# instantiate a model grouped by macroarea and language family
gp_model = gpb.GPModel(
group_data=df[['macroarea', 'Family_ID']].astype(str).to_numpy(),
likelihood="bernoulli_logit",
)
get_results(gp_model) # run the model and print the results
# =====================================================
# Model summary:
# Nb. observations: 1729
# Nb. groups: 6 (Group_1), 175 (Group_2)
# Log-lik AIC BIC
# -957.75 1923.51 1945.33
# -----------------------------------------------------
# Covariance parameters (random effects):
# Param. Std. err.
# Group_1 0.0002 0.0042
# Group_2 6.4533 1.8888
# -----------------------------------------------------
# Linear regression coefficients (fixed effects):
# Param. Std. err. z value P(>|z|)
# Covariate_1 0.1409 0.2730 0.5162 0.6057
# Covariate_2 -0.2606 0.1932 -1.3488 0.1774
# =====================================================Using longitude and latitude in a gaussian process
Now let’s see how we can integrate a gaussian
process term using longitude and latitude. Here
we can simply instantiate the
gp_model with a
gp_coords argument, where we pass
the column names in our dataset containing the
actual coordinates for each datapoint.
gp_model = gpb.GPModel(
gp_coords=df[['longitude', 'latitude']].to_numpy(),
likelihood="bernoulli_logit",
)
get_results(gp_model)
# =====================================================
# Model summary:
# Nb. observations: 1729
# Log-lik AIC BIC
# -876.35 1760.69 1782.51
# -----------------------------------------------------
# Covariance parameters (random effects):
# Param. Std. err.
# GP_var 8.8535 2.6895
# GP_range 20.8785 6.6060
# -----------------------------------------------------
# Linear regression coefficients (fixed effects):
# Param. Std. err. z value P(>|z|)
# Covariate_1 -0.4728 0.7141 -0.6621 0.5079
# Covariate_2 -0.2446 0.2300 -1.0638 0.2874
# =====================================================Combining grouped random effects with the gaussian process
These grouped random effects and the gaussian process coordinates can be combined, as in the following model:
gp_model = gpb.GPModel(
group_data=df[['macroarea', 'Family_ID']].astype(str).to_numpy(),
gp_coords=df[['longitude', 'latitude']].to_numpy(),
likelihood="bernoulli_logit",
)
get_results(gp_model)
# =====================================================
# Model summary:
# Nb. observations: 1729
# Nb. groups: 6 (Group_1), 175 (Group_2)
# Log-lik AIC BIC
# -852.65 1717.29 1750.02
# -----------------------------------------------------
# Covariance parameters (random effects):
# Param. Std. err.
# Group_1 0.5472 0.5459
# Group_2 3.8430 1.4460
# GP_var 2.9193 0.7816
# GP_range 7.9437 1.7111
# -----------------------------------------------------
# Linear regression coefficients (fixed effects):
# Param. Std. err. z value P(>|z|)
# Covariate_1 -0.2748 0.4804 -0.5720 0.5673
# Covariate_2 -0.1067 0.2550 -0.4182 0.6758
# =====================================================Discussion
Here we see that the purported “universal”
relationship between a language with
SOV + noun>genitive profile
(GB133_GB065) and
noun > adj order
(GB193) is unsupported by most
models, and even with a vanilla linear model the
relationship is negative (though significant).
The lack of correlation is evident when we are
able to control for various spatial and temporal
factors.
The benefit of controls in any investigation is that it facilitates the elimination of possibilities. This is crucial when trying to understand how features interact in a system. For typologists that system is language, both at small and large scales. GPS coordinates allow for integration of a more robust set of controls than has been possible with previous approaches, and the fact that we can do this relatively easily means that investigation can proceed apace. As always, feel free to reach out with any comments or questions!