Spatio-temporal dynamics of attacks around deaths of wolves: A statistical assessment of lethal control efficiency in France

This repository contains the supplementary materials (Supplementary_figures.docx, Supplementary_tables.docx) of the manuscript: "Spatio-temporal dynamics of attacks around deaths of wolves: A statistical assessment of lethal control efficiency in France". This repository also provides the R codes and datasets necessary to run the analyses described in the manuscript. 

The R datasets with suffix "_a" have anonymous spatial coordinates to respect confidentiality. Therefore, the preliminary preparation of the data is not provided in the public codes. These datasets, all geolocated and necessary to the analyses, are:

• Attack_sf_a.RData: 19,302 analyzed wolf attacks on sheep
  • ID: unique ID of the attack
  • DATE: date of the attack
  • PASTURE: the related pasture ID from "Pasture_sf_a" where the attack is located
  • STATUS: column resulting from the preparation and the attribution of attacks to pastures (part 2.2.4 of the manuscript); not shown here to respect confidentiality
• Pasture_sf_a.RData: 4987 analyzed pastures grazed by sheep
  • ID: unique ID of the pasture
  • CODE: Official code in the pastoral census
  • FLOCK_SIZE: maximum annual number of sheep grazing in the pasture
  • USED_MONTHS: months for which the pasture is grazed by sheep
• Removal_sf_a.RData:  232 analyzed single wolf removal or groups of wolf removals
  • ID: unique ID of the removal
  • OVERLAP: are they single removal ("non-interacting" in the manuscript => "NO" here), or not ("interacting" in the manuscrit, here "SIMULTANEOUS" for removals occurring during the same operation or "NON-SIMULTANEOUS" if not). 
  • DATE_MIN: date of the single removal or date of the first removal of a group
  • DATE_MAX: date of the single removal or date of the last removal of a group
  • CLASS: administrative type of the removal according to definitions from 2.1 part of the manuscript
  • SEX: sex or sexes of the removed wolves if known
  • AGE: class age of the removed wolves if known
  • BREEDER: breeding status of the removed female wolves, "Yes" for female breeder, "No" for female non-breeder. Males are "No" by default, when necropsied; dead individuals with NA were not found. 
  • SEASON: season of the removal, as defined in part 2.3.4 of the manuscript
  • MASSIF: mountain range attributed to the removal, as defined in part 2.3.4 of the manuscript
• Area_to_exclude_sf_a.RData: one row for each mountain range, corresponding to the area where removal controls of the mountain range could not be sampled, as defined in part 2.3.6 of the manuscript

These datasets were used to run the following analyses codes:

• Code 1 : The file Kernel_wolf_culling_attacks_p.R contains the before-after analyses.
  • We start by delimiting the spatio-temporal buffer for each row of the "Removal_sf_a.RData" dataset. 
    • We identify the attacks from "Attack_sf_a.RData" within each buffer, giving the data frame "Buffer_df" (one row per attack)
    • We select the pastures from "Pasture_sf_a.RData" within each buffer, giving the data frame "Buffer_sf" (one row per removal)
  • We calculate the spatial correction
    • We spatially slice each buffer into 200 rings, giving the data frame "Ring_sf" (one row per ring)
    • We add the total pastoral area of the ring of the attack ("SPATIAL_WEIGHT"), for each attack of each buffer, within Buffer_df ("Buffer_df.RData")
  • We calculate the pastoral correction
    
    • We create the pastoral matrix for each removal, giving a matrix of 200 rows (one for each ring) and 180 columns (one for each day, 90 days before the removal date and 90 day after the removal date), with the total pastoral area in use by sheep for each corresponding cell of the matrix (one element per removal, "Pastoral_matrix_lt.RData")
    • We simulate, for each removal, the random distribution of the attacks from "Buffer_df.RData" according to "Pastoral_matrix_lt.RData". The process is done 100 times (one element per simulation, "Buffer_simulation_lt.RData"). 
  • We estimate the attack intensities
    • We classified the removals into 20 subsets, according to part 2.3.4 of the manuscript ("Variables_lt.RData") (one element per subset)
    • We perform, for each subset, the kernel estimations with the observed attacks ("Kernel_lt.RData"), with the simulated attacks ("Kernel_simulation_lt.RData") and we correct the first kernel computations with the second ("Kernel_controlled_lt.RData") (one element per subset). 
    • We calculate the trend of attack intensities, for each subset, that compares the total attack intensity before and after the removals (part 2.3.5 of the manuscript), giving "Trends_intensities_df.RData". (one row per subset)
    • We calculate the trend of attack intensities, for each subset, along the spatial axis, three times, one for each time analysis scale. This gives "Shift_df" (one row per ring and per time analysis scale.
• Code 2 : The file Control_removals_p.R contains the control-impact analyses.
  • It starts with the simulation of 100 removal control sets ("Control_sf_lt_a.RData") from the real set of removals ("Removal_sf_a.RData"), that is done with the function "Control_fn" (l. 92).
  • The rest of the analyses follows the same process as in the first code "Kernel_wolf_culling_attacks_p.R", in order to apply the before-after analyses to each control set. All objects have the same structure as before, except that they are now a list, with one resulting element per control set. These objects have "control" in their names (not to be confused with "controlled" which refers to the pastoral correction already applied in the first code). 
  • The code is also applied again, from l. 92 to l. 433, this time for the real set of removals (l. 121) - with "Simulated = FALSE" (l. 119). We could not simply use the results from the first code because the set of removals is restricted to removals attributed to mountain ranges only. There are 2 resulting objects: "Kernel_real_lt.RData" (observed real trends) and "Kernel_controlled_real_lt.RData" (real trends corrected for pastoral use). 
  • The part of the code from line 439 to 524 relates to the calculations of the trends (for the real set and the control sets), as in the first code, giving "Trends_intensities_real_df.RData" and "Trends_intensities_control_lt.RData". 
  • The part of the code from line 530 to 588 relates to the calculation of the 95% confidence intervals and the means of the intensity trends for each subset based on the results of the 100 control sets (Trends_intensities_mean_control_df.RData, Trends_intensities_CImin_control_df.RData and Trends_intensities_CImax_control_df.RData). This will be used to test the significativity of the real trends. This comparison is done right after, l. 595-627, and gives the data frame "Trends_comparison_df.RData". 
• Code 3 : The file Figures.R produces part of the figures from the manuscript:
  • "Dataset map": figure 1
  • "Buffer": figure 2 (then pasted in powerpoint)
  • "Kernel construction": figure 5 (then pasted in powerpoint)
  • "Trend distributions": figure 7
  • "Kernels": part of figures 10 and S2
  • "Attack shifts": figure 9 and S1
  • "Significant": figure 8