TACR 3

Introduction

Status and development of European ground squirrel populations in agricultural land of southern Moravia

In Hrusovany and Pavlovice, there are patches that have been divided into habitat categories, such as orchards and vineyards, and categories of vegetation management, such as mowed or ploughed. The density of ground squirrel burrow openings is also recorded for each patch.

Library

# install.packages("pacman")
# writeLines(pacman::p_lib(), "~/Desktop/list_of_R_packages.csv") # to quickly back up packages
# remotes::install_github("ThinkR-open/remedy")

# install.packages("remotes")
# remotes::install_github("statnmap/cartomisc")


pacman::p_load(rio, tidyverse,janitor, data.table, here, stringr, stats, gridExtra, vtable, tableone, skimr,
               ggstats,
               ggstatsplot,
               ggdist,
               gghalves,
               plyr,
               Hmisc,
               RColorBrewer,
               reshape2,
               dplyr,
               cli,
               introdataviz,
               kableExtra 

               
            
         
) # just add needed packages to this line and Pacman will install and load them.

# devtools::install_github("psyteachr/introdataviz")

Data

hru <- import( here::here("data","hrusovany_densities.csv")) %>% 
  janitor::clean_names()%>% 
  mutate_if(is.character, str_to_title) %>% 
  mutate(habitat = ifelse(habitat == "Backyard", "Garden", habitat)) %>% 
  mutate(habitat = ifelse(habitat == "Tree Avenue", "Tree Alley", habitat)) %>% 
  mutate(habitat = ifelse(habitat == "Mowed Lawn", "Short-cut Lawn", habitat)) %>% 
  dplyr::mutate(habitat =  if_else(habitat == "Vineyards", "Vineyard", habitat)) %>% 
  dplyr::mutate(district ="Hrušovany") %>% 
  mutate_if(is.character, as.factor)
  

  summary(hru)

##       district         id                    habitat            management
##  Hrušovany:171   Min.   :  3.0   Orchard         :58   Fallow        : 4  
##                  1st Qu.: 53.5   Vineyard        :49   Mowed         :65  
##                  Median :102.0   Crop Field      :20   No Management :30  
##                  Mean   :101.8   Garden          :13   Ploughed      :69  
##                  3rd Qu.:144.5   Steppe Grassland:10   Ploughed/Mowed: 3  
##                  Max.   :219.0   Shrubland       : 9                      
##                                  (Other)         :12                      
##     size_ha              bo             density      
##  Min.   :0.01464   Min.   :  0.000   Min.   :  0.00  
##  1st Qu.:0.06965   1st Qu.:  0.000   1st Qu.:  0.00  
##  Median :0.13528   Median :  0.000   Median :  0.00  
##  Mean   :0.20014   Mean   :  4.228   Mean   : 21.35  
##  3rd Qu.:0.21022   3rd Qu.:  3.000   3rd Qu.: 27.45  
##  Max.   :1.66754   Max.   :192.000   Max.   :378.10  
## 

  vp <- import( here::here("data","vp_densities.csv"))%>% 
    janitor::clean_names() %>% 
    dplyr::mutate(district = "Pavlovice", .before= "id")%>% 
    mutate_if(is.character, str_to_title) %>% 
    mutate(habitat = ifelse(habitat == "Backyard", "Garden", habitat)) %>% 
    mutate(habitat = ifelse(habitat == "Tree Avenue", "Tree Alley", habitat)) %>% 
    mutate(habitat = ifelse(habitat == "Mowed Lawn", "Short-cut Lawn", habitat)) %>% 
    dplyr::mutate(habitat =  if_else(habitat == "Vineyards", "Vineyard", habitat)) %>% 
    mutate_if(is.character, as.factor)

 summary(vp)

##       district         id                  habitat             management 
##  Pavlovice:840   Min.   :  1.0   Vineyard      :286   Fallow        : 27  
##                  1st Qu.:228.8   Orchard       :213   Grazed        :  4  
##                  Median :458.5   Crop Field    :102   Mowed         :344  
##                  Mean   :461.9   Garden        : 68   No Management :123  
##                  3rd Qu.:688.2   Short-cut Lawn: 63   Ploughed      :307  
##                  Max.   :940.0   Shrubland     : 54   Ploughed/Mowed: 35  
##                                  (Other)       : 54                       
##     size_ha              bo             density      
##  Min.   : 0.0174   Min.   : 0.0000   Min.   : 0.000  
##  1st Qu.: 0.1105   1st Qu.: 0.0000   1st Qu.: 0.000  
##  Median : 0.2215   Median : 0.0000   Median : 0.000  
##  Mean   : 0.5068   Mean   : 0.5667   Mean   : 1.368  
##  3rd Qu.: 0.4696   3rd Qu.: 0.0000   3rd Qu.: 0.000  
##  Max.   :25.4100   Max.   :30.0000   Max.   :77.700  
## 

 colonies <-  rbind(hru, vp) %>% 
   dplyr::mutate(district = as.factor(district))
 
 summary(colonies)

##       district         id                  habitat             management 
##  Hrušovany:171   Min.   :  1.0   Vineyard      :335   Fallow        : 31  
##  Pavlovice:840   1st Qu.:142.0   Orchard       :271   Mowed         :409  
##                  Median :366.0   Crop Field    :122   No Management :153  
##                  Mean   :401.0   Garden        : 81   Ploughed      :376  
##                  3rd Qu.:645.5   Short-cut Lawn: 71   Ploughed/Mowed: 38  
##                  Max.   :940.0   Shrubland     : 63   Grazed        :  4  
##                                  (Other)       : 68                       
##     size_ha               bo             density       
##  Min.   : 0.01464   Min.   :  0.000   Min.   :  0.000  
##  1st Qu.: 0.10368   1st Qu.:  0.000   1st Qu.:  0.000  
##  Median : 0.19892   Median :  0.000   Median :  0.000  
##  Mean   : 0.45492   Mean   :  1.186   Mean   :  4.748  
##  3rd Qu.: 0.41384   3rd Qu.:  0.000   3rd Qu.:  0.000  
##  Max.   :25.41000   Max.   :192.000   Max.   :378.100  
## 

 # Presence only
 
 hru_presence <- hru %>% 
   dplyr::filter(bo>0)
 
 vp_presence <- vp %>% 
      dplyr::filter(bo>0)
 
 colonies_pres <- colonies %>% 
   dplyr::filter(bo>0)
 
 
 plot_hru <- import(here::here("data","plot_size_hru.csv")) %>% 
   clean_names() %>% 
   dplyr::mutate(district = "Hrusovany", .before= "id")
 
 
 plot_pa <- import(here::here("data","plot_size_pa.csv")) %>% 
   clean_names() %>% 
   dplyr::mutate(district = "Pavlovice", .before= "id")
 
  plot_data <-  rbind(plot_hru, plot_pa) %>% 
   dplyr::mutate(district = as.factor(district))

Explore data

Size of agricultural plots

# Create violin plot for Hrusovany
ggplot(plot_hru, aes(x=1, y=size_ha)) + 
  geom_violin() +
  geom_text(aes(label = paste("mean = ",round(mean(size_ha), 3),", n = ",length(size_ha))),
            y=0.05, nudge_y= 0.05, size=3) +
  ggtitle("Hrusovany")+
  theme_classic()

# Create violin plot for Pavlovice


ggplot(plot_pa, aes(x=1, y=size_ha)) + 
  geom_violin() +
  geom_text(aes(label = paste("mean = ",round(mean(size_ha), 3),", n = ",length(size_ha))),
            y= -0.5, nudge_y= 0.05, size=3) +
  ggtitle("Pavlovice") +
  theme_classic()

ggbetweenstats(plot_data, district, size_ha)

Density

table1 <- colonies %>% 
  dplyr::group_by(district) %>% 
  dplyr::summarize(mean_density = mean(density),
            median_density = median(density),
            min_density = min(density),
            max_density = max(density),
              sd_density = sd(density),
            n = n())

# Modify the column names to capitalize the first letter and remove hyphens
colnames(table1) <- str_to_title(gsub("_", " ", colnames(table1)))
  

kable(table1, caption = "Summary statistics per district") %>% 
   kableExtra::kable_styling(font_size = 10,
                latex_options = "striped")

Summary statistics per district

District	Mean Density	Median Density	Min Density	Max Density	Sd Density	N
Hrušovany	21.351462	0	0	378.1	46.744447	171
Pavlovice	1.367738	0	0	77.7	6.045053	840

table1 <- colonies %>% 
  dplyr::group_by(district, habitat) %>% 
  dplyr::summarize(mean_density = mean(density),
            median_density = median(density),
            min_density = min(density),
            max_density = max(density),
                 sd_density = sd(density),
            n = n())

# Modify the column names to capitalize the first letter and remove hyphens
colnames(table1) <- str_to_title(gsub("_", " ", colnames(table1)))

kable(table1, caption = "Summary statistics per district and habitat") %>% 
  kableExtra::kable_styling(font_size = 10,
                latex_options = "striped")

Summary statistics per district and habitat

District	Habitat	Mean Density	Median Density	Min Density	Max Density	Sd Density	N
Hrušovany	Crop Field	18.2500000	1.95	0	164.9	38.146966	20
Hrušovany	Garden	26.9615385	0.00	0	267.4	72.966277	13
Hrušovany	Orchard	26.4672414	10.40	0	378.1	54.191206	58
Hrušovany	Short-cut Lawn	18.6375000	8.95	0	63.5	23.653628	8
Hrušovany	Shrubland	3.3888889	0.00	0	22.2	7.569420	9
Hrušovany	Steppe Grassland	11.2400000	0.00	0	56.2	23.696001	10
Hrušovany	Tall Lawn	0.0000000	0.00	0	0.0	0.000000	3
Hrušovany	Tree Stands	0.0000000	0.00	0	0.0	NA	1
Hrušovany	Vineyard	22.6224490	4.90	0	261.9	44.509583	49
Pavlovice	Crop Field	1.5803922	0.00	0	47.1	6.713722	102
Pavlovice	Garden	0.4323529	0.00	0	12.8	2.079139	68
Pavlovice	Orchard	1.6206573	0.00	0	39.5	5.238513	213
Pavlovice	Short-cut Lawn	3.0111111	0.00	0	77.7	13.019726	63
Pavlovice	Shrubland	0.0000000	0.00	0	0.0	0.000000	54
Pavlovice	Steppe Grassland	0.0000000	0.00	0	0.0	0.000000	15
Pavlovice	Tall Lawn	1.9555556	0.00	0	19.9	5.047565	18
Pavlovice	Tree Stands	0.0000000	0.00	0	0.0	0.000000	7
Pavlovice	Vineyard	1.3377622	0.00	0	51.3	5.501825	286
Pavlovice	Pasture	0.0000000	0.00	0	0.0	0.000000	3
Pavlovice	Ruderal	0.0000000	0.00	0	0.0	0.000000	3
Pavlovice	Tree Alley	0.7000000	0.00	0	4.6	1.614222	8

table1 <- colonies %>% 
  dplyr::group_by(district, management) %>% 
    dplyr::summarize(mean_density = mean(density),
            median_density = median(density),
            min_density = min(density),
            max_density = max(density),
                sd_density = sd(density),
            n = n())

# Modify the column names to capitalize the first letter and remove hyphens
colnames(table1) <- str_to_title(gsub("_", " ", colnames(table1)))

kable(table1, caption = "Summary statistics per district and management type") %>% 
   kableExtra::kable_styling(font_size = 10,
                latex_options = "striped")

Summary statistics per district and management type

District	Management	Mean Density	Median Density	Min Density	Max Density	Sd Density	N
Hrušovany	Fallow	48.9000000	15.35	0	164.9	77.738064	4
Hrušovany	Mowed	17.1800000	4.10	0	117.1	24.739185	65
Hrušovany	No Management	6.0533333	0.00	0	56.2	15.476093	30
Hrušovany	Ploughed	30.8840580	7.10	0	378.1	65.208567	69
Hrušovany	Ploughed/Mowed	8.7333333	10.30	0	15.9	8.064945	3
Pavlovice	Fallow	0.0000000	0.00	0	0.0	0.000000	27
Pavlovice	Mowed	1.5680233	0.00	0	77.7	6.712902	344
Pavlovice	No Management	0.3365854	0.00	0	19.9	2.039418	123
Pavlovice	Ploughed	1.7785016	0.00	0	51.3	6.841090	307
Pavlovice	Ploughed/Mowed	0.6314286	0.00	0	7.2	1.722988	35
Pavlovice	Grazed	0.0000000	0.00	0	0.0	0.000000	4

# Habitat and management
table1 <- colonies %>% 
  dplyr::group_by(habitat, management) %>% 
    dplyr::summarize(mean_density = mean(density),
            median_density = median(density),
            min_density = min(density),
            max_density = max(density),
              sd_density = sd(density),
            n = n())

# Modify the column names to capitalize the first letter and remove hyphens
colnames(table1) <- str_to_title(gsub("_", " ", colnames(table1)))

kable(table1, caption = "Summary statistics per habitat and management type") %>% 
   kableExtra::kable_styling(font_size = 10,
                latex_options = "striped")

Summary statistics per habitat and management type

Habitat	Management	Mean Density	Median Density	Min Density	Max Density	Sd Density	N
Crop Field	Fallow	6.309677	0	0	164.9	29.698646	31
Crop Field	Ploughed	3.632967	0	0	56.7	10.657080	91
Garden	Ploughed	4.748750	0	0	267.4	30.155135	80
Garden	Ploughed/Mowed	0.000000	0	0	0.0	NA	1
Orchard	Mowed	5.385714	0	0	117.1	15.181616	182
Orchard	No Management	1.666667	0	0	33.8	6.900452	24
Orchard	Ploughed	13.652381	0	0	378.1	49.918918	63
Orchard	Ploughed/Mowed	0.000000	0	0	0.0	0.000000	2
Short-cut Lawn	Mowed	4.771831	0	0	77.7	15.193816	71
Shrubland	No Management	0.484127	0	0	22.2	2.970185	63
Steppe Grassland	Mowed	0.000000	0	0	0.0	NA	1
Steppe Grassland	No Management	4.886957	0	0	56.2	16.191448	23
Steppe Grassland	Grazed	0.000000	0	0	0.0	NA	1
Tall Lawn	No Management	1.676190	0	0	19.9	4.706156	21
Tree Stands	Mowed	0.000000	0	0	0.0	0.000000	3
Tree Stands	No Management	0.000000	0	0	0.0	0.000000	5
Vineyard	Mowed	2.302083	0	0	56.0	7.976894	144
Vineyard	No Management	0.350000	0	0	4.9	1.309580	14
Vineyard	Ploughed	7.791549	0	0	261.9	27.972632	142
Vineyard	Ploughed/Mowed	1.380000	0	0	15.9	3.460432	35
Pasture	Grazed	0.000000	0	0	0.0	0.000000	3
Ruderal	No Management	0.000000	0	0	0.0	0.000000	3
Tree Alley	Mowed	0.700000	0	0	4.6	1.614222	8

Presence only

table1 <- colonies_pres %>% 
  dplyr::group_by(district) %>% 
    dplyr::summarize(mean_density = mean(density),
            median_density = median(density),
            min_density = min(density),
            max_density = max(density),
                sd_density = sd(density),
            n = n())

# Modify the column names to capitalize the first letter and remove hyphens
colnames(table1) <- str_to_title(gsub("_", " ", colnames(table1)))
  

kable(table1, caption = "Summary statistics per district - Presence only") %>% 
   kableExtra::kable_styling(font_size = 10,
                latex_options = "striped")

Summary statistics per district - Presence only

District	Mean Density	Median Density	Min Density	Max Density	Sd Density	N
Hrušovany	42.95412	27.5	0.6	378.1	59.01815	85
Pavlovice	10.54037	6.0	0.2	77.7	13.64933	109

table1 <- colonies_pres %>% 
  dplyr::group_by(district, habitat) %>% 
    dplyr::summarize(mean_density = mean(density),
            median_density = median(density),
            min_density = min(density),
            max_density = max(density),
                 sd_density = sd(density),
            n = n())

# Modify the column names to capitalize the first letter and remove hyphens
colnames(table1) <- str_to_title(gsub("_", " ", colnames(table1)))

kable(table1, caption = "Summary statistics per district and habitat - Presence only") %>% 
   kableExtra::kable_styling(font_size = 10,
                latex_options = "striped")

Summary statistics per district and habitat - Presence only

District	Habitat	Mean Density	Median Density	Min Density	Max Density	Sd Density	N
Hrušovany	Crop Field	33.181818	19.30	0.6	164.9	47.115302	11
Hrušovany	Garden	58.416667	18.60	7.0	267.4	102.824033	6
Hrušovany	Orchard	46.518182	33.80	4.1	378.1	65.338290	33
Hrušovany	Short-cut Lawn	37.275000	33.85	17.9	63.5	19.473122	4
Hrušovany	Shrubland	15.250000	15.25	8.3	22.2	9.828784	2
Hrušovany	Steppe Grassland	56.200000	56.20	56.2	56.2	0.000000	2
Hrušovany	Vineyard	41.055556	23.60	2.1	261.9	53.586683	27
Pavlovice	Crop Field	16.120000	11.75	1.0	47.1	15.665589	10
Pavlovice	Garden	9.800000	8.90	7.7	12.8	2.666458	3
Pavlovice	Orchard	9.862857	6.80	0.4	39.5	9.350215	35
Pavlovice	Short-cut Lawn	31.616667	14.15	6.8	77.7	31.885258	6
Pavlovice	Tall Lawn	8.800000	7.10	1.1	19.9	8.004166	4
Pavlovice	Vineyard	7.808163	3.20	0.2	51.3	11.320517	49
Pavlovice	Tree Alley	2.800000	2.80	1.0	4.6	2.545584	2

table1 <- colonies_pres %>% 
  dplyr::group_by(district, management) %>% 
    dplyr::summarize(mean_density = mean(density),
            median_density = median(density),
            min_density = min(density),
            max_density = max(density),
                 sd_density = sd(density),
            n = n())

# Modify the column names to capitalize the first letter and remove hyphens
colnames(table1) <- str_to_title(gsub("_", " ", colnames(table1)))

kable(table1, caption = "Summary statistics per district and management type - Presence only") %>% 
   kableExtra::kable_styling(font_size = 10,
                latex_options = "striped")

Summary statistics per district and management type - Presence only

District	Management	Mean Density	Median Density	Min Density	Max Density	Sd Density	N
Hrušovany	Fallow	65.200000	19.30	11.4	164.9	86.433038	3
Hrušovany	Mowed	32.844118	29.60	2.1	117.1	25.630015	34
Hrušovany	No Management	30.266667	28.00	4.9	56.2	22.577127	6
Hrušovany	Ploughed	53.275000	27.60	0.6	378.1	78.680072	40
Hrušovany	Ploughed/Mowed	13.100000	13.10	10.3	15.9	3.959798	2
Pavlovice	Mowed	10.177359	5.60	0.2	77.7	14.420040	53
Pavlovice	No Management	6.900000	4.55	1.1	19.9	6.865858	6
Pavlovice	Ploughed	13.317073	8.20	0.9	51.3	14.160595	41
Pavlovice	Ploughed/Mowed	2.455556	0.80	0.2	7.2	2.752776	9

# Habitat and management
table1 <- colonies_pres %>% 
  dplyr::group_by(habitat, management) %>% 
    dplyr::summarize(mean_density = mean(density),
            median_density = median(density),
            min_density = min(density),
            max_density = max(density),
                 sd_density = sd(density),
            n = n())

# Modify the column names to capitalize the first letter and remove hyphens
colnames(table1) <- str_to_title(gsub("_", " ", colnames(table1)))

kable(table1, caption = "Summary statistics per habitat and management type - Presence only") %>% 
  kableExtra::kable_styling(font_size = 10,
                latex_options = "striped")

Summary statistics per habitat and management type - Presence only

Habitat	Management	Mean Density	Median Density	Min Density	Max Density	Sd Density	N
Crop Field	Fallow	65.200000	19.30	11.4	164.9	86.433038	3
Crop Field	Ploughed	18.366667	11.75	0.6	56.7	17.741013	18
Garden	Ploughed	42.211111	12.80	7.0	267.4	84.856709	9
Orchard	Mowed	20.855319	12.00	0.4	117.1	24.026180	47
Orchard	No Management	13.333333	3.40	2.8	33.8	17.727192	3
Orchard	Ploughed	47.783333	21.15	2.6	378.1	85.797870	18
Short-cut Lawn	Mowed	33.880000	23.50	6.8	77.7	26.452885	10
Shrubland	No Management	15.250000	15.25	8.3	22.2	9.828784	2
Steppe Grassland	No Management	56.200000	56.20	56.2	56.2	0.000000	2
Tall Lawn	No Management	8.800000	7.10	1.1	19.9	8.004166	4
Vineyard	Mowed	11.839286	6.30	0.2	56.0	14.829434	28
Vineyard	No Management	4.900000	4.90	4.9	4.9	NA	1
Vineyard	Ploughed	30.733333	12.30	0.9	261.9	49.264656	36
Vineyard	Ploughed/Mowed	4.390909	2.50	0.2	15.9	5.115751	11
Tree Alley	Mowed	2.800000	2.80	1.0	4.6	2.545584	2

Size of plots

table1 <- colonies %>% 
  dplyr::group_by(district) %>% 
  dplyr::summarize(mean_area = mean(size_ha),
            median_area = median(size_ha),
            min_area = min(size_ha),
            max_area = max(size_ha),
              sd_area = sd(size_ha),
            n = n())

# Modify the column names to capitalize the first letter and remove hyphens
colnames(table1) <- str_to_title(gsub("_", " ", colnames(table1)))
  

kable(table1, caption = "Summary statistics per district") %>% 
   kableExtra::kable_styling(font_size = 10,
                latex_options = "striped")

Summary statistics per district

District	Mean Area	Median Area	Min Area	Max Area	Sd Area	N
Hrušovany	0.2001361	0.1352804	0.0146362	1.667536	0.2262519	171
Pavlovice	0.5067885	0.2214649	0.0174018	25.410000	1.2011445	840

table1 <- colonies %>% 
  dplyr::group_by(district, habitat) %>% 
  dplyr::summarize(mean_area = mean(size_ha),
            median_area = median(size_ha),
            min_area = min(size_ha),
            max_area = max(size_ha),
              sd_area = sd(size_ha),
            n = n())

# Modify the column names to capitalize the first letter and remove hyphens
colnames(table1) <- str_to_title(gsub("_", " ", colnames(table1)))

kable(table1, caption = "Summary statistics per district and habitat") %>% 
   kableExtra::kable_styling(font_size = 10,
                latex_options = "striped")

Summary statistics per district and habitat

District	Habitat	Mean Area	Median Area	Min Area	Max Area	Sd Area	N
Hrušovany	Crop Field	0.3966076	0.2568766	0.0303225	1.6675363	0.4236616	20
Hrušovany	Garden	0.1342761	0.0935023	0.0246976	0.5349839	0.1392483	13
Hrušovany	Orchard	0.1918580	0.1430234	0.0196987	0.9856231	0.1809191	58
Hrušovany	Short-cut Lawn	0.1435141	0.1334912	0.0202375	0.3857541	0.1141826	8
Hrušovany	Shrubland	0.1517028	0.1209329	0.0361502	0.5738480	0.1640914	9
Hrušovany	Steppe Grassland	0.1150404	0.1429195	0.0477385	0.1487839	0.0425300	10
Hrušovany	Tall Lawn	0.1154959	0.0935520	0.0551552	0.1977806	0.0738014	3
Hrušovany	Tree Stands	0.0206813	0.0206813	0.0206813	0.0206813	NA	1
Hrušovany	Vineyard	0.1915666	0.1309018	0.0146362	1.0428849	0.1963686	49
Pavlovice	Crop Field	0.5412581	0.1865635	0.0197371	11.3900000	1.2850366	102
Pavlovice	Garden	0.1482291	0.0903777	0.0174018	1.3176442	0.1926856	68
Pavlovice	Orchard	0.3539295	0.1951953	0.0289295	3.9477607	0.4824126	213
Pavlovice	Short-cut Lawn	0.2185341	0.1497017	0.0353690	1.3425311	0.2214651	63
Pavlovice	Shrubland	0.6007065	0.2386269	0.0196223	4.5183010	0.8993255	54
Pavlovice	Steppe Grassland	0.7657877	0.3211392	0.0441818	4.2158971	1.0531802	15
Pavlovice	Tall Lawn	0.2794462	0.1555036	0.0536121	1.3159649	0.3380599	18
Pavlovice	Tree Stands	0.4762270	0.2923001	0.1131347	1.1440361	0.4469978	7
Pavlovice	Vineyard	0.7396431	0.3251528	0.0277983	25.4100000	1.7639577	286
Pavlovice	Pasture	0.7046851	0.3873776	0.3479647	1.3787130	0.5840579	3
Pavlovice	Ruderal	0.1945866	0.1747893	0.1040609	0.3049097	0.1018774	3
Pavlovice	Tree Alley	0.5919291	0.3618700	0.1540836	1.9100263	0.5753529	8

table1 <- colonies %>% 
  dplyr::group_by(district, management) %>% 
  dplyr::summarize(mean_area = mean(size_ha),
            median_area = median(size_ha),
            min_area = min(size_ha),
            max_area = max(size_ha),
              sd_area = sd(size_ha),
            n = n())

# Modify the column names to capitalize the first letter and remove hyphens
colnames(table1) <- str_to_title(gsub("_", " ", colnames(table1)))

kable(table1, caption = "Summary statistics per district and management type") %>% 
 kableExtra::kable_styling(font_size = 10,
                latex_options = "striped")

Summary statistics per district and management type

District	Management	Mean Area	Median Area	Min Area	Max Area	Sd Area	N
Hrušovany	Fallow	0.1012883	0.0838573	0.0303225	0.2071161	0.0749903	4
Hrušovany	Mowed	0.1837003	0.1417454	0.0196987	0.9856231	0.1705828	65
Hrušovany	No Management	0.1547822	0.1391000	0.0206813	0.5738480	0.1292096	30
Hrušovany	Ploughed	0.2354955	0.1351917	0.0146362	1.6675363	0.2951959	69
Hrušovany	Ploughed/Mowed	0.3283180	0.3770540	0.0247399	0.5831603	0.2823822	3
Pavlovice	Fallow	0.2059319	0.1197429	0.0244871	1.5442739	0.2902970	27
Pavlovice	Mowed	0.5168195	0.2264078	0.0289295	25.4100000	1.5139634	344
Pavlovice	No Management	0.4819026	0.2313789	0.0196223	4.5183010	0.6720863	123
Pavlovice	Ploughed	0.4106173	0.1832134	0.0174018	11.3900000	0.8485079	307
Pavlovice	Ploughed/Mowed	1.4483644	0.7468093	0.0852919	7.4042496	1.6776319	35
Pavlovice	Grazed	1.5824881	0.8830453	0.3479647	4.2158971	1.8192219	4

# Habitat and management
table1 <- colonies %>% 
  dplyr::group_by(habitat, management) %>% 
  dplyr::summarize(mean_area = mean(size_ha),
            median_area = median(size_ha),
            min_area = min(size_ha),
            max_area = max(size_ha),
              sd_area = sd(size_ha),
            n = n())

# Modify the column names to capitalize the first letter and remove hyphens
colnames(table1) <- str_to_title(gsub("_", " ", colnames(table1)))

kable(table1, caption = "Summary statistics per habitat and management type") %>% 
  kableExtra::kable_styling(font_size = 10,
                latex_options = "striped")

Summary statistics per habitat and management type

Habitat	Management	Mean Area	Median Area	Min Area	Max Area	Sd Area	N
Crop Field	Fallow	0.1924295	0.1040679	0.0244871	1.5442739	0.2736239	31
Crop Field	Ploughed	0.6282985	0.2479741	0.0197371	11.3900000	1.3495035	91
Garden	Ploughed	0.1475053	0.0931697	0.0174018	1.3176442	0.1851267	80
Garden	Ploughed/Mowed	0.0247399	0.0247399	0.0247399	0.0247399	NA	1
Orchard	Mowed	0.3499308	0.1760018	0.0196987	3.9477607	0.5006071	182
Orchard	No Management	0.4442826	0.3518460	0.0906520	2.0177169	0.4076211	24
Orchard	Ploughed	0.1879936	0.1234033	0.0281398	0.7861068	0.1568365	63
Orchard	Ploughed/Mowed	0.1604780	0.1604780	0.0852919	0.2356642	0.1063293	2
Short-cut Lawn	Mowed	0.2100811	0.1440802	0.0202375	1.3425311	0.2128754	71
Shrubland	No Management	0.5365631	0.2000103	0.0196223	4.5183010	0.8484931	63
Steppe Grassland	Mowed	0.2159259	0.2159259	0.2159259	0.2159259	NA	1
Steppe Grassland	No Management	0.3567564	0.1487839	0.0441818	1.5364734	0.4115373	23
Steppe Grassland	Grazed	4.2158971	4.2158971	4.2158971	4.2158971	NA	1
Tall Lawn	No Management	0.2560248	0.1534268	0.0536121	1.3159649	0.3180290	21
Tree Stands	Mowed	0.1895977	0.1633583	0.1131347	0.2923001	0.0924199	3
Tree Stands	No Management	0.5570954	0.3723737	0.0206813	1.1440361	0.5268096	5
Vineyard	Mowed	0.7333548	0.2944272	0.0359544	25.4100000	2.2523277	144
Vineyard	No Management	0.1785764	0.1375409	0.0309899	0.4949014	0.1376545	14
Vineyard	Ploughed	0.4330249	0.2448186	0.0146362	4.5212598	0.5783951	142
Vineyard	Ploughed/Mowed	1.4666289	0.7468093	0.1811133	7.4042496	1.6649259	35
Pasture	Grazed	0.7046851	0.3873776	0.3479647	1.3787130	0.5840579	3
Ruderal	No Management	0.1945866	0.1747893	0.1040609	0.3049097	0.1018774	3
Tree Alley	Mowed	0.5919291	0.3618700	0.1540836	1.9100263	0.5753529	8

Presence only

table1 <- colonies_pres %>% 
  dplyr::group_by(district) %>% 
  dplyr::summarize(mean_area = mean(size_ha),
            median_area = median(size_ha),
            min_area = min(size_ha),
            max_area = max(size_ha),
              sd_area = sd(size_ha),
            n = n())

# Modify the column names to capitalize the first letter and remove hyphens
colnames(table1) <- str_to_title(gsub("_", " ", colnames(table1)))
  

kable(table1, caption = "Summary statistics per district - Presence only") %>% 
  kableExtra::kable_styling(font_size = 10,
                latex_options = "striped")

Summary statistics per district - Presence only

District	Mean Area	Median Area	Min Area	Max Area	Sd Area	N
Hrušovany	0.2303903	0.1430515	0.0146362	1.667536	0.2709657	85
Pavlovice	1.1054466	0.4560539	0.0424719	25.410000	2.6039602	109

table1 <- colonies_pres %>% 
  dplyr::group_by(district, habitat) %>% 
  dplyr::summarize(mean_area = mean(size_ha),
            median_area = median(size_ha),
            min_area = min(size_ha),
            max_area = max(size_ha),
              sd_area = sd(size_ha),
            n = n())

# Modify the column names to capitalize the first letter and remove hyphens
colnames(table1) <- str_to_title(gsub("_", " ", colnames(table1)))

kable(table1, caption = "Summary statistics per district and habitat - Presence only") %>% 
  kableExtra::kable_styling(font_size = 10,
                latex_options = "striped")

Summary statistics per district and habitat - Presence only

District	Habitat	Mean Area	Median Area	Min Area	Max Area	Sd Area	N
Hrušovany	Crop Field	0.4052782	0.2071161	0.0303225	1.6675363	0.5214459	11
Hrušovany	Garden	0.2127082	0.1416330	0.0881456	0.5349839	0.1718984	6
Hrušovany	Orchard	0.2069387	0.1431608	0.0205713	0.9856231	0.2067632	33
Hrušovany	Short-cut Lawn	0.1182010	0.1055692	0.0558886	0.2057768	0.0694956	4
Hrušovany	Shrubland	0.1281067	0.1281067	0.1209329	0.1352804	0.0101452	2
Hrušovany	Steppe Grassland	0.0889138	0.0889138	0.0889138	0.0889138	0.0000000	2
Hrušovany	Vineyard	0.2264090	0.1496861	0.0146362	1.0428849	0.2361435	27
Pavlovice	Crop Field	0.3851903	0.2303097	0.0424719	1.0139817	0.3360204	10
Pavlovice	Garden	0.1069472	0.1120245	0.0783104	0.1305068	0.0264660	3
Pavlovice	Orchard	0.7705831	0.4560539	0.0475209	3.9477607	0.8521115	35
Pavlovice	Short-cut Lawn	0.2814993	0.1946690	0.0985939	0.7891762	0.2591121	6
Pavlovice	Tall Lawn	0.3283990	0.1463513	0.1002790	0.9206143	0.3961123	4
Pavlovice	Vineyard	1.7187884	0.7487816	0.0906884	25.4100000	3.7286426	49
Pavlovice	Tree Alley	1.0636540	1.0636540	0.2172817	1.9100263	1.1969512	2

table1 <- colonies_pres %>% 
  dplyr::group_by(district, management) %>% 
  dplyr::summarize(mean_area = mean(size_ha),
            median_area = median(size_ha),
            min_area = min(size_ha),
            max_area = max(size_ha),
              sd_area = sd(size_ha),
            n = n())

# Modify the column names to capitalize the first letter and remove hyphens
colnames(table1) <- str_to_title(gsub("_", " ", colnames(table1)))

kable(table1, caption = "Summary statistics per district and management type - Presence only") %>% 
 kableExtra::kable_styling(font_size = 10,
                latex_options = "striped")

Summary statistics per district and management type - Presence only

District	Management	Mean Area	Median Area	Min Area	Max Area	Sd Area	N
Hrušovany	Fallow	0.1083297	0.0875504	0.0303225	0.2071161	0.0902099	3
Hrušovany	Mowed	0.2101207	0.1604011	0.0205713	0.9856231	0.1994039	34
Hrušovany	No Management	0.1652889	0.1281067	0.0889138	0.3547129	0.1018198	6
Hrušovany	Ploughed	0.2540533	0.1377032	0.0146362	1.6675363	0.3409933	40
Hrušovany	Ploughed/Mowed	0.4801071	0.4801071	0.3770540	0.5831603	0.1457391	2
Pavlovice	Mowed	1.4487969	0.6213596	0.0636014	25.4100000	3.5779642	53
Pavlovice	No Management	0.3281995	0.2357093	0.1002790	0.9206143	0.3074804	6
Pavlovice	Ploughed	0.5404016	0.3692687	0.0424719	3.3848277	0.6055274	41
Pavlovice	Ploughed/Mowed	2.1757540	1.4495717	0.3736405	5.3673770	1.5082394	9

# Habitat and management
table1 <- colonies_pres %>% 
  dplyr::group_by(habitat, management) %>% 
  dplyr::summarize(mean_area = mean(size_ha),
            median_area = median(size_ha),
            min_area = min(size_ha),
            max_area = max(size_ha),
              sd_area = sd(size_ha),
            n = n())

# Modify the column names to capitalize the first letter and remove hyphens
colnames(table1) <- str_to_title(gsub("_", " ", colnames(table1)))

kable(table1, caption = "Summary statistics per habitat and management type - Presence only") %>% 
  kableExtra::kable_styling(font_size = 10,
                latex_options = "striped")

Summary statistics per habitat and management type - Presence only

Habitat	Management	Mean Area	Median Area	Min Area	Max Area	Sd Area	N
Crop Field	Fallow	0.1083297	0.0875504	0.0303225	0.2071161	0.0902099	3
Crop Field	Ploughed	0.4436097	0.2568766	0.0424719	1.6675363	0.4493147	18
Garden	Ploughed	0.1774546	0.1305068	0.0783104	0.5349839	0.1464228	9
Orchard	Mowed	0.6134802	0.2481750	0.0205713	3.9477607	0.7926764	47
Orchard	No Management	0.3367714	0.3547129	0.2961331	0.3594682	0.0352740	3
Orchard	Ploughed	0.2197501	0.1654932	0.0400045	0.5676248	0.1621767	18
Short-cut Lawn	Mowed	0.2161800	0.1532783	0.0558886	0.7891762	0.2145236	10
Shrubland	No Management	0.1281067	0.1281067	0.1209329	0.1352804	0.0101452	2
Steppe Grassland	No Management	0.0889138	0.0889138	0.0889138	0.0889138	0.0000000	2
Tall Lawn	No Management	0.3283990	0.1463513	0.1002790	0.9206143	0.3961123	4
Vineyard	Mowed	1.8145594	0.4412720	0.1293170	25.4100000	4.8549069	28
Vineyard	No Management	0.2029797	0.2029797	0.2029797	0.2029797	NA	1
Vineyard	Ploughed	0.5216953	0.2856495	0.0146362	3.3848277	0.6592842	36
Vineyard	Ploughed/Mowed	1.8674546	1.3805248	0.3736405	5.3673770	1.5140816	11
Tree Alley	Mowed	1.0636540	1.0636540	0.2172817	1.9100263	1.1969512	2

Boxplots

Habitat type

I would remove the means from the plots and leave the stats for Table 1 (perhaps Table 1 a, b)

# Compute the number of observations for each habitat type
n <- colonies %>% 
  group_by(habitat, district) %>% 
  dplyr::summarise(n = n())

# Find the maximum density value for each district
max_density <- colonies %>% group_by(district) %>% dplyr::summarise(max_density = max(density))

# Join the n and max_density data frames
n <- left_join(n, max_density, by = c("district" = "district"))

ggplot(colonies, aes(x=habitat, y=density)) + 
  geom_boxplot() +
  xlab("") + ylab("Density (burrow openings/ha)") +
  ggtitle("Density of burrow openings by habitat type")+
  facet_grid(. ~ district, scales = "free_x")+
  geom_text(data=n, aes(x=habitat, y=max_density-5, label=paste("n =", n)), vjust=-1.1, size=3.5, position=position_dodge(0.4))+
  coord_flip()+
  theme(legend.position = "none")

Same scale example (I prefer the version with different scales)

# Compute the number of observations for each habitat type
n <- colonies %>% 
  group_by(habitat, district) %>% 
  dplyr::summarise(n = n())

# Find the maximum density value for each district
max_density <- colonies %>% group_by(district) %>% dplyr::summarise(max_density = max(density))

# Join the n and max_density data frames
n <- left_join(n, max_density, by = c("district" = "district"))

ggplot(colonies, aes(x=habitat, y=density)) + 
  geom_boxplot() +
  xlab("") + ylab("Density (burrow openings/ha)") +
  ggtitle("Density of burrow openings by habitat type")+
  facet_grid(. ~ district)+
  geom_text(data=n, aes(x=habitat, y=max_density-5, label=paste("n =", n)), vjust=-1.1, size=3.5, position=position_dodge(0.4))+
  coord_flip()+
  theme(legend.position = "none")

Management

# Compute the number of observations for each habitat type
n <- colonies %>% 
  group_by(management, district) %>% 
  dplyr::summarise(n = n())

# Find the maximum density value for each district
max_density <- colonies %>% group_by(district) %>% dplyr::summarise(max_density = max(density))

# Join the n and max_density data frames
n <- left_join(n, max_density, by = c("district" = "district"))

ggplot(colonies, aes(x=management, y=density)) + 
  geom_boxplot() +
  xlab("") + ylab("Density (burrow openings/ha)") +
  ggtitle("Density of burrow openings by management type")+
  facet_grid(. ~ district, scales = "free_x")+
  geom_text(data=n, aes(x=management, y=max_density-10, label=paste("n =", n)), vjust=-1.1, size=3.5, position=position_dodge(0.4))+
  coord_flip()+
  theme(legend.position = "none")

ggplot(hru, aes(x=habitat, y=density)) + 
  geom_boxplot() +
  xlab("") + ylab("Density (burrow openings/ha)") +
  ggtitle("Density of burrow openings by Habitat in Hrusovany")+
  coord_flip()

ggplot(vp, aes(x=habitat, y=density)) + 
  geom_boxplot() +
  xlab("") + ylab("Density (burrow openings/ha)") +
  ggtitle("Density of burrow openings by Habitat in Velke Pavlovice")+
  coord_flip()

ggplot(hru, aes(x=management, y=density)) + 
  geom_boxplot() +
  xlab("") + ylab("Density (burrow openings/ha)") +
  ggtitle("Density of burrow openings by management type in Hrusovany")+
  coord_flip()

ggplot(vp, aes(x=management, y=density)) + 
  geom_boxplot() +
  xlab("") + ylab("Density (burrow openings/ha)") +
  ggtitle("Density of burrow openings by management type in Velke Pavlovice")+
  coord_flip()

Presence only

For BOs > 0 only. It can be misleading. Needs to be clearly stated.

ggplot(colonies_pres, aes(x=habitat, y=density)) + 
  geom_boxplot() +
  xlab("") + ylab("Density (burrow openings/ha)") +
  ggtitle("Density of burrow openings by Habitat - Presence only")+
  facet_grid(. ~ district, scales = "free_x")+
  coord_flip()

ggplot(colonies_pres, aes(x=management, y=density)) + 
  geom_boxplot() +
  xlab("") + ylab("Density (burrow openings/ha)") +
  ggtitle("Density of burrow openings by Management - Presence only")+
    facet_grid(. ~ district, scales = "free_x")+
  coord_flip()

Both with n

Habitat

# Load the dplyr package


# Compute the number of observations for each habitat and management type
nc_habitat <- colonies_pres %>% 
  group_by(habitat, district) %>% 
  dplyr::summarise(n = n())

# Find the maximum density value for each district
max_density <- colonies_pres %>% group_by(district) %>% dplyr::summarise(max_density = max(density))

# Join the n and max_density data frames
nc_habitat <- left_join(nc_habitat, max_density, by = c("district" = "district"))

# Generate the plot
p <- ggplot(colonies_pres, aes(x=habitat, y=density)) + 
  geom_boxplot() +
  geom_text(data=nc_habitat, aes(x=habitat, y=max_density-5, label=paste("n =", n)), vjust=-1.1, size=3.5, position=position_dodge(0.4)) +  # add n to the facets
  xlab("") + ylab("Density (burrow openings/ha)") +
  ggtitle("Density of burrow openings by habitat type - Presence only")+
  facet_grid(. ~ district, scales = "free_x")+
  coord_flip()+
  theme(legend.position = "none")

# Show the plot
p

Same scale example (I still prefer the version with different scales)

# Compute the number of observations for each habitat and management type
nc_habitat <- colonies_pres %>% 
  group_by(habitat, district) %>% 
  dplyr::summarise(n = n())

# Find the maximum density value for each district
max_density <- colonies_pres %>% group_by(district) %>% dplyr::summarise(max_density = max(density))

# Join the n and max_density data frames
nc_habitat <- left_join(nc_habitat, max_density, by = c("district" = "district"))

# Generate the plot
p <- ggplot(colonies_pres, aes(x=habitat, y=density)) + 
  geom_boxplot() +
  geom_text(data=nc_habitat, aes(x=habitat, y=max_density-5, label=paste("n =", n)), vjust=-1.1, size=3.5, position=position_dodge(0.4)) +  # add n to the facets
  xlab("") + ylab("Density (burrow openings/ha)") +
  ggtitle("Density of burrow openings by habitat type - Presence only")+
  facet_grid(. ~ district)+
  coord_flip()+
  theme(legend.position = "none")

# Show the plot
p

Management

# Compute the number of observations for each management and district type
nc_management <- colonies_pres %>% 
  group_by(management, district) %>% 
  dplyr::summarise(n = n())

# Find the maximum density value for each district
max_density <- colonies_pres %>% group_by(district) %>% dplyr::summarise(max_density = max(density))

# Join the n and max_density data frames
nc_management <- left_join(nc_management, max_density, by = c("district" = "district"))

# Generate the plot
p <- ggplot(colonies_pres, aes(x=management, y=density)) + 
  geom_boxplot() +
  geom_text(data=nc_management, aes(x=management, y=max_density-5, label=paste("n =", n)), vjust=-1.1, size=3.5, position=position_dodge(0.4)) +  # add n to the facets
  xlab("") + ylab("Density (burrow openings/ha)") +
  ggtitle("Density of burrow openings by management type - Presence only")+
  facet_grid(. ~ district, scales = "free_x")+
  coord_flip()+
  theme(legend.position = "none")

# Show the plot
p

ggplot(hru_presence, aes(x=habitat, y=density)) + 
  geom_boxplot() +
  xlab("") + ylab("Density (burrow openings/ha)") +
  ggtitle("Density of burrow openings by Habitat in Hrusovany - Presence only")+
  coord_flip()

ggplot(vp_presence, aes(x=habitat, y=density)) + 
  geom_boxplot() +
  xlab("") + ylab("Density (burrow openings/ha)") +
  ggtitle("Density of burrow openings by Habitat in Velke Pavlovice - Presence only")+
  coord_flip()

ggplot(hru_presence, aes(x=management, y=density)) + 
  geom_boxplot() +
  xlab("") + ylab("Density (burrow openings/ha)") +
  ggtitle("Density of burrow openings by management type in Hrusovany - Presence only")+
  coord_flip()

ggplot(vp_presence, aes(x=management, y=density)) + 
  geom_boxplot() +
  xlab("") + ylab("Density (burrow openings/ha)") +
  ggtitle("Density of burrow openings by management type in Velke Pavlovice - Presence only")+
  coord_flip()

Example color plot

Not so sure about the mean, with so many outliers. I’d leave just the n. The boxplot gives already the median. Color is also optional.

Once you all decide the details, I’ll change it for all the plots.

# Create a color palette for the management types
palette <- c("#0072B2", "#D55E00", "#CC79A7", "#F0E442")
names(palette) <- c("Ploughed", "Ploughed/Mowed", "Mowed", "No Management")

# Compute the number of observations for each management type
n <- vp_presence %>% 
  group_by(management) %>% 
  dplyr::summarise(n = n())


# Modify the plot code to include the color aesthetic and appropriate labels
ggplot(vp_presence, aes(x=management, y=density, fill=management)) + 
  geom_boxplot() +
  xlab("") + ylab("Density (burrow openings/ha)") +
  ggtitle("Density of burrow openings by Management Type in Velke Pavlovice - Presence Only") +
  scale_fill_manual(values=palette) +
  coord_flip()+
    geom_text(data=n, aes(x=management, y=70, label=paste("n =", n)), vjust=-1, size=3.5, position=position_dodge(0.4))+
   labs(fill = "Management")+
  theme(legend.position = "none")

Anovas

And tukey post-hoc tests to identify which specific levels of the variables are significantly different from one another.

# Perform a two-way ANOVA to test the effect of management and habitat on density
aov_results <- aov(density ~ management + habitat + district, data=colonies)

summary(aov_results)

##              Df Sum Sq Mean Sq F value Pr(>F)    
## management    5   4596     919   2.315 0.0419 *  
## habitat      11   5224     475   1.196 0.2851    
## district      1  54646   54646 137.590 <2e-16 ***
## Residuals   993 394388     397                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# Perform post-hoc tests to identify significant differences between levels
tukey_results <- TukeyHSD(aov_results)

tukey_results

##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = density ~ management + habitat + district, data = colonies)
## 
## $management
##                                    diff         lwr       upr     p adj
## Mowed-Fallow                 -2.2605332 -12.8609007  8.339834 0.9904100
## No Management-Fallow         -4.8521611 -16.0599369  6.355615 0.8190226
## Ploughed-Fallow               0.8100034  -9.8230874 11.443094 0.9999334
## Ploughed/Mowed-Fallow        -5.0386248 -18.8103703  8.733121 0.9026602
## Grazed-Fallow                -6.3096774 -36.5412171 23.921862 0.9913146
## No Management-Mowed          -2.5916279  -7.9842238  2.800968 0.7437580
## Ploughed-Mowed                3.0705366  -0.9949856  7.136059 0.2593859
## Ploughed/Mowed-Mowed         -2.7780916 -12.4283157  6.872132 0.9635060
## Grazed-Mowed                 -4.0491443 -32.6395595 24.541271 0.9986120
## Ploughed-No Management        5.6621645   0.2055247 11.118804 0.0367406
## Ploughed/Mowed-No Management -0.1864637 -10.5002056 10.127278 0.9999999
## Grazed-No Management         -1.4575163 -30.2786595 27.363627 0.9999913
## Ploughed/Mowed-Ploughed      -5.8486282 -15.5347860  3.837530 0.5160721
## Grazed-Ploughed              -7.1196809 -35.7222449 21.482883 0.9806772
## Grazed-Ploughed/Mowed        -1.2710526 -31.1826729 28.640568 0.9999963
## 
## $habitat
##                                       diff        lwr       upr     p adj
## Garden-Crop Field                0.2433934  -9.113690  9.600477 1.0000000
## Orchard-Crop Field               5.0261823  -2.091632 12.143996 0.4677268
## Short-cut Lawn-Crop Field        3.3234323  -6.421627 13.068491 0.9939269
## Shrubland-Crop Field             1.6273562  -8.501271 11.755983 0.9999958
## Steppe Grassland-Crop Field      5.5938647  -8.738692 19.926421 0.9816639
## Tall Lawn-Crop Field             2.8194197 -12.604451 18.243290 0.9999842
## Tree Stands-Crop Field           0.1713687 -23.655210 23.997947 1.0000000
## Vineyard-Crop Field              2.0996604  -4.803861  9.003182 0.9977735
## Pasture-Crop Field               2.6007456 -35.552315 40.753806 1.0000000
## Ruderal-Crop Field               1.1432292 -37.009831 39.296290 1.0000000
## Tree Alley-Crop Field           -0.7483987 -24.574977 23.078180 1.0000000
## Orchard-Garden                   4.7827888  -3.484424 13.050001 0.7622488
## Short-cut Lawn-Garden            3.0800388  -7.533618 13.693695 0.9985461
## Shrubland-Garden                 1.3839628  -9.582924 12.350850 0.9999997
## Steppe Grassland-Garden          5.3504713  -9.586252 20.287195 0.9908221
## Tall Lawn-Garden                 2.5760262 -13.410822 18.562874 0.9999957
## Tree Stands-Garden              -0.0720247 -24.266846 24.122796 1.0000000
## Vineyard-Garden                  1.8562670  -6.227181  9.939715 0.9998432
## Pasture-Garden                   2.3573521 -36.026753 40.741457 1.0000000
## Ruderal-Garden                   0.8998358 -37.484270 39.283941 1.0000000
## Tree Alley-Garden               -0.9917921 -25.186613 23.203029 1.0000000
## Short-cut Lawn-Orchard          -1.7027500 -10.406655  7.001155 0.9999684
## Shrubland-Orchard               -3.3988261 -12.530139  5.732487 0.9874332
## Steppe Grassland-Orchard         0.5676825 -13.078341 14.213706 1.0000000
## Tall Lawn-Orchard               -2.2067626 -16.994851 12.581326 0.9999981
## Tree Stands-Orchard             -4.8548135 -28.274840 18.565213 0.9999433
## Vineyard-Orchard                -2.9265219  -8.260413  2.407370 0.8194464
## Pasture-Orchard                 -2.4254367 -40.325935 35.475062 1.0000000
## Ruderal-Orchard                 -3.8829531 -41.783452 34.017546 1.0000000
## Tree Alley-Orchard              -5.7745810 -29.194607 17.645445 0.9996865
## Shrubland-Short-cut Lawn        -1.6960761 -12.995801  9.603649 0.9999979
## Steppe Grassland-Short-cut Lawn  2.2704325 -12.912349 17.453214 0.9999980
## Tall Lawn-Short-cut Lawn        -0.5040126 -16.720994 15.712968 1.0000000
## Tree Stands-Short-cut Lawn      -3.1520635 -27.499559 21.195432 0.9999996
## Vineyard-Short-cut Lawn         -1.2237719  -9.753327  7.305783 0.9999987
## Pasture-Short-cut Lawn          -0.7226867 -39.203210 37.757837 1.0000000
## Ruderal-Short-cut Lawn          -2.1802031 -40.660727 36.300321 1.0000000
## Tree Alley-Short-cut Lawn       -4.0718310 -28.419326 20.275665 0.9999937
## Steppe Grassland-Shrubland       3.9665086 -11.465270 19.398287 0.9995329
## Tall Lawn-Shrubland              1.1920635 -15.258267 17.642394 1.0000000
## Tree Stands-Shrubland           -1.4559875 -25.959527 23.047552 1.0000000
## Vineyard-Shrubland               0.4723042  -8.492974  9.437582 1.0000000
## Pasture-Shrubland                0.9733894 -37.606056 39.552835 1.0000000
## Ruderal-Shrubland               -0.4841270 -39.063573 38.095319 1.0000000
## Tree Alley-Shrubland            -2.3757549 -26.879294 22.127785 1.0000000
## Tall Lawn-Steppe Grassland      -2.7744451 -22.099208 16.550318 0.9999987
## Tree Stands-Steppe Grassland    -5.4224960 -31.941480 21.096488 0.9999505
## Vineyard-Steppe Grassland       -3.4942043 -17.029687 10.041279 0.9995132
## Pasture-Steppe Grassland        -2.9931192 -42.883052 36.896813 1.0000000
## Ruderal-Steppe Grassland        -4.4506355 -44.340568 35.439297 0.9999999
## Tree Alley-Steppe Grassland     -6.3422635 -32.861248 20.176721 0.9997667
## Tree Stands-Tall Lawn           -2.6480509 -29.772397 24.476295 1.0000000
## Vineyard-Tall Lawn              -0.7197593 -15.405906 13.966387 1.0000000
## Pasture-Tall Lawn               -0.2186741 -40.513591 40.076243 1.0000000
## Ruderal-Tall Lawn               -1.6761905 -41.971107 38.618726 1.0000000
## Tree Alley-Tall Lawn            -3.5678184 -30.692164 23.556528 0.9999995
## Vineyard-Tree Stands             1.9282917 -21.427499 25.284083 1.0000000
## Pasture-Tree Stands              2.4293768 -41.768931 46.627684 1.0000000
## Ruderal-Tree Stands              0.9718605 -43.226447 45.170168 1.0000000
## Tree Alley-Tree Stands          -0.9197674 -33.562381 31.722846 1.0000000
## Pasture-Vineyard                 0.5010851 -37.359754 38.361924 1.0000000
## Ruderal-Vineyard                -0.9564312 -38.817270 36.904408 1.0000000
## Tree Alley-Vineyard             -2.8480591 -26.203850 20.507732 0.9999998
## Ruderal-Pasture                 -1.4575163 -54.762681 51.847648 1.0000000
## Tree Alley-Pasture              -3.3491443 -47.547452 40.849163 1.0000000
## Tree Alley-Ruderal              -1.8916279 -46.089935 42.306679 1.0000000
## 
## $district
##                          diff       lwr       upr p adj
## Pavlovice-Hrušovany -19.37638 -22.65736 -16.09541     0

Presence only

And tukey post-hoc tests to identify which specific levels of the variables are significantly different from one another.

# Perform a two-way ANOVA to test the effect of management and habitat on density
aov_results <- aov(density ~ management + habitat + district, data=colonies_pres)

summary(aov_results)

##              Df Sum Sq Mean Sq F value   Pr(>F)    
## management    4  18345    4586   2.817   0.0267 *  
## habitat       8  16107    2013   1.237   0.2801    
## district      1  35371   35371  21.725 6.11e-06 ***
## Residuals   180 293059    1628                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# Perform post-hoc tests to identify significant differences between levels
tukey_results <- TukeyHSD(aov_results)

tukey_results

##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = density ~ management + habitat + district, data = colonies_pres)
## 
## $management
##                                     diff         lwr       upr     p adj
## Mowed-Fallow                 -46.1643678 -111.456790 19.128055 0.2958522
## No Management-Fallow         -46.6166667 -118.388850 25.155517 0.3827034
## Ploughed-Fallow              -32.1506173  -97.523598 33.222363 0.6570096
## Ploughed/Mowed-Fallow        -60.8090909 -133.230810 11.612628 0.1452719
## No Management-Mowed           -0.4522989  -34.691934 33.787336 0.9999996
## Ploughed-Mowed                14.0137505   -3.154035 31.181536 0.1665963
## Ploughed/Mowed-Mowed         -14.6447231  -50.225780 20.936334 0.7882062
## Ploughed-No Management        14.4660494  -19.926956 48.859054 0.7745745
## Ploughed/Mowed-No Management -14.1924242  -60.605327 32.220479 0.9168932
## Ploughed/Mowed-Ploughed      -28.6584736  -64.387143  7.070196 0.1805889
## 
## $habitat
##                                        diff        lwr       upr     p adj
## Garden-Crop Field                21.7469136  -28.72984  72.22366 0.9135625
## Orchard-Crop Field               17.5114558  -14.11806  49.14097 0.7219771
## Short-cut Lawn-Crop Field        27.4295530  -21.24846  76.10757 0.7022079
## Shrubland-Crop Field              9.2518519  -84.50455 103.00825 0.9999974
## Steppe Grassland-Crop Field      50.2018519  -43.55455 143.95825 0.7571806
## Tall Lawn-Crop Field              2.8018519  -66.31628  71.91999 1.0000000
## Vineyard-Crop Field               8.6567798  -22.57747  39.89103 0.9942130
## Tree Alley-Crop Field            -3.6504470  -97.40685  90.10595 1.0000000
## Orchard-Garden                   -4.2354578  -49.17528  40.70436 0.9999982
## Short-cut Lawn-Garden             5.6826394  -52.52999  63.89527 0.9999976
## Shrubland-Garden                -12.4950617 -111.53760  86.54747 0.9999825
## Steppe Grassland-Garden          28.4549383  -70.58760 127.49747 0.9926258
## Tall Lawn-Garden                -18.9450617  -95.07966  57.18954 0.9972494
## Vineyard-Garden                 -13.0901338  -57.75264  31.57237 0.9915753
## Tree Alley-Garden               -25.3973606 -124.43990  73.64517 0.9966074
## Short-cut Lawn-Orchard            9.9180973  -32.99151  52.82770 0.9983680
## Shrubland-Orchard                -8.2596039  -99.15486  82.63565 0.9999987
## Steppe Grassland-Orchard         32.6903961  -58.20486 123.58565 0.9690365
## Tall Lawn-Orchard               -14.7096039  -79.89397  50.47476 0.9986255
## Vineyard-Orchard                 -8.8546760  -30.00328  12.29393 0.9258570
## Tree Alley-Orchard              -21.1619027 -112.05716  69.73335 0.9982816
## Shrubland-Short-cut Lawn        -18.1777011 -116.31572  79.96032 0.9996759
## Steppe Grassland-Short-cut Lawn  22.7722989  -75.36572 120.91032 0.9983219
## Tall Lawn-Short-cut Lawn        -24.6277011  -99.58185  50.32645 0.9822999
## Vineyard-Short-cut Lawn         -18.7727732  -61.39185  23.84631 0.9029041
## Tree Alley-Short-cut Lawn       -31.0800000 -129.21802  67.05802 0.9860084
## Steppe Grassland-Shrubland       40.9500000  -85.74564 167.64564 0.9840556
## Tall Lawn-Shrubland              -6.4500000 -116.17164 103.27164 1.0000000
## Vineyard-Shrubland               -0.5950721  -91.35354  90.16340 1.0000000
## Tree Alley-Shrubland            -12.9022989 -139.59794 113.79334 0.9999967
## Tall Lawn-Steppe Grassland      -47.4000000 -157.12164  62.32164 0.9123100
## Vineyard-Steppe Grassland       -41.5450721 -132.30354  49.21340 0.8818319
## Tree Alley-Steppe Grassland     -53.8522989 -180.54794  72.84334 0.9195673
## Vineyard-Tall Lawn                5.8549279  -59.13856  70.84842 0.9999987
## Tree Alley-Tall Lawn             -6.4522989 -116.17394 103.26934 1.0000000
## Tree Alley-Vineyard             -12.3072268 -103.06570  78.45124 0.9999695
## 
## $district
##                          diff       lwr       upr   p adj
## Pavlovice-Hrušovany -25.65017 -37.17136 -14.12898 1.9e-05