Exploration of dispersal

exploration/scriptR/dispersionCalibration.Rmd

+---
+date: "`r Sys.Date()`"
+author: "P. Lambert, C. Poulet "
+title: "Dispersal calibration for GR3D_NEA"
+output: 
+  officedown::rdocx_document:
+    mapstyles:
+      Normal: ['First Paragraph']
+---
+
+```{r setup, include=FALSE}
+knitr::opts_chunk$set(echo = FALSE, fig.cap = TRUE)
+library(officedown)
+library(officer)
+
+library(ggplot2)
+library(dplyr)
+library(scales)
+library(flextable)
+
+fp <- fp_par(
+  text.align = "center", 
+  padding.bottom = 20, padding.top = 120, 
+  border.bottom = fp_border())
+
+ft <- fp_text(shading.color='#EFEFEF', bold = TRUE)
+```
+
+
+# The kernel function
+
+It corresponds in Rougier et al 2012 to the basin weight between a departure $j_1$ to destination basin $j_2$:
+$$w_{j_1\rightarrow j_2} = \frac {1} {1 +  e ^{\alpha_0 + \alpha_1 \cdot {\frac {( D_{j_1\rightarrow j_2}   - \mu_D)} {\sigma_D} } } }$$
+
+where $D_{j_1\rightarrow j_2}$ is the distance between the departure and destination basins, $\alpha_0$ and $\alpha_1$ are the kernel parameters, $\mu_D$ and $\sigma_D$ are the mean and standard deviation between inter basin distances.
+
+The objective is to find the kernel parameters which correspond to knowledge expert, i.e. 50 % of strayers settle before 19 km, 95% before 111 km.
+Mathematically,  that leads to find $\alpha_0$ and $\alpha_1$ that solve the 2 following equations:
+
+$$0.5 = \frac {\int_0^{19} \frac {1} {1 +  e ^{\alpha_0 + \alpha_1 \cdot {\frac {( x - \mu_D)} {\sigma_D} } } }~dx} {\int_0^{+\infty} \frac {1} {1 +  e ^{\alpha_0 + \alpha_1 \cdot {\frac {( x - \mu_D)} {\sigma_D} } } }~dx}$$
+$$0.95 = \frac {\int_0^{111} \frac {1} {1 +  e ^{\alpha_0 + \alpha_1 \cdot {\frac {( x - \mu_D)} {\sigma_D} } } }~dx} {\int_0^{+\infty} \frac {1} {1 +  e ^{\alpha_0 + \alpha_1 \cdot {\frac {( x - \mu_D)} {\sigma_D} } } }~dx}$$
+
+```{r data_function}
+rm(list = ls())
+
+source("../GR3D_Rdescription/GR3Dfunction.R")
+
+
+# calculate the percentage of settlement before a distance according to a kernel function
+pctSettlementBefore = function(dist, alpha0, alpha1, meanInterDistance, standardDeviationInterDistance){  
+  integralInf <- integrate(logitKernel, lower = 0, upper = Inf,
+                           alpha0 = alpha0, 
+                           alpha1 = alpha1,
+                           meanInterDistance = meanInterDistance,
+                           standardDeviationInterDistance = standardDeviationInterDistance)$value
+  
+  integralX <- sapply(dist, function(x) {integrate(logitKernel, lower = 0, upper = x,
+                                                   alpha0 = alpha0, 
+                                                   alpha1 = alpha1,
+                                                   meanInterDistance = meanInterDistance,
+                                                   standardDeviationInterDistance = standardDeviationInterDistance )$value})
+  
+  return(integralX / integralInf)
+}
+
+
+# targets
+dataTarget = data.frame(dist = c(19, 111), pct = c(0.50,.95))
+
+```
+
+
+# Reference distances with Rougier et al 2015 parametrisation
+
+```{r}
+# values used in Rougier et al 2015
+alpha0 = -2.9
+alpha1 =  19.7
+meanInterDistance = 300
+standardDeviationInterDistance = 978
+
+# reference points in Rougier et al. 2015
+RougierReferenceDistance = sapply(dataTarget$pct, function(target) uniroot(function(x, alpha0, alpha1, meanInterDistance, standardDeviationInterDistance, target) (pctSettlementBefore(x, alpha0, alpha1, meanInterDistance, standardDeviationInterDistance) - target),
+        interval = c(0,500),
+         alpha0 = alpha0, alpha1 = alpha1,
+         meanInterDistance = meanInterDistance,
+         standardDeviationInterDistance = standardDeviationInterDistance,
+        target = target)$root)
+
+```
+
+With  parameter values defined in Rougier et al 2015, `r scales::percent(dataTarget$pct[1])` of strayers settle before `r RougierReferenceDistance[1]` km, and respectively `r scales::percent(dataTarget$pct[2])` before `r RougierReferenceDistance[1]` km.
+
+# Parameters corresponding to the reference distance defined by US experts
+
+```{r optimisation}
+
+# sum of square error for percentage of settlement
+objFn = function(par, dataTarget, meanInterDistance, standardDeviationInterDistance){
+  
+  SSE = sum((pctSettlementBefore(dist =  dataTarget[,1], 
+                                 alpha0 = par[1], 
+                                 alpha1 = par[2],
+                                 meanInterDistance = meanInterDistance,
+                                 standardDeviationInterDistance = standardDeviationInterDistance) -
+               dataTarget[,2])^2)
+  
+  return(SSE)
+}
+
+
+resDispersal_USExpert = optim(objFn, 
+                     par = c(alpha0 = alpha0, alpha1 = alpha1), 
+                     dataTarget = dataTarget, 
+                     meanInterDistance = meanInterDistance,
+                     standardDeviationInterDistance = standardDeviationInterDistance)
+
+
+# data for verification
+dataVerif <- data.frame(dist = 0:500) %>%  mutate(kernel = logitKernel(dist,alpha0 = resDispersal_USExpert$par[1], 
+                                                                       alpha1 = resDispersal_USExpert$par[2],
+                                                                       meanInterDistance = meanInterDistance,
+                                                                       standardDeviationInterDistance = standardDeviationInterDistance ),
+                                                  pctSettlement = pctSettlementBefore(dist, 
+                                                                                      alpha0 = resDispersal_USExpert$par[1], 
+                                                                                      alpha1 = resDispersal_USExpert$par[2],
+                                                                                      meanInterDistance = meanInterDistance,
+                                                                                      standardDeviationInterDistance = standardDeviationInterDistance)) 
+```
+
+
+To have `r dataTarget$pct[1]` and `r dataTarget$pct[2]` of strayer settlement before `r dataTarget$dist[1]` and `r dataTarget$dist[2]`,  kernel parameters $\alpha_0$ and $\alpha_1$ take values `r resDispersal_USExpert$par[1]` and `r resDispersal_USExpert$par[2]`.
+
+With this parametrisation, the curve of settlement before a distance becomes
+
+```{r, include=TRUE}
+dataVerif %>% ggplot(aes( x = dist)) + 
+  geom_path(aes( y = pctSettlement)) +  geom_hline(yintercept = dataTarget$pct, lty = 2) +
+  labs(x = 'distance (km)', y = 'percentage of settlement before a distance')
+```
+
+and the kernel function
+
+```{r, include=TRUE}
+dataVerif %>% ggplot(aes( x = dist)) + 
+  geom_path(aes( y = kernel)) +
+  labs(x = 'distance (km)', y = 'kernel')
+```
+
+
+# Value for sensitivity analysis
+
+ 
+
+```{r}
+ASdata <- data.frame(source = c("US Expert A. sapidissima", 'intermediate', "Rougier et al 2015 A. alosa"),
+           dist_50pct = c(dataTarget$dist[1],  round((RougierReferenceDistance[1] + dataTarget$dist[1])/2 ,1), round(RougierReferenceDistance[1],1)),
+           dist_95pct = c(dataTarget$dist[2],  round((RougierReferenceDistance[2] + dataTarget$dist[2])/2 ,1), round(RougierReferenceDistance[2],1)))
+
+dataTarget_intermediate = data.frame(dist = ASdata %>% filter(source == 'intermediate') %>% select(-source) %>% unlist(use.names = FALSE), 
+                                     pct = dataTarget$pct)
+
+resDispersal_intermediate = optim(objFn, 
+                     par = c(alpha0 = alpha0, alpha1 = alpha1), 
+                     dataTarget = dataTarget_intermediate, 
+                     meanInterDistance = meanInterDistance,
+                     standardDeviationInterDistance = standardDeviationInterDistance)
+
+
+ASdata <- ASdata %>% bind_cols(data.frame(alpha0 = c(resDispersal_USExpert$par[1], resDispersal_intermediate$par[1], alpha0),
+                                alpha1 = c(resDispersal_USExpert$par[2], resDispersal_intermediate$par[2], alpha1)))
+
+
+```
+
+```{r, include=TRUE}
+ASdata %>% 
+  flextable() %>% 
+  autofit()
+```
+
+
+```{r}
+# VERSION ANALYTIQUE QUI NE MARCHE PAS !!!!!
+
+#logitKernelIntegral = function(distance, alpha0, alpha1, meanInterDistance, standardDeviationInterDistance){
+#   
+#   primitive = function(x, alpha0, alpha1, meanInterDistance, standardDeviationInterDistance) {
+#     (-standardDeviationInterDistance * 
+#        log(exp((alpha1 * x / standardDeviationInterDistance) - (alpha1 * meanInterDistance / standardDeviationInterDistance) + alpha0) + 1) / alpha1) +
+#       x +
+#       alpha0 * standardDeviationInterDistance / alpha1 - 
+#       meanInterDistance
+#   }
+#   
+#   return(primitive(distance, alpha0, alpha1, meanInterDistance, standardDeviationInterDistance) - primitive(0, alpha0, alpha1, meanInterDistance, standardDeviationInterDistance)) 
+# } 
+
+# tibble(distance = seq(0,500, 1)) %>% 
+#   mutate(W = logitKernel(distance, 
+#                          alpha0 = alpha0, 
+#                          alpha1 = alpha1,
+#                          meanInterDistance = meanInterDistance,
+#                          standardDeviationInterDistance = standardDeviationInterDistance)) %>% 
+#   mutate(cumW = cumsum(W)) %>% 
+#   mutate(cumW2 = logitKernelIntegral(distance, 
+#                                      alpha0 = alpha0, 
+#                                      alpha1 = alpha1,
+#                                      meanInterDistance = meanInterDistance,
+#                                      standardDeviationInterDistance = standardDeviationInterDistance))
+# 
+# ggplot(aes(x = distance, y = cumW)) + geom_line()
+
+```
+
+
+