---

title: “Gamsaison_latlong” author: “Agathe Dumont” date: “5/24/2021” output: html_document

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) 

# choisir la distribution -----------------------------------------

famille<-"tw"

load models se(trimestre) without the cluster

It’s long but it’s reasonable

#eval=!file.exists("resultats_gam_se_trimestre.95perc")

res2<-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)+te(obs_lon,obs_lat,k=13)", sep="~")) 
    gam <- gam(form, data=data,family = famille) 
    gr<-strat[i]
    sp=gsub("^num_","",espece[i])
    fn=paste(gr,sp,sep=".")
    res2[[paste(gr,sp,sep=".")]] = list(sp=sp,gr=gr,fn=fn,gam=gam)
  }
}

#save(res2,file="~/Desktop/STAGE/resultats_gam_se_trimestre.95perc")

load("resultats_gam_se_trimestre.95perc")

Pvalues

pvalues2<-list()
for (p in res2){
    a<-summary(p$gam)
    b<-c(a$s.table[,'p-value'][[1]],a$s.table[,'p-value'][[2]],a$s.table[,'p-value'][[3]],a$edf[[1]],a$edf[[2]],a$edf[[3]],a$dev.expl)
    pvalues2[[p$fn]]<-b
  }

colnames<-c("Pvalue s(capture) ","Pvalue s(mois)","Pvalue s(spatial)","edf s(capture)","edf s(mois)","edf s(spatial)","deviance expl. ")
pvalues2.table<-t(as.data.frame(pvalues2))
colnames(pvalues2.table)<-colnames
pvalues2.table

##                         Pvalue s(capture)  Pvalue s(mois) Pvalue s(spatial)
## O1G1.Other_bony_fishes        0.000000e+00     0.00000000      0.000000e+00
## O1G2.Other_bony_fishes        0.000000e+00     0.00000000      0.000000e+00
## O2G1.Other_bony_fishes        0.000000e+00     0.00000000      0.000000e+00
## O2G2.Other_bony_fishes        0.000000e+00     0.00000000      0.000000e+00
## O1G1.Silky_shark              0.000000e+00     0.00000000      0.000000e+00
## O1G2.Silky_shark              0.000000e+00     0.00000000      0.000000e+00
## O2G1.Silky_shark              0.000000e+00     0.00000000      0.000000e+00
## O2G2.Silky_shark             -2.140288e-06     0.09012353      5.366183e-07
## O1G1.Common_dolphinfish       0.000000e+00     0.00000000      0.000000e+00
## O1G2.Common_dolphinfish       5.454358e-03     0.00000000      0.000000e+00
## O2G1.Common_dolphinfish       0.000000e+00     0.00000000      0.000000e+00
## O2G2.Common_dolphinfish       0.000000e+00     1.00000639      0.000000e+00
##                         edf s(capture)   edf s(mois) edf s(spatial)
## O1G1.Other_bony_fishes    8.998676e+00  9.997019e+00   1.528212e+02
## O1G2.Other_bony_fishes   -3.522636e+08 -2.500313e+11  -2.022463e+11
## O2G1.Other_bony_fishes    8.990515e+00  9.995958e+00   1.593286e+02
## O2G2.Other_bony_fishes    8.998822e+00  9.906167e+00   9.060048e+01
## O1G1.Silky_shark          8.189983e+00  9.008191e+00   6.440124e+01
## O1G2.Silky_shark          8.791842e+00  7.935461e+00   3.693995e+01
## O2G1.Silky_shark          7.611918e+00  9.213945e+00   1.183819e+02
## O2G2.Silky_shark          7.519901e+00  1.920001e+00   2.688388e+01
## O1G1.Common_dolphinfish   8.902342e+00  9.951751e+00   1.276302e+02
## O1G2.Common_dolphinfish   1.024408e+00  9.391334e+00   4.140644e+01
## O2G1.Common_dolphinfish   8.888255e+00  9.947049e+00   1.468313e+02
## O2G2.Common_dolphinfish   8.604548e+00  7.183159e-04   1.834111e+01
##                         deviance expl. 
## O1G1.Other_bony_fishes        1.0000000
## O1G2.Other_bony_fishes        0.8084719
## O2G1.Other_bony_fishes        1.0000000
## O2G2.Other_bony_fishes        0.9890326
## O1G1.Silky_shark              0.3766025
## O1G2.Silky_shark              0.6530729
## O2G1.Silky_shark              0.1244849
## O2G2.Silky_shark              0.7412516
## O1G1.Common_dolphinfish       0.2200557
## O1G2.Common_dolphinfish       0.9145357
## O2G1.Common_dolphinfish       0.1950071
## O2G2.Common_dolphinfish       0.9782163

Diagnostics

for (p in res2){ 
    print(p$fn)
    par(mfrow=c(2,2))
    try(gam.check(p$gam))
  }

## [1] "O1G1.Other_bony_fishes"

## 
## Method: REML   Optimizer: outer newton
## full convergence after 17 iterations.
## Gradient range [-0.001727455,0.02993617]
## (score 1280320 & scale 1).
## Hessian positive definite, eigenvalue range [3.995457,7260.989].
## Model rank =  188 / 188 
## 
## 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.92  <2e-16 ***
## s(mois)                                 10  10    0.96    0.38    
## te(obs_lon,obs_lat)                    168 153    0.96    0.24    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "O1G2.Other_bony_fishes"

## 
## Method: REML   Optimizer: outer newton
## full convergence after 12 iterations.
## Gradient range [-2.519514e-06,1.216617e-06]
## (score 2300.21 & scale 1).
## Hessian positive definite, eigenvalue range [1.207714,208.0986].
## Model rank =  188 / 188 
## 
## 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.00e+00 -3.52e+08    0.88    0.64    
## s(mois)                                 1.00e+01 -2.50e+11    0.88    0.89    
## te(obs_lon,obs_lat)                     1.68e+02 -2.02e+11    0.74  <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "O2G1.Other_bony_fishes"

## 
## Method: REML   Optimizer: outer newton
## full convergence after 23 iterations.
## Gradient range [-7.35755e-05,2.027488]
## (score 636538 & scale 1).
## Hessian positive definite, eigenvalue range [0.01997831,3589.84].
## Model rank =  188 / 188 
## 
## 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.03 *
## s(mois)                                 10.00  10.00    0.98    0.88  
## te(obs_lon,obs_lat)                    168.00 159.33    0.95    0.18  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "O2G2.Other_bony_fishes"

## 
## Method: REML   Optimizer: outer newton
## full convergence after 9 iterations.
## Gradient range [-4.991576e-06,5.181546e-06]
## (score 573.5645 & scale 1).
## Hessian positive definite, eigenvalue range [1.952279,79.61888].
## Model rank =  188 / 188 
## 
## 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.41    0.46    
## s(mois)                                 10.00   9.91    0.40    0.09 .  
## te(obs_lon,obs_lat)                    168.00  90.60    0.33  <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "O1G1.Silky_shark"

## 
## Method: REML   Optimizer: outer newton
## full convergence after 7 iterations.
## Gradient range [-1.323792e-06,0.0008500532]
## (score 4889.556 & scale 1).
## Hessian positive definite, eigenvalue range [2.379979,1682.142].
## Model rank =  188 / 188 
## 
## 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.19    0.87   0.845  
## s(mois)                                 10.00   9.01    0.84   0.065 .
## te(obs_lon,obs_lat)                    168.00  64.40    0.86   0.545  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "O1G2.Silky_shark"

## 
## Method: REML   Optimizer: outer newton
## full convergence after 8 iterations.
## Gradient range [-8.408506e-07,0.0001429207]
## (score 2055.096 & scale 1).
## Hessian positive definite, eigenvalue range [1.615737,228.8458].
## Model rank =  188 / 188 
## 
## 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.79    0.96   0.575  
## s(mois)                                 10.00   7.94    0.95   0.335  
## te(obs_lon,obs_lat)                    168.00  36.94    0.89   0.035 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "O2G1.Silky_shark"

## 
## Method: REML   Optimizer: outer newton
## full convergence after 9 iterations.
## Gradient range [-3.567433e-06,0.005723907]
## (score 16657.55 & scale 1).
## Hessian positive definite, eigenvalue range [2.080501,3370.01].
## Model rank =  188 / 188 
## 
## 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.61    0.85  <2e-16 ***
## s(mois)                                 10.00   9.21    0.90    0.35    
## te(obs_lon,obs_lat)                    168.00 118.38    0.88    0.04 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "O2G2.Silky_shark"

## 
## Method: REML   Optimizer: outer newton
## full convergence after 10 iterations.
## Gradient range [-1.915054e-06,2.486232e-06]
## (score 222.7481 & scale 1).
## Hessian positive definite, eigenvalue range [0.1510041,23.77446].
## Model rank =  188 / 188 
## 
## 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.52    0.49    0.12    
## s(mois)                                 10.00   1.92    0.48    0.09 .  
## te(obs_lon,obs_lat)                    168.00  26.88    0.37  <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "O1G1.Common_dolphinfish"

## 
## Method: REML   Optimizer: outer newton
## full convergence after 7 iterations.
## Gradient range [0.0005143517,0.0459782]
## (score 30999.1 & scale 1).
## Hessian positive definite, eigenvalue range [3.683676,5004.906].
## Model rank =  188 / 188 
## 
## 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.610  
## s(mois)                                 10.00   9.95    0.92   0.015 *
## te(obs_lon,obs_lat)                    168.00 127.63    0.99   0.930  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "O1G2.Common_dolphinfish"

## 
## Method: REML   Optimizer: outer newton
## full convergence after 11 iterations.
## Gradient range [-0.000136261,-1.029578e-05]
## (score 427.5456 & scale 1).
## Hessian positive definite, eigenvalue range [0.0001362303,52.5642].
## Model rank =  188 / 188 
## 
## 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.26   0.310   
## s(mois)                                 10.00   9.39    0.24   0.005 **
## te(obs_lon,obs_lat)                    168.00  41.41    0.22   0.005 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "O2G1.Common_dolphinfish"

## 
## Method: REML   Optimizer: outer newton
## full convergence after 11 iterations.
## Gradient range [-5.408486e-05,0.08419109]
## (score 54236.95 & scale 1).
## Hessian positive definite, eigenvalue range [3.765559,3595.682].
## Model rank =  188 / 188 
## 
## 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.89    0.93   0.015 *
## s(mois)                                 10.00   9.95    0.97   0.680  
## te(obs_lon,obs_lat)                    168.00 146.83    0.95   0.335  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "O2G2.Common_dolphinfish"

## 
## Method: REML   Optimizer: outer newton
## full convergence after 12 iterations.
## Gradient range [-3.471609e-05,1.760379e-05]
## (score 156.0826 & scale 1).
## Hessian positive definite, eigenvalue range [4.794464e-06,33.05297].
## Model rank =  188 / 188 
## 
## 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.00e+00 8.60e+00    0.20   0.450   
## s(mois)                                1.00e+01 7.18e-04    0.19   0.065 . 
## te(obs_lon,obs_lat)                    1.68e+02 1.83e+01    0.17   0.010 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Plot gam se(saison)

for (p in res2){ 
o=getViz(p$gam)
print(plot(sm(o,2))+ xlab("mois") + ylab("effect mois") + ggtitle(paste(p$fn)))
  }

for (p in res2){ 
o=getViz(p$gam)
print(plot(sm(o,1))+ xlab("capture") + ylab("effect capture") + ggtitle(paste(p$fn)))
  }

for (p in res2){ 
o=getViz(p$gam)
plot<-plot(sm(o,3))+ xlab("longitude") + ylab("latitude") + ggtitle(paste(p$fn)) +l_fitRaster() + l_fitContour()+ annotation_map(map_data("world")) 
print(plot)
  }