BiodiversityR
is an R package for statistical analysis of biodiversity and ecological communities, including species accumulation curves, diversity indices, Renyi profiles, GLMs for analysis of species abundance and presence-absence, distance matrices, Mantel tests, and cluster, constrained and unconstrained ordination analysis.
The package was initially built to provide a graphical user interface and different helper functions for the vegan
community ecology package. As documented in a manual on biodiversity and community ecology analysis, available from this website and this link, most analysis pipelines require a community matrix (typically having sites as rows, species as columns and abundance values as cell values) and an environmental data set (typically providing numerical and categorical variables for the different sites) as inputs. After being launched on CRAN, new methods from vegan
have been integrated in the package as documented in the ChangeLog file.
A major update of BiodiversityR
involved the inclusion of different methods of ensemble suitability modelling, as described here. Most of the functions for ensemble suitability modelling were named with a ensemble.XXXX() naming template.
At the end of 2021, some methods to facilitate plotting results via ggplot
were included in BiodiversityR
. These methods are described in detail in this rpubs series, with some examples given below.
library(BiodiversityR) # also loads vegan
library(ggplot2)
library(ggforce)
library(concaveman)
library(ggrepel)
library(ggsci)
library(dplyr)
library(pvclust)
The data in the examples are data that were used as case study in Data Analysis in Community and Landscape Ecology and also in the Tree Diversity Analysis manual.
Data set dune is a community data set, where variables (columns) typically correspond to different species and data represents abundance of each species. Species names were abbreviated to eight characters, with for example Agrostol representing Agrostis stolonifera.
data(dune)
str(dune)
#> 'data.frame': 20 obs. of 30 variables:
#> $ Achimill: num 1 3 0 0 2 2 2 0 0 4 ...
#> $ Agrostol: num 0 0 4 8 0 0 0 4 3 0 ...
#> $ Airaprae: num 0 0 0 0 0 0 0 0 0 0 ...
#> $ Alopgeni: num 0 2 7 2 0 0 0 5 3 0 ...
#> $ Anthodor: num 0 0 0 0 4 3 2 0 0 4 ...
#> $ Bellpere: num 0 3 2 2 2 0 0 0 0 2 ...
#> $ Bromhord: num 0 4 0 3 2 0 2 0 0 4 ...
#> $ Chenalbu: num 0 0 0 0 0 0 0 0 0 0 ...
#> $ Cirsarve: num 0 0 0 2 0 0 0 0 0 0 ...
#> $ Comapalu: num 0 0 0 0 0 0 0 0 0 0 ...
#> $ Eleopalu: num 0 0 0 0 0 0 0 4 0 0 ...
#> $ Elymrepe: num 4 4 4 4 4 0 0 0 6 0 ...
#> $ Empenigr: num 0 0 0 0 0 0 0 0 0 0 ...
#> $ Hyporadi: num 0 0 0 0 0 0 0 0 0 0 ...
#> $ Juncarti: num 0 0 0 0 0 0 0 4 4 0 ...
#> $ Juncbufo: num 0 0 0 0 0 0 2 0 4 0 ...
#> $ Lolipere: num 7 5 6 5 2 6 6 4 2 6 ...
#> $ Planlanc: num 0 0 0 0 5 5 5 0 0 3 ...
#> $ Poaprat : num 4 4 5 4 2 3 4 4 4 4 ...
#> $ Poatriv : num 2 7 6 5 6 4 5 4 5 4 ...
#> $ Ranuflam: num 0 0 0 0 0 0 0 2 0 0 ...
#> $ Rumeacet: num 0 0 0 0 5 6 3 0 2 0 ...
#> $ Sagiproc: num 0 0 0 5 0 0 0 2 2 0 ...
#> $ Salirepe: num 0 0 0 0 0 0 0 0 0 0 ...
#> $ Scorautu: num 0 5 2 2 3 3 3 3 2 3 ...
#> $ Trifprat: num 0 0 0 0 2 5 2 0 0 0 ...
#> $ Trifrepe: num 0 5 2 1 2 5 2 2 3 6 ...
#> $ Vicilath: num 0 0 0 0 0 0 0 0 0 1 ...
#> $ Bracruta: num 0 0 2 2 2 6 2 2 2 2 ...
#> $ Callcusp: num 0 0 0 0 0 0 0 0 0 0 ...
Data set dune.env is an environmental data set, where variables (columns) correspond to different descriptors (typically continuous and categorical variables) of the sample sites. One of the variables is Management, a categorical variable that describes different management categories, coded as BF (an abbreviation for biological farming), HF (hobby farming), NM (nature conservation management) and SF (standard farming).
data(dune.env)
summary(dune.env)
#> A1 Moisture Management Use Manure
#> Min. : 2.800 1:7 BF:3 Hayfield:7 0:6
#> 1st Qu.: 3.500 2:4 HF:5 Haypastu:8 1:3
#> Median : 4.200 4:2 NM:6 Pasture :5 2:4
#> Mean : 4.850 5:7 SF:6 3:4
#> 3rd Qu.: 5.725 4:3
#> Max. :11.500
For some plotting methods, it is necessary that the environmental data set is attached:
For the examples, we use the constrained ordination method of redundancy analysis.
The ordiplot object is obtained from the result of via function rda. The ordination is done with a community data set that is transformed by the Hellinger method as recommended in this article.
# script generated by the BiodiversityR GUI from the constrained ordination menu
dune.Hellinger <- disttransform(dune, method='hellinger')
Ordination.model1 <- rda(dune.Hellinger ~ Management,
data=dune.env,
scaling="species")
summary(Ordination.model1)
#>
#> Call:
#> rda(formula = dune.Hellinger ~ Management, data = dune.env, scaling = "species")
#>
#> Partitioning of variance:
#> Inertia Proportion
#> Total 0.5559 1.0000
#> Constrained 0.1667 0.2999
#> Unconstrained 0.3892 0.7001
#>
#> Eigenvalues, and their contribution to the variance
#>
#> Importance of components:
#> RDA1 RDA2 RDA3 PC1 PC2 PC3 PC4
#> Eigenvalue 0.09377 0.05304 0.01988 0.1279 0.05597 0.04351 0.03963
#> Proportion Explained 0.16869 0.09542 0.03575 0.2300 0.10069 0.07827 0.07129
#> Cumulative Proportion 0.16869 0.26411 0.29986 0.5299 0.63054 0.70881 0.78010
#> PC5 PC6 PC7 PC8 PC9 PC10 PC11
#> Eigenvalue 0.03080 0.02120 0.01623 0.01374 0.01138 0.009469 0.007651
#> Proportion Explained 0.05541 0.03814 0.02919 0.02471 0.02047 0.017034 0.013764
#> Cumulative Proportion 0.83551 0.87365 0.90284 0.92755 0.94802 0.965056 0.978820
#> PC12 PC13 PC14 PC15 PC16
#> Eigenvalue 0.003957 0.003005 0.002485 0.001670 0.0006571
#> Proportion Explained 0.007117 0.005406 0.004470 0.003004 0.0011821
#> Cumulative Proportion 0.985937 0.991344 0.995814 0.998818 1.0000000
#>
#> Accumulated constrained eigenvalues
#> Importance of components:
#> RDA1 RDA2 RDA3
#> Eigenvalue 0.09377 0.05304 0.01988
#> Proportion Explained 0.56255 0.31822 0.11923
#> Cumulative Proportion 0.56255 0.88077 1.00000
#>
#> 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: 1.80276
#>
#>
#> Species scores
#>
#> RDA1 RDA2 RDA3 PC1 PC2 PC3
#> Achimill 0.036568 0.127978 0.016399 -0.188968 -0.070111 0.175399
#> Agrostol -0.003429 -0.270611 -0.018285 0.361207 0.005321 -0.014204
#> Airaprae -0.127487 -0.011618 -0.005281 -0.120126 0.103133 0.060671
#> Alopgeni 0.235229 -0.190348 0.027102 0.163865 0.129547 -0.078901
#> Anthodor -0.085534 0.115742 -0.079991 -0.240230 0.117203 0.201272
#> Bellpere 0.033345 0.041954 0.083723 -0.081623 -0.079661 -0.076104
#> Bromhord 0.093316 0.120369 0.065335 -0.058227 -0.047198 0.043097
#> Chenalbu 0.016343 -0.027339 0.008563 0.006896 0.028380 0.004920
#> Cirsarve 0.019792 -0.033109 0.010370 -0.009773 -0.008401 -0.025099
#> Comapalu -0.110016 -0.010026 -0.004558 0.091463 -0.045993 0.026311
#> Eleopalu -0.160847 -0.069849 -0.041472 0.352544 -0.115164 0.126956
#> Elymrepe 0.173955 -0.047695 -0.002020 -0.108520 -0.206622 -0.115376
#> Empenigr -0.047886 -0.004364 -0.001984 -0.029817 0.061163 -0.018612
#> Hyporadi -0.130851 0.036162 0.047182 -0.127285 0.148136 0.027482
#> Juncarti -0.057086 -0.003074 -0.101556 0.265930 -0.061713 -0.009618
#> Juncbufo 0.102986 -0.052188 -0.062882 0.033316 0.152206 -0.038314
#> Lolipere 0.260211 0.153503 0.052002 -0.192751 -0.234269 -0.146056
#> Planlanc -0.018339 0.192788 -0.064491 -0.222661 0.025888 0.091459
#> Poaprat 0.183546 0.093202 0.019068 -0.196708 -0.136129 -0.115068
#> Poatriv 0.377233 -0.046485 -0.039312 -0.019700 -0.010085 0.010877
#> Ranuflam -0.113470 -0.073249 -0.022781 0.249877 -0.046445 0.073742
#> Rumeacet 0.118187 0.070051 -0.198589 -0.056786 0.065922 0.042193
#> Sagiproc 0.076848 -0.060054 0.017557 0.035039 0.222856 -0.152162
#> Salirepe -0.197203 -0.017971 -0.008170 -0.014209 0.015031 -0.104243
#> Scorautu -0.146131 0.134240 0.015763 -0.059562 0.117768 -0.088453
#> Trifprat 0.061485 0.069094 -0.135080 -0.066713 -0.001921 0.069306
#> Trifrepe 0.010076 0.151446 0.031054 0.039480 0.082043 -0.081712
#> Vicilath -0.014220 0.076123 0.084100 -0.026628 0.009338 -0.057100
#> Bracruta -0.057834 0.021502 -0.054847 0.083143 0.070717 -0.115855
#> Callcusp -0.107307 -0.059711 0.009214 0.169330 -0.068364 0.107835
#>
#>
#> Site scores (weighted sums of species scores)
#>
#> RDA1 RDA2 RDA3 PC1 PC2 PC3
#> 1 0.5831 0.1826 0.52222 -0.59008 -0.95148 -0.002365
#> 2 0.4752 0.3624 0.90879 0.05773 -0.27927 -0.053465
#> 3 0.5157 -0.3323 0.54177 -0.17847 -0.28527 -0.368932
#> 4 0.4399 -0.3608 0.65164 -0.15066 -0.12951 -0.386916
#> 5 0.2767 0.7135 -0.68027 -0.32785 -0.13648 0.331094
#> 6 0.1630 0.7846 -0.99354 -0.24794 0.06413 0.280023
#> 7 0.3075 0.7908 -0.69144 -0.29554 0.01116 0.283797
#> 8 0.1248 -0.4586 -0.03533 0.58291 -0.02785 -0.296149
#> 9 0.4391 -0.3568 -0.47994 0.28843 0.08904 -0.598766
#> 10 0.1626 0.8889 0.52513 -0.10281 -0.02239 0.398688
#> 11 -0.1034 0.6728 0.63968 0.04508 0.30166 -0.345223
#> 12 0.2387 -0.6802 -0.39191 0.13764 0.90201 -0.091365
#> 13 0.3253 -0.7706 0.08451 0.12875 0.52983 0.091846
#> 14 -0.6531 -0.4778 0.18077 0.45837 -0.26343 0.275022
#> 15 -0.6814 -0.5359 -0.46471 0.55931 -0.24900 0.020745
#> 16 -0.2720 -1.1011 -0.44903 0.65282 -0.06558 0.757733
#> 17 -0.5168 0.5924 -0.02957 -0.74414 0.25120 0.742874
#> 18 -0.3035 0.6161 0.46554 -0.42677 -0.21642 -0.831729
#> 19 -0.7255 0.1808 0.15665 -0.38151 0.78259 -0.238142
#> 20 -0.7960 -0.7107 -0.46098 0.53475 -0.30493 0.031229
#>
#>
#> Site constraints (linear combinations of constraining variables)
#>
#> RDA1 RDA2 RDA3 PC1 PC2 PC3
#> 1 0.3051 -0.51040 0.15987 -0.59008 -0.95148 -0.002365
#> 2 0.1781 0.64135 0.69120 0.05773 -0.27927 -0.053465
#> 3 0.3051 -0.51040 0.15987 -0.17847 -0.28527 -0.368932
#> 4 0.3051 -0.51040 0.15987 -0.15066 -0.12951 -0.386916
#> 5 0.2622 0.29468 -0.57610 -0.32785 -0.13648 0.331094
#> 6 0.2622 0.29468 -0.57610 -0.24794 0.06413 0.280023
#> 7 0.2622 0.29468 -0.57610 -0.29554 0.01116 0.283797
#> 8 0.2622 0.29468 -0.57610 0.58291 -0.02785 -0.296149
#> 9 0.2622 0.29468 -0.57610 0.28843 0.08904 -0.598766
#> 10 0.1781 0.64135 0.69120 -0.10281 -0.02239 0.398688
#> 11 0.1781 0.64135 0.69120 0.04508 0.30166 -0.345223
#> 12 0.3051 -0.51040 0.15987 0.13764 0.90201 -0.091365
#> 13 0.3051 -0.51040 0.15987 0.12875 0.52983 0.091846
#> 14 -0.6127 -0.05584 -0.02538 0.45837 -0.26343 0.275022
#> 15 -0.6127 -0.05584 -0.02538 0.55931 -0.24900 0.020745
#> 16 0.3051 -0.51040 0.15987 0.65282 -0.06558 0.757733
#> 17 -0.6127 -0.05584 -0.02538 -0.74414 0.25120 0.742874
#> 18 -0.6127 -0.05584 -0.02538 -0.42677 -0.21642 -0.831729
#> 19 -0.6127 -0.05584 -0.02538 -0.38151 0.78259 -0.238142
#> 20 -0.6127 -0.05584 -0.02538 0.53475 -0.30493 0.031229
#>
#>
#> Biplot scores for constraining variables
#>
#> RDA1 RDA2 RDA3 PC1 PC2 PC3
#> ManagementHF 0.3756 0.42205 -0.82512 0 0 0
#> ManagementNM -0.9950 -0.09068 -0.04122 0 0 0
#> ManagementSF 0.4955 -0.82890 0.25963 0 0 0
#>
#>
#> Centroids for factor constraints
#>
#> RDA1 RDA2 RDA3 PC1 PC2 PC3
#> ManagementBF 0.1781 0.64135 0.69120 0 0 0
#> ManagementHF 0.2622 0.29468 -0.57610 0 0 0
#> ManagementNM -0.6127 -0.05584 -0.02538 0 0 0
#> ManagementSF 0.3051 -0.51040 0.15987 0 0 0
For the different graphs, the theme for ggplot2 plotting will be modified as follows, resulting in a more empty plotting canvas than the default for ggplot2.
BioR.theme <- theme(
panel.background = element_blank(),
panel.border = element_blank(),
panel.grid = element_blank(),
axis.line = element_line("gray25"),
text = element_text(size = 12),
axis.text = element_text(size = 10, colour = "gray25"),
axis.title = element_text(size = 14, colour = "gray25"),
legend.title = element_text(size = 14),
legend.text = element_text(size = 14),
legend.key = element_blank())
The ordination result needs to be plotted first via ordiplot.
To plot data via ggplot2, information on the locations of sites (circles in the ordiplot) is obtained via function sites.long. Information on species is extracted by function species.long.
sites.long1 <- sites.long(plot1, env.data=dune.env)
head(sites.long1)
#> A1 Moisture Management Use Manure axis1 axis2 labels
#> 1 2.8 1 SF Haypastu 4 0.5831201 0.1825554 1
#> 2 3.5 1 BF Haypastu 2 0.4751761 0.3624230 2
#> 3 4.3 2 SF Haypastu 4 0.5157071 -0.3323011 3
#> 4 4.2 2 SF Haypastu 4 0.4398671 -0.3607994 4
#> 5 6.3 1 HF Hayfield 2 0.2767400 0.7134658 5
#> 6 4.3 1 HF Haypastu 2 0.1630074 0.7845782 6
Information on the labelling of the axes is obtained with function axis.long. This information is extracted from the ordination model and not the ordiplot, hence it is important to select the same axes via argument ‘choices’. The extracted information includes information on the percentage of explained variation. I suggest that you cross-check with the summary of the redundancy analysis, where information on proportion explained is given.
axis.long1 <- axis.long(Ordination.model1, choices=c(1, 2))
axis.long1
#> axis ggplot label
#> 1 1 xlab.label RDA1 (16.9%)
#> 2 2 ylab.label RDA2 (9.5%)
Function centroids.long obtains data on the centroids of each grouping. It is possible to directly add their results as ‘ordispider diagrams’ in ggplot. Here also another layer is added of superellipses obtained from the ggforce::geom_mark_ellipse function.
plotgg1 <- ggplot() +
geom_vline(xintercept = c(0), color = "grey70", linetype = 2) +
geom_hline(yintercept = c(0), color = "grey70", linetype = 2) +
xlab(axis.long1[1, "label"]) +
ylab(axis.long1[2, "label"]) +
scale_x_continuous(sec.axis = dup_axis(labels=NULL, name=NULL)) +
scale_y_continuous(sec.axis = dup_axis(labels=NULL, name=NULL)) +
geom_mark_ellipse(data=sites.long1,
aes(x=axis1, y=axis2, colour=Management,
fill=after_scale(alpha(colour, 0.2))),
expand=0, size=0.2, show.legend=FALSE) +
geom_segment(data=centroids.long(sites.long1, grouping=Management),
aes(x=axis1c, y=axis2c, xend=axis1, yend=axis2, colour=Management),
size=1, show.legend=FALSE) +
geom_point(data=sites.long1,
aes(x=axis1, y=axis2, colour=Management, shape=Management),
size=5) +
BioR.theme +
ggsci::scale_colour_npg() +
coord_fixed(ratio=1)
plotgg1
Function ordiellipse.long reproduces similar ellipses as vegan::ordiellipse.
It is necessary to re-plot the ordiplot. It is also necessary to give the name of the Management variable explicitly as this was not captured in the ordiellipse result.
plot1 <- ordiplot(Ordination.model1, choices=c(1,2))
Management.ellipses <- ordiellipse(plot1,
groups=Management,
display="sites",
kind="sd")
plotgg2 <- ggplot() +
geom_vline(xintercept = c(0), color = "grey70", linetype = 2) +
geom_hline(yintercept = c(0), color = "grey70", linetype = 2) +
xlab(axis.long1[1, "label"]) +
ylab(axis.long1[2, "label"]) +
scale_x_continuous(sec.axis = dup_axis(labels=NULL, name=NULL)) +
scale_y_continuous(sec.axis = dup_axis(labels=NULL, name=NULL)) +
geom_polygon(data=Management.ellipses.long1,
aes(x=axis1, y=axis2,
colour=Management,
fill=after_scale(alpha(colour, 0.2))),
size=0.2, show.legend=FALSE) +
geom_segment(data=centroids.long(sites.long1,
grouping=Management),
aes(x=axis1c, y=axis2c, xend=axis1, yend=axis2,
colour=Management),
size=1, show.legend=FALSE) +
geom_point(data=sites.long1,
aes(x=axis1, y=axis2, colour=Management, shape=Management),
size=5) +
BioR.theme +
ggsci::scale_colour_npg() +
coord_fixed(ratio=1)
plotgg2
Where the previous example examined patterns for a categorical explanatory variable, here we will explore patterns for the continuous variable A1, documenting the thickness of the A1 horizon.
The vegan package includes a method of adding a smooth surface to an ordination diagram. This method is implemented in function ordisurf.
Function ordisurfgrid.long extracts the data to be plotted with ggplot2
plotgg3 <- ggplot() +
geom_contour_filled(data=A1.grid,
aes(x=x, y=y, z=z)) +
geom_vline(xintercept = c(0), color = "grey70", linetype = 2) +
geom_hline(yintercept = c(0), color = "grey70", linetype = 2) +
xlab(axis.long1[1, "label"]) +
ylab(axis.long1[2, "label"]) +
scale_x_continuous(sec.axis = dup_axis(labels=NULL, name=NULL)) +
scale_y_continuous(sec.axis = dup_axis(labels=NULL, name=NULL)) +
geom_point(data=sites.long1,
aes(x=axis1, y=axis2, shape=Management),
colour="red", size=4) +
BioR.theme +
scale_fill_viridis_d() +
labs(fill="A1") +
coord_fixed(ratio=1)
plotgg3
#> Warning: Removed 209 rows containing non-finite values (stat_contour_filled).
First pvclust results need to be obtained
dune.pv <- pvclust(t(dune.Hellinger),
method.hclust="mcquitty",
method.dist="euclidean",
nboot=1000)
#> Bootstrap (r = 0.5)... Done.
#> Bootstrap (r = 0.6)... Done.
#> Bootstrap (r = 0.7)... Done.
#> Bootstrap (r = 0.8)... Done.
#> Bootstrap (r = 0.9)... Done.
#> Bootstrap (r = 1.0)... Done.
#> Bootstrap (r = 1.1)... Done.
#> Bootstrap (r = 1.2)... Done.
#> Bootstrap (r = 1.3)... Done.
#> Bootstrap (r = 1.4)... Done.
plot(dune.pv)
pvrect(dune.pv, alpha=0.89, pv="au")
Function pvclust.long extracts the data both from the pvclust object and a ordicluster object.
plot1 <- ordiplot(Ordination.model1, choices=c(1,2), scaling='species')
cl.data1 <- ordicluster(plot1,
cluster=as.hclust(dune.pv$hclust))
With a similar methodology as available in ordicluster, higher clustering levels in the hierarchy can be removed via variable prune available from the results for nodes and branches.
plotgg4 <- ggplot() +
geom_vline(xintercept = c(0), color = "grey70", linetype = 2) +
geom_hline(yintercept = c(0), color = "grey70", linetype = 2) +
xlab(axis.long1[1, "label"]) +
ylab(axis.long1[2, "label"]) +
scale_x_continuous(sec.axis = dup_axis(labels=NULL, name=NULL)) +
scale_y_continuous(sec.axis = dup_axis(labels=NULL, name=NULL)) +
geom_segment(data=subset(pvlong$segments,
pvlong$segments$prune > 3),
aes(x=x1, y=y1, xend=x2, yend=y2,
colour=au>=0.89,
size=au),
show.legend=TRUE) +
geom_point(data=subset(pvlong$nodes,
pvlong$nodes$prune > 3),
aes(x=x, y=y,
fill=au>=0.89),
shape=21, size=2, colour="black") +
geom_point(data=sites.long1,
aes(x=axis1, y=axis2, shape=Management),
colour="darkolivegreen4", alpha=0.9, size=5) +
geom_text(data=sites.long1,
aes(x=axis1, y=axis2, label=labels)) +
BioR.theme +
ggsci::scale_colour_npg() +
scale_size(range=c(0.3, 2)) +
scale_shape_manual(values=c(15, 16, 17, 18)) +
guides(shape = guide_legend(override.aes = list(linetype = 0))) +
coord_fixed(ratio=1)
plotgg4
sessionInfo()
#> R version 4.0.2 (2020-06-22)
#> Platform: x86_64-w64-mingw32/x64 (64-bit)
#> Running under: Windows 10 x64 (build 19042)
#>
#> 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] pvclust_2.2-0 dplyr_1.0.2 ggsci_2.9
#> [4] ggrepel_0.8.2 concaveman_1.1.0 ggforce_0.3.2
#> [7] ggplot2_3.3.3 BiodiversityR_2.14-2 vegan_2.5-6
#> [10] lattice_0.20-41 permute_0.9-5
#>
#> loaded via a namespace (and not attached):
#> [1] nlme_3.1-148 RColorBrewer_1.1-2 tools_4.0.2
#> [4] backports_1.1.7 R6_2.4.1 rpart_4.1-15
#> [7] Hmisc_4.4-0 nortest_1.0-4 DBI_1.1.0
#> [10] mgcv_1.8-31 colorspace_1.4-1 nnet_7.3-14
#> [13] withr_2.2.0 tidyselect_1.1.0 gridExtra_2.3
#> [16] curl_4.3 compiler_4.0.2 htmlTable_2.0.0
#> [19] isoband_0.2.1 sandwich_2.5-1 labeling_0.3
#> [22] tcltk2_1.2-11 effects_4.1-4 scales_1.1.1
#> [25] checkmate_2.0.0 RcmdrMisc_2.7-2 stringr_1.4.0
#> [28] digest_0.6.25 foreign_0.8-80 relimp_1.0-5
#> [31] minqa_1.2.4 rmarkdown_2.3 rio_0.5.16
#> [34] base64enc_0.1-3 jpeg_0.1-8.1 pkgconfig_2.0.3
#> [37] htmltools_0.5.1.1 Rcmdr_2.7-2 lme4_1.1-23
#> [40] htmlwidgets_1.5.1 rlang_0.4.8 readxl_1.3.1
#> [43] rstudioapi_0.11 farver_2.0.3 generics_0.1.0
#> [46] zoo_1.8-8 acepack_1.4.1 zip_2.0.4
#> [49] car_3.0-8 magrittr_1.5 Formula_1.2-3
#> [52] Matrix_1.2-18 Rcpp_1.0.7 munsell_0.5.0
#> [55] abind_1.4-5 lifecycle_0.2.0 stringi_1.4.6
#> [58] yaml_2.2.1 carData_3.0-4 MASS_7.3-51.6
#> [61] grid_4.0.2 parallel_4.0.2 forcats_0.5.0
#> [64] crayon_1.3.4 haven_2.3.1 splines_4.0.2
#> [67] hms_0.5.3 knitr_1.28 pillar_1.4.4
#> [70] boot_1.3-25 glue_1.4.1 evaluate_0.14
#> [73] mitools_2.4 latticeExtra_0.6-29 data.table_1.12.8
#> [76] png_0.1-7 vctrs_0.3.4 nloptr_1.2.2.1
#> [79] tweenr_1.0.1 cellranger_1.1.0 gtable_0.3.0
#> [82] purrr_0.3.4 polyclip_1.10-0 xfun_0.15
#> [85] openxlsx_4.1.5 survey_4.0 e1071_1.7-3
#> [88] viridisLite_0.3.0 class_7.3-17 survival_3.1-12
#> [91] tibble_3.0.1 cluster_2.1.0 statmod_1.4.34
#> [94] ellipsis_0.3.1