---
title: "From the simulation of forest plantation dynamics to the quantification of bark-stripping damage by ungulates : Analyses"
author: Gauthier Ligot
format: pdf
keep-tex: true
execute:
  echo: false
  warning: false
editor_options: 
  chunk_output_type: inline
output: 
  pdf_document:
    extra_dependencies: ["float"]
---

```{r}
#| message: false
#| warning: false
options(knitr.kable.NA = '')
knitr::opts_chunk$set(fig.pos = "H", out.extra = "")
library(knitr)
library(tidyverse)
library(ggpubr)
library(RColorBrewer)
library(car)
load("tmp/prepared_data.Rdata")
source("configs/myGGplotConfigs.R")
options(dplyr.summarise.inform = FALSE)

treatmentcolors<-brewer.pal(n = 8, name = 'Dark2')[c(1,2,3,8)]
sicolors<-brewer.pal(n=9,name = 'Greens')[c(4,6,7,8,9)]
```

## Results

### Stand characteristics over time

```{r}
#| label: fig-stand
#| fig-cap: Changes in in the main characteristics of the simulated stands
#| fig-width: 10
#| fig-height: 7

si = 27
gdat<-filter(stand,siteIndex == si)

g1<-ggplot(gdat,aes(x=date, y=nha, lty=scenario_stand, col=bpr*100, group = interaction(scenario_stand,bpr,repetition)))+
  geom_line()+
  labs(col = "BSR (%)",lty = "")+
  ylab(expression(paste("Tree density ( ", ha^{-1} ,")")))+xlab("Plantation age (year)")+theme(legend.key=element_rect(fill="white"))

g2<-ggplot(gdat,aes(x=date, y=dg, lty=scenario_stand, col=bpr*100, 
                 group = interaction(scenario_stand,bpr,repetition)))+
  geom_line()+
  labs(col = "BSR (%)",lty = "")+
  ylab("Quadratic mean diameter (cm)")+xlab("Plantation age (year)")+theme(legend.key=element_rect(fill="white"))

g3<-ggplot(gdat,aes(x=date, y=gha, lty=scenario_stand, col=bpr*100, 
                 group = interaction(scenario_stand,bpr,repetition)))+
  geom_line()+
  labs(col = "BSR (%)",lty = "")+
  ylab(expression(paste("Basal area ( ", m^{2} ,"/ha)")))+xlab("Plantation age (year)")+theme(legend.key=element_rect(fill="white"))

g4<-ggplot(gdat,aes(x=date, y=vha, lty=scenario_stand, col=bpr*100, group = interaction(scenario_stand,bpr,repetition)))+
  geom_line()+
  labs(col = "BSR (%)",lty = "")+
  ylab(expression(paste("Standing volume ( ", m^{3}/ha ,")")))+xlab("Plantation age (year)")+theme(legend.key=element_rect(fill="white"))

g5<-ggplot(gdat,aes(x=date, y=rdi, lty=scenario_stand, col=bpr*100, 
                 group = interaction(scenario_stand,bpr,repetition)))+
  geom_line()+
  labs(col = "BSR (%)",lty = "")+
  ylab("Relative density index")+xlab("Plantation age (year)")+theme(legend.key=element_rect(fill="white"))


ggarrange(g1,g2,g3,g4,g5,ncol=3,nrow=2,common.legend = T)

```
Proportion of damaged trees increased with BSR

```{r}
filter(gdat,siteIndex == 27 & date == 25 & scenario_stand == "Unprotected")%>%
  group_by(bpr, scenario_stand)%>%
  summarise(damaged_prop = mean(peeled_prop),
            damaged_prop_sd = sd(peeled_prop))%>%
    kable(format.args = list(scientific = FALSE, digits=4))
```



Mean productivity across the simulations

```{r}
stand%>%
  group_by(siteIndex,bpr,scenario_case,repetition,scenario_stand,date)%>%
  summarize(harvested_vha = max(vha)-min(vha))%>%
  group_by(siteIndex,bpr,scenario_case,scenario_stand,repetition)%>%
  arrange(date)%>%
  mutate(Cumul_harvested_vha = cumsum(harvested_vha))%>%
  summarise(productivity = max(Cumul_harvested_vha)/max(date))%>%
  group_by(siteIndex)%>%
  summarise(productivity_mean = mean(productivity),
            productivity_sd = sd(productivity))%>%
  kable(format.args = list(scientific = FALSE, digits=4))
```

```{r}
#| label: fig-Ndamaged
#| fig-cap: Changes in the mean stand density
#| fig-width: 10

tmp<-stand%>%
  mutate(nha_d = nha_dp + nha_d0)%>%
  group_by(siteIndex,bpr,scenario_stand,date)%>%
  summarise_at(vars(nha_d,nha),.funs = mean)%>%
  gather(key="population",value="nha",nha,nha_d)%>%
  mutate(population = factor(population))%>%
  mutate(bpr_label = factor(paste("BSR =",bpr*100,"%")))%>%
  mutate(population_label = ifelse(population == "nha", "Total", "Damaged"))

g1<-
ggplot(filter(tmp,siteIndex==27 & bpr %in% c(0.02,0.05,0.1)),
       aes(x=date, y=nha, size = population_label, col=scenario_stand, 
           group=interaction(population_label,bpr,scenario_stand)))+
  geom_line()+
  facet_grid(~bpr_label)+
  ylab(expression(paste("Tree density ( ", ha^{-1} ,")")))+xlab("Plantation age (year)")+ 
  scale_color_manual(values = treatmentcolors)+
  scale_size_manual(values=c(1.5, 0.5))+
  guides(lty = guide_legend(order = 2, nrow = 1),
         color = guide_legend(order = 1, nrow = 1)) +
  theme(legend.title = element_blank(),legend.position="top", legend.box="vertical",legend.key=element_rect(fill="white"))

    # ggsave(file="c:/Data_G_current/TG3/1-manuscrits/2-scientific_paper/revision 1/figures/fig1.svg", plot=g1, width=10, height = 5)
```
Statistics of the decay column height
```{r}
results2%>%
  select(siteIndex:scenario,discountRate,rotation,maxnpv,stand2)%>%
  unnest(cols = c(stand2))%>%
  filter(scenario == "Unprotected" & bpr > 0)%>%
  slice(n()-1)%>%
  # group_by(bpr)%>%
  ungroup()%>%
  summarise(m = mean(purge_height),
            min = min(purge_height_min),
            max = max(purge_height_max))%>%
  kable(format.args = list(scientific = FALSE, digits=4))
```

### Losses in timber volume and value

```{r}
#| label: fig-decayed
#| fig-cap: Importance of the decayed wood across the simulations
#| fig-width: 10
#| fig-height: 8

stand2<-results2%>%
  select(siteIndex:scenario,discountRate,rotation,maxnpv,stand2)%>%
  unnest(cols = c(stand2))

g1<-stand2%>%
  slice(n()-1)%>%
  ggplot(aes(x=as_factor(bpr*100),y=peeled_prop*100,fill=scenario_stand,coll=scenario_stand))+
  geom_boxplot()+ggtitle("A.")+
  scale_fill_manual(values = treatmentcolors)+
  ylab("% of trees with damage")+xlab("Bark-stripping rate (%)")+
  theme(legend.position="top",legend.title = element_blank(),legend.key=element_rect(fill="white"))

g3<-stand2%>%
  slice(n()-1)%>%
  ggplot(aes(x=as_factor(bpr*100),y=(vha-vha_eco)/vha*100,fill=scenario_stand,coll=scenario_stand))+
  geom_boxplot()+ggtitle("B.")+
  scale_fill_manual(values = treatmentcolors)+
  ylab("Volume % with decay")+xlab("Bark-stripping rate (%)")+
  theme(legend.position="top",legend.title = element_blank(),legend.key=element_rect(fill="white"))

g2<-stand2%>%
  filter(siteIndex == 27 & scenario == "Unprotected" & bpr == 0.05)%>%
  mutate(purge_height_min = ifelse(purge_height_min == 999, 0, purge_height_min))%>%
  group_by(date, scenario_case)%>%
  summarise(purge_height = mean(purge_height),
            purge_height_min = mean(purge_height_min),
            purge_height_max = mean(purge_height_max))%>%
  ggplot(aes(x=date,y=purge_height))+
  geom_line()+ggtitle("C.")+
  geom_line(aes(y = purge_height_min),col="grey")+
  geom_line(aes(y = purge_height_max),col="grey")+
  ylab("Height of the decay column (m)")+xlab("Plantation age (year)")+
  annotate(geom="text", x=68, y=2.3, label="mean",color="black")+
  annotate(geom="text", x=68, y=4.25, label="max",color="grey")+
  annotate(geom="text", x=68, y=1.2, label="min",color="grey")

ggarrange(g1, g3, g2 ,ncol=2,nrow=2,widths = c(2,2), common.legend = T)


```

```{r}
#| label: fig-losses
#| fig-cap: Undamaged volume and timber value
#| fig-width: 7
#| fig-height: 4
#| 
gdat<-stand%>%
  filter(siteIndex == 27)%>%
  group_by(siteIndex,bpr,scenario_case,repetition,scenario_stand,date)%>%
  summarize(harvested_vha = max(vha)-min(vha),
            harvested_vha_eco = max(vha_eco)-min(vha_eco),
            harvested_value = max(stumpage)-min(stumpage))%>%
  group_by(siteIndex,bpr,scenario_case,scenario_stand,repetition)%>%
  summarise(productivity = sum(harvested_vha)/max(date),
            productivity_eco = sum(harvested_vha_eco)/max(date),
            volume = sum(harvested_vha_eco),
            stumpage = sum(harvested_value))

g1<-ggplot(gdat,aes(x=bpr*100,y=volume,col=scenario_stand,pch=scenario_stand))+
  geom_point()+geom_smooth(method="lm",se=F)+
  xlab("Bark-stripping rate (%)")+ylab("Undamaged timber volume (m³/ha)")+
    scale_color_manual(values = treatmentcolors)+
  theme(legend.title = element_blank(),legend.key=element_rect(fill="white"))#+facet_wrap(~siteIndex)

g2<-ggplot(gdat,aes(x=bpr*100,y=stumpage,col=scenario_stand,pch=scenario_stand))+
  geom_point()+geom_smooth(method="lm",se=F)+
  xlab("Bark-stripping rate (%)")+ylab("Timber value (€/ha)")+#+facet_wrap(~siteIndex)+
    scale_color_manual(values = treatmentcolors)+
  theme(legend.title = element_blank(),legend.key=element_rect(fill="white"))#+facet_wrap(~siteIndex)


ggarrange(g1,g2,nrow = 1,ncol = 2,common.legend = TRUE)
```

### Net present value
#### Optimum rotation length

```{r}
tmp_rotation<-filter(results_maxnpv, bpr == 0 & discountRate == 0.02)%>%
  group_by(siteIndex)%>%
  summarise(mean_rotation = mean(rotation),
            sd_rotation = sd(rotation))

tmp_rotation2<-filter(results_maxnpv, siteIndex == 27 & bpr == 0 & discountRate == 0.02)%>%
  group_by(scenario)%>%
  summarise(mean_rotation = mean(rotation),
            sd_rotation = sd(rotation))

tmp_rotation3<-filter(results_maxnpv, siteIndex == 27 & bpr == 0.10 & discountRate == 0.02)%>%
  group_by(scenario)%>%
  summarise(mean_rotation = mean(rotation),
            sd_rotation = sd(rotation))
```


```{r}
filter(tmp_rotation,siteIndex %in% c(21,33)) %>%
  kable()
```

```{r}
tmp_rotation2%>%
  kable()
```



#### Opportunity cost of bark-stripping damages

Statistics of the opportunity cost

```{r}
opcost<-filter(results_maxnpv, scenario_case == 0)%>%
  group_by(siteIndex,discountRate,repetition)%>%
  arrange(bpr)%>%
  mutate(opportunity_cost = first(maxnpv) - maxnpv)%>%
  ungroup()
summary(opcost$opportunity_cost)
```

```{r}
filter(opcost,siteIndex == 27 & discountRate == 0.02 & bpr %in% c(0,0.1))%>%
  group_by(bpr,scenario)%>%
  summarise(npv = mean(maxnpv),
            cost = mean(opportunity_cost))%>%
  kable()
```


```{r}
model<-lm(opportunity_cost ~ bpr + discountRate:bpr + siteIndex:bpr -1,data=opcost)
smry<-summary(model)
kable(smry$coefficients, format.args = list(scientific = TRUE, digits=3))
```


```{r}
#| label: fig-opcost
#| fig-cap : changes in opportunity cost with BSR  
#| fig-width: 10
#| fig-height: 6
tmp<-filter(results_maxnpv, scenario_case == 0)%>%
  group_by(siteIndex,discountRate,repetition)%>%
  arrange(bpr)%>%
  mutate(opportunity_cost = first(maxnpv) - maxnpv)%>%
  ungroup()%>%
  mutate(r_label = paste0("r ==",discountRate))

ggplot(data=tmp,aes(y=opportunity_cost, x = bpr*100, col=factor(siteIndex), group=siteIndex))+
  geom_point()+
  stat_smooth(method="lm",se=FALSE)+
  facet_grid(~r_label,labeller = "label_parsed")+
  theme(legend.position = "top")+
  labs(col = "Site index")+xlab("Bark-stripping rate (%)")+
  scale_color_manual(values = sicolors)+
  ylab("Opportunity cost (€/ha)")+
  geom_abline(intercept = 6000, col="grey")


```



```{r}
#| label: fig-opcosta
#| fig-cap: Changes in cost annuity with BSR
#| fig-width: 10
#| fig-height: 6
ggplot(data=tmp,aes(y=opportunity_cost*discountRate, x = bpr*100, col=factor(siteIndex), group=siteIndex))+
  geom_point()+
  stat_smooth(method="lm",se=F)+
  facet_grid(~r_label,labeller = "label_parsed")+
  theme(legend.position = "top")+
  labs(col = "Site index")+xlab("Bark-stripping rate (%)")+
  scale_color_manual(values = sicolors)+
  ylab("Opportunity cost annuity (€/ha/an)")
```

#### Net present value across treatments

```{r}
model<-lm(maxnpv ~ scenario + siteIndex + discountRate + bpr + discountRate:bpr + scenario:bpr, data=results_maxnpv)
smry<-summary(model)
# anova(model)
kable(smry$coefficients, format.args = list(scientific = TRUE, digits=4))

# library(xtable)
# xtable(smry$coefficients,digits=3,display=c("s",rep("e",4)))
```

```{r}
#| label: fig-npv
#| fig-cap: NPV across BSR
#| fig-width: 10
#| fig-height: 7
gdat<-filter(results_maxnpv, siteIndex %in% c(21,27,33) & discountRate %in% c(0.01,0.03))%>%
    mutate(r_label = paste0("r ==",discountRate),
         si_label = paste0("SI ==",siteIndex,"~m"))

ggplot(gdat,aes(x=bpr*100, y=maxnpv, col=scenario, pch=scenario))+
  geom_point()+
  stat_smooth(method="lm",se=F)+#geom_point()+
  facet_grid(r_label~si_label,scales="free_y", labeller = "label_parsed")+
  geom_abline(intercept = 0,slope = 0)+
  ylab("Net present value (€/ha)")+xlab("Bark-stripping rate (%)")+
  scale_color_manual(values = treatmentcolors)+
  theme(legend.position = "top", legend.title=element_blank(),legend.key=element_rect(fill="white"))
```

## Discussion

### Model assumption

Tree size at initialization
```{r}
filter(stand,date == 8)%>%
  group_by(siteIndex)%>%
  summarise(dg_m = mean(dg),
            dg_sd = sd(dg))%>%
  kable()
```

### How much does bark-stripping alter timber production and profitability?

Loss at clear-cut :

```{r}
stand2%>%
  slice(n()-1)%>%
  filter(scenario == "Unprotected" & bpr %in% c(0.05,0.1) & siteIndex %in% c(27,33))%>%
  mutate(loss = (vha-vha_eco)/vha*100) %>%
  group_by(bpr,siteIndex)%>%
  summarise(loss_rel = mean(loss),
            peeled_prop = mean(peeled_prop))%>%
  kable()

```

mean cost for the study area

```{r}
filter(opcost, siteIndex %in% c(21,27,33) & discountRate %in% c(0.02) & scenario == "Unprotected" & bpr %in% c(0.04))%>%
  group_by(siteIndex)%>%
  summarise(cost = mean(opportunity_cost),
            cost_rel = mean(opportunity_cost/maxnpv*100),
            annuity = mean(opportunity_cost*discountRate))%>%
  kable()
```

```{r}
filter(results_maxnpv, siteIndex %in% c(21,27,33) & discountRate %in% c(0.04) & scenario == "Unprotected" & bpr %in% c(0,0.1))%>%
  group_by(bpr,siteIndex)%>%
  summarise(npv = mean(maxnpv))%>%
  group_by(siteIndex)%>%
  arrange(desc(npv))%>%
  mutate(pc = (lag(npv)-npv)/lag(npv))%>%
  kable()
  
```

```{r}
filter(results_maxnpv, siteIndex %in% c(27) & discountRate %in% c(0.02) & scenario == "Unprotected" & bpr %in% c(0,0.04))%>%
  group_by(bpr,siteIndex,discountRate)%>%
  summarise(npv = mean(maxnpv))%>%
  group_by(siteIndex,discountRate)%>%
  arrange(desc(npv))%>%
  mutate(cost = (lag(npv)-npv),
         loss = cost / lag(npv) * 100,
         annuity = cost*discountRate)%>%
  kable()
  
```

## Supplementary material

```{r}
#| label: fig-rotation
#| fig-cap: Optimum rotation across simulation
#| fig-width: 10
#| fig-height: 7
gdat<-filter(results_maxnpv, siteIndex %in% c(21,27,33) & discountRate %in% c(0.01,0.02,.03,0.04))%>%
  mutate(r_label = paste0("r ==",discountRate),
         si_label = paste0("SI ==",siteIndex,"~m"))
ggplot(gdat,aes(x=bpr*100, y=rotation, col=as_factor(scenario),pch=as_factor(scenario)))+
  geom_point()+stat_smooth(method="lm",se=F)+
  facet_grid(r_label~si_label,scales = "free", labeller = "label_parsed")+
  geom_abline(intercept = 0,slope = 0)+
  xlab("Bark-stripping rate (%)")+ylab("Optimum rotation length (years)")+
  scale_color_manual(values = treatmentcolors)+
  theme(legend.title = element_blank(),legend.position = "top",legend.key=element_rect(fill="white"))#+facet_wrap(~siteIndex)
```

```{r}
opcost<-filter(results_maxnpv, scenario_case == 0)%>%
  mutate(maxnpv2 = maxnpv * ((1+discountRate)^rotation-1) / ((1+discountRate)^rotation))%>%
  group_by(siteIndex,discountRate,repetition)%>%
  arrange(bpr)%>%
  mutate(opportunity_cost = first(maxnpv2) - maxnpv2)%>%
  ungroup()
model<-lm(opportunity_cost ~ bpr + discountRate:bpr + siteIndex:bpr -1,data=opcost)
smry<-summary(model)
kable(smry$coefficients, format.args = list(scientific = TRUE, digits=3))
```

```{r}
#| label: fig-Ndamaged-si
#| fig-cap: changes in density of damaged trees, the effect of site index.
#| fig-width: 7

tmp<-stand%>%
  mutate(nha_d = nha_dp + nha_d0)%>%
  select(siteIndex,bpr,scenario_stand,date,scenario_case,lofflerModel,repetition,nha,nha_d)%>%
  gather(key="population",value="nha",nha,nha_d)%>%
  mutate(population = factor(population))%>%
  mutate(bpr_label = factor(paste("BSR =",bpr*100,"%")))%>%
  mutate(population_label = ifelse(population == "nha", "Total", "Damaged"))

ggplot(filter(tmp,siteIndex==c(21,33) & bpr %in% c(0.01,0.05) & scenario_case == 2 & repetition == 0),
       aes(x=date, y=nha, lty = as.factor(siteIndex), col=population_label, 
           group=interaction(population_label,scenario_stand,siteIndex)))+
  geom_line()+
  facet_grid(~bpr_label)+
  ylab(expression(paste("Tree density ( ", ha^{-1} ,")")))+xlab("Plantation age (year)")+
    scale_linetype_manual(values=1:2,labels=expression(SI==21~m,SI==33~m))+
  guides(lty = guide_legend(order = 2, nrow = 1),
         color = guide_legend(order = 1, nrow = 1)) +
  scale_color_manual(values = c("darkgrey","black"))+
  theme(legend.position="top", legend.title = element_blank())
```

```{r}
#| label: fig-npvall
#| fig-cap: Changes in NPV across all simulations
#| fig-width: 11
#| fig-height: 11
gdat<-filter(results_maxnpv, siteIndex %in% c(21,24,27,30,33) & discountRate %in% c(0.01,0.02,0.03,0.04))%>%
    mutate(r_label = paste0("r ==",discountRate),
         si_label = paste0("SI ==",siteIndex,"~m"))

ggplot(gdat,aes(x=bpr*100, y=maxnpv, col=scenario, pch=scenario))+
  geom_point()+
  stat_smooth(method="lm",se=F)+#geom_point()+
  facet_grid(r_label~si_label,scales="free_y", labeller = "label_parsed")+
  geom_abline(intercept = 0,slope = 0)+
  ylab("Net present value (€/ha)")+xlab("Bark-stripping rate (%)")+
  scale_color_manual(values = treatmentcolors)+
  theme(legend.position = "top", legend.title=element_blank(),legend.key=element_rect(fill="white"))
```


```{r}
#| label: fig-ct
#| fig-cap: density of protected and damages trees
#| fig-width: 10
#| fig-height: 6

tmp<-filter(stand,scenario_case %in% c(1,2))%>%
  select(scenario_stand,siteIndex:date,nha,nha_dp,nha_0p)%>%
  group_by(scenario_stand,siteIndex,bpr,date)%>%
  summarise_at(vars(starts_with("NHA")),.funs = mean)%>%
  gather(key="population",value="nha",nha,nha_dp,nha_0p)%>%
  mutate(population = factor(population))%>%
  mutate(bpr_label = factor(paste("BPR =",bpr*100,"%")))%>%
  mutate(population_label = ifelse(population == "nha", "Total", ifelse(population == "nha_0p","Protected and not damaged","Protected and damaged")))

ggplot(filter(tmp,siteIndex==27 & bpr %in% c(0.02,0.05,0.07,0.1)),
       aes(x=date, y=nha, lty = population_label,col=population_label,  group = interaction(population_label,bpr)))+
  geom_line()+
  facet_grid(scenario_stand~bpr_label)+
  ylab(expression(paste("Tree density ( ", ha^{-1} ,")")))+xlab("Plantation age (year)")+
  theme(legend.title = element_blank(),legend.key=element_rect(fill="white"),legend.position="top")
```


```{r}
#| label: fig-opcostrel
#| fig-cap: Variability of the opportunity cost
#| fig-width: 4
#| fig-height: 3
tmp2<-filter(results_maxnpv, scenario_case == 0)%>%
  group_by(siteIndex,discountRate,bpr)%>%
  summarise(npv = mean(maxnpv))%>%
  group_by(siteIndex,discountRate)%>%
  arrange(bpr)%>%
  mutate(NPV_ref = first(npv),
         NPV = npv,
         cost = NPV_ref - NPV,
         cost_pc = cost/NPV_ref)%>%
  ungroup()%>%
  mutate(r_label = paste0("r ==",discountRate))


ggplot(data=filter(tmp2,discountRate ==0.02),aes(y=cost_pc, x = bpr*100, col=factor(siteIndex), group=siteIndex))+
  geom_point()+
  stat_smooth(method="lm",se=FALSE)+
  facet_wrap(~r_label,labeller = "label_parsed")+
  theme(legend.position = "top")+
  labs(col = "Site index")+xlab("Bark-stripping rate (%)")+
  scale_color_manual(values = sicolors)+
  ylab("Opportunity cost (%NPV)")+
  geom_abline(intercept = 6000, col="grey")

```





