── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr 1.1.4 ✔ readr 2.1.5
✔ forcats 1.0.0 ✔ stringr 1.5.1
✔ ggplot2 3.5.2 ✔ tibble 3.2.1
✔ lubridate 1.9.4 ✔ tidyr 1.3.1
✔ purrr 1.0.4
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag() masks stats::lag()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
Loading required package: Rcpp
Loading 'brms' package (version 2.22.0). Useful instructions
can be found by typing help('brms'). A more detailed introduction
to the package is available through vignette('brms_overview').
Attaching package: ‘brms’
The following object is masked from ‘package:stats’:
ar
Loading required package: viridisLite
This is loo version 2.8.0
- Online documentation and vignettes at mc-stan.org/loo
- As of v2.0.0 loo defaults to 1 core but we recommend using as many as possible. Use the 'cores' argument or set options(mc.cores = NUM_CORES) for an entire session.
This is bayesplot version 1.11.1
- Online documentation and vignettes at mc-stan.org/bayesplot
- bayesplot theme set to bayesplot::theme_default()
* Does _not_ affect other ggplot2 plots
* See ?bayesplot_theme_set for details on theme setting
Attaching package: ‘bayesplot’
The following object is masked from ‘package:brms’:
rhat
Attaching package: ‘ape’
The following object is masked from ‘package:dplyr’:
where
Loading required package: maps
Attaching package: ‘maps’
The following object is masked from ‘package:viridis’:
unemp
The following object is masked from ‘package:purrr’:
map
Attaching package: ‘Matrix’
The following objects are masked from ‘package:tidyr’:
expand, pack, unpack
This is posterior version 1.6.1
Attaching package: ‘posterior’
The following object is masked from ‘package:bayesplot’:
rhat
The following objects are masked from ‘package:stats’:
mad, sd, var
The following objects are masked from ‘package:base’:
%in%, match
Loading required package: StanHeaders
rstan version 2.32.6 (Stan version 2.32.2)
For execution on a local, multicore CPU with excess RAM we recommend calling
options(mc.cores = parallel::detectCores()).
To avoid recompilation of unchanged Stan programs, we recommend calling
rstan_options(auto_write = TRUE)
For within-chain threading using `reduce_sum()` or `map_rect()` Stan functions,
change `threads_per_chain` option:
rstan_options(threads_per_chain = 1)
Attaching package: ‘rstan’
The following objects are masked from ‘package:posterior’:
ess_bulk, ess_tail
The following object is masked from ‘package:tidyr’:
extract
Attaching package: ‘bridgesampling’
The following object is masked from ‘package:brms’:
bf
Loading required package: nlme
Attaching package: ‘nlme’
The following object is masked from ‘package:dplyr’:
collapse
This is mgcv 1.9-3. For overview type 'help("mgcv-package")'.
Attaching package: ‘mgcv’
The following objects are masked from ‘package:brms’:
s, t2
Attaching package: ‘tidybayes’
The following objects are masked from ‘package:brms’:
dstudent_t, pstudent_t, qstudent_t, rstudent_t
ggtree v3.14.0 Learn more at https://yulab-smu.top/contribution-tree-data/
Please cite:
Guangchuang Yu. Data Integration, Manipulation and Visualization of
Phylogenetic Trees (1st edition). Chapman and Hall/CRC. 2022,
doi:10.1201/9781003279242, ISBN: 9781032233574
Attaching package: ‘ggtree’
The following object is masked from ‘package:nlme’:
collapse
The following object is masked from ‘package:Matrix’:
expand
The following object is masked from ‘package:ape’:
rotate
The following object is masked from ‘package:tidyr’:
expand
Attaching package: ‘cowplot’
The following object is masked from ‘package:patchwork’:
align_plots
The following object is masked from ‘package:ggthemes’:
theme_map
The following object is masked from ‘package:lubridate’:
stamp
Load WALS data
Filter and process data
d =read_csv("../data/affix_adposition.csv")d %>%head()
Rows: 488 Columns: 8
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (4): Glottocode, Name, Macroarea, Family
dbl (4): Longitude, Latitude, Affix, Adposition
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
A tibble: 6 × 8
Glottocode
Name
Macroarea
Longitude
Latitude
Family
Affix
Adposition
<chr>
<chr>
<chr>
<dbl>
<dbl>
<chr>
<dbl>
<dbl>
aari1239
Aari
Africa
36.583333
6.000000
Afro-Asiatic
2
1
abkh1244
Abkhaz
Eurasia
41.000000
43.083333
Northwest Caucasian
4
1
acha1250
Achagua
South America
-72.250000
4.416667
Arawakan
4
1
acol1236
Acholi
Africa
32.666667
3.000000
Eastern Sudanic
4
2
west2632
Acoma
North America
-107.583333
34.916667
Keresan
4
1
adio1239
Adioukrou
Africa
-4.583333
5.416667
Niger-Congo
5
1
load tree and scale its height
tree <-read.tree("../data/affix_adposition.tre")tree$edge.length <- tree$edge.length /median(tree$edge.length)
edge_matrix <- tree$edgeinternal_nodes <-intersect(edge_matrix[, 1], edge_matrix[, 2])# Map: for each node, at what row does it first appear as a child?child_indices <-match(internal_nodes, edge_matrix[,2])# Map: for each node, at what row does it first appear as a parent?parent_indices <-match(internal_nodes, edge_matrix[,1])all(child_indices < parent_indices)
TRUE
arrange data in order of tree tip.labels
d <- d[match(tree$tip.label, d$Glottocode), ]
create numerical features for affix and adposition types, and family affiliation
create tibble with latent values for interior nodes
# Add node labels to the tree based on node numberingtree$node.label <-as.character((length(tree$tip.label) +1):(length(tree$tip.label) + tree$Nnode))# 1. Subset auf Indo-Europeand_ie <- d %>%filter(Family =="Indo-European")glottocodes_ie <- d_ie$Glottocode# 2. Baum prunenie_tree <-drop.tip(tree, setdiff(tree$tip.label, glottocodes_ie), keep.node.label =TRUE)
ie_tree_plot <-ggtree(ie_tree)ie_tree_info <- ie_tree_plot$data# Add column node.labelie_tree_info$node.label <-as.character(ie_tree_info$label)for (i in1:nrow(ie_tree_info)) {if (ie_tree_info$isTip[i]) {# If it's a tip, get the corresponding Glottocode glottocode <- ie_tree_info$label[i]# Find the index of the Glottocode in d_ie index <-which(tree$tip.label == glottocode)# Assign the Affix value to the node label ie_tree_info$node.label[i] <- index } else {# If it's not a tip, copy label ie_tree_info$node.label[i] <-as.integer(ie_tree_info$node.label[i]) }}# trait vectors with original node ID (as in the full tree)node_data <-tibble(node.label =as.character(1:length(z_affix)),z_affix = z_affix + ie_intercepts[1],z_adposition = z_adposition + ie_intercepts[2])# merge with ie_tree_infoie_tree_data <- ie_tree_info %>%left_join(node_data, by ="node.label")ie_tree_data <- ie_tree_data %>%left_join(rename(select(d, Glottocode, affix_ord, adp_bin),label=Glottocode) ) %>%mutate(Affix =factor(affix_ord, levels =names(affix_levels), labels = affix_levels),Adposition =factor(adp_bin, levels =names(adposition_levels), labels = adposition_levels) )
Estimated log Bayes factor (based on medians of log marginal likelihood estimates)
in favor of x1 over x2: -1310.27602
Range of estimates: -1316.56516 to -1305.51547
Interquartile range: 5.96180
Estimated log Bayes factor (based on medians of log marginal likelihood estimates)
in favor of x1 over x2: -925.51298
Range of estimates: -931.28727 to -916.06188
Interquartile range: 6.42169
