GAM_saison
Agathe Dumont
5/31/2021
data & packages
df = readRDS("ecd_obs_agathe.RDS")
df <- df[!is.na(df$obs_ts),]
df[is.na(df)] <- 0
library('mgcv')## Loading required package: nlme## This is mgcv 1.8-34. For overview type 'help("mgcv-package")'.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, qqplotlibrary('ggpubr')## Registered S3 methods overwritten by 'car':
## method from
## influence.merMod lme4
## cooks.distance.influence.merMod lme4
## dfbeta.influence.merMod lme4
## dfbetas.influence.merMod lme4library("lubridate")##
## Attaching package: 'lubridate'## The following objects are masked from 'package:base':
##
## date, intersect, setdiff, unionExtract 95% of captures
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,]Res_saison pour 100% des données de capture
load(file="~/Desktop/STAGE/resultats_gam_saison.100perc")Res_saison pour 95% des données de capture
load(file="~/Desktop/STAGE/resultats_gam_saison.95perc")espèce à considérer
espece<-c("num_Other_bony_fishes","num_Silky_shark","num_Common_dolphinfish")transformation 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<-list("O1G1","O1G2","O2G1","O2G2")
data_list<-list(O1G1=df1,O1G2=df2,O2G1=df3,O2G2=df4) choisir la distribution
famille<-"tw"enregistrer les modèles
res_saison<-list()
for (i in 1:length(espece)){
for (j in 1:length(data_list)){
data <- data_list[[j]]
form <- formula(paste(espece[i],"s(capture_total_target_tunas_corrigee,k=10)+s(mois,bs='cc',k=12)", sep="~"))
gam <- gam(form, data=data,family = famille)
gr=strat[j]
sp=gsub("^num_","",espece[i])
fn=paste(gr,sp,sep=".")
res_saison[[paste(gr,sp,sep=".")]] = list(sp=sp,gr=gr,fn=fn,gam=gam)
}
}les Pvalues
pvalues<-list()
for (v in res_saison){
a<-summary(v$gam)
b<-c(a$s.table[,'p-value'][[1]],a$s.table[,'p-value'][[2]],a$edf[[1]],a$edf[[2]],a$dev.expl)
pvalues[[v$fn]]<-b
}
colnames<-c("Pvalue s(capture) ","Pvalue s(saison)","edf s(capture)","edf s(saison)","deviance expl. ")
pvalues.table<-t(as.data.frame(pvalues))
colnames(pvalues.table)<-colnames
pvalues.table## Pvalue s(capture) Pvalue s(saison) edf s(capture)
## O1G1.Other_bony_fishes 0 0 8.998574
## O1G2.Other_bony_fishes 0 0 8.998063
## O2G1.Other_bony_fishes 0 0 8.991366
## O2G2.Other_bony_fishes 0 0 8.997648
## O1G1.Silky_shark 0 0 7.661930
## O1G2.Silky_shark 0 0 8.852464
## O2G1.Silky_shark 0 0 7.479381
## O2G2.Silky_shark 0 0 8.436158
## O1G1.Common_dolphinfish 0 0 8.901683
## O1G2.Common_dolphinfish 0 0 7.347386
## O2G1.Common_dolphinfish 0 0 8.913729
## O2G2.Common_dolphinfish 0 0 8.970957
## edf s(saison) deviance expl.
## O1G1.Other_bony_fishes 9.997196 1.00000000
## O1G2.Other_bony_fishes 9.977691 0.31853607
## O2G1.Other_bony_fishes 9.995302 1.00000000
## O2G2.Other_bony_fishes 9.807641 0.54901049
## O1G1.Silky_shark 9.736164 0.24089555
## O1G2.Silky_shark 9.702106 0.44612076
## O2G1.Silky_shark 9.342765 0.03808246
## O2G2.Silky_shark 8.075151 0.62543815
## O1G1.Common_dolphinfish 9.957383 0.08701539
## O1G2.Common_dolphinfish 9.445858 0.70058100
## O2G1.Common_dolphinfish 9.959970 0.04691637
## O2G2.Common_dolphinfish 8.392829 0.96338537Réalisation des diagnostiques
for (v in res_saison){
par(mfrow=c(2,2))
try(gam.check(v$gam))
}
##
## Method: REML Optimizer: outer newton
## full convergence after 20 iterations.
## Gradient range [9.107371e-07,4.180816]
## (score 1576236 & scale 1).
## Hessian positive definite, eigenvalue range [3.995208,9104.877].
## Model rank = 20 / 20
##
## 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 9 0.90 <2e-16 ***
## s(mois) 10 10 0.96 0.34
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Method: REML Optimizer: outer newton
## full convergence after 18 iterations.
## Gradient range [-0.01998189,0.009370835]
## (score 4269.912 & scale 1).
## Hessian positive definite, eigenvalue range [0.003186278,208.0463].
## Model rank = 20 / 20
##
## 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 9.00 0.98 0.41
## s(mois) 10.00 9.98 1.00 0.88
##
## Method: REML Optimizer: outer newton
## full convergence after 19 iterations.
## Gradient range [-9.702013e-08,3.299373]
## (score 704702.4 & scale 1).
## Hessian positive definite, eigenvalue range [1.179766,3678.222].
## Model rank = 20 / 20
##
## 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.94 0.040 *
## s(mois) 10.00 10.00 0.94 0.035 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Method: REML Optimizer: outer newton
## full convergence after 12 iterations.
## Gradient range [-2.465634e-05,0.001379617]
## (score 4454.288 & scale 1).
## Hessian positive definite, eigenvalue range [3.767728,78.97735].
## Model rank = 20 / 20
##
## 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 9.00 1.02 0.935
## s(mois) 10.00 9.81 0.92 0.015 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Method: REML Optimizer: outer newton
## full convergence after 6 iterations.
## Gradient range [-0.0008389602,0.001198045]
## (score 5210.386 & scale 1).
## Hessian positive definite, eigenvalue range [2.176351,1921.354].
## Model rank = 20 / 20
##
## 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 7.66 0.88 0.67
## s(mois) 10.00 9.74 0.86 0.10
##
## Method: REML Optimizer: outer newton
## full convergence after 11 iterations.
## Gradient range [-7.100197e-06,0.004164123]
## (score 2682.333 & scale 1).
## Hessian positive definite, eigenvalue range [1.594967,246.5256].
## Model rank = 20 / 20
##
## 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.85 0.97 0.89
## s(mois) 10.00 9.70 0.96 0.67
##
## Method: REML Optimizer: outer newton
## full convergence after 8 iterations.
## Gradient range [-2.018802e-06,0.008425058]
## (score 17319.55 & scale 1).
## Hessian positive definite, eigenvalue range [1.896182,4819.931].
## Model rank = 20 / 20
##
## 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 7.48 0.88 0.01 **
## s(mois) 10.00 9.34 0.91 0.25
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Method: REML Optimizer: outer newton
## full convergence after 25 iterations.
## Gradient range [-7.258224e-05,0.0008652862]
## (score 238.6662 & scale 1).
## Hessian positive definite, eigenvalue range [0.0001027784,31.99959].
## Model rank = 20 / 20
##
## 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.44 0.55 0.075 .
## s(mois) 10.00 8.08 0.58 0.645
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Method: REML Optimizer: outer newton
## full convergence after 8 iterations.
## Gradient range [-3.379468e-07,0.004690382]
## (score 34556.21 & scale 1).
## Hessian positive definite, eigenvalue range [3.666854,9565.964].
## Model rank = 20 / 20
##
## 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.90 0.95 0.475
## s(mois) 10.00 9.96 0.92 0.025 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Method: REML Optimizer: outer newton
## full convergence after 24 iterations.
## Gradient range [-6.053703e-05,0.0005999936]
## (score 597.648 & scale 1).
## Hessian positive definite, eigenvalue range [0.5814478,57.57166].
## Model rank = 20 / 20
##
## 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 7.35 0.82 0.45
## s(mois) 10.00 9.45 0.82 0.16
##
## Method: REML Optimizer: outer newton
## full convergence after 28 iterations.
## Gradient range [3.579359e-13,0.002986933]
## (score 61808.55 & scale 1).
## Hessian positive definite, eigenvalue range [0.002987603,3216.578].
## Model rank = 20 / 20
##
## 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.91 0.93 0.03 *
## s(mois) 10.00 9.96 0.97 0.84
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Method: REML Optimizer: outer newton
## full convergence after 12 iterations.
## Gradient range [-1.729302e-05,0.0002038824]
## (score 189.3396 & scale 1).
## Hessian positive definite, eigenvalue range [0.0001078746,19.99989].
## Model rank = 20 / 20
##
## 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.97 0.58 0.22
## s(mois) 10.00 8.39 0.45 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1gam.check(res_saison\(O1G2.Common_dolphinfish\)gam) Error in while (sum(ind) > 0 && kk < 50) { : valeur manquante là où TRUE / FALSE est requis
plot all the effect
for (v in res_saison){
o<-getViz(v$gam)
print(plot(o, allTerms = T), pages = 1)
}
