Commit ce74c732 authored by patrick.lambert's avatar patrick.lambert
Browse files

Merge branch 'exploration_GR3D_process' of...

Merge branch 'exploration_GR3D_process' of gitlab-ssh.irstea.fr:SimAquaLife/GR3D into SpawnerRunAnalysis
parents 65a76167 185c867a
......@@ -124,3 +124,4 @@ org.*
# qgis file
*.qix
.Rproj.user
/desktop.ini
library(dplyr)
library(tidyr)
library(ggplot2)
distance <- as.matrix(read.csv("../../data/input/northeastamerica/distanceGridNEA.csv", row.names = 1, stringsAsFactors = FALSE))
#distance <- as.matrix(read.csv("../../data/input/atlanticarea/distanceGridAA.csv", row.names = 1, stringsAsFactors = FALSE))
distance <- distance %>%
replace(., col(.) == row(.), NA) %>%
as.data.frame() %>% mutate(destination = row.names(.)) %>%
pivot_longer(cols = -destination, names_to = 'departure', values_to = 'distance')
# true values
meanInterDistance <- mean(distance$distance, na.rm = TRUE )
standardDeviationInterDistance <- sd(distance$distance, na.rm = TRUE )
# asin XML from Rougier's application
meanInterDistance = 300
standardDeviationInterDistance = 978
alpha0 = -2.9
alpha1 = 19.7
distance <- distance %>% mutate(logitW = alpha0 - alpha1 * (distance - meanInterDistance) / standardDeviationInterDistance) %>%
mutate(W = 1/(1 + exp(-logitW)))
#filter(departure =='Pearl', destination == 'Escambia')
distance %>% group_by(departure) %>% summarise(sumW = sum(W, na.rm = TRUE)) %>%
#filter(departure == 'Potomac')
ggplot(aes(x = sumW)) + geom_histogram() +
geom_vline(aes(xintercept = mean(sumW)), color="blue", linetype="dashed", size=1)
distance %>% group_by(departure) %>% summarise(sumW = sum(W, na.rm = TRUE)) %>%
summarise(sumW = mean(sumW))
......@@ -378,10 +378,14 @@ plot(temperature,spawnerSurvivalPostReproductionTempRef(temperature,TrefSpring,
ylab = "Spawner survival in river",
main = "Spawner Survival after reproduction")
lines(temperature, spawnerSurvivalPostReproductionTempRef(temperature,
plot(temperature, spawnerSurvivalPostReproductionTempRef(temperature,
TrefSummer,
1,
0.1), col = "#FF0000")
0.19,
0.20), col = "#FF0000")
lines(temperature, spawnerSurvivalPostReproductionTempRef(temperature,
ToptSpring,
1.,
......
library(tydiverse)
#x50 = Tref
#x5 = 25°C, ie the temperature above which there is no spanwer survival after reproduction (based on Limbrug et al, 2003)
logit2 = function(x, x50, x5){
return( 1/ (1+exp((log(19)/(x5-x50))*(x-x50))))
}
logit2 = function(Triver, Tref,minTempForIteroparity){
return( 1/ (1+exp((log(19)/(minTempForIteroparity-Tref))*(temp-Tref))))
}
#-------------------------------------------
data.frame(temperature = seq(0,40,.1)) %>%
mutate(pctSurvieMax = logit2(temperature, x50=20, x5=25),
survie = pctSurvieMax *.2) %>%
ggplot(aes(x=temperature, y = survie )) +
geom_path()
> log(19)/(25-19.9)
parLogit = c(Tref = 19.9, minTempForIteroparity = 25)
optim_logit2_fn <- function(par,temperature, target, logitPar){
survival = logit2(temp = temperature, Tref = logitPar['Tref'],
minTempForIteroparity = logitPar['minTempForIteroparity']) * par[1]
squareDistance = (mean(survival) - target)^2
return(squareDistance)
}
optim_logit2_fn(par = 0.2,temperature = tempInRiverAvg$summer_river_temperature,
target = 0.1, logitPar = parLogit)
optim_logit2 <- optim(c(coeffb = 0.2),
optim_logit2_fn,
temperature = tempInRiverAvg$summer_river_temperature,
target = 0.1,
logitPar = parLogit,
method = "L-BFGS-B")
ID;River;Latitude;% repeaters
1;St. Johns R.;30;0
2;Ogeechee;32;0
3;Edisto;32,7;0
4;Neuse;35;3
5;James;36,8;28
6;York;37,4;24
7;Potomac;37,9;18
8;Susquehanna;39,9;40
9;Delaware;40,2;2
9;Delaware;40,2;17
10;Hudson 1984;41;72
10;Hudson 1996;41;38
11;Connecticut;41,5;51
11;Connecticut;41,5;63
11;Connecticut;41,5;39
12;St. John;45;71
13;Miramichi;47;59
14;St. Lawrence;49,5;85,7
"";"x"
"tempOptGrow";4,54319875071081
"kOptForFemale";0,574633441952413
"kOptForMale";0,610641589581584
"linfVonBert";75,9100534276917
"lFirstMaturityForMale";40,2428684557835
"lFirstMaturityForFemale";44,8366119284253
"sigmaDeltaLVonBert";1,9415506017532
---
title: "Untitled"
author: "Camille POULET"
date: "24/03/2021"
output: word_document
---
```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
rm(list = ls())
```
```{r library}
library(AlgDesign)
```
#D-optimal design
Select the subset of variables that maximized the determinant of the X'X matrix.
See Atkinson and Donev (1992) for details.
We used the D-optimal design when we have a limited budget and cannot run a completely replicated factorial design.
##Factor scaling
This algorithm deals with both quantitative (continuous) and qualitative (discrete) factors. The levels of quantitative factors are scaled so that the minimum value is -1 and the maximum value is 1.
#Experimental design
4 parameters sets (meanDistStray, pHoming, AlleeEffect, spSurvival)
meanDistStray:
- dist_50 pct {3; 19, 53.6,88.3}
- a1 {3; -40.76, -5.17, -2.9}
- a2 {3; 35.76, 21.25, -19.7}
pHoming:
- phom {2; 0.75, 0.95} Ajouter une valeur intermediaire ?
AlleeEffect:
- S95 {3; 2.4, X, X}
spSurvival:
- alpha {2; 0.58, X}
- beta {2; 0.2, X}
- Spsp {2; 0.1,X}
6-8 parameters, with two or three modalities
```{r sensitivity analysis }
#10 replicates
(3*3*2*2)*10
360*3
1080/60
```
```{r design}
#fractional factorial
factorialPlan = gen.factorial(levels = c(3,2,3,3), nVars=4, varNames = c("meanDistForStraying","pHoming","alleeEffect","survivalAfterSpawning"))
desFactorial = optFederov(~.,factorialPlan)
desFactorial$design
#orthogonal design with Federov algorythm
desOrthogonal = optFederov(~.,factorialPlan,nRepeats=20)
desOrthogonal$design
centralPlan = gen.factorial(levels = c(3,2,3,3), center = TRUE, varNames = c("meanDistForStraying","pHoming","alleeEffect","survivalAfterSpawning"))
desCentral = optFederov(~quad (meanDistForStraying,pHoming,alleeEffect,survivalAfterSpawning), centralPlan,evaluateI = TRUE, nRepeats = 100)
#Latin square
latinPlan <- gen.factorial(levels = c(3,2,3,3))
desLatin <- optFederov(~.,dat,nTrials=25,nRepeats=1000)
cs <- xtabs(~.,desL$design)
{xx<-matrix(0,5,5); for (i in 1:5) xx=xx+cs[1:5,1:5,i]*i;xx}
optBlock
2.4*1500
```
\ No newline at end of file
[ViewState]
Mode=
Vid=
FolderType=Generic
This diff is collapsed.
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment