-
Notifications
You must be signed in to change notification settings - Fork 5
/
Seminar_09_01_VoI.Rmd
249 lines (198 loc) · 7.47 KB
/
Seminar_09_01_VoI.Rmd
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
## Seminar 9.1: Value of Information {#voi_1}
<!-- reference with [Models](#value_of_information_seminar_part2) -->
Welcome to the second part of the Value of Information seminar for the course **Decision Analysis and Forecasting for Agricultural Development**.Feel free to bring up any questions or concerns in the Slack or to [Dr. Cory Whitney](mailto:[email protected]?subject=[Seminar_4]%20Decision%20Analysis%20Lecture) or the course tutor.
### Practical Value of Information example:
Let's use a modified version of sheep apple agroforestry from Seminar 7 (#model_programming_cont). Just for this class let's presume that the farmer has a grassland where they graze sheep. Intervention: implement apple agroforestry. Let's prepare our input estimates for the system.
```{r libraries_voi_part2, eval=FALSE, echo=TRUE}
library(ggplot2)
library(decisionSupport)
```
```{r input_table_voi_test, exercise=TRUE}
input_estimates <- data.frame(
variable = c(
"sheep_income",
"sheep_cost",
"apple_income",
"apple_cost",
"discount_rate"
),
lower = c(3000, 1000, 30000, 15000, 10),
median = NA,
upper = c(5000, 2500, 60000, 30000, 10),
distribution = c("posnorm", "posnorm", "posnorm", "posnorm", "const"),
label = c(
"Income from sheep (euro/year)",
"Cost of sheep (euro/year)",
"Income from apple (euro/year)",
"Cost of apple (euro/year)",
"Discount Rate"
)
)
# show the result
input_estimates
```
```{r input_table_voi_test-solution}
input_estimates <- data.frame(
variable = c(
"sheep_income",
"sheep_cost",
"apple_income",
"apple_cost",
"discount_rate"
),
lower = c(3000, 1000, 30000, 15000, 10),
median = NA,
upper = c(5000, 2500, 60000, 30000, 10),
distribution = c("posnorm", "posnorm", "posnorm", "posnorm", "const"),
label = c(
"Income from sheep (euro/year)",
"Cost of sheep (euro/year)",
"Income from apple (euro/year)",
"Cost of apple (euro/year)",
"Discount Rate"
)
)
# show the result
input_estimates
```
If you want you can write the input estimates into a local file on your machine.
```{r, eval=FALSE, echo=TRUE}
# write the results to a .csv file
df <- data.frame(input_estimates)
write.csv(df, file = "input_table.csv", row.names = FALSE)
```
In the next step, let's build a function to model the system and calculate the NPV of the AF system vs only sheep. Simulate it using the `mcSimulation` function from the `decisionSupport` package with 5000 'numberOfModelRuns'. Using 'plot_distributions' plot NPV of intervention trade off vs no intervention.
```{r apple_sheep_voi_test, exercise=TRUE}
model_function <- function() {
# Estimate the income in a normal season
AF_income <- sheep_income + apple_income
AF_cost <- sheep_cost + apple_cost
# Estimate the final results from the model
AF_final_result <- AF_income - AF_cost
# baseline with sheep only
sheep_only <- sheep_income - sheep_cost
# should I plant trees in the sheep pastures?
Decision_benefit <- AF_final_result - sheep_only
#Calculating NPV
#AF System
AF_NPV <- discount(AF_final_result,
discount_rate = discount_rate,
calculate_NPV = TRUE)#NVP of AF system
NPV_sheep_only <- discount(sheep_only,
discount_rate = discount_rate,
calculate_NPV = TRUE) #NVP of grassland
NPV_decision <- discount(Decision_benefit,
discount_rate = discount_rate,
calculate_NPV = TRUE)
# Generate the list of outputs from the Monte Carlo simulation
return(
list(
NPV_Agroforestry_System = AF_NPV,
NPV_Treeless_System = NPV_sheep_only,
NPV_decision = NPV_decision
)
)
}
```
```{r apple_sheep_voi_test-solution}
model_function <- function() {
# Estimate the income in a normal season
AF_income <- sheep_income + apple_income
AF_cost <- sheep_cost + apple_cost
# Estimate the final results from the model
AF_final_result <- AF_income - AF_cost
sheep_only <- sheep_income - sheep_cost
Decision_benefit <- AF_final_result - sheep_only
#Calculating NPV
#AF System
AF_NPV <- discount(AF_final_result,
discount_rate = discount_rate,
calculate_NPV = TRUE)#NVP of AF system
NPV_sheep_only <- discount(sheep_only,
discount_rate = discount_rate,
calculate_NPV = TRUE) #NVP of grassland
NPV_decision <- discount(Decision_benefit,
discount_rate = discount_rate,
calculate_NPV = TRUE)
# Generate the list of outputs from the Monte Carlo simulation
return(
list(
NPV_Agroforestry_System = AF_NPV,
NPV_Treeless_System = NPV_sheep_only,
NPV_decision = NPV_decision
)
)
}
apple_sheep_mc_simulation <- mcSimulation(
estimate = as.estimate(input_estimates),
model_function = model_function,
numberOfModelRuns = 80,
functionSyntax = "plainNames"
)
plot_distributions(
mcSimulation_object = apple_sheep_mc_simulation,
vars = c("NPV_decision", "NPV_Treeless_System"),
method = 'smooth_simple_overlay',
base_size = 7,
x_axis_name = "Outcome of AF intervention",
scale_x_continuous(
labels = function(x)
x / 100000
),
ggtitle("Net Present Value of the system"),
legend.position = "bottom"
)
```
Now, let's perform sensitivity analysis using Partial Least Squares (PLS) regression using 'plsr' function. Plot the results showing a `cut_off_line` at 1 and `threshold` of 0.5.
```{r apple_sheep_pls_test, exercise=TRUE}
pls_result_AF <- plsr.mcSimulation(
object = apple_sheep_mc_simulation,
resultName = names(apple_sheep_mc_simulation$y)[3],
ncomp = 1
)
```
```{r apple_sheep_pls_test-solution}
pls_result_AF <- plsr.mcSimulation(
object = apple_sheep_mc_simulation,
resultName = names(apple_sheep_mc_simulation$y)[3],
ncomp = 1
)
plot_pls(pls_result_AF,
input_table = input_estimates,
cut_off_line = 1,
threshold = 0.5)
```
Save it on your machine with `ggsave`.
```{r, eval=FALSE, echo=TRUE}
plot_pls(pls_result_AF,
input_table = input_estimates,
cut_off_line = 1,
threshold = 0.5)
ggsave(filename = "PLS_sheep_vs_AppleAF.png",
plot = last_plot(),
width = 5,
height = 3)
```
Now, let's quantify Value of Information by calculating expected value of perfect information (EVPI). Save the results of your MC simulation to a data frame.
```{r assign_results_voi, exercise=TRUE}
# assign the results of x and y to a data frame
df <- data.frame(apple_sheep_mc_simulation$x, apple_sheep_mc_simulation$y[1:3])
```
```{r assign_results_voi-solution}
# assign the results of x and y to a data frame
df <- data.frame(apple_sheep_mc_simulation$x, apple_sheep_mc_simulation$y[1:3])
```
Generate the Variable Importance in the Projection (VIP) scores of the PLS for the input variables based on the apple sheep model.
```{r apple_sheep_evpi_test, exercise=TRUE}
# run evpi on the NPVs (the last few in the table, starting with "NPV_decision")
EVPI <- multi_EVPI(mc = df, first_out_var = "NPV_decision")
# plot the EVPI results for the decision
plot_evpi(EVPIresults = EVPI, decision_vars = "NPV_decision")
```
```{r apple_sheep_evpi_test-solution}
# run evpi on the NPVs (the last few in the table, starting with "NPV_decision")
EVPI <- multi_EVPI(mc = df, first_out_var = "NPV_decision")
# plot the EVPI results for the decision
plot_evpi(EVPIresults = EVPI, decision_vars = "NPV_decision")
```
Share your result with us in Slack and explain the result.