---

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')
# pour le modele gam
library('mgcViz')
library('ggpubr')
library("maps") # pour la carte sur le graph 2D 
library("lubridate")
library(dplyr)
library(cowplot)

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

Diagnostics

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

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)
  }