Qst vs Fst: A Power Analysis
Setup
library(pacman); p_load(ggplot2)
my_theme <- function(){
theme(
legend.background = element_rect(colour = "black", linewidth = 0.5),
legend.key.width = unit(0.5, "cm"),
legend.key.height = unit(0.4, "cm"),
panel.background = element_rect(fill = "white"),
plot.background = element_rect(fill = "white"),
panel.grid.major.x = element_blank(),
panel.grid.minor.x = element_blank(),
panel.grid.major.y = element_line(colour ="gray80",linewidth = 1,linetype="dashed"),
panel.grid.minor.y = element_blank(),
axis.text.x = element_text(colour ="black",size=rel(1.5),hjust=0.5),
axis.text.y = element_text(colour ="black",size=rel(1.5),vjust=0.5),
axis.title.x = element_text(colour ="black",size=rel(1.2)),
axis.title.y = element_text(colour ="black",size=rel(1.2)),
plot.margin=unit(c(1,0.5,0.5,0.5), "cm"),
plot.title=element_text(hjust=0.5,vjust=5,size=rel(1.5),face="bold"),
plot.caption=element_text(hjust=0,size=rel(1.2)),
axis.line.x=element_line(color="gray80",linewidth = 1,linetype="dashed"),
axis.line.y=element_line(color="gray80",linewidth = 1,linetype="dashed"),
axis.ticks.x = element_blank(),
axis.ticks.y = element_blank())}
# Function to calculate critical Qst/Fst values
calculate_critical_values <- function(n_d, alpha) {
# Critical value for the upper tail (Qst > Fst)
chi_square_upper = qchisq(1 - alpha / 2, df = n_d - 1)
delta_upper = chi_square_upper / (n_d - 1)
# Critical value for the lower tail (Qst < Fst)
chi_square_lower = qchisq(alpha / 2, df = n_d - 1)
delta_lower = chi_square_lower / (n_d - 1)
return(c(upper = delta_upper, lower = delta_lower))}
# Function to simulate power
simulate_power <- function(n_d, Fst_bar, Qst_Fst_ratio, num_simulations = 10000, alpha = 0.05) {
critical_values <- calculate_critical_values(n_d, alpha)
# Simulate Qst values under the alternative hypothesis
simulated_Qst = rchisq(num_simulations, df = n_d - 1) * (Fst_bar / (n_d - 1)) * Qst_Fst_ratio
# Calculate power
if (Qst_Fst_ratio > 1) {
# For Qst > Fst, test if significantly higher
power = mean(simulated_Qst > critical_values['upper'] * Fst_bar)
} else {
# For Qst < Fst, test if significantly lower
power = mean(simulated_Qst < critical_values['lower'] * Fst_bar)
}
return(power)}Rationale
Walsh & Lynch (2018, pp. 455-457) described how tests of whether \(Q_{st}\) exceeds or is smaller than \(F_{st}\) have very low power without large numbers of demes or extreme levels of differentiation. Here’s an illustration in the vein of Whitlock & Guillaume (2009), using the parameters mentioned by Walsh & Lynch.
Analysis
set.seed(123)
# Set Qst/Fst ratios to simulate
Qst_Fst_ratios = c(1.15, 1.30, 1.50,
1/1.15, 1/1.30, 1/1.50)
# Set parameters
n_demes_min = 2
n_demes_max = 50
Fst_bar = 0.1 # Average Fst
# Initialize data frame for results
power_results <- expand.grid(n_d = n_demes_min:n_demes_max, Qst_Fst_ratio = Qst_Fst_ratios)
power_results$power <- NA
power_results$critical_value_upper <- NA
power_results$critical_value_lower <- NA
for (i in 1:nrow(power_results)) {
n_d = power_results$n_d[i]
critical_values <- calculate_critical_values(n_d, 0.05)
power_results$power[i] <- simulate_power(n_d, Fst_bar, power_results$Qst_Fst_ratio[i])
power_results$critical_value_upper[i] <- critical_values['upper']
power_results$critical_value_lower[i] <- critical_values['lower']}
# Rename columns for clarity
names(power_results) <- c("Number of Demes", "Qst/Fst Ratio", "Power", "Critical Value (Qst > Fst)", "Critical Value (Qst < Fst)")
power_resultsggplot(power_results, aes(x = `Number of Demes`, y = Power, color = as.factor(`Qst/Fst Ratio`))) +
geom_line(linewidth = 1.5) +
labs(
title = "Power to Detect Selection Over Different Numbers of Demes",
x = "Number of Demes",
y = "Power",
color = "Qst/Fst Ratio") +
my_theme() +
scale_color_manual(values = c("orangered", "#1874cd", "#67238a", "#B33F62", "#F3C677", "darkblue"),
labels = function(x) round(as.numeric(x), 3)) +
theme(
legend.position = c(.10, .5)) +
ylim(0, .8)
Discussion
When it comes to detecting directional or stabilizing selection, the demands of this test can be extreme. In realistic scenarios involving human populations, it is unlikely that this test will be sufficiently powered.
References
Walsh, B., & Lynch, M. (2018). Evolution and Selection of Quantitative Traits. Oxford University Press. https://doi.org/10.1093/oso/9780198830870.001.0001
Whitlock, M. C., & Guillaume, F. (2009). Testing for Spatially Divergent Selection: Comparing QST to FST. Genetics, 183(3), 1055–1063. https://doi.org/10.1534/genetics.108.099812
sessionInfo()## R version 4.3.1 (2023-06-16 ucrt)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 19045)
##
## Matrix products: default
##
##
## locale:
## [1] LC_COLLATE=English_United States.utf8
## [2] LC_CTYPE=English_United States.utf8
## [3] LC_MONETARY=English_United States.utf8
## [4] LC_NUMERIC=C
## [5] LC_TIME=English_United States.utf8
##
## time zone: America/Chicago
## tzcode source: internal
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] ggplot2_3.4.4 pacman_0.5.1
##
## loaded via a namespace (and not attached):
## [1] vctrs_0.6.3 cli_3.6.1 knitr_1.45 rlang_1.1.1
## [5] xfun_0.39 highr_0.10 generics_0.1.3 jsonlite_1.8.7
## [9] labeling_0.4.3 glue_1.6.2 colorspace_2.1-0 htmltools_0.5.7
## [13] sass_0.4.7 fansi_1.0.4 scales_1.2.1 rmarkdown_2.25
## [17] grid_4.3.1 evaluate_0.23 munsell_0.5.0 jquerylib_0.1.4
## [21] tibble_3.2.1 fastmap_1.1.1 yaml_2.3.7 lifecycle_1.0.4
## [25] compiler_4.3.1 dplyr_1.1.2 pkgconfig_2.0.3 rstudioapi_0.15.0
## [29] farver_2.1.1 digest_0.6.33 R6_2.5.1 tidyselect_1.2.0
## [33] utf8_1.2.3 pillar_1.9.0 magrittr_2.0.3 bslib_0.6.0
## [37] withr_2.5.2 tools_4.3.1 gtable_0.3.4 cachem_1.0.8