-
Notifications
You must be signed in to change notification settings - Fork 0
/
investment_scenarios.R
11747 lines (9068 loc) · 356 KB
/
investment_scenarios.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
### This script is for the first part of my final PhD chapter. The analysis is investigating scenarios around investment of conservation funds over time, given increasing threats. Threats here represent increasing human population density, which is represented by increasing user budgets.
### Load libraries ####
library('tidyverse')
library('GMSE')
library('patchwork')
library('viridis')
library('scales')
############## LANDSCAPE DETAILS ##################
# Study period- 50 years
# Landscape dimensions - 100 x 100 cells. which results in 10,000 cells (or ha). With 30 villages, this results in 333.33 ha or 3.33km2 per village. The number of trees is 100,000
# There is no public land. It does not serve any purpose in this study.
# "Resources" are trees. The resources therefore do not move.
# The density of trees in tropical forest landscapes vary hugely. In previous simulations I have assumed 50 stems/ha which is low, but not implausible (e.g. deciduous dipterocarp woodland). A reference for this value can be found here:https://www.jstor.org/stable/44521915. This resulted in >1M trees, which meant that the differences in trees lost between scenarios, and the total number of trees lost, were too low. In the final scenarios, the number of trees was set to 100,000. Note that the trees are distributed randomly across the landscape, and so there will not be exactly 50/cell. This reflects reality.
# Trees in a cell reduce the farmer's yield. The amount a tree reduces yield is governed by an exponential function: yield = (1 - % yield reduction per tree) ^ remaining trees. I want a farmer's yield to be reduced by a significant amount if all trees in a cell are standing. But the trees do not completely eliminate yield. This is a balance between the farmer being able to farm and gain some yield even when there are trees on their cell, but also providing an incentive to cull where possible. I have set each tree on a cell to reduce yield by 8%. See the N1:N1f scenario comparisons which were used to decide on the parameter values for res_consume and tendcrop_yield
# The amount a user can increase their yield by tending crops is governed by tend_crop_yld. I have set this at 0.01 (1%) which means they can increase their yield on a cell by 1% in a single time step if they choose to tend crops. This is lower than the yield gain they would make if they felled some trees. This is set up so that there is an incentive to fell trees and expand their farmland, as it will increase their yield. See the N1:N1f scenario comparisons which were used to decide on the parameter values for res_consume and tendcrop_yield
# For simplicity, I am assuming there is no natural death and no natural birth (forest regeneration). remove_pr is set to 0, and lambda is set to 0, and res_death_type is set to 0 (new value created by Brad that means no natural death at all)
# The carrying capacity of new resources is set to 1 as it has to be a positive number but I want it as low as possible i.e. there is no real recruitment
# Resource carrying capacity is set very high (5,000,000) to reduce density-dependent death. Although res_death_type is set to 0 (no natural death) and so this parameter shouldn't be doing anything.
# The max age of trees is set high - 1000. This is to reduce natural death caused by old age (this should now be obsolete)
# The observation process is set to density-based sampling, with 1 observation per time step. The manager can move in any direction. Currently the manager can see 150 cells, and move 50 cells. We decided to remove any observation error.
# There is no minimum age for resources to be acted upon i.e. all trees in the landscape can be observed/culled
# Agents are permitted to move at the end of each time step. Because land_ownership==TRUE I believe this then only relates to the manager.
# User and manager budgets will vary based on the scenario. But the total amount of budget available to the manager for the whole study period will be the same. The total is 25,000. See the "scenario_details_budgets" spreadsheet
# Group_think == FALSE, and so users act independently
# Only culling is allowed (i.e. cutting down of trees)
# farming is allowed (tend_crops==TRUE, i.e. farmers can increase their yield by tending their crops rather than felling trees)
############## SCENARIOS #####################
# Below are the scenarios that I will run. There are a few "null" scenario, which will explore the basic manager/user dynamics and will potentially reflect the counterfactuals. The scenarios starting with "N.." will not feature in the main text but will be in the supporting information. The scenarios labelled 1 through 5 will be in the main text and will form the main results.
# N1 - Null - Manager and user budgets do not change throughout the study period, and are equal
# N2a - Decreasing Null - Manager budget remains constant, user budget decreases linearly
# N2b - Decreasing Null - User budget remains constant, manager budget decreases linearly
# N3 - Optimistic Null - Manager budget increases linearly over time, user budget remains constant
# 1 - Manager budget remains constant, user budgets increase linearly
# 2 - Manager and user budgets both increase linearly (different slopes and starting values)
# 3 - Sine wave - Manager budget increases and decreases in a predictable/regular way above and below a mean (like a sine wave), user budget increases linearly
# 4 - Random complex wave(s), but with small variation which results in a "core" budget i.e., the manager budget doesn't drop too low
# 5 - Random complex wave(s) - As above, but with high variation i.e., no "core" budget so the manager budget can drop very low, but also go very high
## below there are more scenarios than I have listed above. These were run to work out some of the parameter values (e.g. res_consume and tend_crop_yield)
#### N1 ####
## The original runs of N1 had the landscape set at 200x200 cells, 2 million trees, res_consume set to 0.05, and tend_crop_yield set to 0.02, and 40% public land. This was the set up used to run all the comparisons below between N1a:N1f. All of the comparison plots etc. saved in the N1 comparison folder will be using this set up.
# Then I changed the set up to a landscape of 150x150 cells which results in 22,500 cells (or ha). With 20 villages, this results in 1,125 ha or 11.25km2 per village. This means the number of trees will be 1,125,000. Res_consume will be 0.08 and tend_crop_yield will be 0.01. There will be no public land.
# The final simulations though, have 100x100 cells, 30 villages, and 100,000 trees (see landcape section above). The below code has therefore been changed to produce the final N1
# This null scenario has the manager and user budgets remaining static over the entire study period
N1 <- gmse(
time_max = 50,
land_dim_1 = 100,
land_dim_2 = 100, # landscape is 40,000ha or 400km2
res_movement = 0, # trees don't move
remove_pr = 0, # Assume no death
lambda = 0, # assume no growth
agent_view = 150,
agent_move = 50,
res_birth_K = 1, # must be positive value, but I want it small i.e. no real recruitment
res_death_K = 5000000, # carrying capacity set to way above starting number of resources
res_move_type = 0, # 0=no move,
res_death_type = 0, # no natural death
observe_type = 0, # 0=density-based sampling
times_observe = 1,
obs_move_type = 1, # uniform in any direction
res_min_age = 0, # age of resources before agents record/act on them
res_move_obs = FALSE, # trees don't move
plotting = FALSE,
res_consume = 0.08, # Trees have 8% impact on yield
# all genetic algorithm parameters left to default
move_agents = TRUE,
max_ages = 1000,
minimum_cost = 10,
user_budget = 500,
manager_budget = 500,
usr_budget_rng = 50, # introduce variation around the mean user budget (removes step pattern)
manage_target = 100000,
RESOURCE_ini = 100000,
culling = TRUE,
tend_crops = TRUE,
tend_crop_yld = 0.01, # tending crops increases yield by 1% - less than that of culling trees
stakeholders = 30,
land_ownership = TRUE,
public_land = 0,
manage_freq = 1,
group_think = FALSE
)
# save summary
N1_summary <- as.data.frame(gmse_table(N1))
write.csv(N1_summary, file="outputs/investment/null_scenarios/N1/N1_summary.csv")
# load data
#N1_summary <- read.csv("outputs/investment/null_scenarios/N1/N1_summary.csv")
# custom plots
time_cost_N1 <- ggplot(N1_summary, aes(x=time_step,y=cost_culling))+
geom_line(size=1, color="firebrick3")+
theme_classic()+
theme(axis.text = element_text(size=15),
axis.title = element_text(size=15))+
xlab("Time step")+
ylab("Cost of felling")
time_cull_N1 <- ggplot(N1_summary, aes(x=time_step, y=act_culling))+
geom_line(size=1, color="firebrick3")+
ylim(0,350)+
theme_classic()+
theme(axis.text = element_text(size=15),
axis.title = element_text(size=15))+
xlab("Time step")+
ylab("Number of felling actions")
time_res_N1 <- ggplot(N1_summary, aes(x=time_step, y=resources))+
geom_line(size=1, color="firebrick3")+
ylim(0,100000)+
theme_classic()+
theme(axis.text = element_text(size=15),
axis.title = element_text(size=15))+
xlab("Time step")+
ylab("Remaining trees")
time_yield_N1 <- ggplot(N1_summary, aes(x=time_step, y=crop_yield))+
geom_line(size=1, color="firebrick3")+
ylim(0,1500)+
theme_classic()+
theme(axis.text = element_text(size=15),
axis.title = element_text(size=15))+
xlab("Time step")+
ylab("Yield")
N1_custom_plots <- time_cost_N1 + time_cull_N1 + time_res_N1 + time_yield_N1
N1_custom_plots[[1]] <- N1_custom_plots[[1]] + xlab("")
N1_custom_plots[[2]] <- N1_custom_plots[[2]] + xlab("")
ggsave("outputs/investment/null_scenarios/N1/N1_custom_plots.png",
width = 25, height = 20, units = "cm", dpi=300)
# gmse plots
N1_gmse_plots <- plot_gmse_results(sim_results = N1)
# remove object
#rm(N1)
## N1a ####
# here I have made the yield increases from tending crops and felling trees equal.
N1a <- gmse(
time_max = 50,
land_dim_1 = 200,
land_dim_2 = 200, # landscape is 40,000ha or 400km2
res_movement = 0, # trees don't move
remove_pr = 0, # Assume no death
lambda = 0, # assume no growth
agent_view = 10,
agent_move = 50,
res_birth_K = 1, # must be positive value, but I want it small i.e. no real recruitment
res_death_K = 5000000, # carrying capacity set to way above starting number of resources
res_move_type = 0, # 0=no move,
res_death_type = 0, # no natural death
observe_type = 0, # 0=density-based sampling
times_observe = 1,
obs_move_type = 1, # uniform in any direction
res_min_age = 0, # age of resources before agents record/act on them
res_move_obs = FALSE, # trees don't move
plotting = FALSE,
res_consume = 0.05, # Trees have 5% impact on yield
# all genetic algorithm parameters left to default
move_agents = TRUE,
max_ages = 1000,
minimum_cost = 10,
user_budget = 200,
manager_budget = 200,
usr_budget_rng = 20, # introduce variation around the mean user budget (removes step pattern)
manage_target = 2000000,
RESOURCE_ini = 2000000,
culling = TRUE,
tend_crops = TRUE,
tend_crop_yld = 0.05, # tending crops increases yield by 5% - same as culling trees
stakeholders = 20,
land_ownership = TRUE,
public_land = 0.4,
manage_freq = 1,
group_think = FALSE
)
# As predicted, there is less felling of trees by the users becuase there is no real incentive to do so.
N1a_summary <- as.data.frame(gmse_table(N1a))
#write.csv(N1a_summary, file="outputs/investment/null_scenarios/N1/N1a_summary.csv")
# load
N1a_summary <- read.csv("outputs/investment/null_scenarios/N1/N1a_summary.csv")
# plots
time_cost_N1a <- ggplot(N1a_summary, aes(x=time_step,y=cost_culling))+
geom_line()+
theme_classic()
time_cull_N1a <- ggplot(N1a_summary, aes(x=time_step, y=act_culling))+
geom_line()+
theme_classic()
time_res_N1a <- ggplot(N1a_summary, aes(x=time_step, y=resources))+
geom_line()+
ylim(0,2000000)+
theme_classic()
time_yield_N1a <- ggplot(N1a_summary, aes(x=time_step, y=crop_yield))+
geom_line()+
theme_classic()
N1a_custom_plots <- time_cost_N1a + time_cull_N1a + time_res_N1a + time_yield_N1a
ggsave("outputs/investment/null_scenarios/N1/N1a_custom_plots.png", N1a_custom_plots,
width = 30, height = 20, units="cm", dpi=300)
plot_gmse_results(sim_results = N1a)
### calculate the % yield for each timestep
# Extract total yield for all cells from each timestep
yield_ts <- sapply(1:length(N1a$land), function(x) sum(N1a$land[[x]][,,2]))
# into dataframe
yield.df <- data.frame(time_step = 1:50,
total_yld = yield_ts,
available_yld = 24000)
# add % yield to df
yield.df$perc_yld <- yield.df$total_yld/yield.df$available_yld*100
# add sim number and description
yield.df$sim <- "N1a"
yield.df$desc <- "tend==cull"
## N1b ####
# Here I have increased the res_consume to 0.08 (8%)
N1b <- gmse(
time_max = 50,
land_dim_1 = 200,
land_dim_2 = 200, # landscape is 40,000ha or 400km2
res_movement = 0, # trees don't move
remove_pr = 0, # Assume no death
lambda = 0, # assume no growth
agent_view = 10,
agent_move = 50,
res_birth_K = 1, # must be positive value, but I want it small i.e. no real recruitment
res_death_K = 5000000, # carrying capacity set to way above starting number of resources
res_move_type = 0, # 0=no move,
res_death_type = 0, # no natural death
observe_type = 0, # 0=density-based sampling
times_observe = 1,
obs_move_type = 1, # uniform in any direction
res_min_age = 0, # age of resources before agents record/act on them
res_move_obs = FALSE, # trees don't move
plotting = FALSE,
res_consume = 0.08, # Trees have 8% impact on yield
# all genetic algorithm parameters left to default
move_agents = TRUE,
max_ages = 1000,
minimum_cost = 10,
user_budget = 200,
manager_budget = 200,
usr_budget_rng = 20, # introduce variation around the mean user budget (removes step pattern)
manage_target = 2000000,
RESOURCE_ini = 2000000,
culling = TRUE,
tend_crops = TRUE,
tend_crop_yld = 0.02, # tending crops increases yield by 2% - less than that of culling trees
stakeholders = 20,
land_ownership = TRUE,
public_land = 0.4,
manage_freq = 1,
group_think = FALSE
)
N1b_summary <- as.data.frame(gmse_table(N1b))
#write.csv(N1b_summary, file="outputs/investment/null_scenarios/N1/N1b_summary.csv")
rm(N1b)
# load
#N1b_summary <- read.csv("outputs/investment/null_scenarios/N1/N1b_summary.csv")
## N1c ####
# Following on from Nils' email, I am going to run another scenario with more extreme parameters for res_consume and tend_crop_yld. This is to add further contrast so we can better tease apart the scenarios.
# Here I have increased the res_consume to 0.1 (10%)
N1c <- gmse(
time_max = 50,
land_dim_1 = 200,
land_dim_2 = 200, # landscape is 40,000ha or 400km2
res_movement = 0, # trees don't move
remove_pr = 0, # Assume no death
lambda = 0, # assume no growth
agent_view = 10,
agent_move = 50,
res_birth_K = 1, # must be positive value, but I want it small i.e. no real recruitment
res_death_K = 5000000, # carrying capacity set to way above starting number of resources
res_move_type = 0, # 0=no move,
res_death_type = 0, # no natural death
observe_type = 0, # 0=density-based sampling
times_observe = 1,
obs_move_type = 1, # uniform in any direction
res_min_age = 0, # age of resources before agents record/act on them
res_move_obs = FALSE, # trees don't move
plotting = FALSE,
res_consume = 0.1, # Trees have 10% impact on yield
# all genetic algorithm parameters left to default
move_agents = TRUE,
max_ages = 1000,
minimum_cost = 10,
user_budget = 200,
manager_budget = 200,
usr_budget_rng = 20, # introduce variation around the mean user budget (removes step pattern)
manage_target = 2000000,
RESOURCE_ini = 2000000,
culling = TRUE,
tend_crops = TRUE,
tend_crop_yld = 0.02, # tending crops increases yield by 2% - less than that of culling trees
stakeholders = 20,
land_ownership = TRUE,
public_land = 0.4,
manage_freq = 1,
group_think = FALSE
)
N1c_summary <- as.data.frame(gmse_table(N1c))
write.csv(N1c_summary, file="outputs/investment/null_scenarios/N1/N1c_summary.csv")
rm(N1c)
## N1d ####
# Following on from Nils' email, I am going to run another scenario with more extreme parameters for res_consume and tend_crop_yld. This is to add further contrast so we can better tease apart the scenarios.
# Here I have kept res_consume at 0.08 (8%) but reduced tend_crop_yld to 0.01 (1%)
N1d <- gmse(
time_max = 50,
land_dim_1 = 200,
land_dim_2 = 200, # landscape is 40,000ha or 400km2
res_movement = 0, # trees don't move
remove_pr = 0, # Assume no death
lambda = 0, # assume no growth
agent_view = 10,
agent_move = 50,
res_birth_K = 1, # must be positive value, but I want it small i.e. no real recruitment
res_death_K = 5000000, # carrying capacity set to way above starting number of resources
res_move_type = 0, # 0=no move,
res_death_type = 0, # no natural death
observe_type = 0, # 0=density-based sampling
times_observe = 1,
obs_move_type = 1, # uniform in any direction
res_min_age = 0, # age of resources before agents record/act on them
res_move_obs = FALSE, # trees don't move
plotting = FALSE,
res_consume = 0.08, # Trees have 8% impact on yield
# all genetic algorithm parameters left to default
move_agents = TRUE,
max_ages = 1000,
minimum_cost = 10,
user_budget = 200,
manager_budget = 200,
usr_budget_rng = 20, # introduce variation around the mean user budget (removes step pattern)
manage_target = 2000000,
RESOURCE_ini = 2000000,
culling = TRUE,
tend_crops = TRUE,
tend_crop_yld = 0.01, # tending crops increases yield by 1% - less than that of culling trees
stakeholders = 20,
land_ownership = TRUE,
public_land = 0.4,
manage_freq = 1,
group_think = FALSE
)
N1d_summary <- as.data.frame(gmse_table(N1d))
write.csv(N1d_summary, file="outputs/investment/null_scenarios/N1/N1d_summary.csv")
rm(N1d)
## N1e ####
# After comparing N1:N1d below, there appeared to be a difference between N1:N1a, and N1b:N1d. The latter resulted in users felling as often as possible (see the non-zero plateau in the plot of number of cull actions), whereas the former resulted in the users sometimes choosing to tend crops. N1b:N1d have parameter settings with distance between the values of 0.06, 0.07, and 0.08. N1 and N1a have parameter distances of 0 and 0.03. Therefore, I want to test what happens when the distance between the parameter values are 0.04 and 0.05.
# I am going to keep tend_crop_yld low, at 0.01 and rather change the value of res_consume. This is becuase conceptually, just tending your existing crops is unlikely to have any dramatic increase on your yield beyond what you already produce on that existing land.
# below I will set res_consume to 0.05 and tend_crop_yld 0.01
N1e <- gmse(
time_max = 50,
land_dim_1 = 200,
land_dim_2 = 200, # landscape is 40,000ha or 400km2
res_movement = 0, # trees don't move
remove_pr = 0, # Assume no death
lambda = 0, # assume no growth
agent_view = 10,
agent_move = 50,
res_birth_K = 1, # must be positive value, but I want it small i.e. no real recruitment
res_death_K = 5000000, # carrying capacity set to way above starting number of resources
res_move_type = 0, # 0=no move,
res_death_type = 0, # no natural death
observe_type = 0, # 0=density-based sampling
times_observe = 1,
obs_move_type = 1, # uniform in any direction
res_min_age = 0, # age of resources before agents record/act on them
res_move_obs = FALSE, # trees don't move
plotting = FALSE,
res_consume = 0.05, # Trees have 5% impact on yield
# all genetic algorithm parameters left to default
move_agents = TRUE,
max_ages = 1000,
minimum_cost = 10,
user_budget = 200,
manager_budget = 200,
usr_budget_rng = 20, # introduce variation around the mean user budget (removes step pattern)
manage_target = 2000000,
RESOURCE_ini = 2000000,
culling = TRUE,
tend_crops = TRUE,
tend_crop_yld = 0.01, # tending crops increases yield by 1% - less than that of culling trees
stakeholders = 20,
land_ownership = TRUE,
public_land = 0.4,
manage_freq = 1,
group_think = FALSE
)
N1e_summary <- as.data.frame(gmse_table(N1e))
write.csv(N1e_summary, file="outputs/investment/null_scenarios/N1/N1e_summary.csv")
rm(N1e)
## N1f ####
# see explanation above in N1e
# Here res_consume is 0.06 and tend_crop_yld is 0.01
N1f <- gmse(
time_max = 50,
land_dim_1 = 200,
land_dim_2 = 200, # landscape is 40,000ha or 400km2
res_movement = 0, # trees don't move
remove_pr = 0, # Assume no death
lambda = 0, # assume no growth
agent_view = 10,
agent_move = 50,
res_birth_K = 1, # must be positive value, but I want it small i.e. no real recruitment
res_death_K = 5000000, # carrying capacity set to way above starting number of resources
res_move_type = 0, # 0=no move,
res_death_type = 0, # no natural death
observe_type = 0, # 0=density-based sampling
times_observe = 1,
obs_move_type = 1, # uniform in any direction
res_min_age = 0, # age of resources before agents record/act on them
res_move_obs = FALSE, # trees don't move
plotting = FALSE,
res_consume = 0.06, # Trees have 6% impact on yield
# all genetic algorithm parameters left to default
move_agents = TRUE,
max_ages = 1000,
minimum_cost = 10,
user_budget = 200,
manager_budget = 200,
usr_budget_rng = 20, # introduce variation around the mean user budget (removes step pattern)
manage_target = 2000000,
RESOURCE_ini = 2000000,
culling = TRUE,
tend_crops = TRUE,
tend_crop_yld = 0.01, # tending crops increases yield by 1% - less than that of culling trees
stakeholders = 20,
land_ownership = TRUE,
public_land = 0.4,
manage_freq = 1,
group_think = FALSE
)
N1f_summary <- as.data.frame(gmse_table(N1f))
write.csv(N1f_summary, file="outputs/investment/null_scenarios/N1/N1f_summary.csv")
rm(N1f)
## Comparison N1:N1f ####
# Here I want to compare the differences in yields at each time step between N1, where res_consume is 0.05 and tend_crops_yld is 0.02, N1a, where res_consume and tend_crop_yld are equal, N1b where res_consume is 0.08 and tend_crop_yld is 0.02, N1c where res_consume is 0.1 and tend_crop_yld is 0.02, N1d where res_consume is 0.08 and tend_crop_yld is 0.01, N1e where res_consume is 0.05 and tend_crop_yld is 0.01, and N1f where res_consume is 0.06 and tend_crop_yld is 0.01.
# No need to make the dataframe each time. Load the pre-made dataframe
yield.df <- read.csv("outputs/investment/null_scenarios/N1/N1_comparison_res_versus_tend/yield_df_N1-N1f.csv")
# This was just me working out how to find and extraact the yield from the simulation output object. Extract total yield for all cells from each timestep for all sims. This needs the simulation objects to be in the environment, so won't work unles you run all of the sims again
yield_ts <- sapply(1:length(N1a$land), function(x) sum(N1a$land[[x]][,,2]))
yield_ts_N1 <- sapply(1:length(N1$land), function(x) sum(N1$land[[x]][,,2]))
yield_ts_N1b <- sapply(1:length(N1b$land), function(x) sum(N1b$land[[x]][,,2]))
yield_ts_N1c <- sapply(1:length(N1c$land), function(x) sum(N1c$land[[x]][,,2]))
yield_ts_N1d <- sapply(1:length(N1d$land), function(x) sum(N1d$land[[x]][,,2]))
yield_ts_N1e <- sapply(1:length(N1e$land), function(x) sum(N1e$land[[x]][,,2]))
yield_ts_N1f <- sapply(1:length(N1f$land), function(x) sum(N1f$land[[x]][,,2]))
# load all simulation summaries
N1_summary <- read.csv("outputs/investment/null_scenarios/N1/N1_summary.csv", header=TRUE)
N1a_summary <- read.csv("outputs/investment/null_scenarios/N1/N1a_summary.csv", header=TRUE)
N1b_summary <- read.csv("outputs/investment/null_scenarios/N1/N1b_summary.csv", header=TRUE)
N1c_summary <- read.csv("outputs/investment/null_scenarios/N1/N1c_summary.csv", header=TRUE)
N1d_summary <- read.csv("outputs/investment/null_scenarios/N1/N1d_summary.csv", header=TRUE)
N1e_summary <- read.csv("outputs/investment/null_scenarios/N1/N1e_summary.csv", header=TRUE)
N1f_summary <- read.csv("outputs/investment/null_scenarios/N1/N1f_summary.csv", header=TRUE)
# create dataframe
sim <- c("N1","N1a","N1b", "N1c", "N1d","N1e","N1f")
yield.df <- data.frame(time_step = 1:50,
Simulation = rep(sim, each=50),
available_yld = 24000,
sim_yield = c(N1_summary$crop_yield, N1a_summary$crop_yield,
N1b_summary$crop_yield, N1c_summary$crop_yield,
N1d_summary$crop_yield, N1e_summary$crop_yield,
N1f_summary$crop_yield),
trees = c(N1_summary$resources, N1a_summary$resources, N1b_summary$resources,
N1c_summary$resources, N1d_summary$resources, N1e_summary$resources,
N1f_summary$resources))
# add % yield to df
yield.df$perc_yld <- yield.df$sim_yield/yield.df$available_yld*100
# rename sim columns
yield.df <- yield.df %>% rename(Simulation = sim)
# save dataframe
#write.csv(yield.df, file="outputs/investment/null_scenarios/N1/yield_df_N1-N1f.csv")
# first make a plot with all of the different sims and the parameter values
sim.para <- data.frame(Simulation = c("N1","N1a","N1b", "N1c", "N1d","N1e","N1f"),
res_consume = c(0.05,0.05,0.08,0.1,0.08,0.05,0.06),
tend_crop_yield = c(0.02,0.05,0.02,0.02,0.01,0.01,0.01))
p.sims <- ggplot(sim.para, aes(x=res_consume, y=tend_crop_yield, color=Simulation))+
geom_point(size=7)+
theme_classic()+
theme(axis.title = element_text(size=17),
axis.text = element_text(size=15),
legend.text = element_text(size=15),
legend.title = element_text(size=15))
#ggsave("outputs/investment/null_scenarios/N1/N1_comparison_res_versus_tend/param_values.png", p.sims,
# width = 25, height=20, unit="cm", dpi=300)
# plot
p1 <- ggplot(yield.df, aes(x=time_step, y=trees, group=Simulation, color=Simulation))+
geom_line(size=2)+
theme_classic()+
ylab("Number of trees")+
xlab("Time step")+
theme(legend.text = element_text(size=15),
legend.title = element_text(size=15),
axis.text = element_text(size=15),
axis.title = element_text(size=15))
#ggsave("outputs/investment/null_scenarios/N1/N1_comparison_res_versus_tend/trees_comparison.png", p1,
# width = 25, height=20, unit="cm", dpi=300)
# N1a results in the fewest trees being lost, which makes sense as N1a has no incentive to fell trees (equal parameter values). N1 and N1b:f are all fairly similar in their loss of trees. Interestingly, the simulation with the highest res_consume (N1c) does not end up with the fewest trees, and that must be because tend_crop_yield is higher than some of the others and so users will be more likely to choose to tend crops when costs of felling are very high. The simulation with the most trees lost is N1d, where tend_crop_yield is very low (0.01) and res_consume is quite high (0.08). This is closely followed by N1f which although has a lower res_consume than N1b and N1c, it also has a lower tend_crop_yield. This quite nicely shows the interaction between the two parameters I think. Essentially, I think this shows that small incremental changes in tend_crop_yield are actually more influential than similar increases in res_consume.
p2 <- ggplot(yield.df, aes(x=time_step, y=perc_yld, group=Simulation, color=Simulation))+
geom_line(size=2)+
theme_classic()+
ylab("% Yield")+
xlab("Time step")+
theme(legend.text = element_text(size=15),
legend.title = element_text(size=15),
axis.text = element_text(size=15),
axis.title = element_text(size=15))
ggsave("outputs/investment/null_scenarios/N1/N1_comparison_res_versus_tend/yield_comparison.png", p2,
width = 25, height=20, unit="cm", dpi=300)
# Yield is lowest in N1c, as this simulation has the highest value for res_consume (0.1), followed by N1b and N1d (0.08). N1f then sits on it's own in the middle (0.06). The highest yields are for N1, N1a, and N1e, where res_consume is 0.05 for all. For this last group, we see that N1e is increasing slightly faster, as tend_crop_yld is set lower than N1 and N1a, and so users are more keen to fell trees as tending crops has less value.
## more comparison plots between N1, N1a, N1b, N1c, N1d, N1e, N1f
# combine the simulation summaries for easier plotting
N1_summary$Simulation <- "N1"
N1a_summary$Simulation <- "N1a"
N1b_summary$Simulation <- "N1b"
N1c_summary$Simulation <- "N1c"
N1d_summary$Simulation <- "N1d"
N1e_summary$Simulation <- "N1e"
N1f_summary$Simulation <- "N1f"
Nx_summary <- rbind(N1_summary, N1a_summary, N1b_summary, N1c_summary, N1d_summary, N1e_summary, N1f_summary)
# save
#write.csv(Nx_summary, file="outputs/investment/null_scenarios/N1/Nx_summary.csv")
# load
Nx_summary <- read.csv("outputs/investment/null_scenarios/N1/N1_comparison_res_versus_tend/Nx_summary.csv")
Nx_summary <- Nx_summary %>% rename(Simulation = sim)
# number of cull actions
p_cull <- ggplot(Nx_summary, aes(x=time_step, y=act_culling, group=Simulation, colour=Simulation))+
geom_line(size=1.5)+
xlab("Time step")+
ylab("Number of cull actions")+
theme_classic()+
theme(legend.text = element_text(size=15),
legend.title = element_text(size=15),
axis.text = element_text(size=15),
axis.title = element_text(size=15))
ggsave("outputs/investment/null_scenarios/N1/N1_comparison_res_versus_tend/cull_comparison.png", p_cull,
width = 25, height=20, unit="cm", dpi=300)
# This plot has changed significantly since the introduction of the new res_death_type. The simulations do appear to be broadly split between N1 and N1a, and the rest. N1 and N1a show much more variation in the number of cull actions, with the number of cull actions regulalry dropping to 0. N1a, which has the two parameters set equally at 0.05, we see regular spells of 0 cull actions, where the users are choosing to tend crops (as that produces the same benefits in terms of yield). This happens less frequently with N1, and in N1 there are no occasions when number of cull actions remain at 0 for more than a single time step. This is becuase tend_crop_yield is lower than res_consume, and so it is more beneficial to fell trees. For all of the other simualtions though, there appears to be a minimum number of culls below which they never drop (just over 100). Even N1e, which is very similar to N1 in terms of parameters, never drops below a certain value of cull actions.
p_cost <- ggplot(Nx_summary, aes(x=time_step, y=cost_culling, group=Simulation, colour=Simulation))+
geom_line(size=1.5)+
xlab("Time step")+
ylab("Cost of cull actions")+
theme_classic()+
theme(legend.text = element_text(size=15),
legend.title = element_text(size=15),
axis.text = element_text(size=15),
axis.title = element_text(size=15))
# It's quite difficult to see what is happening in this plot, but as would be expected the manager is adapting the cost of culling regulalry to try and prevent culling, which is being attempted in all scenarios.
N1_Nf_tests <- p.sims / (p1 + p2) / (p_cull + p_cost)
ggsave("outputs/investment/null_scenarios/N1/N1_Nf_param_tests.png", N1_Nf_tests,
width = 30, height = 20, units="cm", dpi=300)
## N1x - reduced landscape ####
# I want to see what the results look like when I remove the public land, reduce the size of the landscape. I am not sure what purpose the public land is serving, and at the moment it is making the landscape larger, and increasing the number of trees, which is making seeing differences in tree loss more difficult.
# I will try a landscape ppf 150 x 150 cells. This results in 22,500 cells (or ha). With 20 villages, this results in 1,125 ha or 11.25km2 per village. This seems plausible. This means the number of trees will be 1,125,000
N1x <- gmse(
time_max = 50,
land_dim_1 = 150,
land_dim_2 = 150, # landscape is 22,500ha or 22.5km2
res_movement = 0, # trees don't move
remove_pr = 0, # Assume no death
lambda = 0, # assume no growth
agent_view = 10,
agent_move = 50,
res_birth_K = 1, # must be positive value, but I want it small i.e. no real recruitment
res_death_K = 5000000, # carrying capacity set to way above starting number of resources
res_move_type = 0, # 0=no move,
res_death_type = 0, # no natural death
observe_type = 0, # 0=density-based sampling
times_observe = 1,
obs_move_type = 1, # uniform in any direction
res_min_age = 0, # age of resources before agents record/act on them
res_move_obs = FALSE, # trees don't move
plotting = FALSE,
res_consume = 0.05, # Trees have 5% impact on yield
# all genetic algorithm parameters left to default
move_agents = TRUE,
max_ages = 1000,
minimum_cost = 10,
user_budget = 200,
manager_budget = 200,
usr_budget_rng = 20, # introduce variation around the mean user budget (removes step pattern)
manage_target = 1125000,
RESOURCE_ini = 1125000,
culling = TRUE,
tend_crops = TRUE,
tend_crop_yld = 0.02, # tending crops increases yield by 2% - less than that of culling trees
stakeholders = 20,
land_ownership = TRUE,
public_land = 0,
manage_freq = 1,
group_think = FALSE
)
plot_gmse_results(sim_results=N1c)
plot_gmse_effort(sim_results = N1c)
N1c_summary <- as.data.frame(gmse_table(N1c))
#write.csv(N1c_summary, file="outputs/investment/null_scenarios/N1/N1c_summary.csv")
# plots
time_cost_N1c <- ggplot(N1c_summary, aes(x=time_step,y=cost_culling))+
geom_line()+
theme_classic()
time_cull_N1c <- ggplot(N1c_summary, aes(x=time_step, y=act_culling))+
geom_line()+
theme_classic()
time_res_N1c <- ggplot(N1c_summary, aes(x=time_step, y=resources))+
geom_line()+
ylim(0,1125000)+
theme_classic()
time_yield_N1c <- ggplot(N1c_summary, aes(x=time_step, y=crop_yield))+
geom_line()+
theme_classic()
time_cost_N1c + time_cull_N1c + time_res_N1c + time_yield_N1c
## N1y - test parameters ####
# following Brad's advice, I am going to test a simulation where res_consume=0, tend_crop_yield=0.01 and user_budget=1. This is just to make sure everything behaves as expected. What I expect is that there will be no cull actions, becuase there is no incentive, nor budget, to cull trees. There will also be no tending of crops because the user budget is lower than the minimum cost
N1y <- gmse(
time_max = 50,
land_dim_1 = 150,
land_dim_2 = 150, # landscape is 22,500ha or 22.5km2
res_movement = 0, # trees don't move
remove_pr = 0, # Assume no death
lambda = 0, # assume no growth
agent_view = 10,
agent_move = 50,
res_birth_K = 1, # must be positive value, but I want it small i.e. no real recruitment
res_death_K = 5000000, # carrying capacity set to way above starting number of resources
res_move_type = 0, # 0=no move,
res_death_type = 0, # no natural death
observe_type = 0, # 0=density-based sampling
times_observe = 1,
obs_move_type = 1, # uniform in any direction
res_min_age = 0, # age of resources before agents record/act on them
res_move_obs = FALSE, # trees don't move
plotting = FALSE,
res_consume = 0, # Trees have 0% impact on yield
# all genetic algorithm parameters left to default
move_agents = TRUE,
max_ages = 1000,
minimum_cost = 10,
user_budget = 1,
manager_budget = 200,
usr_budget_rng = 1, # introduce variation around the mean user budget (removes step pattern)
manage_target = 1125000,
RESOURCE_ini = 1125000,
culling = TRUE,
tend_crops = TRUE,
tend_crop_yld = 0.01, # tending crops increases yield by 1% - less than that of culling trees
stakeholders = 20,
land_ownership = TRUE,
public_land = 0,
manage_freq = 1,
group_think = FALSE
)
N1y_summary <- read.csv("outputs/investment/null_scenarios/N1/N1_comparison_res_versus_tend/N1y_summary.csv")
## N1z - test parameters ####
N1z <- gmse(
time_max = 50,
land_dim_1 = 150,
land_dim_2 = 150, # landscape is 22,500ha or 22.5km2
res_movement = 0, # trees don't move
remove_pr = 0, # Assume no death
lambda = 0, # assume no growth
agent_view = 10,
agent_move = 50,
res_birth_K = 1, # must be positive value, but I want it small i.e. no real recruitment
res_death_K = 5000000, # carrying capacity set to way above starting number of resources
res_move_type = 0, # 0=no move,
res_death_type = 0, # no natural death
observe_type = 0, # 0=density-based sampling
times_observe = 1,
obs_move_type = 1, # uniform in any direction
res_min_age = 0, # age of resources before agents record/act on them
res_move_obs = FALSE, # trees don't move
plotting = FALSE,
res_consume = 0, # Trees have 0% impact on yield
# all genetic algorithm parameters left to default
move_agents = TRUE,
max_ages = 1000,
minimum_cost = 10,
user_budget = 10,
manager_budget = 200,
usr_budget_rng = 1, # introduce variation around the mean user budget (removes step pattern)
manage_target = 1125000,
RESOURCE_ini = 1125000,
culling = TRUE,
tend_crops = TRUE,
tend_crop_yld = 0.01, # tending crops increases yield by 1% - less than that of culling trees
stakeholders = 20,
land_ownership = TRUE,
public_land = 0,
manage_freq = 1,
group_think = FALSE
)
N1z_summary <- as.data.frame(gmse_table(N1z))
write.csv(N1z_summary, file="outputs/investment/null_scenarios/N1/N1z_summary.csv")
#### N2 ####
# The below calls are for the N2a and N2b null scenarios (see details above in "NULL SCENARIOS"). The landscape has been set up to match the landscape in the final scenarios.
# N2c and N2d were used to check the effects of altering the number of cells the manager could see during observation. These can be ignored now as I have set the parameters to ensure perfect observation.
## N2a ####
# user budget decreases linearly, manager budget remains constant. Observation type changed to transect, and agent_view increased so no observation error
UB <- 500
UBR <- 50
N2a_sim_old <- gmse_apply(
res_mod = resource,
obs_mod = observation,
man_mod = manager,
use_mod = user,
get_res = "FUll",
time_max = 50,
land_dim_1 = 100,
land_dim_2 = 100, # landscape is 22,500ha or 22.5km2
res_movement = 0, # trees don't move
remove_pr = 0, # Assume no death
lambda = 0, # assume no growth
agent_view = 150,
agent_move = 50,
res_birth_K = 1, # must be positive value, but I want it small i.e. no real recruitment
res_death_K = 5000000, # carrying capacity set to way above starting number of resources
res_move_type = 0, # 0=no move,
res_death_type = 0, # no natural death
observe_type = 2, # transect
times_observe = 1,
obs_move_type = 1, # uniform in any direction
res_min_age = 0, # age of resources before agents record/act on them
res_move_obs = FALSE, # trees don't move
plotting = FALSE,
res_consume = 0.08, # Trees have 8% impact on yield
# all genetic algorithm parameters left to default
move_agents = TRUE,
max_ages = 1000,
minimum_cost = 10,
user_budget = UB,
manager_budget = 500,
usr_budget_rng = UBR, # introduce variation around the mean user budget (removes step pattern)
manage_target = 100000,
RESOURCE_ini = 100000,
culling = TRUE,
tend_crops = TRUE,
tend_crop_yld = 0.01, # tending crops increases yield by 1% - less than that of culling trees
stakeholders = 30,
land_ownership = TRUE,
public_land = 0,
manage_freq = 1,
group_think = FALSE
)
# matrix for results
N2a <- matrix(data=NA, nrow=50, ncol=6)
# loop the simulation.
for(time_step in 1:50){
sim_new <- gmse_apply(get_res = "Full", old_list = N2a_sim_old, user_budget=UB,
usr_budget_rng = UBR)
N2a[time_step, 1] <- time_step
N2a[time_step, 2] <- sim_new$basic_output$resource_results[1]
N2a[time_step, 3] <- sim_new$basic_output$observation_results[1]
N2a[time_step, 4] <- sim_new$basic_output$manager_results[3]
N2a[time_step, 5] <- sum(sim_new$basic_output$user_results[,3])
N2a[time_step, 6] <- UB
N2a_sim_old <- sim_new
UB <- UB - 3
UBR <- UB/5
print(time_step)
}
colnames(N2a) <- c("Time", "Trees", "Trees_est", "Cull_cost", "Cull_count",
"User_budget")
N2a_summary <- data.frame(N2a)
write.csv(N2a_summary, file = "outputs/investment/null_scenarios/N2/N2a_summary.csv")
## N2b ####
# User budget remains constant, manager budge decreases linearly
MB <- 500
N2b_sim_old <- gmse_apply(
res_mod = resource,
obs_mod = observation,
man_mod = manager,
use_mod = user,
get_res = "FUll",
time_max = 50,
land_dim_1 = 100,
land_dim_2 = 100, # landscape is 22,500ha or 22.5km2
res_movement = 0, # trees don't move
remove_pr = 0, # Assume no death
lambda = 0, # assume no growth
agent_view = 150,
agent_move = 50,
res_birth_K = 1, # must be positive value, but I want it small i.e. no real recruitment
res_death_K = 5000000, # carrying capacity set to way above starting number of resources
res_move_type = 0, # 0=no move,
res_death_type = 0, # no natural death
observe_type = 2, # transect
times_observe = 1,
obs_move_type = 1, # uniform in any direction
res_min_age = 0, # age of resources before agents record/act on them
res_move_obs = FALSE, # trees don't move
plotting = FALSE,
res_consume = 0.08, # Trees have 8% impact on yield
# all genetic algorithm parameters left to default
move_agents = TRUE,
max_ages = 1000,
minimum_cost = 10,
user_budget = 500,