Results

Species richness at TF

Figure 1: TF taxon diversity. B = butterfly; P = pollinator; G = ground dwellers *** p < 0.001; ** p < 0.01; * p < 0.05

A total of 22 species of butterflies were recorded in 82 butterfly counts (surveys with no species name recorded were excluded) at 42 TFs with a mean of 1.76 species per survey. The average number of species recorded per TF was 3.88, with a maximum of 9. The most widespread species were small white (24 TFs), large white (18 TFs), and common blue (8). The most abundant were small white, large white and small heath. Simpson diversity index (SDI) varies with species richness (SR) - TFs with higher SR and SDI have greater butterfly diversity.

Within each taxonomic group, SR and SDI were significantly correlated, but there was little association between groups apart from SDIs for ground dwellers (ants excluded) and pollinators ?@fig-ft-but-sr-1.

SPAC

For butterfly species, 22 butterfly species (gamma diversity) were observed in 51 surveys across 41 TFs. As SACs show (Figure 2), saturation is nearly achieved implying that the vast majority of butterfly species occupying these areas have were spotted in the surveys. For ground-dwellers and pollinators, saturation was reached at 13 and 10 species groups respectively, and for both is achieved within a relatively small number of surveys. [discussion - relatively difficulty of observing butterflies, relatively less ‘effort’ per survey - may need fewer pollinator and ground-dweller surveys to estimate gamma diversity].

##3 Survey level analysis

Model outcomes

# A tibble: 7 × 7
  effect   group          term            estimate std.error statistic  p.value
  <chr>    <chr>          <chr>              <dbl>     <dbl>     <dbl>    <dbl>
1 fixed    <NA>           (Intercept)      1.79       0.0387    46.2    0      
2 fixed    <NA>           month            0.00948    0.0314     0.302  0.763  
3 fixed    <NA>           mean_age_days    0.0316     0.0599     0.528  0.597  
4 fixed    <NA>           y                0.0759     0.0353     2.15   0.0318 
5 fixed    <NA>           mean_ht          0.133      0.0683     1.95   0.0509 
6 fixed    <NA>           shannon         -0.112      0.0424    -2.65   0.00807
7 ran_pars tiny_forest_id sd__(Intercept)  0.179     NA         NA     NA

X.Intercept.	mean_age_days	mean_ht	month	shannon	y	AICc
1.8		0.16		-0.118	0.077	863
1.8	0.031	0.14		-0.113	0.076	865
1.8		0.16	0.0085	-0.118	0.077	865
1.8		0.15		-0.113		866
1.8	0.121			-0.069	0.066	867
1.8	0.032	0.13	0.0095	-0.112	0.076	867
1.8	0.040	0.11		-0.106		867
1.8	0.102				0.067	868
1.8		0.15	0.0080	-0.113		868
1.8	0.116			-0.069		868

# A tibble: 7 × 7
  effect   group          term           estimate std.error statistic    p.value
  <chr>    <chr>          <chr>             <dbl>     <dbl>     <dbl>      <dbl>
1 fixed    <NA>           (Intercept)      1.60      0.0462    34.7    1.64e-263
2 fixed    <NA>           month            0.0265    0.0445     0.595  5.52e-  1
3 fixed    <NA>           mean_age_days   -0.0648    0.0674    -0.961  3.37e-  1
4 fixed    <NA>           y               -0.0682    0.0454    -1.50   1.34e-  1
5 fixed    <NA>           mean_ht          0.169     0.0762     2.21   2.68e-  2
6 fixed    <NA>           shannon         -0.0638    0.0509    -1.25   2.10e-  1
7 ran_pars tiny_forest_id sd__(Intercep…   0.0642   NA         NA     NA

X.Intercept.	mean_age_days	mean_ht	month	shannon	y	AICc
1.6		0.096			-0.067	522
1.6		0.106				522
1.6		0.136		-0.061		523
1.6		0.124		-0.059	-0.066	523
1.6		0.086	0.036		-0.071	524
1.6	-0.048	0.142				524
1.6		0.099	0.028			524
1.6	-0.042	0.128			-0.066	524
1.6	-0.072	0.196		-0.074		524
1.6	-0.065	0.179		-0.071	-0.065	524

# A tibble: 7 × 7
  effect   group          term            estimate std.error statistic   p.value
  <chr>    <chr>          <chr>              <dbl>     <dbl>     <dbl>     <dbl>
1 fixed    <NA>           (Intercept)      8.11e-1    0.0979    8.29    1.13e-16
2 fixed    <NA>           month            2.43e-2    0.102     0.238   8.12e- 1
3 fixed    <NA>           mean_age_days    1.65e-2    0.137     0.120   9.04e- 1
4 fixed    <NA>           y               -6.94e-2    0.105    -0.664   5.07e- 1
5 fixed    <NA>           mean_ht         -2.16e-1    0.150    -1.44    1.50e- 1
6 fixed    <NA>           shannon         -4.11e-3    0.108    -0.0380  9.70e- 1
7 ran_pars tiny_forest_id sd__(Intercept)  5.01e-6   NA        NA      NA

X.Intercept.	mean_age_days	mean_ht	month	shannon	y	AICc
0.81		-0.19				187
0.82						188
0.81		-0.21			-0.071	189
0.82	-0.109					189
0.81	0.028	-0.21				189
0.81		-0.19	0.0175			189
0.81		-0.19		-0.0099		189
0.82				-0.0545		190
0.82					-0.025	190
0.82			0.0043			190

Model outcomes

Descriptive analysis of SR by taxa at TFb

As Figure 12 shows insects, vasculat plants and birds are the most diverse taxa across all habitat type. Insect and plant diversity are greatest in residentail gardens; bird and mamal diversity are siialr betwen gardens, grasslant and arable land.

Summary of environmental predictors

Table 1:
Environmental predictor summary table
Category	Metric	Mean	Sd	Q25	Q75
Blue_infrastructure	still_water	21590	39716	2451	21590
Blue_infrastructure	watercourse	39287	135997	3753	39287
Blue_infrastructure	nearest_water	276	251	82	410
Built	built_area	2086421	653916	1694140	2591003
Climate	tf_temp_2020	10.88	0.76	10.65	11.17
Climate	tf_temp_2022	11.27	0.70	11.02	11.56
Climate	tf_rain_2020	70	20	59	71
Climate	tf_rain_2021	61.5	13.1	54.6	63.1
Climate	tf_rain_2022	55.8	14.7	47.8	57.5
Climate	mean_spring_rain	145	29	130	153
Climate	mean_spring_temp	9.28	0.66	9.02	9.55
Connectivity	nc_y	10.0	3.8	7.0	12.0
Connectivity	lnk_y	136	95	58	205
Connectivity	slc_y	113078	131052	26100	144700
Connectivity	msc_y	21423	53115	5630	21423
Connectivity	ccp	0.48	0.22	0.32	0.66
Connectivity	lcp	0.0046	0.0110	0.0000	0.0046
Connectivity	cpl	2.7	1.1	1.9	3.2
Connectivity	ecs	94889	124745	17650	117182
Connectivity	awf	6042740107	17522014830	173363627	3765155663
Connectivity	iic_y	0.0026	0.0071	0.0000	0.0000
Green_infrastructure	perc_tree	0.0187	0.0245	0.0038	0.0220
Green_infrastructure	perc_garden	0.201	0.089	0.136	0.260
Green_infrastructure	perc_veg	0.41	0.16	0.32	0.50
Green_infrastructure	wood_area_sum	227252	202858	82726	327279
Green_infrastructure	gi_prop	0.55	0.18	0.42	0.68
Green_infrastructure	tree_density_m_2	0.73	0.36	0.49	0.98
Greenness	green_prop	0.30	0.21	0.12	0.44
Greenness	ndvi	0.51	0.10	0.46	0.57
Greenness	savi	0.263	0.071	0.218	0.306
Greenness	gbndvi	0.28	0.12	0.23	0.35
Greenness	evi_2	0.310	0.086	0.256	0.361
Landscape	aggregation_index	92.4	2.1	91.0	94.0
Landscape	modified_simpsons_diversity_index	0.89	0.31	0.69	1.11
Landscape	perimeter_area_fractal_dimension	1.317	0.027	1.303	1.335
Landscape	patch_richness	8.3	1.6	7.0	9.0
Landscape	shape_index	47.1	8.2	41.1	52.4
Linear_features	hedge_len	3870	2898	1575	5692
Linear_features	road_length	40002	11801	33248	47869
Tf_char	tf_age_y	579	204	572	655
Tf_char	total_trees	142	60	100	200
Tf_char	elevation	80	51	32	121
Tf_char	lon	-1.79	1.11	-2.13	-1.16
Tf_char	lat	53.00	1.39	52.36	53.32

Many indicators are highly correlated or cluster as shown in Figure 5 which orders pair-wise correlations of scaled (normalised) metrics by hierarchical cluster. Rainfall metrics cluster and are inversely correlated with temperature indicators. Gardens, built areas and roads form a separate cluster as do connectivity metrics, and vegetation/‘greenness’ indicators.

(a) Non-significant correlation are masked by an X

Gamma diversity

Results

Mean TF butterfly SR was 2.29 (survey sites = 42) and varied between 1 and 8. Taxon richness also varied between 1 and 13 (mean = 6.2, sites = 102) for pollinators and 1 to 10 (mean = 5.03, sites = 81) for ground dwelling invertebrates respectively.

Figure 7: Correlation between TF butterfly SR and insect and butterfly TFb SR

As Figure 7 shows, TF butterfly richness is uncorrelated with TFb butterfly or insect SR.

Intrinsic TF variables were then fitted to butterfly SR using GLMM with a negative binomial error structure and TF as a random effect. The 10 models with the lowest AICc fitted with the dredge function of the MuMIm R package are shown in ?@tbl-butterfly-model (Bartoń (2023)). Each row of the table shows the coefficients of the explanatory variables - the number of surveys (n) is included in all 10 models, with mean age at time of survey (scaled) included in four. Survey month and weighted tree height are included in 2 candidate models.

Model selection can be difficult. It has been common practice in ecology to average models with with an AICc within 7 of the model of the lowest AICc but this has been criticised for leading to faulty inferences about relevant predictors (Galipaud et al. (2014); Walker (2017)).

An alternative is use penalised regression where predictors are removed

Ground dwelling invertebrates

GDM

Generalised dissimilarity models (GDM) were fitted to community matrices generated for each TF site for butterflies, pollinators and ground dwelling invertebrates to test the hypothesis that community composition changed as forest aged. Modelling followed the method set out by Mokany et al. (2022).

Figure 8: Partial dependence plots for butterfly GDM

Figure 9: GDM results for ground-dwellers

Figure 10: GDM results for ground-dwellers

Table 2:
Variable importance
taxon	Important variables	All predictors
Butterflies	Geographic...1	0.00
Butterflies	shannon...2	0.28
Butterflies	mean_ht...3	0.54
Pollinators	Geographic...4	0.20
Pollinators	mean_ht...5	0.58
Pollinators	tf_age...6	0.28
Ground dwellers)	Geographic...7	0.00
Ground dwellers)	shannon...8	0.68
Ground dwellers)	mean_ht...9	0.66
Ground dwellers)	tf_age...10	0.00

Relationship between intrinsic tiny forest variables and tiny forest species, richness

Impacts of characteristics of surrounding area on tiny forest, biodiversity line Abiotic. Biotic

Env
{#fig-glmnet}
taxon	metric	estimate
Butterfly	gi_prop	0.18352018
Butterfly	n	0.03240214
Butterfly	lnk_y	0.00093757
Butterfly	hedge_len	0.00000651
Butterfly	watercourse	0.00000115
Butterfly	wood_area_sum	0.00000012
Butterfly	nc_y	-0.00009932
Butterfly	tf_rain_2021	-0.00303056
Butterfly	perc_tree	-6.97509738
Ground-dwellers	lat	0.02982828
Ground-dwellers	n	0.01232032
Ground-dwellers	tf_age_y	0.00014862
Pollinators	n	0.04499600

Figure 11: Butterfly predictors

Environmnental predictors of taxon richness

The results of glment modelling are shown in Figure 11. For butterflies, green infrastructure, linear features (water courses and hedges), wooded areas, rainfall and connectivity metrics variables were found to be weakly positively associated with butterfly SR, whereas tree cover was strongly negatively predictive. For pollinators and ground-dwellers, SR was associated with the number of surveys, and increased with forest age and latitude for ground-dwellers.

Estimated SR for small areas based on random sampling

Estimated SR for small areas based on SAR

Figure 13: Relationship between log habitat area and log SR

Figure 14: Relationship between log habitat area and log SR

Sampling


  |                                                                            
  |                                                                      |   0%
  |                                                                            
  |                                                                      |   1%
  |                                                                            
  |=                                                                     |   1%
  |                                                                            
  |=                                                                     |   2%
  |                                                                            
  |==                                                                    |   2%
  |                                                                            
  |==                                                                    |   3%
  |                                                                            
  |==                                                                    |   4%
  |                                                                            
  |===                                                                   |   4%
  |                                                                            
  |===                                                                   |   5%
  |                                                                            
  |====                                                                  |   5%
  |                                                                            
  |====                                                                  |   6%
  |                                                                            
  |=====                                                                 |   6%
  |                                                                            
  |=====                                                                 |   7%
  |                                                                            
  |=====                                                                 |   8%
  |                                                                            
  |======                                                                |   8%
  |                                                                            
  |======                                                                |   9%
  |                                                                            
  |=======                                                               |   9%
  |                                                                            
  |=======                                                               |  10%
  |                                                                            
  |=======                                                               |  11%
  |                                                                            
  |========                                                              |  11%
  |                                                                            
  |========                                                              |  12%
  |                                                                            
  |=========                                                             |  12%
  |                                                                            
  |=========                                                             |  13%
  |                                                                            
  |=========                                                             |  14%
  |                                                                            
  |==========                                                            |  14%
  |                                                                            
  |==========                                                            |  15%
  |                                                                            
  |===========                                                           |  15%
  |                                                                            
  |===========                                                           |  16%
  |                                                                            
  |============                                                          |  16%
  |                                                                            
  |============                                                          |  17%
  |                                                                            
  |============                                                          |  18%
  |                                                                            
  |=============                                                         |  18%
  |                                                                            
  |=============                                                         |  19%
  |                                                                            
  |==============                                                        |  19%
  |                                                                            
  |==============                                                        |  20%
  |                                                                            
  |==============                                                        |  21%
  |                                                                            
  |===============                                                       |  21%
  |                                                                            
  |===============                                                       |  22%
  |                                                                            
  |================                                                      |  22%
  |                                                                            
  |================                                                      |  23%
  |                                                                            
  |================                                                      |  24%
  |                                                                            
  |=================                                                     |  24%
  |                                                                            
  |=================                                                     |  25%
  |                                                                            
  |==================                                                    |  25%
  |                                                                            
  |==================                                                    |  26%
  |                                                                            
  |===================                                                   |  26%
  |                                                                            
  |===================                                                   |  27%
  |                                                                            
  |===================                                                   |  28%
  |                                                                            
  |====================                                                  |  28%
  |                                                                            
  |====================                                                  |  29%
  |                                                                            
  |=====================                                                 |  29%
  |                                                                            
  |=====================                                                 |  30%
  |                                                                            
  |=====================                                                 |  31%
  |                                                                            
  |======================                                                |  31%
  |                                                                            
  |======================                                                |  32%
  |                                                                            
  |=======================                                               |  32%
  |                                                                            
  |=======================                                               |  33%
  |                                                                            
  |=======================                                               |  34%
  |                                                                            
  |========================                                              |  34%
  |                                                                            
  |========================                                              |  35%
  |                                                                            
  |=========================                                             |  35%
  |                                                                            
  |=========================                                             |  36%
  |                                                                            
  |==========================                                            |  36%
  |                                                                            
  |==========================                                            |  37%
  |                                                                            
  |==========================                                            |  38%
  |                                                                            
  |===========================                                           |  38%
  |                                                                            
  |===========================                                           |  39%
  |                                                                            
  |============================                                          |  39%
  |                                                                            
  |============================                                          |  40%
  |                                                                            
  |============================                                          |  41%
  |                                                                            
  |=============================                                         |  41%
  |                                                                            
  |=============================                                         |  42%
  |                                                                            
  |==============================                                        |  42%
  |                                                                            
  |==============================                                        |  43%
  |                                                                            
  |==============================                                        |  44%
  |                                                                            
  |===============================                                       |  44%
  |                                                                            
  |===============================                                       |  45%
  |                                                                            
  |================================                                      |  45%
  |                                                                            
  |================================                                      |  46%
  |                                                                            
  |=================================                                     |  46%
  |                                                                            
  |=================================                                     |  47%
  |                                                                            
  |=================================                                     |  48%
  |                                                                            
  |==================================                                    |  48%
  |                                                                            
  |==================================                                    |  49%
  |                                                                            
  |===================================                                   |  49%
  |                                                                            
  |===================================                                   |  50%
  |                                                                            
  |===================================                                   |  51%
  |                                                                            
  |====================================                                  |  51%
  |                                                                            
  |====================================                                  |  52%
  |                                                                            
  |=====================================                                 |  52%
  |                                                                            
  |=====================================                                 |  53%
  |                                                                            
  |=====================================                                 |  54%
  |                                                                            
  |======================================                                |  54%
  |                                                                            
  |======================================                                |  55%
  |                                                                            
  |=======================================                               |  55%
  |                                                                            
  |=======================================                               |  56%
  |                                                                            
  |========================================                              |  56%
  |                                                                            
  |========================================                              |  57%
  |                                                                            
  |========================================                              |  58%
  |                                                                            
  |=========================================                             |  58%
  |                                                                            
  |=========================================                             |  59%
  |                                                                            
  |==========================================                            |  59%
  |                                                                            
  |==========================================                            |  60%
  |                                                                            
  |==========================================                            |  61%
  |                                                                            
  |===========================================                           |  61%
  |                                                                            
  |===========================================                           |  62%
  |                                                                            
  |============================================                          |  62%
  |                                                                            
  |============================================                          |  63%
  |                                                                            
  |============================================                          |  64%
  |                                                                            
  |=============================================                         |  64%
  |                                                                            
  |=============================================                         |  65%
  |                                                                            
  |==============================================                        |  65%
  |                                                                            
  |==============================================                        |  66%
  |                                                                            
  |===============================================                       |  66%
  |                                                                            
  |===============================================                       |  67%
  |                                                                            
  |===============================================                       |  68%
  |                                                                            
  |================================================                      |  68%
  |                                                                            
  |================================================                      |  69%
  |                                                                            
  |=================================================                     |  69%
  |                                                                            
  |=================================================                     |  70%
  |                                                                            
  |=================================================                     |  71%
  |                                                                            
  |==================================================                    |  71%
  |                                                                            
  |==================================================                    |  72%
  |                                                                            
  |===================================================                   |  72%
  |                                                                            
  |===================================================                   |  73%
  |                                                                            
  |===================================================                   |  74%
  |                                                                            
  |====================================================                  |  74%
  |                                                                            
  |====================================================                  |  75%
  |                                                                            
  |=====================================================                 |  75%
  |                                                                            
  |=====================================================                 |  76%
  |                                                                            
  |======================================================                |  76%
  |                                                                            
  |======================================================                |  77%
  |                                                                            
  |======================================================                |  78%
  |                                                                            
  |=======================================================               |  78%
  |                                                                            
  |=======================================================               |  79%
  |                                                                            
  |========================================================              |  79%
  |                                                                            
  |========================================================              |  80%
  |                                                                            
  |========================================================              |  81%
  |                                                                            
  |=========================================================             |  81%
  |                                                                            
  |=========================================================             |  82%
  |                                                                            
  |==========================================================            |  82%
  |                                                                            
  |==========================================================            |  83%
  |                                                                            
  |==========================================================            |  84%
  |                                                                            
  |===========================================================           |  84%
  |                                                                            
  |===========================================================           |  85%
  |                                                                            
  |============================================================          |  85%
  |                                                                            
  |============================================================          |  86%
  |                                                                            
  |=============================================================         |  86%
  |                                                                            
  |=============================================================         |  87%
  |                                                                            
  |=============================================================         |  88%
  |                                                                            
  |==============================================================        |  88%
  |                                                                            
  |==============================================================        |  89%
  |                                                                            
  |===============================================================       |  89%
  |                                                                            
  |===============================================================       |  90%
  |                                                                            
  |===============================================================       |  91%
  |                                                                            
  |================================================================      |  91%
  |                                                                            
  |================================================================      |  92%
  |                                                                            
  |=================================================================     |  92%
  |                                                                            
  |=================================================================     |  93%
  |                                                                            
  |=================================================================     |  94%
  |                                                                            
  |==================================================================    |  94%
  |                                                                            
  |==================================================================    |  95%
  |                                                                            
  |===================================================================   |  95%
  |                                                                            
  |===================================================================   |  96%
  |                                                                            
  |====================================================================  |  96%
  |                                                                            
  |====================================================================  |  97%
  |                                                                            
  |====================================================================  |  98%
  |                                                                            
  |===================================================================== |  98%
  |                                                                            
  |===================================================================== |  99%
  |                                                                            
  |======================================================================|  99%
  |                                                                            
  |======================================================================| 100%

Butterfly richness estimated by sampling

Based on 300 (of 10,000) randomly sample 4Ha pseudo-controls, median butterfly richness was , with a percentile range of Butterfly diversity for TFs was 22 which is at the 95th centile of the distribution. Butterfly richness was highest in suburban areas approaching near saturation at 34 species. Richness in arable areas, improved grassland and deciduous woodland were lower but based on smaller samples.For a single 4Ha site the highest SR was in fen habitat (15 species), with arable (12), suburban (11), grassland (10) and woodland (8).

References

Bartoń, K. (2023) “MuMIn: Multi-model inference.” Available at: https://CRAN.R-project.org/package=MuMIn.

Galipaud, M. et al. (2014) “Ecologists overestimate the importance of predictor variables in model averaging: a plea for cautious interpretations,” Methods in Ecology and Evolution. Edited by R.B. O’Hara, 5(10), pp. 983–991. doi:10.1111/2041-210x.12251.

Mokany, K. et al. (2022) “A working guide to harnessing generalized dissimilarity modelling for biodiversity analysis and conservation assessment,” Global Ecology and Biogeography. Edited by V. Bahn, 31(4), pp. 802–821. doi:10.1111/geb.13459.

Walker, J.A. (2017) “Model-averaged regression coefficients have a straightforward interpretation using causal conditioning.” Available at: http://dx.doi.org/10.1101/133785.

--- title: "Results" format: html: code-fold: true code-tools: true fig-width: 8 fig-height: 6 toc: true toc-location: right docx: reference-doc: reference.docx fig-width: 8 fig-height: 8 pdf: default execute: echo: false warning: false message: false cache: refresh editor: visual bibliography: ../bibliography/Paperpile - Dissertation - Jan 06.bib csl: "https://github.com/citation-style-language/styles/blob/master/harvard-cite-them-right-11th-edition.csl?raw=TRUE" --- ```{r} #| label: setup #| include: false #| cache: false source("startup.R") #source("../R/glmnet_taxon.R") theme_thesis <- theme_set(ggthemes::theme_base() + theme(plot.title.position = "plot", text = element_text(size = 10), plot.background = element_blank())) env_pred <- env_pred_dt needs(vegan, ggcorrplot, glmnet) glmnet_fit <- function(df, taxon = "Pollinators"){ set.seed(1234) x <- df |> select(-Species) |> as.matrix() y <- df |> select(Species) |> as.matrix() cv.fit <- cv.glmnet(x, y, family = "poisson") best_l <- cv.fit$lambda.min coeffs <- coef(cv.fit, s = best_l) coef_v <- as.matrix(coeffs) |> data.frame() |> rownames_to_column("metric") |> filter(s1 != 0) coef_plot <- coef_v |> filter(metric != "(Intercept)") |> ggplot() + geom_col(aes(reorder(metric, s1), s1 )) + geom_hline(yintercept = 0, colour = "red") + labs(title = "Retained non-zero terms in glmnet model", subtitle = taxon, x = "", y = "Estimate") return(list(plot = coef_plot, df = coef_v)) } conflicted::conflict_prefer("filter", "dplyr") conflicted::conflict_prefer("select", "dplyr") options(na.action = "na.fail") ceh_lookup <- data.frame( land_cover = c("deciduous woodland", "coniferous woodland", "arable", "improved grassland", "neutral grassland", "calcareous grassland", "acid grassland", "fen", "heather", "heather grassland", "bog", "inland rock", "saltwater", "fresh water", "supralittoral rock", "suralittoral sediment", "littoral rock", "littoral sediment", "saltmarsh", "urban", "suburban"), lc_id = 1:21, id = c(1, 2, 3, 4, 5, 5, 5, 5, 6, 6, 6, 6, 7, 8, 9, 9, 9, 9, 9, 10, 10) ) ceh_agg_lu <- data.frame(aggregate_class = c("broadleaf woodland", "coniferous woodland", "arable", "improved grass", "semi-natural grass", "heath and bog", "satlwater", "freshwater", "coastal", "built and gardens"), id = 1:10) ceh_lu <- left_join(ceh_lookup, ceh_agg_lu) ``` ### Species richness at TF ```{r} #| label: fig-ft-but-sr #| fig-cap: "TF taxon diversity. B = butterfly; P = pollinator; G = ground dwellers\n*** p < 0.001; ** p < 0.01; * p < 0.05" #| fig-height: 6 butterfly <- butterfly |> filter(!is.na(qty)) butterfly_spp <- butterfly |> group_by(butterfly_species) |> dplyr::filter(!is.na(qty)) |> reframe(tf = n_distinct(tiny_forest_id), tot = sum(qty)) |> arrange(-tf) butterfly_surveys <- n_distinct(butterfly$survey_id) butterfly_sp_surv <- butterfly |> group_by(survey_id) |> reframe(n_sp = n_distinct(butterfly_species)) butterfly_ssm <- butterfly |> filter(!is.na(qty)) |> group_by(tiny_forest_id, butterfly_species, qty) |> reframe(sumq = sum(qty)) |> select(-qty) |> pivot_wider(names_from = butterfly_species, values_from = sumq, values_fn = \(x) sum(x), values_fill = 0) butterfly_sr <- butterfly_ssm |> vegan::specnumber() butterfly_gamma <- butterfly_ssm[, -1] |> poolaccum() butterfly_simp <- butterfly_ssm |> vegan::diversity("simpson") butterfly_div <- bind_cols(site = butterfly_ssm$tiny_forest_id, `Species richness` = butterfly_sr, `Simpsons diversity index` = butterfly_simp) # butterfly_div |> # pivot_longer(names_to = "Metric", values_to = "Values", cols = 2:3) |> # ggplot() + # geom_histogram(aes(Values), fill = "grey90", colour = "black") + # facet_wrap(~Metric, scales = "free") # # butterfly_div |> # ggplot() + # geom_point(aes(`Species richness`,`Simpsons diversity index` )) + # geom_smooth(aes(`Species richness`,`Simpsons diversity index` ), method = "lm") poll_spp <- pollinators |> group_by(species) |> filter(!is.na(qty)) |> reframe(tf = n_distinct(tiny_forest_id), tot = sum(qty)) |> arrange(-tf) poll_ssm <- pollinators |> filter(!is.na(qty)) |> group_by(tiny_forest_id, species, qty) |> reframe(sumq = sum(qty)) |> select(-qty) |> pivot_wider(names_from = species, values_from = sumq, values_fn = \(x) sum(x), values_fill = 0) poll_sr <- poll_ssm |> vegan::specnumber() poll_simp <- poll_ssm |> vegan::diversity("simpson") poll_div <- bind_cols(site = poll_ssm$tiny_forest_id, `Species richness` = poll_sr, `Simpsons diversity index` = poll_simp) # poll_div |> # pivot_longer(names_to = "Metric", values_to = "Values", cols = 2:3) |> # ggplot() + # geom_histogram(aes(Values), fill = "grey90", colour = "black") + # facet_wrap(~Metric, scales = "free") # # poll_div |> # ggplot() + # geom_point(aes(`Species richness`,`Simpsons diversity index` )) + # geom_smooth(aes(`Species richness`,`Simpsons diversity index` ), method = "lm") gd_spp <- ground_dwellers |> group_by(species) |> filter(!is.na(qty)) |> reframe(tf = n_distinct(tiny_forest_id), tot = sum(qty)) |> arrange(-tf) gd_ssm <- ground_dwellers |> filter(!is.na(qty)) |> group_by(tiny_forest_id, species, qty) |> reframe(sumq = sum(qty)) |> select(-qty) |> pivot_wider(names_from = species, values_from = sumq, values_fn = \(x) sum(x), values_fill = 0) gd_sr <- gd_ssm[, -1] |> vegan::specnumber() gd_simp <- gd_ssm[, -1] |> vegan::diversity("simpson") gd_simp_no_ant <- gd_ssm[, -c(1:2)] |> vegan::diversity("simpson") gd_div <- bind_cols(site = gd_ssm$tiny_forest_id, `Species richness` = gd_sr, `Simpsons diversity index` = gd_simp, SDI_no_ants = gd_simp_no_ant) # gd_div |> # pivot_longer(names_to = "Metric", values_to = "Values", cols = 2:4) |> # ggplot() + # geom_histogram(aes(Values), fill = "grey90", colour = "black") + # facet_wrap(~Metric, scales = "free") # # gd_div |> # ggplot() + # geom_point(aes(`Species richness`, SDI_no_ants )) + # geom_smooth(aes(`Species richness`, SDI_no_ants ), method = "lm") taxon_div <<- butterfly_div |> inner_join(poll_div, by = "site") |> left_join(gd_div, by = "site") taxon_div |> rename(`B SR` = 2, `B SDI` = 3, `P SR` = 4, `P SDI` = 5, `G SR` = 6, `G SDI` = 7, `G SDI (no ants)` = 8) |> ggpairs(2:8, progress = FALSE, lower = list(continuous = "smooth")) ``` A total of `r max(butterfly_gamma$S)` species of butterflies were recorded in `r butterfly_surveys` butterfly counts (surveys with no species name recorded were excluded) at `r nrow(butterfly_ssm)` TFs with a mean of `r mean(butterfly_sp_surv$n_sp)` species per survey. The average number of species recorded per TF was `r mean(butterfly_sr)`, with a maximum of `r max(butterfly_sr)`. The most widespread species were small white (`r butterfly_spp[1,2]` TFs), large white (`r butterfly_spp[2,2]` TFs), and common blue (`r butterfly_spp[3,2]`). The most abundant were small white, large white and small heath. Simpson diversity index (SDI) varies with species richness (SR) - TFs with higher SR and SDI have greater butterfly diversity. Within each taxonomic group, SR and SDI were significantly correlated, but there was little association between groups apart from SDIs for ground dwellers (ants excluded) and pollinators @fig-ft-but-sr-1. ### SPAC ```{r} #| label: fig-gamma #| fig-cap: "Species accumulation curves derived from EW surveys" #| fig-subcap: #| - "Butterflies" #| - "Pollinators" #| - "Ground-dwellers" #| layout_ncol: 3 srb_pool <- butterfly |> count(tiny_forest_id, month, butterfly_species) |> pivot_wider(names_from = butterfly_species, values_from = n, values_fill = 0) |> dplyr::select(-c(1:2)) |> vegan::poolaccum() srpol_pool <- pollinators |> count(tiny_forest_id, month, species) |> pivot_wider(names_from = species, values_from = n, values_fill = 0) |> dplyr::select(-c(1:2)) |> vegan::poolaccum() srgd_pool <- ground_dwellers |> count(tiny_forest_id, month, species) |> pivot_wider(names_from = species, values_from = n, values_fill = 0) |> dplyr::select(-c(1:2)) |> vegan::poolaccum() # srb <- srb |> # bind_cols(survey_counts) |> # rename(sr = 1) # # srp <- srp |> # bind_cols(survey_counts) |> # rename(sr = 1) srb_pool$means |> data.frame() |> ggplot(aes(N, S)) + geom_line() + labs(x = "Surveys", y = "SR", title = "SAC") srpol_pool$means |> data.frame() |> ggplot(aes(N, S)) + geom_line() + labs(x = "Surveys", y = "SR", title = "SAC") srgd_pool$means |> data.frame() |> ggplot(aes(N, S)) + geom_line() + labs(x = "Surveys", y = "SR", title = "SAC") ``` For butterfly species, 22 butterfly species (gamma diversity) were observed in 51 surveys across 41 TFs. As SACs show (@fig-gamma), saturation is nearly achieved implying that the vast majority of butterfly species occupying these areas have were spotted in the surveys. For ground-dwellers and pollinators, saturation was reached at 13 and 10 species groups respectively, and for both is achieved within a relatively small number of surveys. \[discussion - relatively difficulty of observing butterflies, relatively less 'effort' per survey - may need fewer pollinator and ground-dweller surveys to estimate gamma diversity\]. ##3 Survey level analysis ```{r} #| label: taxon-sr-survey-level butterfly_1 <- butterfly |> filter(!is.na(qty)) |> mutate(event_date = dmy(event_date)) |> inner_join(tfb_coords, by = c("tiny_forest_id" = "tf_id")) |> inner_join(tree_shannon) |> mutate(age = as.numeric(event_date - plant_date)) |> dplyr::select(tiny_forest_id, survey_id, qty, month, x, y, shannon, mean_ht, age, butterfly_species) survey_counts <- butterfly_1 |> group_by(tiny_forest_id, month) |> reframe(n = n(), mean_age = mean(age), Y = mean(y), shannon = mean(shannon), mean_ht = mean(mean_ht)) srb <- butterfly_1 |> count(tiny_forest_id, month, butterfly_species) |> pivot_wider(names_from = butterfly_species, values_from = n, values_fill = 0) |> nest_by(tiny_forest_id, month) |> mutate(sr = list(vegan::specnumber(data)), sdi = list(vegan::diversity(data, "simpson"))) |> select(-data) |> unnest(sdi) |> unnest(sr) |> inner_join(tree_shannon) |> inner_join(tfb_coords, by = c("tiny_forest_id" = "tf_id")) srb_pool <- butterfly_1 |> count(tiny_forest_id, month, butterfly_species) |> pivot_wider(names_from = butterfly_species, values_from = n, values_fill = 0) |> dplyr::select(-c(1:2)) |> vegan::poolaccum() # # srb <- srb |> # bind_cols(survey_counts) |> # rename(sr = 1) ####################### ####################### pollinators_1 <- pollinators |> filter(!is.na(qty)) |> mutate(event_date = dmy(event_date)) |> inner_join(tfb_coords, by = c("tiny_forest_id" = "tf_id")) |> inner_join(tree_shannon) |> mutate(age = as.numeric(event_date - plant_date)) |> dplyr::select(tiny_forest_id, survey_id, qty, month, x, y, shannon, mean_ht, age, species) survey_counts_p <- pollinators_1 |> group_by(tiny_forest_id, month) |> reframe(n = n(), mean_age = mean(age), Y = mean(y), shannon = mean(shannon), mean_ht = mean(mean_ht)) srp <- pollinators_1 |> count(tiny_forest_id, month, species) |> pivot_wider(names_from = species, values_from = n, values_fill = 0) |> nest_by(tiny_forest_id, month) |> mutate(sr = list(vegan::specnumber(data)), sdi = list(vegan::diversity(data, "simpson"))) |> select(-data) |> unnest(sdi) |> unnest(sr) |> inner_join(tree_shannon) |> inner_join(tfb_coords, by = c("tiny_forest_id" = "tf_id")) srp_pool <- pollinators_1 |> count(tiny_forest_id, month, species) |> pivot_wider(names_from = species, values_from = n, values_fill = 0) |> dplyr::select(-c(1:2)) |> vegan::poolaccum() # srb <- srb |> # bind_cols(survey_counts) |> # rename(sr = 1) ############################ ############################ gd_1 <- ground_dwellers |> filter(!is.na(qty)) |> mutate(event_date = dmy(event_date)) |> inner_join(tfb_coords, by = c("tiny_forest_id" = "tf_id")) |> inner_join(tree_shannon) |> mutate(age = as.numeric(event_date - plant_date)) |> dplyr::select(tiny_forest_id, survey_id, qty, month, x, y, shannon, mean_ht, age, species) survey_counts_g <- gd_1 |> group_by(tiny_forest_id, month) |> reframe(n = n(), mean_age = mean(age), Y = mean(y), shannon = mean(shannon), mean_ht = mean(mean_ht)) srg <- gd_1 |> count(tiny_forest_id, month, species) |> pivot_wider(names_from = species, values_from = n, values_fill = 0) |> nest_by(tiny_forest_id, month) |> mutate(sr = list(vegan::specnumber(data)), sdi = list(vegan::diversity(data, "simpson"))) |> select(-data) |> unnest(sdi) |> unnest(sr) |> inner_join(tree_shannon) |> inner_join(tfb_coords, by = c("tiny_forest_id" = "tf_id")) ``` ### Model outcomes ```{r} #| label: models-glmer-neg-bin #| tbl-cap: "Model outcomes" srg_scale <- srg |> ungroup() |> select(tiny_forest_id, Species = sr, month, shannon, mean_ht, y, mean_age_days) |> mutate_at(vars(3:7), scale) gd_mod <- glmer.nb(Species ~ month + mean_age_days + y + mean_ht + shannon + (1|tiny_forest_id), data = srg_scale, family = "poisson") gd_gam <- gam(Species ~ s(month) + s(mean_age_days) + s(y) + s(mean_ht) + s(shannon) + s(tiny_forest_id, bs ="re"), data = srg_scale, family = negbin(20)) gd_df <- srg_scale |> ungroup() |> select(Species, month, shannon, mean_ht, y, mean_age_days) gd_glmnet <- glmnet_fit(gd_df, taxon = "Ground dwellers") gd_glmnet$plot tidy(gd_mod) plot(gd_gam, shift = coef(gd_gam)[1]) gd_models <- dredge(gd_mod) gd_models |> slice(1:10) |> data.frame() |> dplyr::select(-c(df, logLik, delta, weight)) |> flextable() ########### pollinators ########### srp_scale <- srp |> ungroup() |> select(tiny_forest_id, Species = sr, month, shannon, mean_ht, y, mean_age_days) |> mutate_at(vars(3:7), scale) p_mod <- glmer.nb(Species ~ month + mean_age_days + y + mean_ht + shannon + (1|tiny_forest_id), data = srp_scale, family = "poisson") p_gam <- gam(Species ~ month + s(mean_age_days) + s(y) + s(mean_ht) + s(shannon) + s(tiny_forest_id, bs ="re"), data = srp_scale, family = negbin(20)) p_df <- srp_scale |> ungroup() |> select(Species, month, shannon, mean_ht, y, mean_age_days) p_glmnet <- glmnet_fit(p_df, taxon = "Pollinators") p_glmnet$plot tidy(p_mod) plot(p_gam, shift = coef(p_gam)[1]) p_models <- dredge(p_mod) p_models |> slice(1:10) |> data.frame() |> dplyr::select(-c(df, logLik, delta, weight)) |> flextable() ############# ############ srb_scale <- srb |> ungroup() |> select(tiny_forest_id, Species = sr, month, shannon, mean_ht, y, mean_age_days) |> mutate_at(vars(3:7), scale) b_mod <- glmer.nb(Species ~ month + mean_age_days + y + mean_ht + shannon + (1|tiny_forest_id), data = srb_scale, family = "poisson") b_gam <- gam(Species ~ month + s(mean_age_days) + s(y) + s(mean_ht) + s(shannon) + s(tiny_forest_id, bs ="re"), data = srb_scale, family = negbin(20)) b_df <- srb_scale |> ungroup() |> select(Species, month, shannon, mean_ht, y, mean_age_days) b_glmnet <- glmnet_fit(b_df, taxon = "Butterflies") b_glmnet$plot tidy(b_mod) plot(b_gam, shift = coef(gd_gam)[1]) b_models <- dredge(b_mod) b_models |> slice(1:10) |> data.frame() |> dplyr::select(-c(df, logLik, delta, weight)) |> flextable() ``` #### Descriptive analysis of SR by taxa at TFb ```{r} #| label: fig-nbn-sr #| fig-cap: "Species richness from NBN data" #| fig-subcap: #| - "TFb SR distribution" #| - "SR by taxon class" #| layout-ncol: 1 #| column: page nbn_d <- read_sf("../data/spatial/nbn_df_sf.gpkg") nbn_d_tf <- nbn_d |> st_intersection(tfb) nbn_d_tf <- nbn_d_tf |> st_drop_geometry() |> group_by(occurrenceID) |> slice(1) |> ungroup() nbn_d_tf_ssm <- nbn_d_tf |> count(tf_id, classs, species) |> pivot_wider(names_from = species, values_from = n, values_fill = 0) nbn_d_tf_acc <- nbn_d_tf_ssm |> ungroup() |> dplyr::select(-c(2)) |> nest_by(tf_id) |> mutate(sr_nbn = list(vegan::specpool(data))) |> dplyr::select(-data) |> unnest(sr_nbn) nbn_d_tf_acc |> ggplot() + geom_histogram(aes(Species), fill = "white", colour = "black") + labs(y = "Number of TFbs", x = "Species richness") nbn_d_tf_class <- nbn_d_tf_ssm |> ungroup() |> dplyr::select(-c(1)) |> nest_by(classs) |> mutate(sr_nbn = list(vegan::specpool(data))) |> dplyr::select(-data) |> unnest(sr_nbn) nbn_d_tf_class |> dplyr::filter(Species > 5) |> ggplot() + geom_point(aes(reorder(classs, Species), Species)) + geom_segment(aes(x = classs, xend = classs, y = 0, yend = Species), lty = "dotted") + coord_flip() + labs(x = "", title = "Species richness by class", subtitle = "Note: log10 scale", caption = "Source: NBN Atlas") + scale_y_log10() ``` ```{r} #| label: fig-sr-lc-class #| fig-cap: "Species richness by taxon class and land cover" #| fig-subcap: #| - "Class by land cover" #| - "Land cover by class" #| layout-nrow: 2 #| column: page nbn <- read_sf("../data/spatial/nbn_tf_lc.gpkg") #source("startup.R") # nbn <- setDT(nbn) nbn_dt_wide <- nbn[, .N, by = .(species, classs, dscrptn.x)] |> pivot_wider(names_from = species, values_from = N, values_fill = 0) nbn_sr_class <- nbn_dt_wide |> nest_by(dscrptn.x, classs) |> mutate(sr = list(vegan::specpool(data))) |> dplyr::select(-data) |> unnest("sr") nbn_sr_class |> write_csv("../data/nbn_sr_class.csv") # nbn_sr_class |> # kableExtra::kable() nbn_sr_class |> dplyr::filter(!is.na(classs), Species > 9, !str_detect(dscrptn.x, "Scattered|Orch|Qua|Rough|Marsh|Mixed|^Scr")) |> top_n(5, Species) |> ggplot() + geom_point(aes(reorder(classs, Species), Species)) + geom_segment(aes(x = classs, y = 0, xend = classs, yend = Species), lty = "dotted") + coord_flip() + facet_wrap(~ dscrptn.x, ncol = 3) + labs(x = "") + scale_y_log10() nbn_sr_class |> dplyr::filter(!is.na(classs), Species > 4, str_detect(classs, "Aves|Mag|Mamm|Rept|Insec|Bryop"), !str_detect(dscrptn.x, "Scattered|Orch|Qua|Rough|Marsh|Mixed|^Scr")) |> top_n(10, Species) |> ggplot() + geom_point(aes(reorder(dscrptn.x, Species), Species)) + geom_segment(aes(x = dscrptn.x, y = 0, xend = dscrptn.x, yend = Species), lty = "dotted") + coord_flip() + facet_wrap(~ classs, ncol = 3) + labs(x = "") + scale_y_log10() ``` As @fig-sr-lc-class-1 shows insects, vasculat plants and birds are the most diverse taxa across all habitat type. Insect and plant diversity are greatest in residentail gardens; bird and mamal diversity are siialr betwen gardens, grasslant and arable land. ## Summary of environmental predictors ```{r} #| label: tbl-env-pred #| tbl-cap: "Environmental predictor summary table" env_table <- summary_df |> dplyr::select(-5) |> flextable::flextable() env_table <- fontsize(env_table, size = 11, part = "header") env_table <- fontsize(env_table, size = 10, part = "body") env_table ``` Many indicators are highly correlated or cluster as shown in @fig-corr_env which orders pair-wise correlations of scaled (normalised) metrics by hierarchical cluster. Rainfall metrics cluster and are inversely correlated with temperature indicators. Gardens, built areas and roads form a separate cluster as do connectivity metrics, and vegetation/'greenness' indicators. ```{r} #| label: fig-corr_env #| fig-cap: "Correlogram of environmental predictors" #| fig-subcap: "Non-significant correlation are masked by an X" #| fig-height: 9 #| fig-width: 9 env_pred_scale <- env_pred |> mutate_all(scale) env_pred_scale <- env_pred_scale |> dplyr::select(-site) p.mat <- cor_pmat(env_pred_scale) env_pred_scale |> cor() |> ggcorrplot(hc.order = TRUE, outline = "white", type = "lower", show.diag = NULL, p.mat = p.mat, ggtheme = ggthemes::theme_base(), colors = c("red", "white", "darkgreen"), tl.cex = 10) + theme(plot.background = element_blank(), legend.text = element_text(size = 8), legend.title = element_text(size = 10), axis.text.x = element_text(angle = 90, hjust = 1), panel.background = element_rect(colour = "white")) ``` ```{r} #| label: fig-bio_env #| fig-cap: "Correlogram of bio predictors" #| fig-subcap: "Non-significant correlation are masked by an X" #| fig-height: 9 #| fig-width: 9 biotic <- fread("../data/biotic_pred.csv") biotic_scale <- scale(biotic[, -1]) #summary(biotic) p.mat_bio <- cor_pmat(biotic_scale) biotic_scale |> cor() |> ggcorrplot(hc.order = TRUE, outline = "white", type = "lower", show.diag = NULL, p.mat = p.mat_bio, ggtheme = ggthemes::theme_base(), colors = c("red", "white", "darkgreen"), tl.cex = 10) + theme(plot.background = element_blank(), legend.text = element_text(size = 8), legend.title = element_text(size = 10), axis.text.x = element_text(angle = 90, hjust = 1), panel.background = element_rect(colour = "white")) # b_cols <- colnames(biotic)[-1] # get column names of "env_pred_dt" excluding the first column # # b_means <- biotic[, map(.SD, mean), .SDcols = b_cols] # calculate the mean of each column in "env_pred_dt" # b_sd <- biotic[, lapply(.SD, sd), .SDcols = b_cols] # calculate the standard deviation of each column in "env_pred_dt" # b_iqr <- biotic[, map(.SD, IQR), .SDcols = b_cols] # calculate the interquartile range of each column in "env_pred_dt" # b_qnt <- biotic[, map(.SD, \(x) quantile(x, probs = c(0.25, 0.5, 0.75, 1))), .SDcols = b_cols] # calculate the 25th and 75th percentiles of each column in "env_pred_dt" # # summary_bio_df <- bind_rows(b_means, b_sd, b_iqr, b_qnt) |> setDF() # combine the mean, standard deviation, interquartile range, and percentiles into a single data frame # rownames(summary_bio_df) <- c('mean', 'sd', 'iqr', 'q25', 'median', 'q75', 'q100') # # summary_bio_df <- summary_bio_df |> # pipe the summary statistics data frame # rownames_to_column('metric') |> # convert the row names to a separate column named "metric" # relocate(metric) |> # move the "metric" column to the front of the data frame # t() |> # transpose the data frame # data.frame() # # colnames(summary_bio_df) <- summary_bio_df |> slice(1) # set the column names of the summary statistics data frame to be the first row # summary_bio_df <- summary_bio_df |> slice(-1) # remove the first row from the summary statistics data frame # summary_bio_df <- rownames_to_column(summary_bio_df, 'metric') # # colnames(summary_bio_df) <- str_to_title(colnames(summary_bio_df)) # # summary_bio_df |> # mutate_at(vars(2:8), as.numeric) |> # dplyr::filter(Q25 > 0.5) |> # arrange(-Mean) |> # flextable() sr_db <- fread('../data/richness_database.csv') |> dplyr::select(classs, Species) ``` ```{r} #| label: ground-dweller-sr ground_dwellers_1 <- ground_dwellers |> filter(!is.na(qty)) |> mutate(event_date = dmy(event_date)) |> inner_join(tfb_coords, by = c("tiny_forest_id" = "tf_id")) |> #inner_join(lc, by = c("tiny_forest_id" = "site")) |> inner_join(tree_shannon) |> mutate(age = as.numeric(event_date - plant_date)) |> dplyr::select(tiny_forest_id, survey_id, qty, month, x, y, shannon, mean_ht, age, species) gd_survey_counts <- ground_dwellers_1 |> group_by(tiny_forest_id, month) |> reframe(n = n(), mean_age = mean(age), Y = mean(y), shannon = mean(shannon), mean_ht = mean(mean_ht)) srgd <- ground_dwellers_1 |> count(tiny_forest_id, month, species) |> pivot_wider(names_from = species, values_from = n, values_fill = 0) |> dplyr::select(-c(1:2)) |> vegan::specnumber() srgd_pool <- ground_dwellers_1 |> count(tiny_forest_id, month, species) |> pivot_wider(names_from = species, values_from = n, values_fill = 0) |> dplyr::select(-c(1:2)) |> vegan::poolaccum() srgd <- srgd |> bind_cols(gd_survey_counts) |> rename(sr = 1) ``` ```{r} #| label: pollinator-sr poll_1 <- pollinators |> filter(!is.na(qty)) |> mutate(event_date = dmy(event_date)) |> inner_join(tfb_coords, by = c("tiny_forest_id" = "tf_id")) |> #inner_join(lc, by = c("tiny_forest_id" = "site")) |> inner_join(tree_shannon) |> mutate(age = as.numeric(event_date - plant_date)) |> dplyr::select(tiny_forest_id, survey_id, qty, month, x, y, shannon, mean_ht, age, species) poll_survey_counts <- poll_1 |> group_by(tiny_forest_id, month) |> reframe(n = n(), mean_age = mean(age), Y = mean(y), shannon = mean(shannon), mean_ht = mean(mean_ht)) srpol <- poll_1 |> count(tiny_forest_id, month, species) |> pivot_wider(names_from = species, values_from = n, values_fill = 0) |> dplyr::select(-c(1:2)) |> vegan::specnumber() srpol_pool <- poll_1 |> count(tiny_forest_id, month, species) |> pivot_wider(names_from = species, values_from = n, values_fill = 0) |> dplyr::select(-c(1:2)) |> vegan::poolaccum() srpol <- srpol |> bind_cols(poll_survey_counts) |> rename(sr = 1) ``` ## Gamma diversity Results ```{r} #| label: correlations gsr_cor <- srgd |> select(tiny_forest_id, sr) |> inner_join( srpol |> dplyr::select(tiny_forest_id, sr), by = "tiny_forest_id") |> inner_join( srb |> dplyr::select(tiny_forest_id, sr), by = "tiny_forest_id") |> group_by(tiny_forest_id) |> summarise(mean_s = mean(sr, na.rm = TRUE), mean_sx = mean(sr.x, na.rm = TRUE), mean_sy = mean(sr.y, na.rm = TRUE)) |> select(-1) |> cor() ``` ```{r} #| label: richness-database srdb <- fread("../data/richness_database.csv") mean_sr <- srdb[str_detect(class, "Ground|Pollin|Butterfli"), .(mean_sr = mean(Species), max_sr = max(Species), min_sr = min(Species), sites = n_distinct(tf_id)), , by = .(class) ] ``` Mean TF butterfly SR was `r mean(srb$sr)` (survey sites = `r mean_sr[3, 5]`) and varied between `r min(srb$sr)` and `r max(srb$sr)`. Taxon richness also varied between `r min(srgd$sr)` and `r max(srgd$sr)` (mean = `r mean(srgd$sr)`, sites = `r mean_sr[1,5]`) for pollinators and `r min(srpol$sr)` to `r max(srpol$sr)` (mean = `r mean(srpol$sr)`, sites = `r mean_sr[2,5]`) for ground dwelling invertebrates respectively. ```{r} #| label: fig-cor-but #| fig-cap: "Correlation between TF butterfly SR and insect and butterfly TFb SR" srdb[str_detect(class, "Butterfl|Insect"), .(tf_id, class, Species)] |> pivot_wider(names_from = class, values_from = Species, values_fn = \(x) mean(x)) |> filter(!is.na(Butterflies), Butterflies < 10) |> ggpairs(columns = 2:5, lower = list(continuous = "smooth"), progress = FALSE) ``` As @fig-cor-but shows, TF butterfly richness is uncorrelated with TFb butterfly or insect SR. Intrinsic TF variables were then fitted to butterfly SR using GLMM with a negative binomial error structure and TF as a random effect. The 10 models with the lowest AICc fitted with the `dredge` function of the `MuMIm` R package are shown in @tbl-butterfly-model (@MuMIn). Each row of the table shows the coefficients of the explanatory variables - the number of surveys (*n*) is included in all 10 models, with mean age at time of survey (scaled) included in four. Survey month and weighted tree height are included in 2 candidate models. Model selection can be difficult. It has been common practice in ecology to average models with with an AICc within 7 of the model of the lowest AICc but this has been criticised for leading to faulty inferences about relevant predictors (@galipaud2014; @walker2017). An alternative is use penalised regression where predictors are removed ```{r} #| label: tbl-butterfly-model #| tbl-cap: "Top 10 butterfly models" #| eval: false ## model options(na.action = "na.fail", digits = 2) null_glm <- glm(sr ~ 1, data = srb, family = poisson ) null_mod <- glmer.nb(sr ~ 1 + (1|tiny_forest_id), data = srb) b_mod <- glmer.nb(sr ~ scale(month) + n + scale(mean_age) + scale(Y) + scale(mean_ht) + scale(shannon) + (1|tiny_forest_id) , data = srb) dredge(b_mod) |> slice(1:10) |> data.frame() |> dplyr::select(-c(df, logLik, delta, weight)) |> flextable() ``` ```{r} #| label: fig-butterfly-model #| fig-cap: "Butterfly GAM" #| fig-width: 8 #| fig-height: 6 #| eval: false b_mod_g <- gam(sr ~ month + n + s(mean_age) + s(Y) + s(shannon) + s(mean_ht) + s(tiny_forest_id, bs = "re"), data = srb, family = negbin(1)) b <- mgcViz::getViz(b_mod_g) plot(b_mod_g, pages = 1, all.terms = F, shift = coef(b_mod_g)[1], shade = F) ``` ### Ground dwelling invertebrates ```{r} #| label: tbl-gd-model #| tbl-cap: "GD models" #| eval: false # srgd$sr |> # histogram() ## model options(na.action = "na.fail", digits = 2) null_glm_gd <- glm(sr ~ 1, data = srgd, family = poisson ) #summary(null_glm_gd) null_mod_gd <- glmer.nb(sr ~ 1 + (1|tiny_forest_id), data = srgd) #summary(null_mod_gd) gd_mod <- glmer.nb(sr ~ scale(month) + n + scale(mean_age) + scale(Y) + scale(mean_ht) + scale(shannon) + (1|tiny_forest_id), data = srgd) gd_models <- dredge(gd_mod) gd_models |> slice(1:10) |> data.frame() |> dplyr::select(-c(df, logLik, delta, weight)) |> flextable() ``` ```{r} #| label: fig-gd-model #| fig-cap: "Ground dwellers GAM" #| fig-height: 6 #| fig-width: 8 #| eval: false g_mod_g <- gam(sr ~ month + n + s(mean_age) + s(Y) + s(mean_ht) + s(shannon) + + s(tiny_forest_id, bs = "re"), data = srgd, family = negbin(1)) #summary(g_mod_g) gg <- mgcViz::getViz(g_mod_g) plot(g_mod_g, pages = 1, all.terms = F, shift = coef(g_mod_g)[1], shade = F) ``` ```{r} #| label: tbl-poll-model #| tbl-cap: "Pollinator models" #| eval: false # srpol$sr |> # histogram() #summary(srpol) ## model options(na.action = "na.fail", digits = 2) null_glm_p <- glm(sr ~ 1, data = srpol, family = poisson ) #summary(null_glm_p) null_mod_p <- glmer.nb(sr ~ 1 + (1|tiny_forest_id), data = srpol) #summary(null_mod_p) p_mod <- glmer.nb(sr ~ scale(month) + n + scale(mean_age) + scale(Y) + scale(mean_ht) + scale(shannon) + (1|tiny_forest_id), data = srpol) dredge(p_mod) |> slice(10) |> data.frame() |> dplyr::select(-c(df, logLik, delta, weight)) |> flextable() ``` ```{r} #| label: fig-p-model #| fig-cap: "Pollinators GAM" #| fig-width: 8 #| fig-height: 6 #| eval: false p_mod_g <- gam(sr ~ month + n + s(mean_age) + s(Y) + s(mean_ht) + s(shannon) + + s(tiny_forest_id, bs = "re"), data = srpol, family = negbin(1)) #summary(p_mod_g) gp <- mgcViz::getViz(p_mod_g) plot(p_mod_g, pages = 1, all.terms = F, shift = coef(g_mod_g)[1], shade = F) #print(plot(gp, allTerms = T), pages = 1) ``` ### GDM Generalised dissimilarity models (GDM) were fitted to community matrices generated for each TF site for butterflies, pollinators and ground dwelling invertebrates to test the hypothesis that community composition changed as forest aged. Modelling followed the method set out by @mokany2022. ```{r} #| label: fig-gdm-butterflies #| fig-cap: "Partial dependence plots for butterfly GDM" #| fig-height: 6 #| fig-width: 8 needs(gdm) butterfly_spp <- butterfly_1[["butterfly_species"]] |> unique() |> enframe() butterfly_gdm <- butterfly_1 |> dplyr::select(species = butterfly_species, site = tiny_forest_id, x, y) |> left_join(butterfly_spp, by = c("species" = "value")) |> relocate(name) |> dplyr::select(-species) |> #mutate(name = factor(name)) |> arrange(name) |> relocate(site) |> mutate(site = as.character(site)) butterfly_env <- butterfly |> filter(!is.na(qty)) |> mutate(event_date = dmy(event_date)) |> inner_join(tfb_coords, by = c("tiny_forest_id" = "tf_id")) |> #inner_join(lc, by = c("tiny_forest_id" = "site")) |> inner_join(tree_shannon) |> mutate(age = as.numeric(event_date - plant_date)) |> dplyr::select(site = tiny_forest_id, x, y, shannon, mean_ht, tf_age) |> mutate(site = as.character(site)) |> distinct() butterfly_gdm <- butterfly_gdm |> filter(site %in% unique(butterfly_env$site)) #summary(butterfly_gdm) sitePairTab <- formatsitepair(bioData = butterfly_gdm, bioFormat = 2 , XColumn="x", YColumn="y", sppColumn="name", siteColumn="site", predData= butterfly_env) gd_b <- gdm::gdm(sitePairTab, geo = TRUE) b_varImm <- gdm::gdm.varImp(sitePairTab, geo = TRUE) plot(gd_b, plot.layout = c(3,2)) ``` ```{r} #| label: gdm-pollinators #| fig-height: 10 needs(gdm) poll_spp <- poll_1[["species"]] |> unique() |> enframe() poll_gdm <- poll_1 |> dplyr::select(species = species, site = tiny_forest_id, x, y, qty) |> left_join(poll_spp, by = c("species" = "value")) |> relocate(name) |> dplyr::select(-species) |> #mutate(name = factor(name)) |> arrange(name) |> relocate(site) |> mutate(site = as.character(site)) poll_env <- pollinators |> filter(!is.na(qty)) |> mutate(event_date = dmy(event_date)) |> inner_join(tfb_coords, by = c("tiny_forest_id" = "tf_id")) |> #inner_join(lc, by = c("tiny_forest_id" = "site")) |> inner_join(tree_shannon) |> mutate(age = as.numeric(event_date - plant_date)) |> dplyr::select(site = tiny_forest_id, x, y, shannon, mean_ht, tf_age) |> mutate(site = as.character(site)) |> distinct() poll_gdm <- poll_gdm |> filter(site %in% unique(butterfly_env$site)) #summary(poll_gdm) gdm_pol <- formatsitepair(bioData = poll_gdm, bioFormat = 2 , XColumn="x", YColumn="y", sppColumn="name", siteColumn="site", abundColumn = "qty", predData= poll_env) gd_pol <- gdm::gdm(gdm_pol, geo = TRUE) pol_varImm <- gdm::gdm.varImp(gdm_pol, geo = TRUE) plot(gd_pol, plot.layout = c(3,2)) ``` ```{r} #| label: gdm-gd #| eval: true needs(gdm) gd_spp <- ground_dwellers_1[["species"]] |> unique() |> enframe() gd_gdm <- ground_dwellers_1 |> dplyr::select(species = species, site = tiny_forest_id, x, y, qty) |> left_join(gd_spp, by = c("species" = "value")) |> relocate(name) |> dplyr::select(-species) |> #mutate(name = factor(name)) |> arrange(name) |> relocate(site) |> mutate(site = as.character(site)) gd_env <- ground_dwellers |> filter(!is.na(qty)) |> mutate(event_date = dmy(event_date)) |> inner_join(tfb_coords, by = c("tiny_forest_id" = "tf_id")) |> #inner_join(lc, by = c("tiny_forest_id" = "site")) |> inner_join(tree_shannon) |> mutate(age = as.numeric(event_date - plant_date)) |> dplyr::select(site = tiny_forest_id, x, y, shannon, mean_ht, tf_age) |> mutate(site = as.character(site)) |> distinct() gd_gdm <- gd_gdm |> filter(site %in% unique(gd_env$site)) ``` ```{r} #| label: fig-gd-gdm-1 #| fig-cap: "GDM results for ground-dwellers" #| eval: true gdm_gd <- formatsitepair(bioData = gd_gdm, bioFormat = 2 , XColumn="x", YColumn="y", sppColumn="name", siteColumn="site", abundColumn = "qty", predData= gd_env) gd_gd <- gdm::gdm(gdm_gd, geo = TRUE) gd_varImm <- gdm::gdm.varImp(gdm_gd, geo = TRUE, parallel = TRUE, cores = 6) plot(gd_gd, layout = c(3, 2)) ``` ```{r} #| label: tbl-importance-taxon #| tbl-cap: "Variable importance" bind_rows( b_varImm$`Predictor p-values` |> mutate(taxon = "Butterflies"), pol_varImm$`Predictor p-values` |> mutate(taxon = "Pollinators"), gd_varImm$`Predictor p-values` |> mutate(taxon = "Ground dwellers)" ) ) |> rownames_to_column("Important variables") |> relocate(taxon) |> flextable() ``` Relationship between intrinsic tiny forest variables and tiny forest species, richness Impacts of characteristics of surrounding area on tiny forest, biodiversity line Abiotic. Biotic ```{r} #| label: taxon richness butterfly_ssm <- butterfly |> filter(!is.na(qty)) |> group_by(tiny_forest_id, survey_id, species = butterfly_species) |> count() |> dplyr::select(-2) |> pivot_wider(names_from = species, values_from = c("n"), values_fill = 0) butterfly_richness <- butterfly_ssm |> select(-1) |> ungroup() |> nest_by(tiny_forest_id) |> mutate(sr = list(vegan::specpool(data))) |> dplyr::select(-data) |> unnest(sr) pollinator_ssm <- pollinators |> filter(!is.na(qty)) |> group_by(tiny_forest_id, survey_id, species) |> count() |> dplyr::select(-2) |> pivot_wider(names_from = species, values_from = c("n"), values_fill = 0) pollinator_richness <- pollinator_ssm |> select(-1) |> ungroup() |> nest_by(tiny_forest_id) |> mutate(sr = list(vegan::specpool(data))) |> dplyr::select(-data) |> unnest(sr) gd_ssm <- ground_dwellers |> filter(!is.na(qty)) |> group_by(tiny_forest_id, survey_id, species) |> count() |> dplyr::select(-2) |> pivot_wider(names_from = species, values_from = c("n"), values_fill = 0) gd_richness <- gd_ssm |> select(-1) |> ungroup() |> nest_by(tiny_forest_id) |> mutate(sr = list(vegan::specpool(data))) |> dplyr::select(-data) |> unnest(sr) ``` ```{r} #| label: fig-glmnet #| tbl-cap: "Env" #| fig-cap: "Environmnental predictors of taxon richness" #| fig-subcap: #| - "Butterfly predictors" #| - "Pollinator predictors" #| - "Ground dweller predictors" #| layout-ncol: 3 #| layout-nrow: 1 #| column: page needs(glmnet) set.seed(123) butterfly_richness_env <- butterfly_richness |> rename(site = tiny_forest_id) |> inner_join(env_pred, by = "site") |> distinct() |> dplyr::select(1:2, 10:55) p_richness_env <- pollinator_richness |> rename(site = tiny_forest_id) |> inner_join(env_pred, by = "site") |> distinct() |> dplyr::select(1:2, 10:55) gd_richness_env <- gd_richness |> rename(site = tiny_forest_id) |> inner_join(env_pred, by = "site") |> dplyr::select(1:2, 10:55) glmnet_fit <- function(df, taxon = "Pollinators"){ set.seed(1234) x <- df |> select(-Species) |> as.matrix() y <- df |> select(Species) |> as.matrix() cv.fit <- cv.glmnet(x, y, family = "poisson") best_l <- cv.fit$lambda.min coeffs <- coef(cv.fit, s = best_l) coef_v <- as.matrix(coeffs) |> data.frame() |> rownames_to_column("metric") |> filter(s1 != 0) coef_plot <- coef_v |> filter(metric != "(Intercept)") |> ggplot() + geom_col(aes(reorder(metric, s1), s1 )) + geom_hline(yintercept = 0, colour = "red") + labs(title = "Retained non-zero terms in glmnet model", subtitle = taxon, x = "", y = "Estimate") return(list(plot = coef_plot, df = coef_v)) } b_fit <- glmnet_fit(butterfly_richness_env[, -1], taxon = "Butterflies") p_fit <- glmnet_fit(p_richness_env[, -1], taxon = "Pollinators") gd_fit <- glmnet_fit(gd_richness_env[, -1], taxon = "Ground dwellers") b_fit$df |> mutate(taxon = "Butterfly") |> bind_rows(p_fit$df |> mutate(taxon = "Pollinators"), gd_fit$df |> mutate(taxon = "Ground-dwellers") ) |> filter(metric != "(Intercept)") |> relocate(taxon) |> arrange(taxon, -s1) |> rename(estimate = s1) |> flextable::flextable() ``` The results of `glment` modelling are shown in @fig-glmnet. For butterflies, green infrastructure, linear features (water courses and hedges), wooded areas, rainfall and connectivity metrics variables were found to be weakly positively associated with butterfly SR, whereas tree cover was strongly negatively predictive. For pollinators and ground-dwellers, SR was associated with the number of surveys, and increased with forest age and latitude for ground-dwellers. ```{r} #| label: pollin-mod #| eval: false p_richness_env <- pollinator_richness |> rename(site = tiny_forest_id) |> inner_join(env_pred, by = "site") |> distinct() |> dplyr::select(Species, elevation:last_col()) p_richness_env |> cor() |> corrplot(order = "hclust", tl.col = "black", tl.cex = 0.6) p_richness_mod <- glm(Species ~ ., data = p_richness_env, family = poisson) tidy(p_richness_mod) #p_richness_mod_1 <- glm(Species ~ hedge_len + lsm_l_shdi, data = p_richness_env, family = poisson) # p_mod_sel <- dredge(p_richness_mod) # # data.frame(.pred = exp(predict(p_richness_mod_1, p_richness_env)), .obs = p_richness_env$Species) |> # mutate(rmse = sqrt(mean((.pred - .obs)^2))) ``` ```{r} #| eval: false xp <- p_richness_env |> ungroup() |> select(-Species ) |> as.matrix() yp <- p_richness_env |> ungroup() |> select(Species) |> as.matrix() p.cv <- cv.glmnet(xp, yp, family = "poisson") plot(p.cv) #p.cv$lambda.min #coef(p.cv, s = "lambda.min") p_glm <- glmnet(xp, yp, family = poisson, lambda = p.cv$lambda.min) best_mod_p <- coef(p_glm) as.data.frame(as.matrix(best_mod_p)) |> filter(s0 > 0) |> flextable() #glance(butterfly_richness_mod) rmse_poll <- data.frame(.pred = exp(predict(p_glm, xp)), .obs = p_richness_env$Species) |> mutate(rmse = sqrt(mean((s0 - .obs)^2))) rmse_poll |> ggplot() + geom_point(aes(s0, .obs)) + geom_smooth(aes(s0, .obs), method = "lm") ``` ```{r} #| label: fig-ground-dwellers-mod #| fig-cap: "Correlogram of ground dweller environmental predictors" #| eval: false gd_richness_env <- gd_richness |> rename(site = tiny_forest_id) |> inner_join(env_pred, by = "site") |> distinct() |> dplyr::select(Species, elevation:last_col()) #summary(gd_richness_env) gd_richness_env |> cor() |> corrplot(order = "hclust", tl.col = "black", tl.cex = 0.5, method = "ellipse") gd_richness_mod <- glm(Species ~ ., data = gd_richness_env, family = poisson) tidy(gd_richness_mod) # g_richness_mod_1 <- glm(Species ~ tf_age.y + site_tree_ht.x + hedge_len + lsm_l_shdi + mean_spring_rain, data = gd_richness_env, family = poisson) # #tidy(g_richness_mod_1, conf.int = TRUE) # dredge(g_richness_mod_1) |> # slice(1:10) |> # kableExtra::kable() # ``` ```{r} #| label: glmnet-gd #| eval: false # data.frame(.pred = exp(predict(g_richness_mod_1, gd_richness_env)), .obs = gd_richness_env$Species) |> # mutate(rmse = sqrt(mean((.pred - .obs)^2))) xg <- gd_richness_env |> ungroup() |> select(-site) |> as.matrix() yg <- gd_richness_env |> ungroup() |> select(Species) |> as.matrix() g.cv <- cv.glmnet(xg, yg, family = "poisson") plot(g.cv) g.cv$lambda.min coef(g.cv, s = "lambda.min") g_glm <- glmnet(xg, yg, family = poisson, lambda = g.cv$lambda.min) best_mod_g <- coef(g_glm) as.data.frame(as.matrix(best_mod_g)) |> filter(s0 > 0) |> flextable() rmse_gd <- data.frame(.pred = exp(predict(g_glm, xg)), .obs = gd_richness_env$Species) |> mutate(rmse = sqrt(mean((s0 - .obs)^2))) rmse_poll |> ggplot() + geom_point(aes(s0, .obs)) + geom_smooth(aes(s0, .obs), method = "lm") ``` ```{r} #| label: nbn nbn <- read_sf("../data/spatial/nbn_tf_lc.gpkg") nbn_dt <- setDT(nbn) nbn_dt_wide <- nbn_dt[, .N, by = .(species, classs, dscrptn.x)] |> pivot_wider(names_from = species, values_from = N, values_fill = 0) nbn_sr_class <- nbn_dt_wide |> nest_by(dscrptn.x, classs) |> mutate(sr = list(vegan::specpool(data))) |> dplyr::select(-data) |> unnest("sr") ``` ```{r} #| label: fig-sr-lc-class-1 #| fig-cap: "Species richness by taxon class and land cover" #| fig-subcap: #| - "Class by land cover" #| - "Land cover by class" #| layout-nrow: 2 #| column: page nbn_sr_class |> filter(!is.na(classs), Species > 4, !str_detect(dscrptn.x, "Scattered|Orch|Qua|Rough|Marsh|Mixed|^Scr")) |> ggplot() + geom_point(aes(reorder(classs, Species), Species)) + geom_segment(aes(x = classs, y = 0, xend = classs, yend = Species), lty = "dotted") + coord_flip() + facet_wrap(~ dscrptn.x, ncol = 3) + labs(x = "") + scale_y_log10() nbn_sr_class |> filter(!is.na(classs), Species > 4, str_detect(classs, "Aves|Mag|Mamm|Rept|Insec|Bryop"), !str_detect(dscrptn.x, "Scattered|Orch|Qua|Rough|Marsh|Mixed|^Scr")) |> ggplot() + geom_point(aes(reorder(dscrptn.x, Species), Species)) + geom_segment(aes(x = dscrptn.x, y = 0, xend = dscrptn.x, yend = Species), lty = "dotted") + coord_flip() + facet_wrap(~ classs, ncol = 3) + labs(x = "") + scale_y_log10() ``` ### Estimated SR for small areas based on random sampling ### Estimated SR for small areas based on SAR ```{r} os <- read_sf("/Users/julianflowers/Library/CloudStorage/GoogleDrive-julian.flowers12@gmail.com/My Drive/dissertation/data/os_data.gpkg") nbn_data_dt <- read_sf("../data/spatial/nbn_sf_tf_dedup.gpkg") nbn_os <- nbn_data_dt |> st_join(os, st_intersects) ``` ```{r} #| label: fig-sp-area #| fig-cap: "Relationship between log habitat area and log SR" nbn_os_ssm <- nbn_os |> st_drop_geometry() |> group_by(dscrptn, osid, species) |> summarise(spp = n(), sum_area = sum(p_area)) |> pivot_wider(names_from = species, values_from = c("spp", "sum_area"), values_fill = 0) nbn_sr_lc <- nbn_os_ssm |> ungroup() |> dplyr::select(dscrptn, contains("spp")) |> nest_by(dscrptn) |> mutate(sr = list(vegan::specpool(data))) |> dplyr::select(-data) |> unnest(sr) nbn_plot_sr <- nbn_os |> st_drop_geometry() |> group_by(dscrptn) |> reframe(sum_area = sum(p_area)) |> left_join(nbn_sr_lc) |> ggplot() + geom_point(aes(sum_area , Species)) + geom_smooth(aes(sum_area, Species), method = "lm") + geom_vline(xintercept = 40000, lty = "dotted") + geom_hline(yintercept = 10, lty = "dotted") + scale_x_continuous(trans = 'log10') + scale_y_continuous(trans = 'log10', n.breaks = 10) nbn_plot_chao <- nbn_os |> st_drop_geometry() |> group_by(dscrptn) |> reframe(sum_area = sum(p_area)) |> left_join(nbn_sr_lc) |> ggplot() + geom_point(aes(sum_area , chao)) + stat_smooth(aes(sum_area, chao), method = "lm") + geom_vline(xintercept = 40000, lty = "dotted") + geom_hline(yintercept = 20, lty = "dotted") + scale_x_continuous(trans = 'log10') + scale_y_continuous(trans = 'log10', n.breaks = 10) nbn_plot_sr nbn_plot_chao ``` ```{r} nbn_sr_osid <- nbn_os_ssm |> ungroup() |> select(osid, contains("spp")) |> specnumber() |> bind_cols(nbn_os_ssm$osid) |> rename(spp =...1, osid = ...2) nbn_os |> st_drop_geometry() |> count(osid, p_area) |> left_join(nbn_sr_osid) |> inner_join(os) |> filter(!str_detect(dscrptn, "Constr|Rough Grass|Scattered|Quarry|Marsh|Tennis")) |> ggplot() + geom_point(aes(p_area, spp)) + geom_smooth(aes(p_area, spp), method = "lm", se = FALSE) + scale_y_log10() + scale_x_log10() + facet_wrap(~ dscrptn) ``` ### Sampling ```{r} needs(stars, sf, terra, exactextractr) # raster ceh <- rast("../data/spatial/gblcm10m2021.tif") ## lc values for TFb tfbp <- tfb |> st_union() |> st_cast("MULTIPOINT") |> st_cast("POINT") ## combine buffers into a single polygon and sample 10000 points and create 112m buffers around each tfb_samples <- tfb |> st_union() |> # st_cast("MULTIPOINT") |> # st_cast("POINT") |> st_sample(10000) |> st_buffer(sqrt(40000 / pi)) tfb_samples_id <- tfb_samples |> bind_cols(sampleID = 1:10000) ## check sample distribution - should be normal distn of samples per TF and they should be randomly sampled across the whole area nn <- st_join(tfb_samples |> st_as_sf(), tfb, left = FALSE) |> count(tf_id) med <- quantile(nn$n, probs = c(.025, .5, .975)) nn |> ggplot() + geom_density(aes(n), fill = "grey90") + geom_vline(xintercept = med[2], colour = "red") + geom_vline(xintercept = med[1], colour = "red", lty = "dotted") + geom_vline(xintercept = med[3], colour = "red", lty = "dotted") ## extract landcover point values for each sample ceh_tf <- exact_extract(ceh, tfb_samples, 'mode') tfb_lc <- ceh_tf |> cbind(tfb_samples_id) ``` ```{r} sub_smpl <- tfb_lc |> filter(mode.gblcm10m2021_1 %in% c(1, 2, 3, 4, 8, 14, 20, 21)) |> rename(lc_id = 1) |> left_join(ceh_lu) ``` ```{r} #| label: fig-sac-4Ha #| fig-cap: "Species accumulation curve for randomly sampled 4Ha sites" #| fig-subcap: #| - "All species" #| - "Butterflies" tot_test_tot <- st_join(nbn_data_tf_sf |> st_transform(27700), sub_smpl |> st_as_sf(), left = FALSE) tot_pool <- tot_test_tot |> setDT() tot_pool <- tot_pool[, .N, keyby = .(sampleID, land_cover, species)] per_sample_sn <- tot_pool |> pivot_wider(names_from = species, values_from = N, values_fill = 0) |> select(-land_cover) |> nest_by(sampleID) |> mutate(srr = list(vegan::specpool(data))) |> select(-data) |> unnest("srr") per_lc_sn <- tot_pool |> pivot_wider(names_from = species, values_from = N, values_fill = 0) |> #select(-`NA`) |> nest_by(land_cover) |> mutate(srr = list(vegan::specaccum(data))) |> select(-data) quants <- quantile(per_sample_sn$Species, probs = c(0, 0.025, .25, .5, .75, .975, 1)) meanSR <- mean(per_sample_sn$Species) samples <- map(1:7, \(x) per_lc_sn$srr[[x]]$sites) sr <- map(1:7, \(x) per_lc_sn$srr[[x]]$richness) #sd - map(1:7, \(x) per_lc_sn$srr[[x]]$sd) sr_df <- map_dfr(1:7, \(x) bind_cols(site = samples[[x]], sr = sr[[x]]) |> add_row(site = 0, sr = 0, .before = 1) |> mutate(lc = per_lc_sn$land_cover[x])) sr_df |> ggplot() + geom_line(aes(site, sr, group = lc, colour = factor(lc))) ``` ```{r, eval=FALSE} tot_test_tot_but <- st_join(nbn_data_tf_sf |> filter(family %in% c("Papilionidae", "Pieridae", "Hesperiidae", "Rionidae", "Nymphalidae", "Lycaenidae")) |> st_transform(27700), sub_smpl |> st_as_sf(), left = FALSE) tot_pool_bu <- tot_test_tot_but |> setDT() tot_pool_bu <- tot_pool_bu[, .N, keyby = .(sampleID, land_cover, species)] per_sample_sn_bu <- tot_pool_bu |> pivot_wider(names_from = species, values_from = N, values_fill = 0) |> select(-land_cover) |> nest_by(sampleID) |> mutate(srr = list(vegan::specpool(data))) |> select(-data) |> unnest("srr") per_sample_sn_bu_sppool <- tot_pool_bu |> pivot_wider(names_from = species, values_from = N, values_fill = 0) |> select(-land_cover) |> poolaccum() max_SR <- data.frame(per_sample_sn_bu_sppool$means) |> filter(S == max(S)) per_lc_sn_bu_sp <- tot_pool_bu |> pivot_wider(names_from = species, values_from = N, values_fill = 0) |> #select(-`NA`) |> nest_by(land_cover) |> mutate(srr = list(vegan::specaccum(data))) |> select(-data) quants_bu <- quantile(per_sample_sn_bu$Species, probs = c(0, 0.025, .25, .5, .75, .975, 1)) meanSR_bu <- mean(per_sample_sn_bu$Species) ggplot(data = per_sample_sn_bu) + geom_density(aes(jack1)) + geom_vline(xintercept = quants_bu[4], colour = "red") + geom_vline(xintercept = quants_bu[1], colour = "red", lty = "dotted") + geom_vline(xintercept = quants_bu[6], colour = "red", lty = "dotted") + geom_vline(xintercept = 22) + geom_vline(xintercept = 37) + scale_x_continuous(trans = "log10") ``` ```{r, eval=FALSE} samples_bu <- map(c(1:2, 4:7), \(x) per_lc_sn_bu_sp$srr[[x]]$sites) sr_bu <- map(1:7, \(x) per_lc_sn_bu_sp$srr[[x]]$richness) #sd - map(1:7, \(x) per_lc_sn$srr[[x]]$sd) ``` ```{r, eval=FALSE} sr_df_bu <- map_dfr(1:7, \(x) bind_cols(site = samples_bu[[x]], sr = sr_bu[[x]]) |> add_row(site = 0, sr = 0, .before = 1) |> mutate(lc = per_lc_sn$land_cover[x])) sr_df_bu |> ggplot() + geom_line(aes(site, sr, group = lc, colour = factor(lc)), lwd = 1.5) + scale_color_viridis_d(option = "turbo") sr_df_bu |> filter(site == 1) ``` ### Butterfly richness estimated by sampling Based on 300 (of 10,000) randomly sample 4Ha pseudo-controls, median butterfly richness was , with a percentile range of Butterfly diversity for TFs was 22 which is at the 95th centile of the distribution. Butterfly richness was highest in suburban areas approaching near saturation at 34 species. Richness in arable areas, improved grassland and deciduous woodland were lower but based on smaller samples.For a single 4Ha site the highest SR was in fen habitat (15 species), with arable (12), suburban (11), grassland (10) and woodland (8).