Description
Predict and Visualize Population-Level Changes in Allele Frequencies in Response to Climate Change.
Description
Methods (<doi:10.7717/peerj.11534>) are provided of calibrating and predicting shifts in allele frequencies through redundancy analysis ('vegan::rda()') and generalized additive models ('mgcv::gam()'). Visualization functions for predicted changes in allele frequencies include 'shift.dot.ggplot()', 'shift.pie.ggplot()', 'shift.moon.ggplot()', 'shift.waffle.ggplot()' and 'shift.surf.ggplot()' that are made with input data sets that are prepared by helper functions for each visualization method. Examples in the documentation show how to prepare animated climate change graphics through a time series with the 'gganimate' package. Function 'amova.rda()' shows how Analysis of Molecular Variance can be directly conducted with the results from redundancy analysis.
README.md
AlleleShift
AlleleShift is an R package to predict and visualize population-level changes in allele frequencies.
The following manuscript provides details on the methods used by the package: Kindt R. 2021 (in prep.) AlleleShift: An R package to predict and visualize population-level changes in allele frequencies in response to climate change.
Installation
To install from CRAN, use:
install.packages('AlleleShift')
The package can be installed from github via:
# install.packages('devtools')
devtools::install_github('RoelandKindt/AlleleShift')
Packages needed
library(BiodiversityR) # also loads vegan
library(poppr) # also loads adegenet
library(ggplot2)
library(ggsci)
library(ggforce)
library(dplyr)
library(ggrepel)
library(patchwork)
library(GGally)
library(mgcv)
library(ggmap)
library(gggibbous)
library(gganimate)
library(AlleleShift)
Predicting shifts in allele frequencies
The following script is also available from the documentation of the ‘count.model’ function.
Input data ‘Poptri.baseline.env’ and ‘Poptri.future.env’ document bioclimatic conditions for populations in baseline and future climates. ‘Poptri.genind’ is information on allele counts of individuals in the ‘adegenet::genind’ format (note that packages ‘adegenet’ and ‘poppr’ offer various methods of importing data from other applications).
# 1. Reduce the number of explanatory variables
data(Poptri.baseline.env)
data(Poptri.future.env)
VIF.select <- VIF.subset(Poptri.baseline.env,
keep=c("MAT", "CMI"),
cor.plot=TRUE)
#> Step 1: Keeping these vars:
#> [1] "MAT" "CMI" "TD" "MSP" "AHM" "SHM" "DD.0" "bFFP"
#> [9] "PAS" "Eref" "CMD" "cmiJJA" "PPT_sm" "Tmax07"
#>
#> Variance inflation (package: car)
#> MAT TD MSP AHM SHM DD.0 bFFP
#> 1956.71402 55.35151 325.04715 258.27807 199.03579 383.71475 208.48116
#> PAS Eref CMD CMI cmiJJA PPT_sm Tmax07
#> 418.39278 509.80399 702.97284 729.58048 2108.47844 1221.60927 404.91923
#>
#> VIF directly calculated from linear model with focal numeric variable as response
#> cmiJJA MAT PPT_sm CMI CMD Eref PAS
#> 2108.47844 1956.71402 1221.60927 729.58048 702.97284 509.80399 418.39278
#> Tmax07 DD.0 MSP AHM bFFP SHM TD
#> 404.91923 383.71475 325.04715 258.27807 208.48116 199.03579 55.35151
#>
#> Variance inflation (package: car)
#> MAT TD MSP AHM SHM DD.0 bFFP
#> 1152.51188 51.00442 280.39905 196.06676 189.87780 289.82940 177.01912
#> PAS Eref CMD CMI PPT_sm Tmax07
#> 373.46930 302.59709 652.89020 548.59673 527.99218 186.61090
#>
#> VIF directly calculated from linear model with focal numeric variable as response
#> MAT CMD CMI PPT_sm PAS Eref DD.0
#> 1152.51188 652.89020 548.59673 527.99218 373.46930 302.59709 289.82940
#> MSP AHM SHM Tmax07 bFFP TD
#> 280.39905 196.06676 189.87780 186.61090 177.01912 51.00442
#>
#> Variance inflation (package: car)
#> MAT TD MSP AHM SHM DD.0 bFFP PAS
#> 775.21530 48.40917 275.40371 179.50626 90.11375 260.05137 146.29625 368.40367
#> Eref CMI PPT_sm Tmax07
#> 102.83305 545.11114 394.89879 180.19703
#>
#> VIF directly calculated from linear model with focal numeric variable as response
#> MAT CMI PPT_sm PAS MSP DD.0 Tmax07 AHM
#> 775.21530 545.11114 394.89879 368.40367 275.40371 260.05137 180.19703 179.50626
#> bFFP Eref SHM TD
#> 146.29625 102.83305 90.11375 48.40917
#>
#> Variance inflation (package: car)
#> MAT TD MSP AHM SHM DD.0 bFFP PAS
#> 528.61146 32.01349 145.93063 78.09567 53.38943 222.17330 145.00874 188.09374
#> Eref CMI Tmax07
#> 102.40225 306.50400 111.10074
#>
#> VIF directly calculated from linear model with focal numeric variable as response
#> MAT CMI DD.0 PAS MSP bFFP Tmax07 Eref
#> 528.61146 306.50400 222.17330 188.09374 145.93063 145.00874 111.10074 102.40225
#> AHM SHM TD
#> 78.09567 53.38943 32.01349
#>
#> Variance inflation (package: car)
#> MAT TD MSP AHM SHM bFFP PAS Eref
#> 291.46937 30.38197 64.59840 46.78232 37.94467 45.37584 45.57930 99.66232
#> CMI Tmax07
#> 93.58775 109.62229
#>
#> VIF directly calculated from linear model with focal numeric variable as response
#> MAT Tmax07 Eref CMI MSP AHM PAS bFFP
#> 291.46937 109.62229 99.66232 93.58775 64.59840 46.78232 45.57930 45.37584
#> SHM TD
#> 37.94467 30.38197
#>
#> Variance inflation (package: car)
#> MAT TD MSP AHM SHM bFFP PAS
#> 237.300281 4.513064 62.232526 46.486704 37.871152 44.616210 44.930713
#> Eref CMI
#> 33.411744 85.667431
#>
#> VIF directly calculated from linear model with focal numeric variable as response
#> MAT CMI MSP AHM PAS bFFP SHM
#> 237.300281 85.667431 62.232526 46.486704 44.930713 44.616210 37.871152
#> Eref TD
#> 33.411744 4.513064
#>
#> Variance inflation (package: car)
#> MAT TD AHM SHM bFFP PAS Eref
#> 218.240423 4.170754 25.872364 5.680412 43.733182 39.989330 32.625175
#> CMI
#> 27.162467
#>
#> VIF directly calculated from linear model with focal numeric variable as response
#> MAT bFFP PAS Eref CMI AHM SHM
#> 218.240423 43.733182 39.989330 32.625175 27.162467 25.872364 5.680412
#> TD
#> 4.170754
#>
#> Variance inflation (package: car)
#> MAT TD AHM SHM PAS Eref CMI
#> 56.168045 4.154401 24.524555 5.474406 31.877700 13.105090 22.368212
#>
#> VIF directly calculated from linear model with focal numeric variable as response
#> MAT PAS AHM CMI Eref SHM TD
#> 56.168045 31.877700 24.524555 22.368212 13.105090 5.474406 4.154401
#>
#> Variance inflation (package: car)
#> MAT TD AHM SHM Eref CMI
#> 8.713339 4.003336 19.299959 5.465614 9.025702 14.431931
#>
#> VIF directly calculated from linear model with focal numeric variable as response
#> AHM CMI Eref MAT SHM TD
#> 19.299959 14.431931 9.025702 8.713339 5.465614 4.003336
#>
#> Summary of VIF selection process:
#> MAT TD MSP AHM SHM DD.0 bFFP
#> step_1 1956.714025 55.351511 325.04715 258.27807 199.035789 383.7148 208.48116
#> step_2 1152.511883 51.004421 280.39905 196.06676 189.877801 289.8294 177.01912
#> step_3 775.215296 48.409171 275.40371 179.50626 90.113750 260.0514 146.29625
#> step_4 528.611459 32.013493 145.93063 78.09567 53.389431 222.1733 145.00874
#> step_5 291.469373 30.381965 64.59840 46.78232 37.944668 NA 45.37584
#> step_6 237.300281 4.513064 62.23253 46.48670 37.871152 NA 44.61621
#> step_7 218.240423 4.170754 NA 25.87236 5.680412 NA 43.73318
#> step_8 56.168045 4.154401 NA 24.52456 5.474406 NA NA
#> step_9 8.713339 4.003336 NA 19.29996 5.465614 NA NA
#> PAS Eref CMD CMI cmiJJA PPT_sm Tmax07
#> step_1 418.39278 509.803994 702.9728 729.58048 2108.478 1221.6093 404.9192
#> step_2 373.46930 302.597094 652.8902 548.59673 NA 527.9922 186.6109
#> step_3 368.40367 102.833047 NA 545.11114 NA 394.8988 180.1970
#> step_4 188.09374 102.402249 NA 306.50400 NA NA 111.1007
#> step_5 45.57930 99.662321 NA 93.58775 NA NA 109.6223
#> step_6 44.93071 33.411744 NA 85.66743 NA NA NA
#> step_7 39.98933 32.625175 NA 27.16247 NA NA NA
#> step_8 31.87770 13.105090 NA 22.36821 NA NA NA
#> step_9 NA 9.025702 NA 14.43193 NA NA NA
#>
#> Final selection of variables:
#> [1] "AHM" "CMI" "Eref" "MAT" "SHM" "TD"
VIF.select$VIF$vars.included
#> [1] "AHM" "CMI" "Eref" "MAT" "SHM" "TD"
baseline.env <- Poptri.baseline.env[, VIF.select$VIF$vars.included]
summary(baseline.env)
#> AHM CMI Eref MAT
#> Min. : 5.50 Min. : 29.0 Min. :455.0 Min. : 1.690
#> 1st Qu.: 8.45 1st Qu.: 55.5 1st Qu.:602.0 1st Qu.: 5.620
#> Median :10.80 Median : 97.0 Median :675.0 Median : 8.980
#> Mean :12.06 Mean :109.5 Mean :680.9 Mean : 7.903
#> 3rd Qu.:16.80 3rd Qu.:150.5 3rd Qu.:773.5 3rd Qu.:10.045
#> Max. :19.50 Max. :232.0 Max. :920.0 Max. :11.290
#> SHM TD
#> Min. : 23.50 Min. :13.90
#> 1st Qu.: 38.45 1st Qu.:14.25
#> Median : 50.20 Median :15.30
#> Mean : 56.53 Mean :16.11
#> 3rd Qu.: 74.40 3rd Qu.:17.60
#> Max. :102.90 Max. :20.10
future.env <- Poptri.future.env[, VIF.select$VIF$vars.included]
# 2. Create the genpop object
data(Poptri.genind)
Poptri.genpop <- adegenet::genind2genpop(Poptri.genind)
#>
#> Converting data from a genind to a genpop object...
#>
#> ...done.
# Get to know the populations and the alleles
poppr::poppr(Poptri.genind)
#> Pop N MLG eMLG SE H G lambda E.5 Hexp Ia
#> 1 Chilliwack 38 14 5.67 1.244 1.972 3.97 0.748 0.480 0.168 0.29948
#> 2 Columbia 180 31 6.69 1.276 2.634 8.17 0.878 0.555 0.219 0.03214
#> 3 Dean 7 4 4.00 0.000 1.352 3.77 0.735 0.967 0.145 -0.16667
#> 4 Harrison 22 10 5.87 1.093 1.855 4.10 0.756 0.575 0.191 0.06387
#> 5 Homathko 9 7 7.00 0.000 1.831 5.40 0.815 0.840 0.182 -0.17703
#> 6 Kitimat 8 4 4.00 0.000 1.321 3.56 0.719 0.930 0.127 -0.13245
#> 7 Klinaklini 8 5 5.00 0.000 1.494 4.00 0.750 0.868 0.155 -0.41248
#> 8 Nisqually 14 7 5.51 0.830 1.567 3.38 0.704 0.627 0.157 0.00845
#> 9 Nooksack 26 11 6.44 1.072 2.124 6.76 0.852 0.782 0.211 0.05709
#> 10 Olympic_Penninsula 5 3 3.00 0.000 0.950 2.27 0.560 0.802 0.191 -0.25000
#> 11 Puyallup 186 29 6.03 1.303 2.396 5.81 0.828 0.482 0.178 -0.02222
#> 12 Salmon 11 7 6.64 0.481 1.846 5.76 0.826 0.892 0.193 -0.19978
#> 13 Skagit 163 20 5.41 1.268 2.102 4.72 0.788 0.518 0.168 -0.01194
#> 14 Skwawka 14 6 5.29 0.656 1.631 4.45 0.776 0.841 0.160 -0.04804
#> 15 Skykomish 156 31 6.97 1.249 2.730 9.26 0.892 0.576 0.233 0.07931
#> 16 Squamish 14 10 7.86 0.802 2.206 8.17 0.878 0.887 0.205 -0.28884
#> 17 Tahoe 12 9 7.82 0.649 2.095 7.20 0.861 0.870 0.293 -0.11481
#> 18 Vancouver I. East 12 4 3.50 0.584 0.837 1.71 0.417 0.546 0.107 0.48361
#> 19 Willamette 17 10 6.65 1.016 2.007 5.45 0.817 0.691 0.246 0.11902
#> 20 Total 902 59 6.38 1.323 2.663 6.93 0.856 0.445 0.200 0.03621
#> rbarD File
#> 1 0.08166 Poptri.genind
#> 2 0.00883 Poptri.genind
#> 3 -0.16667 Poptri.genind
#> 4 0.02067 Poptri.genind
#> 5 -0.04770 Poptri.genind
#> 6 -0.14222 Poptri.genind
#> 7 -0.21444 Poptri.genind
#> 8 0.00445 Poptri.genind
#> 9 0.01646 Poptri.genind
#> 10 -0.25000 Poptri.genind
#> 11 -0.00606 Poptri.genind
#> 12 -0.07663 Poptri.genind
#> 13 -0.00373 Poptri.genind
#> 14 -0.01849 Poptri.genind
#> 15 0.02249 Poptri.genind
#> 16 -0.08273 Poptri.genind
#> 17 -0.03060 Poptri.genind
#> 18 0.27081 Poptri.genind
#> 19 0.04139 Poptri.genind
#> 20 0.00998 Poptri.genind
adegenet::makefreq(Poptri.genpop)
#>
#> Finding allelic frequencies from a genpop object...
#>
#> ...done.
#> X01_10838495.1 X01_10838495.2 X01_16628872.1 X01_16628872.2
#> Chilliwack 0.10526316 0.8947368 0.19736842 0.8026316
#> Columbia 0.23055556 0.7694444 0.20000000 0.8000000
#> Dean 0.21428571 0.7857143 0.00000000 1.0000000
#> Harrison 0.15909091 0.8409091 0.25000000 0.7500000
#> Homathko 0.05555556 0.9444444 0.22222222 0.7777778
#> Kitimat 0.25000000 0.7500000 0.12500000 0.8750000
#> Klinaklini 0.06250000 0.9375000 0.43750000 0.5625000
#> Nisqually 0.10714286 0.8928571 0.25000000 0.7500000
#> Nooksack 0.17307692 0.8269231 0.26923077 0.7307692
#> Olympic_Penninsula 0.50000000 0.5000000 0.10000000 0.9000000
#> Puyallup 0.17204301 0.8279570 0.16935484 0.8306452
#> Salmon 0.04545455 0.9545455 0.27272727 0.7272727
#> Skagit 0.23619632 0.7638037 0.18098160 0.8190184
#> Skwawka 0.03571429 0.9642857 0.32142857 0.6785714
#> Skykomish 0.20192308 0.7980769 0.31410256 0.6858974
#> Squamish 0.28571429 0.7142857 0.25000000 0.7500000
#> Tahoe 0.62500000 0.3750000 0.08333333 0.9166667
#> Vancouver I. East 0.16666667 0.8333333 0.00000000 1.0000000
#> Willamette 0.35294118 0.6470588 0.20588235 0.7941176
#> X01_23799134.1 X01_23799134.2 X01_38900650.1 X01_38900650.2
#> Chilliwack 0.06578947 0.9342105 0.05263158 0.9473684
#> Columbia 0.10833333 0.8916667 0.06944444 0.9305556
#> Dean 0.00000000 1.0000000 0.21428571 0.7857143
#> Harrison 0.11363636 0.8863636 0.02272727 0.9772727
#> Homathko 0.05555556 0.9444444 0.05555556 0.9444444
#> Kitimat 0.00000000 1.0000000 0.00000000 1.0000000
#> Klinaklini 0.00000000 1.0000000 0.00000000 1.0000000
#> Nisqually 0.00000000 1.0000000 0.00000000 1.0000000
#> Nooksack 0.01923077 0.9807692 0.07692308 0.9230769
#> Olympic_Penninsula 0.00000000 1.0000000 0.00000000 1.0000000
#> Puyallup 0.04569892 0.9543011 0.05107527 0.9489247
#> Salmon 0.00000000 1.0000000 0.04545455 0.9545455
#> Skagit 0.03987730 0.9601227 0.01840491 0.9815951
#> Skwawka 0.07142857 0.9285714 0.00000000 1.0000000
#> Skykomish 0.11217949 0.8878205 0.03525641 0.9647436
#> Squamish 0.03571429 0.9642857 0.03571429 0.9642857
#> Tahoe 0.04166667 0.9583333 0.41666667 0.5833333
#> Vancouver I. East 0.00000000 1.0000000 0.08333333 0.9166667
#> Willamette 0.00000000 1.0000000 0.14705882 0.8529412
#> X01_41284594.1 X01_41284594.2
#> Chilliwack 0.05263158 0.9473684
#> Columbia 0.05000000 0.9500000
#> Dean 0.00000000 1.0000000
#> Harrison 0.02272727 0.9772727
#> Homathko 0.11111111 0.8888889
#> Kitimat 0.00000000 1.0000000
#> Klinaklini 0.06250000 0.9375000
#> Nisqually 0.10714286 0.8928571
#> Nooksack 0.09615385 0.9038462
#> Olympic_Penninsula 0.10000000 0.9000000
#> Puyallup 0.07526882 0.9247312
#> Salmon 0.22727273 0.7727273
#> Skagit 0.03374233 0.9662577
#> Skwawka 0.07142857 0.9285714
#> Skykomish 0.07692308 0.9230769
#> Squamish 0.03571429 0.9642857
#> Tahoe 0.12500000 0.8750000
#> Vancouver I. East 0.04166667 0.9583333
#> Willamette 0.08823529 0.9117647
# 3. Calibrate the models
# Note that the ordistep procedure is not needed
Poptri.count.model <- count.model(Poptri.genpop,
env.data=baseline.env,
ordistep=TRUE)
#>
#> Finding allelic frequencies from a genpop object...
#>
#> ...done.
#>
#>
#> Call:
#> rda(formula = gen.comm.b ~ AHM + CMI + Eref + MAT + SHM + TD, data = env.data)
#>
#> Partitioning of variance:
#> Inertia Proportion
#> Total 13675 1.0000
#> Constrained 6841 0.5002
#> Unconstrained 6834 0.4998
#>
#> Eigenvalues, and their contribution to the variance
#>
#> Importance of components:
#> RDA1 RDA2 RDA3 RDA4 RDA5
#> Eigenvalue 5203.4791 893.09032 489.16133 249.89783 5.3852777
#> Proportion Explained 0.3805 0.06531 0.03577 0.01827 0.0003938
#> Cumulative Proportion 0.3805 0.44580 0.48157 0.49985 0.5002410
#> PC1 PC2 PC3 PC4 PC5
#> Eigenvalue 4597.4331 1443.5788 458.75766 243.85111 90.80236
#> Proportion Explained 0.3362 0.1056 0.03355 0.01783 0.00664
#> Cumulative Proportion 0.8364 0.9420 0.97553 0.99336 1.00000
#>
#> Accumulated constrained eigenvalues
#> Importance of components:
#> RDA1 RDA2 RDA3 RDA4 RDA5
#> Eigenvalue 5203.4791 893.0903 489.1613 249.89783 5.3852777
#> Proportion Explained 0.7606 0.1305 0.0715 0.03653 0.0007872
#> Cumulative Proportion 0.7606 0.8912 0.9627 0.99921 1.0000000
#>
#> Scaling 2 for species and site scores
#> * Species are scaled proportional to eigenvalues
#> * Sites are unscaled: weighted dispersion equal on all dimensions
#> * General scaling constant of scores: 22.27427
#>
#>
#> Species scores
#>
#> RDA1 RDA2 RDA3 RDA4 RDA5 PC1
#> X01_10838495.1 7.7466 -0.7259 -1.69114 0.187490 -0.01266 6.28376
#> X01_10838495.2 -7.7466 0.7259 1.69114 -0.187490 0.01266 -6.28376
#> X01_16628872.1 -2.2584 2.5545 -1.43674 0.084325 -0.17373 -5.83911
#> X01_16628872.2 2.2584 -2.5545 1.43674 -0.084325 0.17373 5.83911
#> X01_23799134.1 -0.4671 1.9706 -0.83990 -0.001968 0.25745 0.03723
#> X01_23799134.2 0.4671 -1.9706 0.83990 0.001968 -0.25745 -0.03723
#> X01_38900650.1 5.2463 2.2731 1.80100 0.256469 -0.03033 2.62821
#> X01_38900650.2 -5.2463 -2.2731 -1.80100 -0.256469 0.03033 -2.62821
#> X01_41284594.1 1.2400 0.3130 0.01204 -2.103589 -0.01203 -1.70511
#> X01_41284594.2 -1.2400 -0.3130 -0.01204 2.103589 0.01203 1.70511
#>
#>
#> Site scores (weighted sums of species scores)
#>
#> RDA1 RDA2 RDA3 RDA4 RDA5 PC1
#> Chilliwack -3.6364 1.6841 5.3973 1.4271 39.264 2.5785
#> Columbia 0.4218 2.5336 -3.8742 4.1743 71.373 0.4335
#> Dean 4.6551 -6.6279 24.0513 14.2372 83.823 8.0787
#> Harrison -3.3140 4.2633 -5.5469 7.1447 53.583 -2.0779
#> Homathko -5.0101 4.0804 8.0351 -9.1976 12.167 -4.4760
#> Kitimat 0.2023 -11.0757 -2.1893 11.1417 24.158 6.4103
#> Klinaklini -7.9981 10.8293 -7.6871 -0.3698 -176.410 -12.6881
#> Nisqually -4.7481 -0.5922 0.4640 -8.6869 -56.562 -7.2177
#> Nooksack -1.4050 4.2429 -0.2436 -4.0265 -61.991 1.2854
#> Olympic_Penninsula 8.5460 -15.9908 -18.3575 -2.4474 24.067 4.1194
#> Puyallup -1.2328 -1.8999 2.9514 -1.6765 34.025 1.7945
#> Salmon -5.3069 5.0051 6.9613 -29.2596 -79.379 -6.6457
#> Skagit -0.2454 -4.5544 -4.5800 5.8729 23.210 5.9223
#> Skwawka -7.8707 7.7636 -1.3073 -3.1379 -28.712 -4.6465
#> Skykomish -2.0503 8.2869 -11.4201 -0.8444 3.752 0.6351
#> Squamish 1.0273 -0.6021 -10.7934 7.1369 -31.007 1.5425
#> Tahoe 21.1711 4.7608 3.7894 3.8075 20.644 4.3325
#> Vancouver I. East 0.6973 -12.3697 17.5610 3.5557 99.335 1.0511
#> Willamette 6.0970 0.2625 -3.2114 1.1486 -55.342 -0.4317
#>
#>
#> Site constraints (linear combinations of constraining variables)
#>
#> RDA1 RDA2 RDA3 RDA4 RDA5 PC1
#> Chilliwack -5.3411 2.94631 -5.73690 -1.2669 3.6318 2.5785
#> Columbia 1.0348 1.48239 -4.35834 -1.9144 -1.5967 0.4335
#> Dean -0.7121 3.05536 13.59960 5.7383 -2.3972 8.0787
#> Harrison -0.6115 4.57450 -6.67046 12.7199 1.7520 -2.0779
#> Homathko -0.2101 4.51571 2.95659 0.1637 3.7693 -4.4760
#> Kitimat -5.7357 -7.72450 0.08812 2.8441 6.3705 6.4103
#> Klinaklini 0.9827 -2.48137 3.56265 6.7122 -7.5810 -12.6881
#> Nisqually 1.1218 -6.23684 -1.27479 0.7219 -0.9320 -7.2177
#> Nooksack -3.6818 3.17251 -4.89865 -2.9837 -3.8484 1.2854
#> Olympic_Penninsula 5.1563 -8.68730 -3.34262 -4.1673 0.8657 4.1194
#> Puyallup -2.4480 -0.01102 -0.54619 0.1670 -4.9542 1.7945
#> Salmon -1.1061 -0.60968 4.52845 -8.3433 -0.9208 -6.6457
#> Skagit -5.9406 -1.76614 0.27739 1.5322 -6.8831 5.9223
#> Skwawka -4.4370 1.80687 6.44400 -4.1364 14.0686 -4.6465
#> Skykomish -1.9900 10.49790 -2.02179 -9.1465 -5.7696 0.6351
#> Squamish -1.1546 0.32081 1.41105 4.4112 0.9792 1.5425
#> Tahoe 16.7611 4.67546 3.53694 -0.3315 1.2278 4.3325
#> Vancouver I. East 2.2338 -10.24318 1.07117 -4.5355 -3.0133 1.0511
#> Willamette 6.0779 0.71222 -8.62622 1.8151 5.2313 -0.4317
#>
#>
#> Biplot scores for constraining variables
#>
#> RDA1 RDA2 RDA3 RDA4 RDA5 PC1
#> AHM 0.13353 -0.29604 -0.5719 -0.16042 -0.4093 0
#> CMI -0.31497 0.20319 0.5713 0.04335 0.5768 0
#> Eref 0.29012 0.03393 -0.7039 -0.35686 -0.1588 0
#> MAT -0.07412 -0.21901 -0.7120 -0.54165 -0.1552 0
#> SHM 0.63769 -0.32011 -0.5075 -0.10019 -0.1988 0
#> TD -0.18759 0.19007 0.3043 0.67274 0.4208 0
#>
#> Permutation test for rda under reduced model
#> Permutation: free
#> Number of permutations: 99
#>
#> Model: rda(formula = gen.comm.b ~ AHM + CMI + Eref + MAT + SHM + TD, data = env.data)
#> Df Variance F Pr(>F)
#> Model 6 6841.0 2.0019 0.08 .
#> Residual 12 6834.4
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> Permutation test for rda under reduced model
#> Terms added sequentially (first to last)
#> Permutation: free
#> Number of permutations: 99
#>
#> Model: rda(formula = gen.comm.b ~ AHM + CMI + Eref + MAT + SHM + TD, data = env.data)
#> Df Variance F Pr(>F)
#> AHM 1 338.4 0.5942 0.58
#> CMI 1 1803.1 3.1659 0.07 .
#> Eref 1 1054.0 1.8507 0.14
#> MAT 1 872.8 1.5325 0.20
#> SHM 1 2259.1 3.9666 0.03 *
#> TD 1 513.5 0.9017 0.35
#> Residual 12 6834.4
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> Permutation test for rda under reduced model
#> Marginal effects of terms
#> Permutation: free
#> Number of permutations: 99
#>
#> Model: rda(formula = gen.comm.b ~ AHM + CMI + Eref + MAT + SHM + TD, data = env.data)
#> Df Variance F Pr(>F)
#> AHM 1 680.0 1.1940 0.32
#> CMI 1 135.1 0.2373 0.86
#> Eref 1 762.5 1.3388 0.27
#> MAT 1 538.0 0.9447 0.36
#> SHM 1 2141.2 3.7595 0.04 *
#> TD 1 513.5 0.9017 0.23
#> Residual 12 6834.4
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> Start: gen.comm.b ~ 1
#>
#> Df AIC F Pr(>F)
#> + SHM 1 180.36 3.5025 0.040 *
#> + Eref 1 182.90 0.9356 0.315
#> + CMI 1 182.90 0.9379 0.405
#> + TD 1 183.39 0.4790 0.630
#> + MAT 1 183.36 0.5027 0.670
#> + AHM 1 183.44 0.4313 0.670
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> Step: gen.comm.b ~ SHM
#>
#> Df AIC F Pr(>F)
#> - SHM 1 181.92 3.5025 0.045 *
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> Df AIC F Pr(>F)
#> + MAT 1 178.52 3.5833 0.030 *
#> + AHM 1 177.28 4.9039 0.040 *
#> + Eref 1 180.08 2.0396 0.105
#> + CMI 1 180.18 1.9444 0.165
#> + TD 1 181.03 1.1544 0.330
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> Step: gen.comm.b ~ SHM + MAT
#>
#> Df AIC F Pr(>F)
#> - MAT 1 180.36 3.5833 0.035 *
#> - SHM 1 183.36 6.9396 0.010 **
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> Df AIC F Pr(>F)
#> + AHM 1 178.56 1.6263 0.215
#> + Eref 1 178.70 1.5047 0.220
#> + CMI 1 179.68 0.6724 0.540
#> + TD 1 179.79 0.5843 0.620
#>
#> Reduced model:
#> Permutation test for rda under reduced model
#> Permutation: free
#> Number of permutations: 99
#>
#> Model: rda(formula = gen.comm.b ~ SHM + MAT, data = env.data)
#> Df Variance F Pr(>F)
#> Model 2 4411.0 3.809 0.02 *
#> Residual 16 9264.4
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> Permutation test for rda under reduced model
#> Terms added sequentially (first to last)
#> Permutation: free
#> Number of permutations: 99
#>
#> Model: rda(formula = gen.comm.b ~ SHM + MAT, data = env.data)
#> Df Variance F Pr(>F)
#> SHM 1 2336.2 4.0348 0.02 *
#> MAT 1 2074.8 3.5833 0.05 *
#> Residual 16 9264.4
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> Permutation test for rda under reduced model
#> Marginal effects of terms
#> Permutation: free
#> Number of permutations: 99
#>
#> Model: rda(formula = gen.comm.b ~ SHM + MAT, data = env.data)
#> Df Variance F Pr(>F)
#> SHM 1 4018.2 6.9396 0.01 **
#> MAT 1 2074.8 3.5833 0.05 *
#> Residual 16 9264.4
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Poptri.pred.baseline <- count.pred(Poptri.count.model, env.data=baseline.env)
head(Poptri.pred.baseline)
#> Pop Pop.label Pop.index N Allele Allele.freq A B
#> P1 Chilliwack 1 1 76 X01_10838495.1 0.10526316 8 68
#> P2 Columbia 2 2 360 X01_10838495.1 0.23055556 83 277
#> P3 Dean 3 3 14 X01_10838495.1 0.21428571 3 11
#> P4 Harrison 4 4 44 X01_10838495.1 0.15909091 7 37
#> P5 Homathko 5 5 18 X01_10838495.1 0.05555556 1 17
#> P6 Kitimat 6 6 16 X01_10838495.1 0.25000000 4 12
#> Ap Bp N.e1 Freq.e1
#> P1 8.950722 67.04928 76 0.1177727
#> P2 88.915206 271.08479 360 0.2469867
#> P3 1.816999 12.18300 14 0.1297857
#> P4 9.875440 34.12456 44 0.2244418
#> P5 3.289572 14.71043 18 0.1827540
#> P6 1.694003 14.30600 16 0.1058752
Poptri.freq.model <- freq.model(Poptri.pred.baseline)
#>
#> Family: binomial
#> Link function: logit
#>
#> Formula:
#> cbind(A, B) ~ s(Freq.e1)
#>
#> Parametric coefficients:
#> Estimate Std. Error z value Pr(>|z|)
#> (Intercept) -2.20412 0.04298 -51.28 <2e-16 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> Approximate significance of smooth terms:
#> edf Ref.df Chi.sq p-value
#> s(Freq.e1) 7.368 8.094 437.4 <2e-16 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> R-sq.(adj) = 0.734 Deviance explained = 71.5%
#> UBRE = 1.2969 Scale est. = 1 n = 95
Poptri.freq.baseline <- freq.pred(Poptri.freq.model,
count.predicted=Poptri.pred.baseline)
head(Poptri.freq.baseline)
#> Pop Pop.label Pop.index N Allele Allele.freq A B
#> P1 Chilliwack 1 1 76 X01_10838495.1 0.10526316 8 68
#> P2 Columbia 2 2 360 X01_10838495.1 0.23055556 83 277
#> P3 Dean 3 3 14 X01_10838495.1 0.21428571 3 11
#> P4 Harrison 4 4 44 X01_10838495.1 0.15909091 7 37
#> P5 Homathko 5 5 18 X01_10838495.1 0.05555556 1 17
#> P6 Kitimat 6 6 16 X01_10838495.1 0.25000000 4 12
#> Ap Bp N.e1 Freq.e1 Freq.e2 LCL UCL increasing
#> P1 8.950722 67.04928 76 0.1177727 0.1144837 0.09803780 0.1309296 TRUE
#> P2 88.915206 271.08479 360 0.2469867 0.2079033 0.18915823 0.2266484 FALSE
#> P3 1.816999 12.18300 14 0.1297857 0.1290953 0.11110733 0.1470833 FALSE
#> P4 9.875440 34.12456 44 0.2244418 0.1898990 0.17229199 0.2075060 TRUE
#> P5 3.289572 14.71043 18 0.1827540 0.2008387 0.17374709 0.2279304 TRUE
#> P6 1.694003 14.30600 16 0.1058752 0.1050042 0.09154483 0.1184635 FALSE
# 4. Check how well the models predict baseline allele frequencies
# Populations are split in those with R2 > 0.50 and those with R2 < 0.50
# Better populations
plotA1 <- freq.ggplot(Poptri.freq.baseline,
plot.best=TRUE,
ylim=c(0.0, 0.8))
#> Pop Pop.label N GAM.rsq
#> 11 Puyallup 11 186 0.99
#> 17 Tahoe 17 12 0.99
#> 10 Olympic_Penninsula 10 5 0.98
#> 16 Squamish 16 14 0.98
#> 1 Chilliwack 1 38 0.95
#> 19 Willamette 19 17 0.91
#> 4 Harrison 4 22 0.90
#> 14 Skwawka 14 14 0.90
#> 15 Skykomish 15 156 0.89
#> 9 Nooksack 9 26 0.84
#> 2 Columbia 2 180 0.82
#> 18 Vancouver I. East 18 12 0.59
#> 13 Skagit 13 163 0.55
#> 8 Nisqually 8 14 0.54
#> 7 Klinaklini 7 8 0.42
#> 6 Kitimat 6 8 0.41
#> 12 Salmon 12 11 0.34
#> 5 Homathko 5 9 0.17
#> 3 Dean 3 7 0.15
#> [1] "Chilliwack" "Columbia" "Harrison"
#> [4] "Nisqually" "Nooksack" "Olympic_Penninsula"
#> [7] "Puyallup" "Skagit" "Skwawka"
#> [10] "Skykomish" "Squamish" "Tahoe"
#> [13] "Vancouver I. East" "Willamette"
plotA1
# Populations with low R2
manual.colour.values1 <- ggsci::pal_npg()(5)
plotB1 <- freq.ggplot(Poptri.freq.baseline,
plot.best=FALSE,
manual.colour.values=manual.colour.values1,
xlim=c(0, 0.5),
ylim=c(0, 0.25))
#> Pop Pop.label N GAM.rsq
#> 11 Puyallup 11 186 0.99
#> 17 Tahoe 17 12 0.99
#> 10 Olympic_Penninsula 10 5 0.98
#> 16 Squamish 16 14 0.98
#> 1 Chilliwack 1 38 0.95
#> 19 Willamette 19 17 0.91
#> 4 Harrison 4 22 0.90
#> 14 Skwawka 14 14 0.90
#> 15 Skykomish 15 156 0.89
#> 9 Nooksack 9 26 0.84
#> 2 Columbia 2 180 0.82
#> 18 Vancouver I. East 18 12 0.59
#> 13 Skagit 13 163 0.55
#> 8 Nisqually 8 14 0.54
#> 7 Klinaklini 7 8 0.42
#> 6 Kitimat 6 8 0.41
#> 12 Salmon 12 11 0.34
#> 5 Homathko 5 9 0.17
#> 3 Dean 3 7 0.15
#> [1] "Dean" "Homathko" "Kitimat" "Klinaklini" "Salmon"
plotB1
# Colouring by alleles
plotA2 <- freq.ggplot(Poptri.freq.baseline,
colour.Pop=FALSE,
plot.best=TRUE,
ylim=c(0.0, 0.8),
manual.colour.values=manual.colour.values1)
#> Pop Pop.label N GAM.rsq
#> 11 Puyallup 11 186 0.99
#> 17 Tahoe 17 12 0.99
#> 10 Olympic_Penninsula 10 5 0.98
#> 16 Squamish 16 14 0.98
#> 1 Chilliwack 1 38 0.95
#> 19 Willamette 19 17 0.91
#> 4 Harrison 4 22 0.90
#> 14 Skwawka 14 14 0.90
#> 15 Skykomish 15 156 0.89
#> 9 Nooksack 9 26 0.84
#> 2 Columbia 2 180 0.82
#> 18 Vancouver I. East 18 12 0.59
#> 13 Skagit 13 163 0.55
#> 8 Nisqually 8 14 0.54
#> 7 Klinaklini 7 8 0.42
#> 6 Kitimat 6 8 0.41
#> 12 Salmon 12 11 0.34
#> 5 Homathko 5 9 0.17
#> 3 Dean 3 7 0.15
#> [1] "Chilliwack" "Columbia" "Harrison"
#> [4] "Nisqually" "Nooksack" "Olympic_Penninsula"
#> [7] "Puyallup" "Skagit" "Skwawka"
#> [10] "Skykomish" "Squamish" "Tahoe"
#> [13] "Vancouver I. East" "Willamette"
plotA2
plotB2 <- freq.ggplot(Poptri.freq.baseline,
colour.Pop=FALSE,
plot.best=FALSE,
manual.colour.values=manual.colour.values1,
xlim=c(0, 0.5),
ylim=c(0, 0.25))
#> Pop Pop.label N GAM.rsq
#> 11 Puyallup 11 186 0.99
#> 17 Tahoe 17 12 0.99
#> 10 Olympic_Penninsula 10 5 0.98
#> 16 Squamish 16 14 0.98
#> 1 Chilliwack 1 38 0.95
#> 19 Willamette 19 17 0.91
#> 4 Harrison 4 22 0.90
#> 14 Skwawka 14 14 0.90
#> 15 Skykomish 15 156 0.89
#> 9 Nooksack 9 26 0.84
#> 2 Columbia 2 180 0.82
#> 18 Vancouver I. East 18 12 0.59
#> 13 Skagit 13 163 0.55
#> 8 Nisqually 8 14 0.54
#> 7 Klinaklini 7 8 0.42
#> 6 Kitimat 6 8 0.41
#> 12 Salmon 12 11 0.34
#> 5 Homathko 5 9 0.17
#> 3 Dean 3 7 0.15
#> [1] "Dean" "Homathko" "Kitimat" "Klinaklini" "Salmon"
plotB2
# Note that you can also compare data with:
poppr::poppr(Poptri.genind)
#> Pop N MLG eMLG SE H G lambda E.5 Hexp Ia
#> 1 Chilliwack 38 14 5.67 1.244 1.972 3.97 0.748 0.480 0.168 0.29948
#> 2 Columbia 180 31 6.69 1.276 2.634 8.17 0.878 0.555 0.219 0.03214
#> 3 Dean 7 4 4.00 0.000 1.352 3.77 0.735 0.967 0.145 -0.16667
#> 4 Harrison 22 10 5.87 1.093 1.855 4.10 0.756 0.575 0.191 0.06387
#> 5 Homathko 9 7 7.00 0.000 1.831 5.40 0.815 0.840 0.182 -0.17703
#> 6 Kitimat 8 4 4.00 0.000 1.321 3.56 0.719 0.930 0.127 -0.13245
#> 7 Klinaklini 8 5 5.00 0.000 1.494 4.00 0.750 0.868 0.155 -0.41248
#> 8 Nisqually 14 7 5.51 0.830 1.567 3.38 0.704 0.627 0.157 0.00845
#> 9 Nooksack 26 11 6.44 1.072 2.124 6.76 0.852 0.782 0.211 0.05709
#> 10 Olympic_Penninsula 5 3 3.00 0.000 0.950 2.27 0.560 0.802 0.191 -0.25000
#> 11 Puyallup 186 29 6.03 1.303 2.396 5.81 0.828 0.482 0.178 -0.02222
#> 12 Salmon 11 7 6.64 0.481 1.846 5.76 0.826 0.892 0.193 -0.19978
#> 13 Skagit 163 20 5.41 1.268 2.102 4.72 0.788 0.518 0.168 -0.01194
#> 14 Skwawka 14 6 5.29 0.656 1.631 4.45 0.776 0.841 0.160 -0.04804
#> 15 Skykomish 156 31 6.97 1.249 2.730 9.26 0.892 0.576 0.233 0.07931
#> 16 Squamish 14 10 7.86 0.802 2.206 8.17 0.878 0.887 0.205 -0.28884
#> 17 Tahoe 12 9 7.82 0.649 2.095 7.20 0.861 0.870 0.293 -0.11481
#> 18 Vancouver I. East 12 4 3.50 0.584 0.837 1.71 0.417 0.546 0.107 0.48361
#> 19 Willamette 17 10 6.65 1.016 2.007 5.45 0.817 0.691 0.246 0.11902
#> 20 Total 902 59 6.38 1.323 2.663 6.93 0.856 0.445 0.200 0.03621
#> rbarD File
#> 1 0.08166 Poptri.genind
#> 2 0.00883 Poptri.genind
#> 3 -0.16667 Poptri.genind
#> 4 0.02067 Poptri.genind
#> 5 -0.04770 Poptri.genind
#> 6 -0.14222 Poptri.genind
#> 7 -0.21444 Poptri.genind
#> 8 0.00445 Poptri.genind
#> 9 0.01646 Poptri.genind
#> 10 -0.25000 Poptri.genind
#> 11 -0.00606 Poptri.genind
#> 12 -0.07663 Poptri.genind
#> 13 -0.00373 Poptri.genind
#> 14 -0.01849 Poptri.genind
#> 15 0.02249 Poptri.genind
#> 16 -0.08273 Poptri.genind
#> 17 -0.03060 Poptri.genind
#> 18 0.27081 Poptri.genind
#> 19 0.04139 Poptri.genind
#> 20 0.00998 Poptri.genind
adegenet::makefreq(Poptri.genpop)
#>
#> Finding allelic frequencies from a genpop object...
#>
#> ...done.
#> X01_10838495.1 X01_10838495.2 X01_16628872.1 X01_16628872.2
#> Chilliwack 0.10526316 0.8947368 0.19736842 0.8026316
#> Columbia 0.23055556 0.7694444 0.20000000 0.8000000
#> Dean 0.21428571 0.7857143 0.00000000 1.0000000
#> Harrison 0.15909091 0.8409091 0.25000000 0.7500000
#> Homathko 0.05555556 0.9444444 0.22222222 0.7777778
#> Kitimat 0.25000000 0.7500000 0.12500000 0.8750000
#> Klinaklini 0.06250000 0.9375000 0.43750000 0.5625000
#> Nisqually 0.10714286 0.8928571 0.25000000 0.7500000
#> Nooksack 0.17307692 0.8269231 0.26923077 0.7307692
#> Olympic_Penninsula 0.50000000 0.5000000 0.10000000 0.9000000
#> Puyallup 0.17204301 0.8279570 0.16935484 0.8306452
#> Salmon 0.04545455 0.9545455 0.27272727 0.7272727
#> Skagit 0.23619632 0.7638037 0.18098160 0.8190184
#> Skwawka 0.03571429 0.9642857 0.32142857 0.6785714
#> Skykomish 0.20192308 0.7980769 0.31410256 0.6858974
#> Squamish 0.28571429 0.7142857 0.25000000 0.7500000
#> Tahoe 0.62500000 0.3750000 0.08333333 0.9166667
#> Vancouver I. East 0.16666667 0.8333333 0.00000000 1.0000000
#> Willamette 0.35294118 0.6470588 0.20588235 0.7941176
#> X01_23799134.1 X01_23799134.2 X01_38900650.1 X01_38900650.2
#> Chilliwack 0.06578947 0.9342105 0.05263158 0.9473684
#> Columbia 0.10833333 0.8916667 0.06944444 0.9305556
#> Dean 0.00000000 1.0000000 0.21428571 0.7857143
#> Harrison 0.11363636 0.8863636 0.02272727 0.9772727
#> Homathko 0.05555556 0.9444444 0.05555556 0.9444444
#> Kitimat 0.00000000 1.0000000 0.00000000 1.0000000
#> Klinaklini 0.00000000 1.0000000 0.00000000 1.0000000
#> Nisqually 0.00000000 1.0000000 0.00000000 1.0000000
#> Nooksack 0.01923077 0.9807692 0.07692308 0.9230769
#> Olympic_Penninsula 0.00000000 1.0000000 0.00000000 1.0000000
#> Puyallup 0.04569892 0.9543011 0.05107527 0.9489247
#> Salmon 0.00000000 1.0000000 0.04545455 0.9545455
#> Skagit 0.03987730 0.9601227 0.01840491 0.9815951
#> Skwawka 0.07142857 0.9285714 0.00000000 1.0000000
#> Skykomish 0.11217949 0.8878205 0.03525641 0.9647436
#> Squamish 0.03571429 0.9642857 0.03571429 0.9642857
#> Tahoe 0.04166667 0.9583333 0.41666667 0.5833333
#> Vancouver I. East 0.00000000 1.0000000 0.08333333 0.9166667
#> Willamette 0.00000000 1.0000000 0.14705882 0.8529412
#> X01_41284594.1 X01_41284594.2
#> Chilliwack 0.05263158 0.9473684
#> Columbia 0.05000000 0.9500000
#> Dean 0.00000000 1.0000000
#> Harrison 0.02272727 0.9772727
#> Homathko 0.11111111 0.8888889
#> Kitimat 0.00000000 1.0000000
#> Klinaklini 0.06250000 0.9375000
#> Nisqually 0.10714286 0.8928571
#> Nooksack 0.09615385 0.9038462
#> Olympic_Penninsula 0.10000000 0.9000000
#> Puyallup 0.07526882 0.9247312
#> Salmon 0.22727273 0.7727273
#> Skagit 0.03374233 0.9662577
#> Skwawka 0.07142857 0.9285714
#> Skykomish 0.07692308 0.9230769
#> Squamish 0.03571429 0.9642857
#> Tahoe 0.12500000 0.8750000
#> Vancouver I. East 0.04166667 0.9583333
#> Willamette 0.08823529 0.9117647
# 5. Predict future allele frequencies
Poptri.pred.future <- count.pred(Poptri.count.model, env.data=future.env)
head(Poptri.pred.future)
#> Pop Pop.label Pop.index N Allele Allele.freq A B
#> P1 Chilliwack 1 1 76 X01_10838495.1 0.10526316 8 68
#> P2 Columbia 2 2 360 X01_10838495.1 0.23055556 83 277
#> P3 Dean 3 3 14 X01_10838495.1 0.21428571 3 11
#> P4 Harrison 4 4 44 X01_10838495.1 0.15909091 7 37
#> P5 Homathko 5 5 18 X01_10838495.1 0.05555556 1 17
#> P6 Kitimat 6 6 16 X01_10838495.1 0.25000000 4 12
#> Ap Bp N.e1 Freq.e1
#> P1 14.795377 61.20462 76 0.1946760
#> P2 180.233278 179.76672 360 0.5006480
#> P3 2.158522 11.84148 14 0.1541801
#> P4 14.307597 29.69240 44 0.3251727
#> P5 4.778511 13.22149 18 0.2654728
#> P6 1.939753 14.06025 16 0.1212346
Poptri.freq.future <- freq.pred(Poptri.freq.model,
count.predicted=Poptri.pred.future)
# The key results are variables 'Allele.freq' representing the baseline allele frequencies
# and variables 'Freq.e2', the predicted frequency for the future/ past climate.
# Variable 'Freq.e1' is the predicted allele frequency in step 1
head(Poptri.freq.future)
#> Pop Pop.label Pop.index N Allele Allele.freq A B
#> P1 Chilliwack 1 1 76 X01_10838495.1 0.10526316 8 68
#> P2 Columbia 2 2 360 X01_10838495.1 0.23055556 83 277
#> P3 Dean 3 3 14 X01_10838495.1 0.21428571 3 11
#> P4 Harrison 4 4 44 X01_10838495.1 0.15909091 7 37
#> P5 Homathko 5 5 18 X01_10838495.1 0.05555556 1 17
#> P6 Kitimat 6 6 16 X01_10838495.1 0.25000000 4 12
#> Ap Bp N.e1 Freq.e1 Freq.e2 LCL UCL increasing
#> P1 14.795377 61.20462 76 0.1946760 0.1991648 0.1744943 0.2238354 TRUE
#> P2 180.233278 179.76672 360 0.5006480 0.5742628 0.4337850 0.7147405 TRUE
#> P3 2.158522 11.84148 14 0.1541801 0.1706418 0.1510378 0.1902458 FALSE
#> P4 14.307597 29.69240 44 0.3251727 0.3542175 0.3008643 0.4075707 TRUE
#> P5 4.778511 13.22149 18 0.2654728 0.2399923 0.2146111 0.2653735 TRUE
#> P6 1.939753 14.06025 16 0.1212346 0.1181220 0.1010202 0.1352239 FALSE
Visualizations
AlleleShift::population.shift
data(Poptri.baseline.env)
data(Poptri.future.env)
VIF.select <- VIF.subset(Poptri.baseline.env,
keep=c("MAT", "CMI"),
cor.plot=FALSE)
#> Step 1: Keeping these vars:
#> [1] "MAT" "CMI" "TD" "MSP" "AHM" "SHM" "DD.0" "bFFP"
#> [9] "PAS" "Eref" "CMD" "cmiJJA" "PPT_sm" "Tmax07"
#>
#> Variance inflation (package: car)
#> MAT TD MSP AHM SHM DD.0 bFFP
#> 1956.71402 55.35151 325.04715 258.27807 199.03579 383.71475 208.48116
#> PAS Eref CMD CMI cmiJJA PPT_sm Tmax07
#> 418.39278 509.80399 702.97284 729.58048 2108.47844 1221.60927 404.91923
#>
#> VIF directly calculated from linear model with focal numeric variable as response
#> cmiJJA MAT PPT_sm CMI CMD Eref PAS
#> 2108.47844 1956.71402 1221.60927 729.58048 702.97284 509.80399 418.39278
#> Tmax07 DD.0 MSP AHM bFFP SHM TD
#> 404.91923 383.71475 325.04715 258.27807 208.48116 199.03579 55.35151
#>
#> Variance inflation (package: car)
#> MAT TD MSP AHM SHM DD.0 bFFP
#> 1152.51188 51.00442 280.39905 196.06676 189.87780 289.82940 177.01912
#> PAS Eref CMD CMI PPT_sm Tmax07
#> 373.46930 302.59709 652.89020 548.59673 527.99218 186.61090
#>
#> VIF directly calculated from linear model with focal numeric variable as response
#> MAT CMD CMI PPT_sm PAS Eref DD.0
#> 1152.51188 652.89020 548.59673 527.99218 373.46930 302.59709 289.82940
#> MSP AHM SHM Tmax07 bFFP TD
#> 280.39905 196.06676 189.87780 186.61090 177.01912 51.00442
#>
#> Variance inflation (package: car)
#> MAT TD MSP AHM SHM DD.0 bFFP PAS
#> 775.21530 48.40917 275.40371 179.50626 90.11375 260.05137 146.29625 368.40367
#> Eref CMI PPT_sm Tmax07
#> 102.83305 545.11114 394.89879 180.19703
#>
#> VIF directly calculated from linear model with focal numeric variable as response
#> MAT CMI PPT_sm PAS MSP DD.0 Tmax07 AHM
#> 775.21530 545.11114 394.89879 368.40367 275.40371 260.05137 180.19703 179.50626
#> bFFP Eref SHM TD
#> 146.29625 102.83305 90.11375 48.40917
#>
#> Variance inflation (package: car)
#> MAT TD MSP AHM SHM DD.0 bFFP PAS
#> 528.61146 32.01349 145.93063 78.09567 53.38943 222.17330 145.00874 188.09374
#> Eref CMI Tmax07
#> 102.40225 306.50400 111.10074
#>
#> VIF directly calculated from linear model with focal numeric variable as response
#> MAT CMI DD.0 PAS MSP bFFP Tmax07 Eref
#> 528.61146 306.50400 222.17330 188.09374 145.93063 145.00874 111.10074 102.40225
#> AHM SHM TD
#> 78.09567 53.38943 32.01349
#>
#> Variance inflation (package: car)
#> MAT TD MSP AHM SHM bFFP PAS Eref
#> 291.46937 30.38197 64.59840 46.78232 37.94467 45.37584 45.57930 99.66232
#> CMI Tmax07
#> 93.58775 109.62229
#>
#> VIF directly calculated from linear model with focal numeric variable as response
#> MAT Tmax07 Eref CMI MSP AHM PAS bFFP
#> 291.46937 109.62229 99.66232 93.58775 64.59840 46.78232 45.57930 45.37584
#> SHM TD
#> 37.94467 30.38197
#>
#> Variance inflation (package: car)
#> MAT TD MSP AHM SHM bFFP PAS
#> 237.300281 4.513064 62.232526 46.486704 37.871152 44.616210 44.930713
#> Eref CMI
#> 33.411744 85.667431
#>
#> VIF directly calculated from linear model with focal numeric variable as response
#> MAT CMI MSP AHM PAS bFFP SHM
#> 237.300281 85.667431 62.232526 46.486704 44.930713 44.616210 37.871152
#> Eref TD
#> 33.411744 4.513064
#>
#> Variance inflation (package: car)
#> MAT TD AHM SHM bFFP PAS Eref
#> 218.240423 4.170754 25.872364 5.680412 43.733182 39.989330 32.625175
#> CMI
#> 27.162467
#>
#> VIF directly calculated from linear model with focal numeric variable as response
#> MAT bFFP PAS Eref CMI AHM SHM
#> 218.240423 43.733182 39.989330 32.625175 27.162467 25.872364 5.680412
#> TD
#> 4.170754
#>
#> Variance inflation (package: car)
#> MAT TD AHM SHM PAS Eref CMI
#> 56.168045 4.154401 24.524555 5.474406 31.877700 13.105090 22.368212
#>
#> VIF directly calculated from linear model with focal numeric variable as response
#> MAT PAS AHM CMI Eref SHM TD
#> 56.168045 31.877700 24.524555 22.368212 13.105090 5.474406 4.154401
#>
#> Variance inflation (package: car)
#> MAT TD AHM SHM Eref CMI
#> 8.713339 4.003336 19.299959 5.465614 9.025702 14.431931
#>
#> VIF directly calculated from linear model with focal numeric variable as response
#> AHM CMI Eref MAT SHM TD
#> 19.299959 14.431931 9.025702 8.713339 5.465614 4.003336
#>
#> Summary of VIF selection process:
#> MAT TD MSP AHM SHM DD.0 bFFP
#> step_1 1956.714025 55.351511 325.04715 258.27807 199.035789 383.7148 208.48116
#> step_2 1152.511883 51.004421 280.39905 196.06676 189.877801 289.8294 177.01912
#> step_3 775.215296 48.409171 275.40371 179.50626 90.113750 260.0514 146.29625
#> step_4 528.611459 32.013493 145.93063 78.09567 53.389431 222.1733 145.00874
#> step_5 291.469373 30.381965 64.59840 46.78232 37.944668 NA 45.37584
#> step_6 237.300281 4.513064 62.23253 46.48670 37.871152 NA 44.61621
#> step_7 218.240423 4.170754 NA 25.87236 5.680412 NA 43.73318
#> step_8 56.168045 4.154401 NA 24.52456 5.474406 NA NA
#> step_9 8.713339 4.003336 NA 19.29996 5.465614 NA NA
#> PAS Eref CMD CMI cmiJJA PPT_sm Tmax07
#> step_1 418.39278 509.803994 702.9728 729.58048 2108.478 1221.6093 404.9192
#> step_2 373.46930 302.597094 652.8902 548.59673 NA 527.9922 186.6109
#> step_3 368.40367 102.833047 NA 545.11114 NA 394.8988 180.1970
#> step_4 188.09374 102.402249 NA 306.50400 NA NA 111.1007
#> step_5 45.57930 99.662321 NA 93.58775 NA NA 109.6223
#> step_6 44.93071 33.411744 NA 85.66743 NA NA NA
#> step_7 39.98933 32.625175 NA 27.16247 NA NA NA
#> step_8 31.87770 13.105090 NA 22.36821 NA NA NA
#> step_9 NA 9.025702 NA 14.43193 NA NA NA
#>
#> Final selection of variables:
#> [1] "AHM" "CMI" "Eref" "MAT" "SHM" "TD"
VIF.select$vars.included
#> [1] "AHM" "CMI" "Eref" "MAT" "SHM" "TD"
baseline.env <- Poptri.baseline.env[, VIF.select$vars.included]
future.env <- Poptri.future.env[, VIF.select$vars.included]
plotA <- population.shift(baseline.env,
future.env,
option="PCA")
#> OK
plotA
plotB <- population.shift(baseline.env,
future.env,
option="RDA")
#> OK
plotB
AlleleShift::shift.dot.ggplot
# The data can be obtained via the count.model and freq.model calibrations.
# These procedures are not repeated here.
data(Poptri.freq.future)
ggdot1 <- shift.dot.ggplot(Poptri.freq.future)
ggdot1
AlleleShift::shift.pie.ggplot
# The data can be obtained via the count.model and freq.model calibrations.
# These procedures are not repeated here.
data(Poptri.freq.baseline)
data(Poptri.freq.future)
Poptri.baseline.pie <- pie.baker(Poptri.freq.baseline, r0=0.1,
sort.index="Latitude.index")
Poptri.future.pie <- pie.baker(Poptri.freq.future, r0=0.1,
freq.focus="Freq.e2",
sort.index="Latitude.index",
ypos=1)
ggpie1 <- shift.pie.ggplot(Poptri.baseline.pie,
Poptri.future.pie)
ggpie1
AlleleShift::shift.moon.ggplot
# The data can be obtained via the count.model and freq.model calibrations.
# These procedures are not repeated here.
data(Poptri.freq.baseline)
data(Poptri.freq.future)
Poptri.baseline.moon <- moon.waxer(Poptri.freq.baseline,
sort.index="Latitude.index")
Poptri.future.moon <- moon.waxer(Poptri.freq.future,
sort.index="Latitude.index",
freq.focus="Freq.e2",
ypos=1)
ggmoon1 <- shift.moon.ggplot(Poptri.baseline.moon,
Poptri.future.moon)
ggmoon1
AlleleShift::shift.waffle.ggplot
# The data can be obtained via the count.model and freq.model calibrations.
# These procedures are not repeated here.
data(Poptri.freq.baseline)
data(Poptri.freq.future)
Poptri.future.waffle <- waffle.baker(Poptri.freq.future,
sort.index="Latitude.index")
ggwaffle1 <- shift.waffle.ggplot(Poptri.future.waffle)
ggwaffle1
AlleleShift::shift.surf.ggplot
# The data can be obtained via the count.model and freq.model calibrations.
# These procedures are not repeated here.
data(Poptri.freq.baseline)
data(Poptri.freq.future)
# Plots for the first allele
# Symbols and colours indicate future change (green, ^ = future increase)
# Symbol size reflects the frequency in the climate shown
# Baseline climate
plotA <- shift.surf.ggplot(Poptri.freq.future,
xcoord="Long", ycoord="Lat",
Allele.focus=unique(Poptri.freq.future$Allele)[1],
freq.focus="Allele.freq")
plotA
# Future/past climate
plotB <- shift.surf.ggplot(Poptri.freq.future,
xcoord="Long", ycoord="Lat",
Allele.focus=unique(Poptri.freq.future$Allele)[1],
freq.focus="Freq.e2")
plotB
# Plots for the fifth allele
# Baseline climate
plotC <- shift.surf.ggplot(Poptri.freq.future,
xcoord="Long", ycoord="Lat",
Allele.focus=unique(Poptri.freq.future$Allele)[5],
freq.focus="Allele.freq")
plotC
plotD <- shift.surf.ggplot(Poptri.freq.future,
xcoord="Long", ycoord="Lat",
Allele.focus=unique(Poptri.freq.future$Allele)[5],
freq.focus="Freq.e2")
plotD
Animations
With graphs that are generated via ggplot2, it is straightforward to create animations. The following example is also avaiable from the documentation of ‘AlleleShift::shift.dot.ggplot’
# create an animation
library(ggplot2)
library(gganimate)
library(gifski)
#> Warning: package 'gifski' was built under R version 4.0.3
# The data is an interpolation and extrapolation between the baseline and future climate.
# For actual application, interpolate between climate data from available sources
data(Poptri.1985to2085)
ggdot.all <- ggplot(data=Poptri.1985to2085, group=Decade) +
scale_y_continuous(limits=c(-0.1, 1.1),
breaks=c(0.0, 0.25, 0.5, 0.75, 1.0)) +
geom_errorbar(aes(x=Pop, ymin=LCL, ymax=UCL),
colour="grey30", width=0.8, show.legend=FALSE) +
geom_segment(aes(x=Pop, y=Allele.freq, xend=Pop, yend=Freq.e2, colour=increasing),
size=1.2) +
geom_point(aes(x=Pop, y=Allele.freq),
colour="black", size=10, alpha=0.7) +
geom_point(aes(x=Pop, y=Freq.e2),
colour="dodgerblue3", size=10, alpha=0.7) +
coord_flip() +
xlab(element_blank()) +
ylab("Allele frequencies") +
theme(panel.grid.minor = element_blank()) +
labs(colour="Future change in allele frequencies") +
scale_colour_manual(values=c("firebrick3", "chartreuse4"),
labels=c("decreasing", "increasing")) +
theme(axis.text.x=element_text(angle=90, vjust=0.5, size=10)) +
theme(legend.position="top") +
facet_grid( ~ Allele, scales="free")
ggdot.all
ggdot.anim <- ggdot.all +
transition_states(as.factor(Decade), transition_length = 10, state_length = 100) +
labs(title = "Decade: {closest_state}s")
Show the animation
animate(ggdot.anim, fps=5, width=1280, height=720)
sessionInfo()
#> R version 4.0.2 (2020-06-22)
#> Platform: x86_64-w64-mingw32/x64 (64-bit)
#> Running under: Windows 10 x64 (build 18363)
#>
#> Matrix products: default
#>
#> locale:
#> [1] LC_COLLATE=English_United Kingdom.1252
#> [2] LC_CTYPE=English_United Kingdom.1252
#> [3] LC_MONETARY=English_United Kingdom.1252
#> [4] LC_NUMERIC=C
#> [5] LC_TIME=English_United Kingdom.1252
#>
#> attached base packages:
#> [1] tcltk stats graphics grDevices utils datasets methods
#> [8] base
#>
#> other attached packages:
#> [1] gifski_0.8.6 AlleleShift_0.9 gganimate_1.0.7
#> [4] gggibbous_0.1.0 ggmap_3.0.0 mgcv_1.8-31
#> [7] nlme_3.1-148 GGally_2.0.0 patchwork_1.1.1
#> [10] ggrepel_0.8.2 dplyr_1.0.2 ggforce_0.3.2
#> [13] ggsci_2.9 ggplot2_3.3.2 poppr_2.8.6
#> [16] adegenet_2.1.3 ade4_1.7-16 BiodiversityR_2.12-3
#> [19] vegan_2.5-6 lattice_0.20-41 permute_0.9-5
#>
#> loaded via a namespace (and not attached):
#> [1] readxl_1.3.1 backports_1.1.7 Hmisc_4.4-0
#> [4] fastmatch_1.1-0 plyr_1.8.6 igraph_1.2.6
#> [7] sp_1.4-2 splines_4.0.2 digest_0.6.25
#> [10] htmltools_0.5.0 gdata_2.18.0 relimp_1.0-5
#> [13] magrittr_1.5 checkmate_2.0.0 cluster_2.1.0
#> [16] openxlsx_4.1.5 tcltk2_1.2-11 gmodels_2.18.1
#> [19] sandwich_2.5-1 prettyunits_1.1.1 jpeg_0.1-8.1
#> [22] colorspace_1.4-1 mitools_2.4 haven_2.3.1
#> [25] xfun_0.15 crayon_1.3.4 lme4_1.1-23
#> [28] survival_3.1-12 zoo_1.8-8 phangorn_2.5.5
#> [31] ape_5.4 glue_1.4.1 polyclip_1.10-0
#> [34] gtable_0.3.0 seqinr_4.2-4 polysat_1.7-4
#> [37] car_3.0-8 abind_1.4-5 scales_1.1.1
#> [40] DBI_1.1.0 Rcpp_1.0.4.6 isoband_0.2.1
#> [43] viridisLite_0.3.0 xtable_1.8-4 progress_1.2.2
#> [46] spData_0.3.8 htmlTable_2.0.0 units_0.6-7
#> [49] foreign_0.8-80 spdep_1.1-5 Formula_1.2-3
#> [52] survey_4.0 httr_1.4.2 htmlwidgets_1.5.1
#> [55] RColorBrewer_1.1-2 acepack_1.4.1 ellipsis_0.3.1
#> [58] reshape_0.8.8 pkgconfig_2.0.3 farver_2.0.3
#> [61] nnet_7.3-14 deldir_0.2-3 labeling_0.3
#> [64] tidyselect_1.1.0 rlang_0.4.8 reshape2_1.4.4
#> [67] later_1.1.0.1 munsell_0.5.0 cellranger_1.1.0
#> [70] tools_4.0.2 generics_0.1.0 evaluate_0.14
#> [73] stringr_1.4.0 fastmap_1.0.1 yaml_2.2.1
#> [76] knitr_1.28 zip_2.0.4 purrr_0.3.4
#> [79] RgoogleMaps_1.4.5.3 mime_0.9 compiler_4.0.2
#> [82] rstudioapi_0.11 curl_4.3 png_0.1-7
#> [85] e1071_1.7-3 tibble_3.0.1 statmod_1.4.34
#> [88] tweenr_1.0.1 stringi_1.4.6 forcats_0.5.0
#> [91] Matrix_1.2-18 classInt_0.4-3 nloptr_1.2.2.1
#> [94] vctrs_0.3.4 effects_4.1-4 RcmdrMisc_2.7-0
#> [97] pillar_1.4.4 LearnBayes_2.15.1 lifecycle_0.2.0
#> [100] bitops_1.0-6 data.table_1.12.8 raster_3.4-5
#> [103] httpuv_1.5.4 R6_2.4.1 latticeExtra_0.6-29
#> [106] promises_1.1.1 KernSmooth_2.23-17 gridExtra_2.3
#> [109] rio_0.5.16 codetools_0.2-16 boot_1.3-25
#> [112] MASS_7.3-51.6 gtools_3.8.2 rjson_0.2.20
#> [115] withr_2.2.0 nortest_1.0-4 Rcmdr_2.6-2
#> [118] pegas_0.14 expm_0.999-5 parallel_4.0.2
#> [121] hms_0.5.3 quadprog_1.5-8 grid_4.0.2
#> [124] rpart_4.1-15 tidyr_1.1.2 coda_0.19-3
#> [127] class_7.3-17 minqa_1.2.4 rmarkdown_2.3
#> [130] carData_3.0-4 sf_0.9-6 shiny_1.4.0.2
#> [133] base64enc_0.1-3