Analysis conducted by:

Dr. Yoni Gavish & Ron Chen

Data ; Science ; Environment

Statistical and Ecological Solutions

   email: 

   phone: +972-58-4479847

General information on the current analysis:

Project - HaMaarag, State of Nature report, 2023

Group - Birds

Scale - National

Species group - Synanthrope species

Factors - Time, Unit

Response - Abundance

Take home messages:

- Abundance decreased non-significantly (p=0.05805) by 1.457% per year (1.015 fold per year) and by 12.371% (1.141 fold) over the entire period:

  • Predicted value 2012: 12.89 Individuals per sample

  • Predicted value 2021: 11.29 Individuals per sample

A few general comments:

1. The aim of the analysis is to explore the effect of time and/or proximity to settlements and/or proximity to agriculture on the species richness, total abundance, geometric mean of abundance and community structure of birds in Israel

2. Definition of non rare species:

  • For species richness:

    • All species were included in the analysis

  • For total abundance and geometric mean of abundance (satisfy both criteria):

    • Appeared in at least 50 samples over all campaigns

    • At least 10 individuals over all campaigns

  • For community structure (satisfy both criteria):

    • Appeared in at least 10 samples over all campaigns

    • At least 10 individuals over all campaigns

3. The starting model (maximum complexity) for all analyses was kept constant at:

  • Proximity to agricultural + temporal contrasts:
-  Response ~ agriculture * time + unit + cos_td_rad + sin_td_rad + (1|site/point)
  • Proximity to settlements + temporal contrast:
-  Response ~ settlements * time + unit + cos_td_rad + sin_td_rad + (1|site/point)
  • Temporal contrast:
-  Response ~ time + unit + cos_td_rad + sin_td_rad + (1|site/point)

4. Accounting for spatial structure and repeated measures:

  • Each sampling point was visited several times in different years/campaigns

  • Each sampling point is nested within sites

  • Site and point are thus included as a nested random factor in the model

- For community analysis we used (1|point) instead of (1|site/point) as random factor

5. Accounting for activity levels (random factors):

  • “cos_td_rad”, “sin_td_rad” - measure the date difference from 21 June, as the cos or sin of the difference in days in radians

  • Wind, Precipitation, Temperature and Cloud are not included in the analysis although they may affect activity level since the data is too sparse (many NAs)

  • The hour of sampling was not recorded in T0 - the two related variables c(“cos_hsun”, “sin_hsun”) are ignored in the analysis

6. Model complexity:

  • This analysis was restricted to monotonic increase/decrease

  • More flexible models that allow some curvature may unveil different patterns and trends

1. Prepare the working environment

2. Load data and source functions

3. Arrange the Traits table

  • Add to the traits table the species code, order, family, IUCN category and group, and residency status in Israel

  • Table shows number of NAs, class, mean value (for integers and numeric) and number of factor levels (for character and factors)

Column NAs Class Mean Levels
SciName 0 factor NA 189
SPECIES_CODE 0 character NA 189
HebName 0 factor NA 189
PRIMARY_COM_NAME 0 factor NA 189
Order 0 character NA 20
Family 0 character NA 53
Synanthrope 0 integer 0.180 NA
AlienInvasive 0 integer 0.021 NA
NativeInvasive 0 integer 0.021 NA
synanthrope_or_invasive 0 logical NA NA
batha_specialist 70 logical NA NA
oblig_insectivore 32 logical NA NA
IUCN_AllCat 43 character NA 9
IUCN 0 character NA 3
Resident 0 integer 0.434 NA
Resident_IL 0 factor NA 2

4. Arrange and filter the Plots table

  • Table shows number of NAs, class, mean value (for integers and numeric) and number of factor levels (for character and factors)
Column NAs Class Mean Levels
unit 0 factor NA 8
subunit 926 factor NA 4
site 0 factor NA 70
year 0 integer 2016.665 NA
year_ct 0 numeric 4.665 NA
settlements 791 factor NA 3
agriculture 1169 factor NA 3
habitat 1145 factor NA 5
dunes 1559 factor NA 3
land_use 1455 factor NA 4
point_name 0 factor NA 574
date 0 POSIXct NA NA
date 0 POSIXt NA NA
datetime 441 POSIXct NA NA
datetime 441 POSIXt NA NA
td_sc 0 numeric 0.129 NA
td_rad 0 numeric -0.783 NA
cos_td_rad 0 numeric 0.629 NA
sin_td_rad 0 numeric -0.635 NA
timediff_Jun21 0 difftime NA NA
monitors_name 1023 factor NA 14
wind 1449 ordered NA 5
wind 1449 factor NA 5
precipitation 1390 ordered NA 3
precipitation 1390 factor NA 3
temperature 1433 ordered NA 5
temperature 1433 factor NA 5
clouds 1519 ordered NA 5
clouds 1519 factor NA 5
h_from_sunrise 441 difftime NA NA
cos_hsun 441 numeric 0.817 NA
sin_hsun 441 numeric 0.442 NA
pilot 0 logical NA NA
campaign 0 character NA 5
Unit_Label 0 character NA 8
Region 0 character NA 3

5. Clean Data and identify non-rare species

  • Filter columns less relevant to the analysis

  • Remove all observations with 0 abundance (legacy from former filtering)

  • Add features from the trait table

  • Keep only Synanthrope species

  • Extract general properties of each column in Data : BasicInfo.R

  • Identify species that meet the two criteria for non-rare : Get_InciAbunRare.R

  • Table shows number of NAs, class, mean value (for integers and numeric) and number of factor levels (for character and factors)

Column NAs Class Mean Levels
SciName 0 character NA 26
year 0 integer 2016.874 NA
point_name 0 character NA 543
campaign 0 character NA 5
unit 0 character NA 8
site 0 character NA 70
year_ct 0 integer 4.874 NA
count_under_250 0 integer 3.045 NA
agriculture 6002 character NA 3
settlements 3388 character NA 3
cos_td_rad 0 numeric 0.644 NA
sin_td_rad 0 numeric -0.644 NA
SPECIES_CODE 0 data.table NA NA
SPECIES_CODE 0 data.frame NA NA
HebName 0 data.table NA NA
HebName 0 data.frame NA NA
PRIMARY_COM_NAME 0 data.table NA NA
PRIMARY_COM_NAME 0 data.frame NA NA
Order 0 data.table NA NA
Order 0 data.frame NA NA
Family 0 data.table NA NA
Family 0 data.frame NA NA
Synanthrope 0 data.table NA NA
Synanthrope 0 data.frame NA NA
AlienInvasive 0 data.table NA NA
AlienInvasive 0 data.frame NA NA
NativeInvasive 0 data.table NA NA
NativeInvasive 0 data.frame NA NA
synanthrope_or_invasive 0 data.table NA NA
synanthrope_or_invasive 0 data.frame NA NA
batha_specialist 526 data.table NA NA
batha_specialist 526 data.frame NA NA
oblig_insectivore 0 data.table NA NA
oblig_insectivore 0 data.frame NA NA
IUCN_AllCat 0 data.table NA NA
IUCN_AllCat 0 data.frame NA NA
IUCN 0 data.table NA NA
IUCN 0 data.frame NA NA
Resident_IL 0 data.table NA NA
Resident_IL 0 data.frame NA NA

A total of 21 species satisfied the non-rare criteria, including:

Streptopelia decaocto ; Turdus merula ; Pycnonotus xanthopygos ; Parus major ; Prinia gracilis ; Passer domesticus ; Corvus cornix ; Cinnyris osea ; Garrulus glandarius ; Spilopelia senegalensis ; Columba livia ; Streptopelia turtur ; Falco tinnunculus ; Acridotheres tristis ; Cecropis daurica ; Dendrocopos syriacus ; Psittacula krameri ; Corvus monedula ; Vanellus spinosus ; Upupa epops ; Hirundo rustica

6. Create a dataframe with the relevant variables for analysis

  • For Richness, we use all species.

  • For abundance, mean abundance (geometric) and community structure we filter rare species.

  • Create a data-table with a row for each combination of ‘year’ and ‘point_name’

    • First two columns are ‘year’ and ‘point_name’

    • Column for each species

  • Add empty plots (plots with no observations)

  • Filter Rare species

  • Calculate total abundance

  • Calculate the geometric mean abundance (GMA)

    • For each row

    • For species with at least one individual (Note)

  • Add the relevant predictor variables

Note:

when calculating a geometric mean (GMA), any species with 0 abundance would result with a geometric mean of 0 (since the geometric mean relies on the product). Since most samples include at least one species with 0 abundance, we based the geometric mean only on the species with at least one individual. For empty plots- NaN is returned and these rows should be removed when modeling GMA.

  • Table shows number of NAs, class, mean value (for integers and numeric) and number of factor levels (for character and factors)
Column NAs Class Mean Levels
year 0 numeric 2016.667 NA
point_name 0 character NA 574
Streptopelia decaocto 0 numeric 2.747 NA
Turdus merula 0 numeric 1.779 NA
Pycnonotus xanthopygos 0 numeric 1.089 NA
Parus major 0 numeric 0.835 NA
Prinia gracilis 0 numeric 0.715 NA
Passer domesticus 0 numeric 2.560 NA
Corvus cornix 0 numeric 0.572 NA
Cinnyris osea 0 numeric 0.496 NA
Garrulus glandarius 0 numeric 0.582 NA
Spilopelia senegalensis 0 numeric 0.523 NA
Columba livia 0 numeric 1.130 NA
Streptopelia turtur 0 numeric 0.254 NA
Falco tinnunculus 0 numeric 0.142 NA
Acridotheres tristis 0 numeric 0.291 NA
Cecropis daurica 0 numeric 0.427 NA
Dendrocopos syriacus 0 numeric 0.077 NA
Psittacula krameri 0 numeric 0.105 NA
Corvus monedula 0 numeric 0.278 NA
Vanellus spinosus 0 numeric 0.119 NA
Upupa epops 0 numeric 0.052 NA
Hirundo rustica 0 numeric 0.144 NA
TotalAbun 0 numeric 14.918 NA
GMA 161 numeric 2.643 NA
site 0 factor NA 70
campaign 0 character NA 5
year_ct 0 numeric 4.667 NA
agriculture 1168 factor NA 3
settlements 791 factor NA 3
unit 0 factor NA 8
Unit_Label 0 character NA 8
Region 0 character NA 3
cos_td_rad 0 numeric 0.629 NA
sin_td_rad 0 numeric -0.636 NA
year_ct_Sq 0 numeric 29.703 NA
Response 0 numeric 14.918 NA

The data contains 161 empty plots.

Of these, 36 plots were empty throughout the entire monitoring period.

These plots represent 67 plot sampling events.

Of these plots, 16 plots span more than a single year.

7. Explore the Abundance data frame

  • Visualize the relationship of the response to various explanatory variables

Main conclusions:

  • Considerable variability in total abundance between units, with highest values in the Mediterranean maquis unit and lowest in the inland sands unit

  • Considerable variability between and within sites

Notes:

  • unit Batha vegetation refers to Herbaceous and Dwarf-Shrub Vegetation

  • unit Transition Zone refers to Mediterranean-Desert Transition Zone

  • unit Loess Areas refers to Loess Covered Areas in the Northern Negev

samples with abundance > 110 (not shown in figures)
761 762 763 875 1389 1408
year 2017 2017 2017 2018 2021 2021
point_name Givat Yeshayahu Near 1 Givat Yeshayahu Near 2 Givat Yeshayahu Near 3 Givot Bar Bedouin Agriculture 2 Kerem Maharal Near 3 Nir Etzion Far 11
Unit_Label Mediterranean Maquis Mediterranean Maquis Mediterranean Maquis Loess Areas Mediterranean Maquis Mediterranean Maquis
Region Mediterranean Mediterranean Mediterranean Transition Mediterranean Mediterranean
site Givat Yeshayahu Givat Yeshayahu Givat Yeshayahu Givot Bar Kerem Maharal Nir Etzion
TotalAbun 119 125 197 115 228 112
Number of samples
Year Arid South Negev Highlands Inland Sands Loess Areas Transition Zone Planted Conifer Forest Mediterranean Maquis Batha Vegetation
2012 30 30 0 24 30 0 89 0
2014 30 40 0 27 29 45 0 60
2015 0 0 0 0 0 45 89 0
2016 30 40 0 29 28 0 0 60
2017 0 0 5 0 0 45 88 0
2018 30 39 0 30 29 0 0 60
2019 0 0 5 0 0 45 89 0
2020 30 39 0 30 29 0 0 60
2021 0 0 5 0 0 45 89 0

8. Model Abundance - Idenitify Family and Fixed factors

We fit generalized linear mixed-effects model (glmmTMB) with:

  • The response is: TotalAbun (internally, “Response”)

  • The fixed factors are: year_ct ; unit ; cos_td_rad ; sin_td_rad

  • The groping factors are: site ; point_name - with the later nested within the first (each site contains several plots)

Objectives:

  • Identify the best combination of fixed factors

  • Decide on the family (poisson or negative binomial)

Workflow:

  • Identify the best model for Poisson:

    • Fit the most complex model with glmmTMB

    • Simplify the model with MuMIn::dredge(rank = “AIC”) - removal of fixed terms

    • Select the best model

  • Identify the best model for Negative binomial:

    • Fit the most complex model with glmmTMB

    • Simplify the model with MuMIn::dredge(rank = “AIC”) - removal of fixed terms

    • Select the best model

  • Select the optimal model in terms of family and fixed factors.

    • Use model selection to compare the best poisson and best negative binomial

    • Diagnose the best model with glmmTMB:diagnose and look at the sigma (for NB)

Simplification of the poisson models

  • Best model kept all terms

  • Second best model DeltaAIC > 2

  • Keep first model

cond((Int)) disp((Int)) cond(cos_td_rad) cond(sin_td_rad) cond(Unit_Label) cond(year_ct) df logLik AIC delta weight
16 0.570
1.013 -1.126
-0.017 13 -7411.108 14848.22 0.000 9.999928e-01
8 0.560
1.025 -1.046
NA 12 -7423.948 14871.90 23.679 7.214279e-06
12 0.964
1.024 -1.132 NA -0.016 6 -7446.260 14904.52 56.304 5.938473e-13
4 0.943
1.035 -1.059 NA NA 5 -7457.110 14924.22 76.004 3.132167e-17
15 1.805
NA -0.338
-0.021 12 -7607.051 15238.10 389.886 2.174363e-85
7 1.808
NA -0.233
NA 11 -7625.206 15272.41 424.196 7.710186e-93
11 2.135
NA -0.335 NA -0.019 5 -7645.516 15301.03 452.816 4.702645e-99
3 2.122
NA -0.237 NA NA 4 -7661.403 15330.81 482.589 1.611158e-105
14 1.855
0.110 NA
0.008 12 -7663.366 15350.73 502.515 7.590160e-110
6 1.904
0.073 NA
NA 11 -7666.079 15354.16 505.942 1.368211e-110
5 1.956
NA NA
NA 10 -7669.876 15359.75 511.537 8.340557e-112
13 1.951
NA NA
0.002 11 -7669.703 15361.41 513.190 3.650166e-112
10 2.154
0.114 NA NA 0.009 5 -7700.423 15410.84 562.629 6.708107e-123
2 2.217
0.070 NA NA NA 4 -7704.267 15416.53 568.317 3.904441e-124
1 2.263
NA NA NA NA 3 -7707.776 15421.55 573.336 3.173890e-125
9 2.252
NA NA NA 0.003 4 -7707.258 15422.51 574.299 1.961412e-125

Simplification of the negative binomial models

  • Best model kept all the terms

  • Second model DeltaAIC > 2

  • Keep first model

cond((Int)) disp((Int)) cond(cos_td_rad) cond(sin_td_rad) cond(Unit_Label) cond(year_ct) df logLik AIC delta weight
16 0.845
1.004 -1.010
-0.029 14 -5547.217 11122.43 0.000 9.760425e-01
8 0.797
1.037 -0.895
NA 13 -5551.924 11129.85 7.414 2.395746e-02
12 1.175
1.043 -1.006 NA -0.025 7 -5576.316 11166.63 44.198 2.465886e-10
4 1.116
1.073 -0.908 NA NA 6 -5579.826 11171.65 49.219 2.003412e-11
15 2.041
NA -0.296
-0.035 13 -5576.651 11179.30 56.868 4.372321e-13
6 1.947
0.231 NA
NA 12 -5580.500 11185.00 62.566 2.532097e-14
14 1.987
0.203 NA
-0.006 13 -5580.295 11186.59 64.157 1.142495e-14
13 2.167
NA NA
-0.016 12 -5582.921 11189.84 67.410 2.247244e-15
7 2.029
NA -0.128
NA 12 -5583.277 11190.55 68.121 1.574878e-15
5 2.112
NA NA
NA 11 -5584.964 11191.93 69.494 7.926085e-16
2 2.209
0.259 NA NA NA 5 -5608.441 11226.88 104.449 2.035717e-23
11 2.348
NA -0.259 NA -0.032 6 -5607.456 11226.91 104.478 2.006022e-23
10 2.230
0.245 NA NA -0.003 6 -5608.389 11228.78 106.346 7.885083e-24
9 2.441
NA NA NA -0.015 5 -5612.163 11234.33 111.893 4.922017e-25
3 2.316
NA -0.106 NA NA 5 -5612.896 11235.79 113.359 2.365806e-25
1 2.380
NA NA NA NA 4 -5614.021 11236.04 113.609 2.086967e-25

Model selection

  • The Negative binomial model (AIC weight = 1) outperformed the Poisson model
cond((Int)) disp((Int)) cond(cos_td_rad) cond(sin_td_rad) cond(Unit_Label) cond(year_ct) family df logLik AIC delta weight
Best_NB 0.845
1.004 -1.010
-0.029 nbinom2(log) 14 -5547.217 11122.43 0.000 1
Best_Po 0.570
1.013 -1.126
-0.017 poisson(log) 13 -7411.108 14848.22 3725.783 0

Diagnostics and sigma

  • Diagnostic of best negative binomial model reveals minor issues
## Unusually large Z-statistics (|x|>5):
## 
##                     cos_td_rad                     sin_td_rad 
##                       7.856566                      -8.341510 
## Unit_LabelMediterranean Maquis                  d~(Intercept) 
##                       5.823163                      16.098282 
##      theta_1|point_name:site.1 
##                      -6.286655 
## 
## Large Z-statistics (estimate/std err) suggest a *possible* failure of
## the Wald approximation - often also associated with parameters that are
## at or near the edge of their range (e.g. random-effects standard
## deviations approaching 0).  (Alternately, they may simply represent
## very well-estimated parameters; intercepts of non-centered models may
## fall in this category.) While the Wald p-values and standard errors
## listed in summary() may be unreliable, profile confidence intervals
## (see ?confint.glmmTMB) and likelihood ratio test p-values derived by
## comparing models (e.g. ?drop1) are probably still OK.  (Note that the
## LRT is conservative when the null value is on the boundary, e.g. a
## variance or zero-inflation value of 0 (Self and Liang 1987; Stram and
## Lee 1994; Goldman and Whelan 2000); in simple cases these p-values are
## approximately twice as large as they should be.)
  • The dispersion parameter of the best negative binomial model was >1, suggesting the negative binomial model is better
## [1] 2.555403

Main conclusions:

  • Continue with the best negative binomial model
## Formula:          
## Response ~ cos_td_rad + sin_td_rad + Unit_Label + year_ct + (1 |  
##     site/point_name)
## Data: MainData
##       AIC       BIC    logLik  df.resid 
## 11122.433 11197.250 -5547.217      1533 
## Random-effects (co)variances:
## 
## Conditional model:
##  Groups          Name        Std.Dev.
##  point_name:site (Intercept) 0.7187  
##  site            (Intercept) 0.1378  
## 
## Number of obs: 1547 / Conditional model: point_name:site, 558; site, 70
## 
## Dispersion parameter for nbinom2 family (): 2.56 
## 
## Fixed Effects:
## 
## Conditional model:
##                      (Intercept)                        cos_td_rad  
##                          0.84459                           1.00418  
##                       sin_td_rad         Unit_LabelNegev Highlands  
##                         -1.00956                          -0.08824  
##           Unit_LabelInland Sands             Unit_LabelLoess Areas  
##                         -0.57540                           0.12194  
##        Unit_LabelTransition Zone  Unit_LabelPlanted Conifer Forest  
##                          0.32653                           0.34132  
##   Unit_LabelMediterranean Maquis        Unit_LabelBatha Vegetation  
##                          0.92402                           0.54454  
##                          year_ct  
##                         -0.02876
## glmmTMB::glmmTMB(formula = Response ~ cos_td_rad + sin_td_rad + 
##     Unit_Label + year_ct + (1 | site/point_name), data = MainData, 
##     family = "nbinom2", ziformula = ~0, dispformula = ~1, na.action = na.fail)

9. Model Abundance - Identfiy best family variant of the negative binomial model

The glmmTMB function support multiple variants of many families and various link function

  • For poisson: ‘poisson’ vs. ‘compois’

  • For negative binomial: ‘nbinom1’ vs. ‘nbinom2’

  • For Gamma: Different link functions (log vs. identity vs. inverse)

Objectives:

  • Identify the best model family/link variant based on the selected fixed factors

Workflow:

  • Fit an alternative model with family == ‘nbinom1’

  • Simplify the alternative model with MuMIn::dredge()

  • Select the best alternative model

  • Use model selection to compare the best original and alternative nodels (MuMIn::model.sel)

  • Select the optimal model in terms of family variant

  • Diagnose the best model with glmmTMB:diagnose

Simplification of the negative binomial models (nbinom1)

  • Best model kept all terms

  • Second best model DeltaAIC > 2

  • Keep the first

cond((Int)) disp((Int)) cond(cos_td_rad) cond(sin_td_rad) cond(Unit_Label) cond(year_ct) df logLik AIC delta weight
16 1.091
0.915 -0.874
-0.015 14 -5537.805 11103.61 0.000 6.883004e-01
8 1.082
0.925 -0.804
NA 13 -5539.597 11105.19 1.584 3.116996e-01
12 1.333
0.964 -0.871 NA -0.012 7 -5570.795 11155.59 51.980 3.551256e-12
4 1.313
0.973 -0.817 NA NA 6 -5571.897 11155.79 52.185 3.205972e-12
6 2.149
0.184 NA
NA 12 -5567.266 11158.53 54.923 8.154498e-13
14 2.108
0.213 NA
0.006 13 -5566.978 11159.95 56.346 4.003827e-13
15 2.209
NA -0.192
-0.018 13 -5567.537 11161.07 57.464 2.288806e-13
7 2.211
NA -0.096
NA 12 -5570.128 11164.25 60.646 4.663647e-14
5 2.277
NA NA
NA 11 -5571.506 11165.01 61.403 3.193567e-14
13 2.296
NA NA
-0.006 12 -5571.145 11166.29 62.681 1.685622e-14
2 2.322
0.217 NA NA NA 5 -5600.782 11211.56 107.954 2.488698e-24
10 2.260
0.259 NA NA 0.008 6 -5600.190 11212.38 108.771 1.653825e-24
11 2.427
NA -0.164 NA -0.016 6 -5603.507 11219.01 115.405 5.996606e-26
1 2.462
NA NA NA NA 4 -5606.468 11220.94 117.326 2.294517e-26
3 2.412
NA -0.081 NA NA 5 -5605.490 11220.98 117.370 2.245302e-26
9 2.484
NA NA NA -0.005 5 -5606.147 11222.29 118.684 1.163969e-26

Model selection

  • The ‘nbinom1’ model outperformed the ‘nbinom2’ model.
cond((Int)) disp((Int)) cond(cos_td_rad) cond(sin_td_rad) cond(Unit_Label) cond(year_ct) family df logLik AIC delta weight
Best_Opt1 1.091
0.915 -0.874
-0.015 nbinom1(log) 14 -5537.805 11103.61 0.000 9.999183e-01
BestModel 0.845
1.004 -1.010
-0.029 nbinom2(log) 14 -5547.217 11122.43 18.824 8.173733e-05

Diagnostics and sigma

  • Diagnostic of best model reveals minor issues
## Unusually large Z-statistics (|x|>5):
## 
##                    (Intercept)                     cos_td_rad 
##                       5.735627                       7.404251 
##                     sin_td_rad Unit_LabelMediterranean Maquis 
##                      -7.381818                       6.194520 
##                  d~(Intercept)      theta_1|point_name:site.1 
##                      32.165651                     -10.517364 
## 
## Large Z-statistics (estimate/std err) suggest a *possible* failure of
## the Wald approximation - often also associated with parameters that are
## at or near the edge of their range (e.g. random-effects standard
## deviations approaching 0).  (Alternately, they may simply represent
## very well-estimated parameters; intercepts of non-centered models may
## fall in this category.) While the Wald p-values and standard errors
## listed in summary() may be unreliable, profile confidence intervals
## (see ?confint.glmmTMB) and likelihood ratio test p-values derived by
## comparing models (e.g. ?drop1) are probably still OK.  (Note that the
## LRT is conservative when the null value is on the boundary, e.g. a
## variance or zero-inflation value of 0 (Self and Liang 1987; Stram and
## Lee 1994; Goldman and Whelan 2000); in simple cases these p-values are
## approximately twice as large as they should be.)

Main conclusions:

  • Continue with the ‘nbinom2’ model
## Formula:          
## Response ~ cos_td_rad + sin_td_rad + Unit_Label + year_ct + (1 |  
##     site/point_name)
## Data: MainData
##       AIC       BIC    logLik  df.resid 
## 11103.610 11178.427 -5537.805      1533 
## Random-effects (co)variances:
## 
## Conditional model:
##  Groups          Name        Std.Dev.
##  point_name:site (Intercept) 0.583207
##  site            (Intercept) 0.000541
## 
## Number of obs: 1547 / Conditional model: point_name:site, 558; site, 70
## 
## Dispersion parameter for nbinom1 family (): 6.31 
## 
## Fixed Effects:
## 
## Conditional model:
##                      (Intercept)                        cos_td_rad  
##                          1.09116                           0.91463  
##                       sin_td_rad         Unit_LabelNegev Highlands  
##                         -0.87442                          -0.12401  
##           Unit_LabelInland Sands             Unit_LabelLoess Areas  
##                         -0.63510                           0.09470  
##        Unit_LabelTransition Zone  Unit_LabelPlanted Conifer Forest  
##                          0.29626                           0.24469  
##   Unit_LabelMediterranean Maquis        Unit_LabelBatha Vegetation  
##                          0.75946                           0.41806  
##                          year_ct  
##                         -0.01467
## glmmTMB::glmmTMB(formula = Response ~ cos_td_rad + sin_td_rad + 
##     Unit_Label + year_ct + (1 | site/point_name), data = MainData, 
##     family = "nbinom1", ziformula = ~0, dispformula = ~1, na.action = na.fail)

10. Diagnose the best model with simulated residuals (DHARMa)

General info

To diagnose the model, we simulate residuals with the DHARMa package. For more details why DHARMa should be used with GLMM, see:

  1. https://cran.r-project.org/web/packages/DHARMa/vignettes/DHARMa.html

To interpret the residuals:

  • Scaled residual value of 0.5 means that half of the simulated data are higher than the observed value, and half of them lower

  • Scaled residual value of 0.99 would mean that nearly all simulated data are lower than the observed value.

  • The minimum/maximum values for the residuals are 0 and 1

  • For a correctly specified model we would expect asymptotically

    • A uniform (flat) distribution of the scaled residuals

    • Uniformity in y direction if we plot against any predictor

Dispersion:

  • Compare the variance of the observed raw residuals (red line) against the variance of the simulated residuals (histogram)

    • Dispersion > 1 indicates overdispersion

    • Dispersion < 1 indicates underdispersion

  • Observed variance (red line) did not differ from the simulated variances (histogram)

  • Dispersion > 1 but n.s.

  • No evidence for over or under dispersion

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 1.2008, p-value = 0.1
## alternative hypothesis: two.sided

Outliers:

  • Understanding outliers:

    • When all simulated values for the observation are higher or smaller than the observation itself

    • Tests if the number of Outliers is larger or smaller than expected

  • There are 7 outliers in the simulated residuals

  • Test is n.s.

  • No evidence for excessive number of outliers

## 
##  DHARMa outlier test based on exact binomial test with approximate
##  expectations
## 
## data:  SimulatedRes
## outliers at both margin(s) = 10, observations = 1547, p-value = 0.1504
## alternative hypothesis: true probability of success is not equal to 0.003992016
## 95 percent confidence interval:
##  0.003104022 0.011855527
## sample estimates:
## frequency of outliers (expected: 0.00399201596806387 ) 
##                                            0.006464124

Zero inflation:

  • The expected distribution of zeros (histogram) and the observed value (red line)

    • ratioObsSim < 1 : the observed data has less zeros than expected

    • ratioObsSim > 1 : the observed data has more zeros than expected

  • No evidence for zero inflation:

    • Red line within histogram

    • Test is n.s

## 
##  DHARMa zero-inflation test via comparison to expected zeros with
##  simulation under H0 = fitted model
## 
## data:  simulationOutput
## ratioObsSim = 1.0586, p-value = 0.58
## alternative hypothesis: two.sided

Quantiles:

  • For numeric:

    • Solid black / red lines: qunatile regressions (0.25, 0.5, 0.75) of the observed residuals against predicted value

      • Significant regressions in red

    • Points - simulated residuals

    • Dashed horizontal lines - the expected value for the 0.25, 0.5, and 0.75 quantiles (simulated residuals have a uniform distribution)

  • For categorical:

    • Within-group uniformity of categorical predictors [asymptotic one-sample Kolmogorov-Smirnov test]

    • Gray: the residual within the group are distributed according to the model assumptions (dashed lines)

    • Red: the residual within the group deviate from the model assumptions (dashed lines)

    • homogeneity of variances of categorical predictors [Levene’s Test]

    • Non-significant (gray): group variances are constant

    • Significant (red): group variances are not constant

  • Results:

    • The 0.25 and 0.5 quantile regressions of the observed residuals against predicted values were significant, but lines are mostly horizontal and deviation occurs mainly in the extremes

    • The 0.75 quantile regression of the observed residuals against predicted values was non-significant

## 
##  Test for location of quantiles via qgam
## 
## data:  simulationOutput
## p-value = 0.06703
## alternative hypothesis: both

year:

  • See categorical above for details (years treated as categorical)

  • Residuals in 2017 deviated from expected

  • Residuals in all other years did not deviate from expected

  • Variances of different years differ significantly from one another

Unit:

  • See categorical above for details

  • Residuals of 1 out of 8 units deviated from expected:

    • Negev Highlands

  • Variances of the units differ significantly from one another

## $uniformity
## $uniformity$details
## catPred: Arid South
## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  dd[x, ]
## D = 0.097438, p-value = 0.1159
## alternative hypothesis: two-sided
## 
## ------------------------------------------------------------ 
## catPred: Negev Highlands
## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  dd[x, ]
## D = 0.11969, p-value = 0.00916
## alternative hypothesis: two-sided
## 
## ------------------------------------------------------------ 
## catPred: Inland Sands
## 
##  Exact one-sample Kolmogorov-Smirnov test
## 
## data:  dd[x, ]
## D = 0.1662, p-value = 0.7425
## alternative hypothesis: two-sided
## 
## ------------------------------------------------------------ 
## catPred: Loess Areas
## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  dd[x, ]
## D = 0.067955, p-value = 0.5376
## alternative hypothesis: two-sided
## 
## ------------------------------------------------------------ 
## catPred: Transition Zone
## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  dd[x, ]
## D = 0.15062, p-value = 0.002779
## alternative hypothesis: two-sided
## 
## ------------------------------------------------------------ 
## catPred: Planted Conifer Forest
## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  dd[x, ]
## D = 0.15011, p-value = 7.892e-05
## alternative hypothesis: two-sided
## 
## ------------------------------------------------------------ 
## catPred: Mediterranean Maquis
## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  dd[x, ]
## D = 0.075977, p-value = 0.01188
## alternative hypothesis: two-sided
## 
## ------------------------------------------------------------ 
## catPred: Batha Vegetation
## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  dd[x, ]
## D = 0.063859, p-value = 0.2817
## alternative hypothesis: two-sided
## 
## 
## $uniformity$p.value
## [1] 1.158670e-01 9.160053e-03 7.425363e-01 5.375670e-01 2.779292e-03
## [6] 7.891615e-05 1.188075e-02 2.816840e-01
## 
## $uniformity$p.value.cor
## [1] 0.4634679459 0.0549603186 1.0000000000 1.0000000000 0.0194550409
## [6] 0.0006313292 0.0594037734 0.8450521044
## 
## 
## $homogeneity
## Levene's Test for Homogeneity of Variance (center = median)
##         Df F value    Pr(>F)    
## group    7  24.437 < 2.2e-16 ***
##       1539                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Uniformity of residuals:

  • Kolmogorov-Smirnov Test (KS) of the uniformity of simulated residuals

    • Residuals did not deviate from expected

  • Dispersion test n.s. (see above)

  • Outliers test n.s. (see above)

  • qq plot -

    • Quantiles of the observed residuals vs. qunatiles of the simulated residuals

    • The expected value is the 1:1 line (red)

  • qqline is close to the expected line

## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  simulationOutput$scaledResiduals
## D = 0.033974, p-value = 0.05624
## alternative hypothesis: two-sided

Main conclusions:

  • The 0.25 and 0.5 quantile regressions of the observed residuals against predicted values were significant, but lines are mostly horizontal and deviation occurs mainly in the extremes

  • Some issues with uniformity of residuals within groups (overall, years, units)

  • Other than that, no other issues

  • Continue to compute marginal effects

11. Compute marginal effects for figures

  • Marginal effects are extracted with ggeffects:ggeffect() that calls effects::Effect() internally

  • Marginal effects are only calculated and plotted for terms that appear in the final model:

    • year

    • unit

  • When calculating change for years:

    • FoldIncrease: Value_year / Value_year_0

    • PercentIncrease: 100 * [(Value_year - Value_year_0) / Value_year_0]

    • FoldDecline: Value_year_0 / Value_year

    • PercentDecline: 100 * [(Value_year_0 - Value_year) / Value_year_0]

    • Delta: Value_year - Value_year_0

  • When calculating change for Agriculture / Settlements / Units:

    • Fold: Value_level / Value_min

    • PercentDiffernce: 100 * [(Value_level - Value_min) / Value_min]

    • Delta: Value_level - Value_min

13. Figures for years

14. Figures for unit

15. Calculate the folds and percent change

  • Calculate the folds and percent change - contrast (agriculture / settlements):

    • If Near > Far (coefficient > 0):

      • Contrast_Type: “Increase”

      • Contrast_Fold: exp(coefficient)

      • Contrast_Percent: 100 * (exp(coefficient) - 1)

    • If Near < Far (coefficient > 0):

      • Contrast_Type: “Decrease”

      • Contrast_Fold: 1 / exp(coefficient)

      • Contrast_Percent: 100 * (1 - exp(coefficient))

    • If Near = Far (coefficient = 0):

      • Contrast_Type: “NoChange”

      • Contrast_Fold: 1

      • Contrast_Percent: 0

  • Calculate the folds and percent change - temporal trend (years):

    • Num_Years: Number of years

    • If LastYear > FirstYear (coefficient > 0):

      • Year_Type: “Increase”

      • YearFold_1Year: exp(coefficient)

      • YearPercent_1Year: 100 * (exp(coefficient) - 1)

      • YearFold_NumYear: exp(coefficient * Num_Years)

      • YearPercent_NumYear: 100 * (exp(coefficien * Num_Years) - 1)

    • If LastYear < FirstYear (coefficient > 0):

      • Year_Type: “Decrease”

      • YearFold_1Year: 1 / exp(coefficient)

      • YearPercent_1Year: 100 * (1 - exp(coefficient))

      • YearFold_NumYear: 1 / exp(coefficient * Num_Years)

      • YearPercent_NumYear: 100 * (1 - exp(coefficient * Num_Years))

    • If LastYear = FirstYear (coefficient = 0):

      • Year_Type: “NoChange”

      • YearFold_1Year: 1

      • YearPercent_1Year: 0

      • YearFold_NumYear: 1

      • YearPercent_NumYear: 0

16. Extract information and write to files

  • Summary(BestModel) is written into:

    • SoN23_BirdsNational_Time_Abundance_Synanthrope_ModSummary.txt

  • Fixed effects from the best model are written into:

    • SoN23_BirdsNational_Time_Abundance_Synanthrope_Coeff.csv

  • General information on the model is written into:

    • SoN23_BirdsNational_Time_Abundance_Synanthrope_ModelInfo.csv

  • Marginal effects of years / units are written to:

    • SoN23_BirdsNational_Time_Abundance_Synanthrope_MargEffe_Years.csv

    • SoN23_BirdsNational_Time_Abundance_Synanthrope_MargEffe_Units.csv

  • A figure is written (png, svg, pdf) for each of the marginal effect figures plotted in this document

    • Names are the predictors (interactions) + indication if the points are observations, partial residuals or none:

      • “AgriUnit_Observation” : agriculture * unit ; observations

      • “SettYearUnit_PartResi” : settlements * year * unit ; partial residuals

      • “AgriYear_None” : agriculture * year, without points

    • Names with “CapY_100” means that the maximal value of Y was truncated at 100 and some samples are not plotted

    • Names with “_Heb” mean that all text elements are in Hebrew

17. Summary of main findings

Data structure:

  • The data contains 1547 samples in which 21 species passed as non-rare (and thus included in the analysis)

  • Total abundances per sample ranges from 0 to 228 individuals

Visual inspection of Response:

  • No apparent trends in total abundance with time

  • Considerable variability in total abundance between units, with highest values in the Mediterranean maquis unit and lowest in the inland sands unit

Model fitting (glmmTMB):

  • Negative binomial outperformed poisson

  • nbinom1 outperformed nbinom2

  • Diagnostics with DHARMa:

    • The 0.25 and 0.5 quantile regressions of the observed residuals against predicted values were significant, but lines are mostly horizontal and deviation occurs mainly in the extremes

    • Some issues with uniformity of residuals within groups (overall, years, units)

Best model call:

## glmmTMB::glmmTMB(formula = Response ~ cos_td_rad + sin_td_rad + 
##     Unit_Label + year_ct + (1 | site/point_name), data = MainData, 
##     family = "nbinom1", ziformula = ~0, dispformula = ~1, na.action = na.fail)

Best model summary:

##  Family: nbinom1  ( log )
## Formula:          
## Response ~ cos_td_rad + sin_td_rad + Unit_Label + year_ct + (1 |  
##     site/point_name)
## Data: MainData
## 
##      AIC      BIC   logLik deviance df.resid 
##  11103.6  11178.4  -5537.8  11075.6     1533 
## 
## Random effects:
## 
## Conditional model:
##  Groups          Name        Variance  Std.Dev.
##  point_name:site (Intercept) 3.401e-01 0.583207
##  site            (Intercept) 2.927e-07 0.000541
## Number of obs: 1547, groups:  point_name:site, 558; site, 70
## 
## Dispersion parameter for nbinom1 family (): 6.31 
## 
## Conditional model:
##                                   Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                       1.091160   0.190243   5.736 9.72e-09 ***
## cos_td_rad                        0.914628   0.123527   7.404 1.32e-13 ***
## sin_td_rad                       -0.874425   0.118457  -7.382 1.56e-13 ***
## Unit_LabelNegev Highlands        -0.124012   0.145277  -0.854  0.39332    
## Unit_LabelInland Sands           -0.635103   0.388749  -1.634  0.10232    
## Unit_LabelLoess Areas             0.094698   0.159607   0.593  0.55296    
## Unit_LabelTransition Zone         0.296265   0.149851   1.977  0.04804 *  
## Unit_LabelPlanted Conifer Forest  0.244692   0.130702   1.872  0.06119 .  
## Unit_LabelMediterranean Maquis    0.759457   0.122601   6.195 5.85e-10 ***
## Unit_LabelBatha Vegetation        0.418060   0.141216   2.960  0.00307 ** 
## year_ct                          -0.014674   0.007742  -1.895  0.05805 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Best model coefficients:

Effect Estimate Std. Error z value Pr(>|z|) 2.5 % 97.5 %
(Intercept) 1.091 0.190 5.736 0.000 0.718 1.464
cos_td_rad 0.915 0.124 7.404 0.000 0.673 1.157
sin_td_rad -0.874 0.118 -7.382 0.000 -1.107 -0.642
Unit_LabelNegev Highlands -0.124 0.145 -0.854 0.393 -0.409 0.161
Unit_LabelInland Sands -0.635 0.389 -1.634 0.102 -1.397 0.127
Unit_LabelLoess Areas 0.095 0.160 0.593 0.553 -0.218 0.408
Unit_LabelTransition Zone 0.296 0.150 1.977 0.048 0.003 0.590
Unit_LabelPlanted Conifer Forest 0.245 0.131 1.872 0.061 -0.011 0.501
Unit_LabelMediterranean Maquis 0.759 0.123 6.195 0.000 0.519 1.000
Unit_LabelBatha Vegetation 0.418 0.141 2.960 0.003 0.141 0.695
year_ct -0.015 0.008 -1.895 0.058 -0.030 0.001

Effect plots - year:

  • Slightly decreasing, almost horizontal

  • Partial residuals look better

x predicted std.error conf.low conf.high group FoldIncrease PercentIncrease FoldDecline PercentDecline Delta
0 12.885 0.048 11.737 14.145 1 1.000 0.000 1.000 0.000 0.000
2 12.512 0.039 11.599 13.498 1 0.971 -2.892 1.030 2.892 -0.373
3 12.330 0.036 11.492 13.230 1 0.957 -4.307 1.045 4.307 -0.555
4 12.151 0.035 11.351 13.007 1 0.943 -5.701 1.060 5.701 -0.735
5 11.974 0.035 11.175 12.830 1 0.929 -7.074 1.076 7.074 -0.912
6 11.799 0.037 10.966 12.696 1 0.916 -8.428 1.092 8.428 -1.086
7 11.627 0.041 10.732 12.597 1 0.902 -9.762 1.108 9.762 -1.258
8 11.458 0.045 10.481 12.526 1 0.889 -11.076 1.125 11.076 -1.427
9 11.291 0.051 10.221 12.474 1 0.876 -12.372 1.141 12.372 -1.594

Effect plots - unit:

  • Partial residuals looks better

  • Highest in Herbaceous and Dwarf-Shrub Vegetation (Batha vegetation in the figures) and Mediterranean Maquis

  • Lowest in the Inland Sands unit

  • Intermediate in all other units

x predicted std.error conf.low conf.high group Fold PercentDiffernce Delta
Arid South 8.625 0.110 6.947 10.708 1 1.887 88.722 4.055
Negev Highlands 7.619 0.099 6.278 9.246 1 1.667 66.711 3.049
Inland Sands 4.570 0.374 2.196 9.513 1 1.000 0.000 0.000
Loess Areas 9.482 0.116 7.545 11.915 1 2.075 107.467 4.911
Transition Zone 11.599 0.105 9.434 14.261 1 2.538 153.798 7.029
Planted Conifer Forest 11.016 0.075 9.516 12.753 1 2.410 141.041 6.446
Mediterranean Maquis 18.432 0.058 16.437 20.670 1 4.033 303.320 13.862
Batha Vegetation 13.101 0.089 11.003 15.600 1 2.867 186.671 8.531

Results summary:

Taxa Birds
Scale National
SpeGroup Resident
DataGroup Time
Response Abundance
Call glmmTMB::glmmTMB(formula = Response ~ cos_td_rad + sin_td_rad + Unit_Label + year_ct + (1 | site/point_name), data = MainData, family = “nbinom1”, ziformula = ~0, dispformula = ~1, na.action = na.fail)
N_Obs 1547
df 14
df.residuals 1533
logLIk -5537.805
AIC 11103.61
BIC 11178.43
deviance 11075.61
sigma 6.307987
ContrastFar_Marginal NA
ContrastNear_Marginal NA
Contrast_Type NA
Contrast_Fold NA
Num_Years 9
Year_Start_Marginal 12.88508
Year_End_Marginal 11.29098
Year_Type Decrease
YearFold_1Year 1.014782
YearFold_NumYear 1.141184
YearPercent_1Year 1.45669
YearPercent_NumYear 12.37171

Session info:

## R version 4.2.3 (2023-03-15 ucrt)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 22621)
## 
## Matrix products: default
## 
## locale:
## [1] LC_COLLATE=he_IL.UTF-8  LC_CTYPE=he_IL.UTF-8    LC_MONETARY=he_IL.UTF-8
## [4] LC_NUMERIC=C            LC_TIME=he_IL.UTF-8    
## 
## attached base packages:
## [1] grid      splines   stats     graphics  grDevices utils     datasets 
## [8] methods   base     
## 
## other attached packages:
##  [1] extrafont_0.19    gridExtra_2.3     ggeffects_1.2.1   ggnewscale_0.4.9 
##  [5] gghalves_0.1.4    GGally_2.1.2      data.table_1.14.8 effects_4.2-3    
##  [9] carData_3.0-5     DHARMa_0.4.6      MuMIn_1.47.5      glmmTMB_1.1.8    
## [13] EnvStats_2.8.1    kableExtra_1.3.4  knitr_1.42        lubridate_1.9.2  
## [17] forcats_1.0.0     stringr_1.5.0     dplyr_1.1.1       purrr_1.0.1      
## [21] readr_2.1.4       tidyr_1.3.0       tibble_3.2.1      ggplot2_3.4.2    
## [25] tidyverse_2.0.0   rmarkdown_2.21   
## 
## loaded via a namespace (and not attached):
##   [1] TH.data_1.1-2      minqa_1.2.5        colorspace_2.1-1  
##   [4] ellipsis_0.3.2     sjlabelled_1.2.0   snakecase_0.11.0  
##   [7] estimability_1.4.1 rstudioapi_0.14    farver_2.1.1      
##  [10] fansi_1.0.4        mvtnorm_1.2-0      xml2_1.3.3        
##  [13] codetools_0.2-19   doParallel_1.0.17  cachem_1.0.7      
##  [16] jsonlite_1.8.4     nloptr_2.0.3       Cairo_1.6-0       
##  [19] Rttf2pt1_1.3.12    shiny_1.7.4        compiler_4.2.3    
##  [22] httr_1.4.5         emmeans_1.8.6      Matrix_1.5-5      
##  [25] fastmap_1.1.1      survey_4.2         cli_3.6.1         
##  [28] later_1.3.0        htmltools_0.5.5    tools_4.2.3       
##  [31] coda_0.19-4        gtable_0.3.3       glue_1.6.2        
##  [34] Rcpp_1.0.10        jquerylib_0.1.4    vctrs_0.6.1       
##  [37] svglite_2.1.1      nlme_3.1-162       extrafontdb_1.0   
##  [40] iterators_1.0.14   insight_0.19.1     xfun_0.38         
##  [43] rbibutils_2.2.16   lme4_1.1-32        rvest_1.0.3       
##  [46] mime_0.12          timechange_0.2.0   lifecycle_1.0.3   
##  [49] gap.datasets_0.0.6 MASS_7.3-58.3      zoo_1.8-11        
##  [52] scales_1.2.1       ragg_1.2.5         promises_1.2.0.1  
##  [55] hms_1.1.3          parallel_4.2.3     sandwich_3.1-0    
##  [58] TMB_1.9.6          RColorBrewer_1.1-3 yaml_2.3.7        
##  [61] sass_0.4.5         reshape_0.8.9      stringi_1.7.12    
##  [64] highr_0.10         gap_1.5-3          foreach_1.5.2     
##  [67] boot_1.3-28.1      qgam_1.3.4         Rdpack_2.6        
##  [70] rlang_1.1.0        pkgconfig_2.0.3    systemfonts_1.0.4 
##  [73] evaluate_0.20      lattice_0.21-8     labeling_0.4.2    
##  [76] tidyselect_1.2.0   plyr_1.8.8         magrittr_2.0.3    
##  [79] R6_2.5.1           generics_0.1.3     multcomp_1.4-23   
##  [82] DBI_1.1.3          pillar_1.9.0       withr_2.5.0       
##  [85] mgcv_1.8-42        survival_3.5-5     nnet_7.3-18       
##  [88] utf8_1.2.3         tzdb_0.3.0         digest_0.6.31     
##  [91] webshot_0.5.4      xtable_1.8-6       httpuv_1.6.9      
##  [94] numDeriv_2022.9-1  textshaping_0.3.6  stats4_4.2.3      
##  [97] munsell_0.5.0      viridisLite_0.4.1  bslib_0.4.2       
## [100] mitools_2.4