Discriminant analysis of principal components will give us some more insights into possible clustering within the population.
We would like to se if there are effects of our treatments, variety and fungicide programmes.
Another interesting thing to look at is year to year variation of the pathogen population structure.
list.of.packages <-
c(
"tidyverse",
"devtools",
"here",
"readxl",
"poppr",
"egg",
"parallel",
"conflicted"
)
new.packages <-
list.of.packages[!(list.of.packages %in% installed.packages()[, "Package"])]
#Download packages that are not already present in the library
if (length(new.packages))
install.packages(new.packages)
packages_load <-
lapply(list.of.packages, require, character.only = TRUE)
#Print warning if there is a problem with installing/loading some of packages
if (any(as.numeric(packages_load) == 0)) {
warning(paste("Package/s: ", paste(list.of.packages[packages_load != TRUE], sep = ", "), "not loaded!"))
} else {
print("All packages were successfully loaded.")
}
## [1] "All packages were successfully loaded."
rm(list.of.packages, new.packages, packages_load)
#if instal is not working try
#install.packages("package_name", repos = c(CRAN="https://cran.r-project.org/"))
path <- here::here("data", "gen data", "final", "pop.csv")
monpop <- read.genalex(path, ploidy = 3)
splitStrata(monpop) <- ~variety/genotype/year
monpop
##
## This is a genclone object
## -------------------------
## Genotype information:
##
## 29 original multilocus genotypes
## 1287 triploid individuals
## 12 codominant loci
##
## Population information:
##
## 3 strata - variety, genotype, year
## 52 populations defined -
## KE_6A1_2016, KE_8A1_2016, SE_6A1_2016, ..., KE_36A2_2019, SM_13A2_2019, SM_6A1_2019
Very good tutorials are avaialble at adegent page and Grünwald lab web page.
Seting the paralel processing. This argument should be removed if the code is to be reproduced on machine which can not devote muliple core.
We will also assign a number of cores to be used. We will set it to be one less than the maximum no. of cores, so we dont ‘sufficate’ the computer.
#Determine the OS to assign apropriate character string for parallel processing
parallel_proc <-
ifelse(Sys.info()['sysname'] == "Windows", "snow", "multicore")
cpus <- ifelse(detectCores()>1, c(detectCores()-1), 1)
cpus
## [1] 3
setPop(monpop) <- ~variety
set.seed(999)
system.time(pramx_var <- xvalDapc(tab(monpop, NA.method = "mean"), pop(monpop),
n.pca = 5:20, n.rep = 3000,
parallel = parallel_proc,
ncpus = cpus))
## user system elapsed
## 1.25 5.69 901.80
#Reorder levels of the factor for plotting
pramx_var$DAPC$grp <-
factor(pramx_var$DAPC$grp, levels = c("KE", "BQ", "RO", "SE", "CL", "SM"))
saveRDS(pramx_var, file = here::here("results", "gen","dapc", "dapc_var.RDS"))
pramx_var[-1]
## $`Median and Confidence Interval for Random Chance`
## 2.5% 50% 97.5%
## 0.1467338 0.1661992 0.1886241
##
## $`Mean Successful Assignment by Number of PCs of PCA`
## 5 6 7 8 9 10 11
## 0.2727234 0.2728408 0.3548089 0.3339103 0.3713358 0.3708384 0.3535601
## 12 13 14 15 16 17 18
## 0.3537474 0.3523758 0.3521737 0.3551396 0.3539692 0.3492212 0.3481205
## 19 20
## 0.3483753 0.3431488
##
## $`Number of PCs Achieving Highest Mean Success`
## [1] "9"
##
## $`Root Mean Squared Error by Number of PCs of PCA`
## 5 6 7 8 9 10 11
## 0.7278316 0.7280140 0.6474076 0.6672821 0.6302772 0.6308389 0.6488356
## 12 13 14 15 16 17 18
## 0.6487290 0.6501647 0.6503099 0.6473630 0.6484855 0.6533128 0.6543361
## 19 20
## 0.6542152 0.6592532
##
## $`Number of PCs Achieving Lowest MSE`
## [1] "9"
##
## $DAPC
## #################################################
## # Discriminant Analysis of Principal Components #
## #################################################
## class: dapc
## $call: dapc.data.frame(x = as.data.frame(x), grp = ..1, n.pca = ..2,
## n.da = ..3)
##
## $n.pca: 9 first PCs of PCA used
## $n.da: 5 discriminant functions saved
## $var (proportion of conserved variance): 0.988
##
## $eig (eigenvalues): 173.9 25.19 12.21 1.484 0.7428 vector length content
## 1 $eig 5 eigenvalues
## 2 $grp 1287 prior group assignment
## 3 $prior 6 prior group probabilities
## 4 $assign 1287 posterior group assignment
## 5 $pca.cent 60 centring vector of PCA
## 6 $pca.norm 60 scaling vector of PCA
## 7 $pca.eig 23 eigenvalues of PCA
##
## data.frame nrow ncol
## 1 $tab 1287 9
## 2 $means 6 9
## 3 $loadings 9 5
## 4 $ind.coord 1287 5
## 5 $grp.coord 6 5
## 6 $posterior 1287 6
## 7 $pca.loadings 60 9
## 8 $var.contr 60 5
## content
## 1 retained PCs of PCA
## 2 group means
## 3 loadings of variables
## 4 coordinates of individuals (principal components)
## 5 coordinates of groups
## 6 posterior membership probabilities
## 7 PCA loadings of original variables
## 8 contribution of original variables
pramx_var[-1]$`Number of PCs Achieving Highest Mean Success`
## [1] "9"
scatter(
pramx_var$DAPC,
legend = TRUE,
cex = 3,
# pch=18:23,
posi.leg = "topleft",
scree.da = FALSE,
clab =2,
cleg = 2.1,
cell = 1, #zoom
xax = 1,
yax = 2,
inset.solid = 1,
cstar=0.5,
lwd=.1
)
Prepare the data in genalex format once again. We can define only thre stratums, so we will prepare another file which will be re-used for same purpose as above.
fin <- readRDS(file = here::here("data", "gen data", "final", "gendata.rds") )
fin <-
unite(fin, "Ind", idd, year, remove = F)
fin <-
unite(fin, "Pop", treatment, Genotype, year)
popcol <- match("Pop",names(fin))
idcol <- match("Ind",names(fin))
gen <-
fin[ , c(idcol,popcol,match("D13",names(fin)):match("X__24",names(fin)))]
path <- here::here("data", "gen data", "final", "pop2.csv")
data.frame(ssr = c(12,"PInf_Ireland"),
samples = nrow(gen),
pop = length(unique(fin$Pop)))
write.table(data.frame(ssr = 12, samples = nrow(gen), pop = length(unique(fin$Pop))),
path,
sep=",", col.names=FALSE, row.names = FALSE)
write.table( "PInf_Ireland",
path,
sep = ",", row.names = FALSE,
col.names = FALSE, append = T)
names(gen)[grep("X__", names(gen))] <- ""
write.table(gen,
path,
sep = ",",
row.names = FALSE,
col.names = !file.exists("myDF.csv"), append = T)
So now we have created another file in genalex format, pop2.csv.
path <- here::here("data", "gen data", "final", "pop2.csv")
monpop <- read.genalex(path, ploidy = 3)
splitStrata(monpop) <- ~treatment/genotype/year
monpop
##
## This is a genclone object
## -------------------------
## Genotype information:
##
## 29 original multilocus genotypes
## 1287 triploid individuals
## 12 codominant loci
##
## Population information:
##
## 3 strata - treatment, genotype, year
## 61 populations defined -
## 0_6A1_2016, 0_8A1_2016, 100_6A1_2016, ..., IR_6A1_2019, 50_6A1_2019, MIR_36A2_2019
setPop(monpop) <- ~treatment
set.seed(999)
system.time(pramx_trt <- xvalDapc(tab(monpop, NA.method = "mean"), pop(monpop),
n.pca = 5:20,
n.rep = 3000,
parallel = parallel_proc,
ncpus = cpus))
## user system elapsed
## 1.50 7.36 1041.64
pramx_trt$DAPC$grp <-
factor(pramx_trt$DAPC$grp, levels = c("0", "100", "50", "IR", "BM", "MIR"))
saveRDS(pramx_trt, file = here::here("results", "gen","dapc", "dapc_trt.RDS"))
pramx_trt[-1]
## $`Median and Confidence Interval for Random Chance`
## 2.5% 50% 97.5%
## 0.1451753 0.1667619 0.1879970
##
## $`Mean Successful Assignment by Number of PCs of PCA`
## 5 6 7 8 9 10 11
## 0.1617465 0.1789775 0.1876726 0.1775959 0.1817178 0.1964336 0.1939121
## 12 13 14 15 16 17 18
## 0.1928035 0.1929071 0.1915203 0.1911883 0.1887293 0.1946328 0.1935559
## 19 20
## 0.1935899 0.1919025
##
## $`Number of PCs Achieving Highest Mean Success`
## [1] "10"
##
## $`Root Mean Squared Error by Number of PCs of PCA`
## 5 6 7 8 9 10 11
## 0.8383041 0.8212036 0.8126004 0.8226078 0.8185173 0.8038847 0.8063984
## 12 13 14 15 16 17 18
## 0.8075041 0.8073959 0.8087795 0.8091192 0.8115764 0.8057035 0.8067724
## 19 20
## 0.8067440 0.8084247
##
## $`Number of PCs Achieving Lowest MSE`
## [1] "10"
##
## $DAPC
## #################################################
## # Discriminant Analysis of Principal Components #
## #################################################
## class: dapc
## $call: dapc.data.frame(x = as.data.frame(x), grp = ..1, n.pca = ..2,
## n.da = ..3)
##
## $n.pca: 10 first PCs of PCA used
## $n.da: 5 discriminant functions saved
## $var (proportion of conserved variance): 0.991
##
## $eig (eigenvalues): 15.81 9.464 3.786 1.946 1.354 vector length content
## 1 $eig 5 eigenvalues
## 2 $grp 1287 prior group assignment
## 3 $prior 6 prior group probabilities
## 4 $assign 1287 posterior group assignment
## 5 $pca.cent 60 centring vector of PCA
## 6 $pca.norm 60 scaling vector of PCA
## 7 $pca.eig 23 eigenvalues of PCA
##
## data.frame nrow ncol
## 1 $tab 1287 10
## 2 $means 6 10
## 3 $loadings 10 5
## 4 $ind.coord 1287 5
## 5 $grp.coord 6 5
## 6 $posterior 1287 6
## 7 $pca.loadings 60 10
## 8 $var.contr 60 5
## content
## 1 retained PCs of PCA
## 2 group means
## 3 loadings of variables
## 4 coordinates of individuals (principal components)
## 5 coordinates of groups
## 6 posterior membership probabilities
## 7 PCA loadings of original variables
## 8 contribution of original variables
scatter(
pramx_trt$DAPC,
legend = TRUE,
cex = 3,
# pch=18:23,
posi.leg = "topleft",
scree.da = FALSE,
clab =2,
cleg = 2.1,
cell = 1, #zoom
xax = 1,
yax = 2,
inset.solid = 1,
cstar=0.5,
lwd=.1
)
setPop(monpop) <- ~year
set.seed(999)
system.time(pramx_yr <- xvalDapc(tab(monpop, NA.method = "mean"), pop(monpop),
n.pca = 5:20,
n.rep = 3000,
parallel = parallel_proc,
ncpus = cpus))
## user system elapsed
## 1.52 5.57 868.18
pramx_yr$DAPC$grp <-
factor(pramx_yr$DAPC$grp, levels = c("2016", "2017", "2018", "2019"))
saveRDS(pramx_yr, file = here::here("results", "gen","dapc", "dapc_yr.RDS"))
pramx_yr[-1]
## $`Median and Confidence Interval for Random Chance`
## 2.5% 50% 97.5%
## 0.2257941 0.2507039 0.2773436
##
## $`Mean Successful Assignment by Number of PCs of PCA`
## 5 6 7 8 9 10 11
## 0.5812706 0.6430049 0.6420921 0.6559491 0.6551287 0.6686413 0.6685358
## 12 13 14 15 16 17 18
## 0.6687640 0.6752763 0.6753330 0.6759402 0.6759702 0.6777685 0.6763974
## 19 20
## 0.6774251 0.6774868
##
## $`Number of PCs Achieving Highest Mean Success`
## [1] "17"
##
## $`Root Mean Squared Error by Number of PCs of PCA`
## 5 6 7 8 9 10 11
## 0.4224605 0.3594239 0.3605585 0.3465186 0.3474015 0.3338579 0.3340508
## 12 13 14 15 16 17 18
## 0.3339987 0.3273422 0.3272862 0.3266113 0.3266299 0.3247128 0.3260877
## 19 20
## 0.3252855 0.3252106
##
## $`Number of PCs Achieving Lowest MSE`
## [1] "17"
##
## $DAPC
## #################################################
## # Discriminant Analysis of Principal Components #
## #################################################
## class: dapc
## $call: dapc.data.frame(x = as.data.frame(x), grp = ..1, n.pca = ..2,
## n.da = ..3)
##
## $n.pca: 17 first PCs of PCA used
## $n.da: 3 discriminant functions saved
## $var (proportion of conserved variance): 0.999
##
## $eig (eigenvalues): 1247 372.5 251.7 vector length content
## 1 $eig 3 eigenvalues
## 2 $grp 1287 prior group assignment
## 3 $prior 4 prior group probabilities
## 4 $assign 1287 posterior group assignment
## 5 $pca.cent 60 centring vector of PCA
## 6 $pca.norm 60 scaling vector of PCA
## 7 $pca.eig 23 eigenvalues of PCA
##
## data.frame nrow ncol
## 1 $tab 1287 17
## 2 $means 4 17
## 3 $loadings 17 3
## 4 $ind.coord 1287 3
## 5 $grp.coord 4 3
## 6 $posterior 1287 4
## 7 $pca.loadings 60 17
## 8 $var.contr 60 3
## content
## 1 retained PCs of PCA
## 2 group means
## 3 loadings of variables
## 4 coordinates of individuals (principal components)
## 5 coordinates of groups
## 6 posterior membership probabilities
## 7 PCA loadings of original variables
## 8 contribution of original variables
scatter(
pramx_yr$DAPC,
legend = TRUE,
cex = 3,
# pch=18:23,
posi.leg = "topleft",
scree.da = FALSE,
clab =1.8,
cleg = 2,
cell = 1, #zoom
xax = 1,
yax = 2,
inset.solid = 1,
cstar=0.5,
lwd=.1
)
session_info()
## - Session info ----------------------------------------------------------
## setting value
## version R version 3.6.3 (2020-02-29)
## os Windows 10 x64
## system x86_64, mingw32
## ui RTerm
## language (EN)
## collate English_United States.1252
## ctype English_United States.1252
## tz Europe/London
## date 2020-03-17
##
## - Packages --------------------------------------------------------------
## package * version date lib source
## ade4 * 1.7-15 2020-02-13 [1] CRAN (R 3.6.3)
## adegenet * 2.1.2 2020-01-20 [1] CRAN (R 3.6.3)
## ape 5.3 2019-03-17 [1] CRAN (R 3.6.1)
## assertthat 0.2.1 2019-03-21 [1] CRAN (R 3.6.1)
## backports 1.1.5 2019-10-02 [1] CRAN (R 3.6.1)
## boot 1.3-24 2019-12-20 [2] CRAN (R 3.6.3)
## broom 0.5.5 2020-02-29 [1] CRAN (R 3.6.3)
## callr 3.4.2 2020-02-12 [1] CRAN (R 3.6.3)
## cellranger 1.1.0 2016-07-27 [1] CRAN (R 3.6.1)
## class 7.3-15 2019-01-01 [2] CRAN (R 3.6.3)
## classInt 0.4-2 2019-10-17 [1] CRAN (R 3.6.1)
## cli 2.0.2 2020-02-28 [1] CRAN (R 3.6.3)
## cluster 2.1.0 2019-06-19 [2] CRAN (R 3.6.3)
## coda 0.19-3 2019-07-05 [1] CRAN (R 3.6.1)
## codetools 0.2-16 2018-12-24 [2] CRAN (R 3.6.3)
## colorspace 1.4-1 2019-03-18 [1] CRAN (R 3.6.1)
## conflicted * 1.0.4 2019-06-21 [1] CRAN (R 3.6.1)
## crayon 1.3.4 2017-09-16 [1] CRAN (R 3.6.1)
## DBI 1.1.0 2019-12-15 [1] CRAN (R 3.6.3)
## deldir 0.1-25 2020-02-03 [1] CRAN (R 3.6.2)
## desc 1.2.0 2018-05-01 [1] CRAN (R 3.6.1)
## devtools * 2.2.2 2020-02-17 [1] CRAN (R 3.6.3)
## digest 0.6.25 2020-02-23 [1] CRAN (R 3.6.3)
## dplyr * 0.8.5 2020-03-07 [1] CRAN (R 3.6.3)
## e1071 1.7-3 2019-11-26 [1] CRAN (R 3.6.3)
## egg * 0.4.5 2019-07-13 [1] CRAN (R 3.6.1)
## ellipsis 0.3.0 2019-09-20 [1] CRAN (R 3.6.1)
## evaluate 0.14 2019-05-28 [1] CRAN (R 3.6.1)
## expm 0.999-4 2019-03-21 [1] CRAN (R 3.6.1)
## fansi 0.4.1 2020-01-08 [1] CRAN (R 3.6.3)
## fastmap 1.0.1 2019-10-08 [1] CRAN (R 3.6.1)
## fastmatch 1.1-0 2017-01-28 [1] CRAN (R 3.6.0)
## forcats * 0.4.0 2019-02-17 [1] CRAN (R 3.6.1)
## fs 1.3.1 2019-05-06 [1] CRAN (R 3.6.1)
## gdata 2.18.0 2017-06-06 [1] CRAN (R 3.6.0)
## generics 0.0.2 2018-11-29 [1] CRAN (R 3.6.1)
## ggplot2 * 3.2.1 2019-08-10 [1] CRAN (R 3.6.1)
## glue 1.3.1 2019-03-12 [1] CRAN (R 3.6.3)
## gmodels 2.18.1 2018-06-25 [1] CRAN (R 3.6.1)
## gridExtra * 2.3 2017-09-09 [1] CRAN (R 3.6.1)
## gtable 0.3.0 2019-03-25 [1] CRAN (R 3.6.1)
## gtools 3.8.1 2018-06-26 [1] CRAN (R 3.6.0)
## haven 2.1.1 2019-07-04 [1] CRAN (R 3.6.1)
## here * 0.1 2017-05-28 [1] CRAN (R 3.6.1)
## hms 0.5.2 2019-10-30 [1] CRAN (R 3.6.1)
## htmltools 0.4.0 2019-10-04 [1] CRAN (R 3.6.1)
## httpuv 1.5.2 2019-09-11 [1] CRAN (R 3.6.1)
## httr 1.4.1 2019-08-05 [1] CRAN (R 3.6.1)
## igraph 1.2.4.1 2019-04-22 [1] CRAN (R 3.6.1)
## jsonlite 1.6.1 2020-02-02 [1] CRAN (R 3.6.3)
## KernSmooth 2.23-16 2019-10-15 [2] CRAN (R 3.6.3)
## knitr 1.25 2019-09-18 [1] CRAN (R 3.6.1)
## later 1.0.0 2019-10-04 [1] CRAN (R 3.6.1)
## lattice 0.20-38 2018-11-04 [2] CRAN (R 3.6.3)
## lazyeval 0.2.2 2019-03-15 [1] CRAN (R 3.6.1)
## LearnBayes 2.15.1 2018-03-18 [1] CRAN (R 3.6.0)
## lifecycle 0.1.0 2019-08-01 [1] CRAN (R 3.6.1)
## lubridate 1.7.4 2018-04-11 [1] CRAN (R 3.6.1)
## magrittr 1.5 2014-11-22 [1] CRAN (R 3.6.1)
## MASS 7.3-51.5 2019-12-20 [2] CRAN (R 3.6.3)
## Matrix 1.2-18 2019-11-27 [2] CRAN (R 3.6.3)
## memoise 1.1.0 2017-04-21 [1] CRAN (R 3.6.1)
## mgcv 1.8-31 2019-11-09 [2] CRAN (R 3.6.3)
## mime 0.7 2019-06-11 [1] CRAN (R 3.6.0)
## modelr 0.1.5 2019-08-08 [1] CRAN (R 3.6.1)
## munsell 0.5.0 2018-06-12 [1] CRAN (R 3.6.1)
## nlme 3.1-144 2020-02-06 [2] CRAN (R 3.6.3)
## pegas 0.12 2019-10-05 [1] CRAN (R 3.6.1)
## permute 0.9-5 2019-03-12 [1] CRAN (R 3.6.1)
## phangorn 2.5.5 2019-06-19 [1] CRAN (R 3.6.1)
## pillar 1.4.2 2019-06-29 [1] CRAN (R 3.6.1)
## pkgbuild 1.0.6 2019-10-09 [1] CRAN (R 3.6.1)
## pkgconfig 2.0.3 2019-09-22 [1] CRAN (R 3.6.1)
## pkgload 1.0.2 2018-10-29 [1] CRAN (R 3.6.1)
## plyr 1.8.4 2016-06-08 [1] CRAN (R 3.6.1)
## polysat 1.7-4 2019-03-06 [1] CRAN (R 3.6.1)
## poppr * 2.8.3 2019-06-18 [1] CRAN (R 3.6.1)
## prettyunits 1.0.2 2015-07-13 [1] CRAN (R 3.6.1)
## processx 3.4.1 2019-07-18 [1] CRAN (R 3.6.1)
## promises 1.1.0 2019-10-04 [1] CRAN (R 3.6.1)
## ps 1.3.0 2018-12-21 [1] CRAN (R 3.6.1)
## purrr * 0.3.3 2019-10-18 [1] CRAN (R 3.6.1)
## quadprog 1.5-7 2019-05-06 [1] CRAN (R 3.6.0)
## R6 2.4.1 2019-11-12 [1] CRAN (R 3.6.3)
## raster 3.0-7 2019-09-24 [1] CRAN (R 3.6.1)
## Rcpp 1.0.2 2019-07-25 [1] CRAN (R 3.6.1)
## readr * 1.3.1 2018-12-21 [1] CRAN (R 3.6.1)
## readxl * 1.3.1 2019-03-13 [1] CRAN (R 3.6.1)
## remotes 2.1.1 2020-02-15 [1] CRAN (R 3.6.3)
## reshape2 1.4.3 2017-12-11 [1] CRAN (R 3.6.1)
## rlang 0.4.5 2020-03-01 [1] CRAN (R 3.6.3)
## rmarkdown 1.16 2019-10-01 [1] CRAN (R 3.6.1)
## rprojroot 1.3-2 2018-01-03 [1] CRAN (R 3.6.1)
## rstudioapi 0.11 2020-02-07 [1] CRAN (R 3.6.3)
## rvest 0.3.4 2019-05-15 [1] CRAN (R 3.6.1)
## scales 1.0.0 2018-08-09 [1] CRAN (R 3.6.1)
## seqinr 3.6-1 2019-09-07 [1] CRAN (R 3.6.1)
## sessioninfo 1.1.1 2018-11-05 [1] CRAN (R 3.6.1)
## sf 0.8-0 2019-09-17 [1] CRAN (R 3.6.1)
## shiny 1.4.0 2019-10-10 [1] CRAN (R 3.6.1)
## sp 1.3-1 2018-06-05 [1] CRAN (R 3.6.1)
## spData 0.3.2 2019-09-19 [1] CRAN (R 3.6.1)
## spdep 1.1-3 2019-09-18 [1] CRAN (R 3.6.1)
## stringi 1.4.6 2020-02-17 [1] CRAN (R 3.6.2)
## stringr * 1.4.0 2019-02-10 [1] CRAN (R 3.6.1)
## testthat 2.2.1 2019-07-25 [1] CRAN (R 3.6.1)
## tibble * 2.1.3 2019-06-06 [1] CRAN (R 3.6.1)
## tidyr * 1.0.0 2019-09-11 [1] CRAN (R 3.6.1)
## tidyselect 0.2.5 2018-10-11 [1] CRAN (R 3.6.1)
## tidyverse * 1.2.1 2017-11-14 [1] CRAN (R 3.6.1)
## units 0.6-5 2019-10-08 [1] CRAN (R 3.6.1)
## usethis * 1.5.1 2019-07-04 [1] CRAN (R 3.6.1)
## vctrs 0.2.0 2019-07-05 [1] CRAN (R 3.6.1)
## vegan 2.5-6 2019-09-01 [1] CRAN (R 3.6.1)
## withr 2.1.2 2018-03-15 [1] CRAN (R 3.6.1)
## xfun 0.10 2019-10-01 [1] CRAN (R 3.6.1)
## xml2 1.2.2 2019-08-09 [1] CRAN (R 3.6.1)
## xtable 1.8-4 2019-04-21 [1] CRAN (R 3.6.1)
## yaml 2.2.0 2018-07-25 [1] CRAN (R 3.6.0)
## zeallot 0.1.0 2018-01-28 [1] CRAN (R 3.6.1)
##
## [1] C:/Users/mlade/Documents/R/win-library/3.6
## [2] C:/Program Files/R/R-3.6.3/library