data & packages
# Définir environement de travail
df = readRDS("ecd_obs_agathe.RDS")
# Charger les données
df <- df[!is.na(df$obs_ts),]
# Tableau sans les "vrais" NA, ie les calais sans données d'observateurs
df[is.na(df)] <- 0
# avec les 0 ( changer les " faux NA " par des 0 )
library('mgcv')
## Loading required package: nlme
## This is mgcv 1.8-34. For overview type 'help("mgcv-package")'.
# pour le modele gam
library('mgcViz')
## Loading required package: qgam
## Loading required package: ggplot2
## Loading required package: rgl
## This build of rgl does not include OpenGL functions. Use
## rglwidget() to display results, e.g. via options(rgl.printRglwidget = TRUE).
## Registered S3 method overwritten by 'GGally':
## method from
## +.gg ggplot2
## Registered S3 method overwritten by 'mgcViz':
## method from
## +.gg GGally
##
## Attaching package: 'mgcViz'
## The following objects are masked from 'package:stats':
##
## qqline, qqnorm, qqplot
library('ggpubr')
## Registered S3 methods overwritten by 'car':
## method from
## influence.merMod lme4
## cooks.distance.influence.merMod lme4
## dfbeta.influence.merMod lme4
## dfbetas.influence.merMod lme4
library("maps") # pour la carte sur le graph 2D
library("lubridate")
##
## Attaching package: 'lubridate'
## The following objects are masked from 'package:base':
##
## date, intersect, setdiff, union
library(dplyr)
##
## Attaching package: 'dplyr'
## The following object is masked from 'package:nlme':
##
## collapse
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
library(cowplot)
##
## Attaching package: 'cowplot'
## The following object is masked from 'package:lubridate':
##
## stamp
## The following object is masked from 'package:ggpubr':
##
## get_legend
Extract 95% OU 100% of captures and load file
cc = sort(df$capture_total_target_tunas_corrigee)
xx = cumsum(cc)
I = min(which(xx > 0.95*max(xx)))
cc_max = cc[I]
df = df[df$capture_total_target_tunas_corrigee < cc_max,]
Species to consider + Stratification+ Legende+ Distribution
# espèce à considérer -----------------------------------------------------
espece<-c("num_Other_bony_fishes","num_Silky_shark","num_Common_dolphinfish")
# transformation & création des variables -------------
df$num_Other_bony_fishes=df$`num_Other bony fishes`
df$num_Common_dolphinfish=df$`num_Common dolphinfish`
df$num_Silky_shark=df$`num_Silky shark`
df$mois<-month(df$d_act)
# stratification ---------------
df$group = paste0("O",df$ocean,"G",df$code_assoc_groupe)
# Création d'une variable groupe
# O1G1, O1G2, O2G1, O2G2
# O1: Atlantique 02: Indien
# G1: Floating object G2: Free swiming school
## Par océan et banc
df1<-df[df$group =="O1G1",]
# df1 contient tous les calais qui sont sont O1G1
df2<-df[df$group =="O1G2",]
# df2 contient tous les calais qui sont O1G2
df3<-df[df$group =="O2G1",]
# df3 contient tous les calais qui sont O2G1
df4<-df[df$group =="O2G2",]
# df4 contient tous les calais qui sont O2G2
# légende -----------------------------------------------
strat<-c("O1G1","O1G2","O2G1","O2G2")
data_list<-list(O1G1=df1,O1G2=df2,O2G1=df3,O2G2=df4)
load models by=trimestre with the cluster
Fonction to load the results from the cluster
load_rdata_to_list <- function(file)
{
env <- new.env()
nm <- load(file, env)
as.list(env)[nm]
}
load the results from the cluster
fns = dir("gam_all_factors_tweedie_bigger_k/95percent","RData$")
xx = do.call(rbind,strsplit(fns,".",fixed=TRUE))
grps=xx[,1]
sps=gsub("^num_","",xx[,2])
res = list()
for (k in seq_along(fns)) {
fn=fns[k]
x = load_rdata_to_list(file.path("gam_all_factors_tweedie_bigger_k","95percent",fn))
gr = grps[k]
sp = sps[k]
res[[paste(gr,sp,sep=".")]] = list(sp=sp,gr=gr,fn=fn,gam=x$gam)
}
saveRDS(res,"resultats_gam_by_trimestre_tweedie.95perc.RDS")
res = readRDS("resultats_gam_by_trimestre.95perc.RDS")
res = readRDS("resultats_gam_by_trimestre.100perc.RDS")
res = readRDS("resultats_gam_by_trimestre_tweedie.95perc.RDS")
Pvalues
pvalues<-list()
for (v in res){
a<-summary(v$gam)
b<-c(a$s.table[,'p-value'][[1]],a$s.table[,'p-value'][[2]],a$s.table[,'p-value'][[3]],a$s.table[,'p-value'][[4]],a$s.table[,'p-value'][[5]],a$edf[[1]],a$edf[[2]],a$edf[[3]],a$edf[[4]],a$edf[[5]],a$dev.expl)
pvalues[[v$fn]]<-b
}
colnames<-c("Pvalue s(capture) ","Pvalue trim 1","Pvalue trim 2","Pvalue trim 3","Pvalue trim 4","edf s(capture)","edf trim 1","edf trim 2","edf trim 3","edf trim 4","deviance expl. ")
pvalues.table<-t(as.data.frame(pvalues))
colnames(pvalues.table)<-colnames
pvalues.table
## Pvalue s(capture) Pvalue trim 1
## O1G1.num_Common_dolphinfish.RData 0.000000e+00 0.000000e+00
## O1G1.num_Silky_shark.RData 0.000000e+00 0.000000e+00
## O1G2.num_Common_dolphinfish.RData 3.607580e-06 3.779866e-04
## O1G2.num_Other_bony_fishes.RData 0.000000e+00 0.000000e+00
## O1G2.num_Silky_shark.RData 0.000000e+00 1.562330e-02
## O2G1.num_Common_dolphinfish.RData 0.000000e+00 0.000000e+00
## O2G1.num_Other_bony_fishes.RData 0.000000e+00 0.000000e+00
## O2G1.num_Silky_shark.RData 0.000000e+00 0.000000e+00
## O2G2.num_Common_dolphinfish.RData 0.000000e+00 0.000000e+00
## O2G2.num_Other_bony_fishes.RData 0.000000e+00 0.000000e+00
## O2G2.num_Silky_shark.RData 1.748851e-04 4.429159e-05
## Pvalue trim 2 Pvalue trim 3 Pvalue trim 4
## O1G1.num_Common_dolphinfish.RData 0.000000e+00 0.0000000 0.0000000000
## O1G1.num_Silky_shark.RData 0.000000e+00 0.0000000 0.0000000000
## O1G2.num_Common_dolphinfish.RData 0.000000e+00 0.0000000 0.0000000000
## O1G2.num_Other_bony_fishes.RData 2.345917e-150 0.0000000 0.0000000000
## O1G2.num_Silky_shark.RData 0.000000e+00 0.0000000 0.0000000000
## O2G1.num_Common_dolphinfish.RData 0.000000e+00 0.0000000 0.0000000000
## O2G1.num_Other_bony_fishes.RData 0.000000e+00 0.0000000 0.0000000000
## O2G1.num_Silky_shark.RData 0.000000e+00 0.0000000 0.0000000000
## O2G2.num_Common_dolphinfish.RData 0.000000e+00 0.0000000 0.0000000000
## O2G2.num_Other_bony_fishes.RData 0.000000e+00 0.0000000 0.0000000000
## O2G2.num_Silky_shark.RData 3.880141e-02 0.3154147 0.0008148199
## edf s(capture) edf trim 1 edf trim 2
## O1G1.num_Common_dolphinfish.RData 8.766676 88.63795 78.872862
## O1G1.num_Silky_shark.RData 5.867743 20.55815 26.888249
## O1G2.num_Common_dolphinfish.RData 1.017513 16.01513 15.256713
## O1G2.num_Other_bony_fishes.RData 9.000514 67.17326 52.231305
## O1G2.num_Silky_shark.RData 8.778859 11.05967 23.283866
## O2G1.num_Common_dolphinfish.RData 8.898762 132.67701 125.526028
## O2G1.num_Other_bony_fishes.RData 8.988825 154.70276 152.028189
## O2G1.num_Silky_shark.RData 4.314473 96.46966 79.021452
## O2G2.num_Common_dolphinfish.RData 9.000000 -357.62743 3.133591
## O2G2.num_Other_bony_fishes.RData 8.995224 26.50747 35.684851
## O2G2.num_Silky_shark.RData 1.075783 22.72786 12.227037
## edf trim 3 edf trim 4 deviance expl.
## O1G1.num_Common_dolphinfish.RData 71.607727 121.185805 0.3705957
## O1G1.num_Silky_shark.RData 28.423969 24.428163 0.3031578
## O1G2.num_Common_dolphinfish.RData 27.336749 3.121850 0.9073282
## O1G2.num_Other_bony_fishes.RData 150.882566 47.604903 0.8512399
## O1G2.num_Silky_shark.RData 30.864023 10.533253 0.7052324
## O2G1.num_Common_dolphinfish.RData 126.836178 140.647688 0.3609168
## O2G1.num_Other_bony_fishes.RData 159.343662 160.103911 1.0000000
## O2G1.num_Silky_shark.RData 90.345499 114.574210 0.2318903
## O2G2.num_Common_dolphinfish.RData 497.542975 3.011305 0.9794741
## O2G2.num_Other_bony_fishes.RData 14.402668 37.095994 0.9887014
## O2G2.num_Silky_shark.RData 7.087074 14.724926 0.8204418
write.csv(pvalues.table,file = "tableau des pvalues.csv")
Diagnostic
for (v in res){
print(v$fn)
par(mfrow=c(2,2))
try(gam.check(v$gam))
}
## [1] "O1G1.num_Common_dolphinfish.RData"
##
## Method: REML Optimizer: outer newton
## full convergence after 8 iterations.
## Gradient range [-1.806193e-05,0.004081038]
## (score 26949.67 & scale 1).
## Hessian positive definite, eigenvalue range [3.370876,3379.998].
## Model rank = 682 / 682
##
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
##
## k' edf k-index p-value
## s(capture_total_target_tunas_corrigee) 9.00 8.77 0.95 0.68
## te(obs_lon,obs_lat):as.factor(trimestre)1 168.00 88.64 0.98 0.97
## te(obs_lon,obs_lat):as.factor(trimestre)2 168.00 78.87 0.99 1.00
## te(obs_lon,obs_lat):as.factor(trimestre)3 168.00 71.61 0.98 0.96
## te(obs_lon,obs_lat):as.factor(trimestre)4 168.00 121.19 0.99 1.00
## [1] "O1G1.num_Silky_shark.RData"
##
## Method: REML Optimizer: outer newton
## full convergence after 7 iterations.
## Gradient range [-0.0007654101,0.001118902]
## (score 4692.888 & scale 1).
## Hessian positive definite, eigenvalue range [0.4261497,1791.87].
## Model rank = 682 / 682
##
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
##
## k' edf k-index p-value
## s(capture_total_target_tunas_corrigee) 9.00 5.87 0.86 0.890
## te(obs_lon,obs_lat):as.factor(trimestre)1 168.00 20.56 0.80 <2e-16 ***
## te(obs_lon,obs_lat):as.factor(trimestre)2 168.00 26.89 0.80 0.005 **
## te(obs_lon,obs_lat):as.factor(trimestre)3 168.00 28.42 0.81 0.035 *
## te(obs_lon,obs_lat):as.factor(trimestre)4 168.00 24.43 0.81 0.005 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "O1G2.num_Common_dolphinfish.RData"
##
## Method: REML Optimizer: outer newton
## full convergence after 13 iterations.
## Gradient range [-0.0001175687,0.0003380411]
## (score 367.4419 & scale 1).
## Hessian positive definite, eigenvalue range [7.062788e-06,87.34817].
## Model rank = 682 / 682
##
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
##
## k' edf k-index p-value
## s(capture_total_target_tunas_corrigee) 9.00 1.02 0.31 0.190
## te(obs_lon,obs_lat):as.factor(trimestre)1 168.00 16.02 0.28 0.010 **
## te(obs_lon,obs_lat):as.factor(trimestre)2 168.00 15.26 0.27 <2e-16 ***
## te(obs_lon,obs_lat):as.factor(trimestre)3 168.00 27.34 0.28 0.005 **
## te(obs_lon,obs_lat):as.factor(trimestre)4 168.00 3.12 0.27 0.005 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "O1G2.num_Other_bony_fishes.RData"
## Error in while (sum(ind) > 0 && kk < 50) { :
## valeur manquante là où TRUE / FALSE est requis
## [1] "O1G2.num_Silky_shark.RData"
##
## Method: REML Optimizer: outer newton
## full convergence after 15 iterations.
## Gradient range [-6.235185e-05,1.535316e-05]
## (score 1918.173 & scale 1).
## Hessian positive definite, eigenvalue range [0.5272671,340.1929].
## Model rank = 682 / 682
##
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
##
## k' edf k-index p-value
## s(capture_total_target_tunas_corrigee) 9.00 8.78 0.96 0.280
## te(obs_lon,obs_lat):as.factor(trimestre)1 168.00 11.06 0.90 0.015 *
## te(obs_lon,obs_lat):as.factor(trimestre)2 168.00 23.28 0.86 <2e-16 ***
## te(obs_lon,obs_lat):as.factor(trimestre)3 168.00 30.86 0.87 <2e-16 ***
## te(obs_lon,obs_lat):as.factor(trimestre)4 168.00 10.53 0.85 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "O2G1.num_Common_dolphinfish.RData"
##
## Method: REML Optimizer: outer newton
## full convergence after 12 iterations.
## Gradient range [-8.211114e-05,0.06204618]
## (score 46267.65 & scale 1).
## Hessian positive definite, eigenvalue range [3.844242,3593.743].
## Model rank = 682 / 682
##
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
##
## k' edf k-index p-value
## s(capture_total_target_tunas_corrigee) 9.0 8.9 0.96 0.32
## te(obs_lon,obs_lat):as.factor(trimestre)1 168.0 132.7 0.97 0.69
## te(obs_lon,obs_lat):as.factor(trimestre)2 168.0 125.5 0.96 0.58
## te(obs_lon,obs_lat):as.factor(trimestre)3 168.0 126.8 0.97 0.59
## te(obs_lon,obs_lat):as.factor(trimestre)4 168.0 140.7 0.97 0.66
## [1] "O2G1.num_Other_bony_fishes.RData"
##
## Method: REML Optimizer: outer newton
## full convergence after 16 iterations.
## Gradient range [-0.0008004842,0.009866675]
## (score 531568.7 & scale 1).
## Hessian positive definite, eigenvalue range [3.968763,3619.525].
## Model rank = 682 / 682
##
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
##
## k' edf k-index p-value
## s(capture_total_target_tunas_corrigee) 9.00 8.99 0.96 0.035 *
## te(obs_lon,obs_lat):as.factor(trimestre)1 168.00 154.70 0.96 0.185
## te(obs_lon,obs_lat):as.factor(trimestre)2 168.00 152.03 0.97 0.495
## te(obs_lon,obs_lat):as.factor(trimestre)3 168.00 159.34 0.97 0.275
## te(obs_lon,obs_lat):as.factor(trimestre)4 168.00 160.10 0.98 0.545
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "O2G1.num_Silky_shark.RData"
##
## Method: REML Optimizer: outer newton
## full convergence after 7 iterations.
## Gradient range [-0.01249161,0.02659204]
## (score 15913.31 & scale 1).
## Hessian positive definite, eigenvalue range [0.931433,3955.474].
## Model rank = 682 / 682
##
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
##
## k' edf k-index p-value
## s(capture_total_target_tunas_corrigee) 9.00 4.31 0.87 <2e-16 ***
## te(obs_lon,obs_lat):as.factor(trimestre)1 168.00 96.47 0.90 0.27
## te(obs_lon,obs_lat):as.factor(trimestre)2 168.00 79.02 0.90 0.42
## te(obs_lon,obs_lat):as.factor(trimestre)3 168.00 90.35 0.90 0.34
## te(obs_lon,obs_lat):as.factor(trimestre)4 168.00 114.57 0.90 0.35
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "O2G2.num_Common_dolphinfish.RData"
## Error in while (sum(ind) > 0 && kk < 50) { :
## valeur manquante là où TRUE / FALSE est requis
## [1] "O2G2.num_Other_bony_fishes.RData"
##
## Method: REML Optimizer: outer newton
## full convergence after 12 iterations.
## Gradient range [-9.877641e-05,0.0001059462]
## (score 528.0641 & scale 1).
## Hessian positive definite, eigenvalue range [0.05304179,100.2689].
## Model rank = 682 / 682
##
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
##
## k' edf k-index p-value
## s(capture_total_target_tunas_corrigee) 9.0 9.0 0.45 0.39
## te(obs_lon,obs_lat):as.factor(trimestre)1 168.0 26.5 0.36 <2e-16 ***
## te(obs_lon,obs_lat):as.factor(trimestre)2 168.0 35.7 0.35 <2e-16 ***
## te(obs_lon,obs_lat):as.factor(trimestre)3 168.0 14.4 0.36 <2e-16 ***
## te(obs_lon,obs_lat):as.factor(trimestre)4 168.0 37.1 0.38 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "O2G2.num_Silky_shark.RData"
##
## Method: REML Optimizer: outer newton
## step failed after 10 iterations.
## Gradient range [-0.3465182,0.29709]
## (score 156.5845 & scale 1).
## Hessian positive definite, eigenvalue range [0.01402357,23.21733].
## Model rank = 682 / 682
##
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
##
## k' edf k-index p-value
## s(capture_total_target_tunas_corrigee) 9.00 1.08 0.59 0.03 *
## te(obs_lon,obs_lat):as.factor(trimestre)1 168.00 22.73 0.36 <2e-16 ***
## te(obs_lon,obs_lat):as.factor(trimestre)2 168.00 12.23 0.36 <2e-16 ***
## te(obs_lon,obs_lat):as.factor(trimestre)3 168.00 7.09 0.37 <2e-16 ***
## te(obs_lon,obs_lat):as.factor(trimestre)4 168.00 14.72 0.36 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot gam by = saison
for (v in res){
o=getViz(v$gam)
print(paste(v$sp,v$gr))
trim1<-plot(sm(o,2))+xlab("Longitude") + ylab("Latitude") + ggtitle("trimestre 1") + l_fitRaster() + l_fitContour() +annotation_map(map_data("world"))
trim2<-plot(sm(o,3))+xlab("Longitude") + ylab("Latitude") + ggtitle("trimestre 2") + l_fitRaster() + l_fitContour() +annotation_map(map_data("world"))
trim3<-plot(sm(o,4))+xlab("Longitude") + ylab("Latitude") + ggtitle("trimestre 3") + l_fitRaster() + l_fitContour() +annotation_map(map_data("world"))
trim4<-plot(sm(o,5))+xlab("Longitude") + ylab("Latitude") + ggtitle("trimestre 4") + l_fitRaster() + l_fitContour() +annotation_map(map_data("world"))
gridPrint(trim1,trim2,trim3,trim4,ncol=2)
}
## [1] "Common_dolphinfish O1G1"
## [1] "Silky_shark O1G1"
## [1] "Common_dolphinfish O1G2"
## [1] "Other_bony_fishes O1G2"
## [1] "Silky_shark O1G2"
## [1] "Common_dolphinfish O2G1"
## [1] "Other_bony_fishes O2G1"
## [1] "Silky_shark O2G1"
## [1] "Common_dolphinfish O2G2"
## [1] "Other_bony_fishes O2G2"
## [1] "Silky_shark O2G2"
for (v in res){
o=getViz(v$gam)
print(plot(sm(o,1))+xlab("capture_total_target_tunas_corrigee") + ylab("effect capture")+ggtitle(v$fn))
}
