-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathseg_sais.R
213 lines (166 loc) · 8.51 KB
/
seg_sais.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
#param <- read.table("/Users/aureliechagnon-lafortune/Desktop/MAITRISE/DATA/FINAL_DATA/param_plus.txt")
param <- read.table("DATA/FINAL_DATA/param_plus.txt")
library(dplyr)
library(segmented)
library(ggplot2)
#library(esquisse)
library(MuMIn)
library(png)
library(grid)
options(scipen = 999)
seg_sais <- function(param, type){
data <- param
# Enlever données manquantes ----
data <- dplyr::filter(data, !is.na(Snowmelt_day))
data<- dplyr::filter(data, !is.na(somme_sais_raw))
# data<- dplyr::filter(data, !Field_site == "Herschel")
# data <- dplyr::filter(data, Field_site == "Bylot"|Field_site == "Canning"|Field_site == "Colville"|Field_site == "Hochstetter"|Field_site == "Ikpikpuk"|Field_site == "Medusa"|Field_site == "Nome"|Field_site == "Southampton"|Field_site == "Utqiagvik"|Field_site == "Zackenberg")
#ajout colonne w ----
field_sites <- unique(data$Field_site)
m_data <- data.frame()
for (i in 1:length(field_sites)){
fs <- field_sites[i]
tab1<- data[which(data$Field_site==fs),]
tab1$w <- 1/length(tab1$Period_Year)
m_data <- rbind(m_data, tab1)
}
# m_data <- m_data[-c(28,29),]
# Modèles >1 variable ----
m_nul <- lm(somme_sais_raw~1, data=m_data, weights=w)
m_base <- lm(somme_sais_raw~dj196, data = m_data, weights = w)
seg_base<- segmented(m_base, seg.Z = ~dj196, data = m_data, weights = w)
# m_quad<-lm(somme_sais_raw~dj196+I(dj196^2), data = m_data, weights = w)
#m_qsnow<-lm(somme_sais_raw~dj196+I(dj196^2)+Snowmelt_day, data = m_data, weights = w)
m_lsnow <- lm(somme_sais_raw ~ dj196 + Snowmelt_day, data = m_data, weights = w)
seg_lsnow <- segmented(m_lsnow, seg.Z = ~dj196, data=m_data, weights = w )
m_lpr <- lm(somme_sais_raw ~ dj196 + pr196, data = m_data, weights = w)
#m_qpr <- lm(somme_sais_raw~dj196 +I(dj196^2) + pr196, data = m_data, weights = w)
seg_lpr <- segmented(m_lpr, seg.Z = ~dj196, data=m_data, weights = w )
m_lrad <-lm(somme_sais_raw ~ dj196 + rjuin5_juil15, data = m_data, weights = w)
#m_qrad <-lm(somme_sais_raw~dj196 +I(dj196^2)+rMoy_juin_juillet, data = m_data, weights = w)
seg_lrad <- segmented(m_lrad, seg.Z = ~ dj196, data=m_data, weights = w )
m_snow <- lm(somme_sais_raw ~ Snowmelt_day, data = m_data, weights = w)
models_sais <- list(mnul=m_nul, m_base=m_base, m_lrad=m_lrad, m_lpr=m_lpr, seg_base=seg_base, m_lsnow=m_lsnow, m_snow=m_snow, seg_lsnow=seg_lsnow, seg_lpr=seg_lpr, seg_lrad=seg_lrad)
mods_simple_sais <- model.sel(m_nul, m_base, m_lrad, m_lpr, seg_base, m_lsnow, m_snow, seg_lsnow, seg_lpr, seg_lrad)
write.csv(file="scripts/model_selection_sais.csv", mods_simple_sais)
# print(xtable(mods_simple_sais[,c(10:14)]))
# print(xtable(mods_simple_sais[,-c(5,8:10)]))
#réimporter data pour remettre toutes les données enlevées à cause de neige
data <- param
data<- dplyr::filter(data, !is.na(somme_sais_raw))
# data<- dplyr::filter(data, !Field_site == "Herschel")
field_sites <- unique(data$Field_site)
m_data <- data.frame()
for (i in 1:length(field_sites)){
fs <- field_sites[i]
tab1<- data[which(data$Field_site==fs),]
tab1$w <- 1/length(tab1$Period_Year)
m_data <- rbind(m_data, tab1) }
m_base <- lm(somme_sais_raw~dj196, data = m_data, weights = w)
seg_base<- segmented(m_base, seg.Z = ~dj196)
pt <- seg_base$psi[2]
sum <- summary(seg_base)
conf <- confint(seg_base)
int <- intercept(seg_base)
sl <- slope(seg_base)
# tab <- as.data.frame(slope(seg_base)$dj196)
# xtable(tab)
prd_dj <- data.frame(dj196=seq(from=range(m_data$dj196)[1], to=pt, length.out=200 ))
err <- predict(seg_base, newdata=prd_dj, re.form=NA, se.fit=T, interval = "confidence")
t <- as.data.frame(err$fit)
prd_dj$lci <- t$lwr
prd_dj$fit <- t$fit
prd_dj$uci <- t$upr
prd_dj2 <- data.frame(dj196=seq(from=pt, to=range(m_data$dj196)[2], length.out=200 ))
err <- predict(seg_base, newdata=prd_dj2, re.form=NA, se.fit=T, interval = "confidence")
t <- as.data.frame(err$fit)
prd_dj2$lci <- t$lwr
prd_dj2$fit <- t$fit
prd_dj2$uci <- t$upr
# gestion couleur ----
pal<-colorRamps::matlab.like(19)
pal[1]<-"#000000"
pal[2]<-"#5221B4"
pal[16]<-"#f24d06"
pal[17]<-"#e55709"
pal[18]<- "#b2154c"
pal[19]<-"#7a0e0e"
shape <- rep(c(21, 23, 24),7)
#icone<-readPNG("/Users/aureliechagnon-lafortune/Desktop/MAITRISE/visuel/icone_seas_bio.png")
icone<-readPNG("visuel/icone_seas_bio.png")
i <- rasterGrob(icone, interpolate=TRUE)
ymax<- range(data$somme_sais_raw, na.rm=T)[2]
ymin<-range(data$somme_sais_raw, na.rm=T)[1]
y_loc <- ymax-((ymax-ymin)/3)
if( type == "sep") {
data<-param
data<- data%>%
mutate(Field_site = forcats::fct_reorder(Field_site, TMOY))
# png(filename="/Users/aureliechagnon-lafortune/Desktop/MAITRISE/memoire/figures/saison_seg1.png",
# units="in",
# width=7,
# height=4,
# pointsize=12,
# res=300)
g <-
ggplot() +
theme_bw() + theme(
plot.background = element_blank()
,panel.grid.major = element_blank()
,panel.grid.minor = element_blank()
,panel.border = element_rect(color = "black")
,plot.margin = margin(0, 0, 0.5, 0.2, "cm")
) +theme(axis.title.x = element_blank(),legend.position="none", axis.line = element_line(color = 'black'), axis.ticks.length = unit(-0.15, "cm"), axis.text.y = element_text(margin=unit(c(0.35,0.25,0.35,0.25), "cm")), axis.text.x = element_text(margin=unit(c(-0.30,-0.3,-0.3,-0.3), "cm")))+
geom_ribbon(data=prd_dj, aes(x=dj196, ymin=lci, ymax=uci), fill="gray86")+
geom_smooth(data=prd_dj, aes(x=dj196, y= fit), stat = "identity", col= "gray35", linetype="solid")+
geom_ribbon(data=prd_dj2, aes(x=dj196, ymin=lci, ymax=uci), fill="gray86")+
geom_smooth(data=prd_dj2, aes(x=dj196, y= fit), stat = "identity", col= "gray35", linetype="dashed")+
geom_point(data = data, aes(x = dj196, y = somme_sais_raw, fill= Field_site, shape=Field_site), alpha=0.55, size =2)+scale_shape_manual(limits = param$Field_site, values=shape)+ scale_fill_manual(limits = param$Field_site, values=pal, name="Study site")+
guides(col = guide_legend(ncol = 2))+theme(axis.title.y = element_text(margin = unit(c(0,0,0,0), "cm")))+
labs(x= "Cumulative thaw degree days on July 15th", y= "Seasonal biomass (mg)")+ theme(legend.title = element_blank(),legend.text = element_text( size = 8), legend.key.size = unit(0.5, "cm"))+ annotate("text", x = 415, y=4400, label= "c)", size=4)+ annotation_custom(i, xmin=275, xmax=395, ymin=y_loc)
g
# print(g1)
# dev.off()
}else{
data<- param
tabmoy <- data.frame()
field_sites <- unique(data$Field_site)
for (i in 1:length(field_sites)){
Field_site <- field_sites[i]
tab1<- data[which(data$Field_site==Field_site),]
sais <- mean(tab1$somme_sais_raw)
sais_sd <- sd(tab1$somme_sais_raw)
dj196 <- mean(tab1$dj196)
dj196_sd<- sd(tab1$dj196)
TMOY<- unique(tab1$TMOY)
tab2 <- data.frame(Field_site, sais, dj196, TMOY, sais_sd, dj196_sd)
tabmoy <- rbind(tabmoy, tab2)
}
#graph moyennes ----
tabmoy<- tabmoy%>%
mutate(Field_site = forcats::fct_reorder(Field_site, TMOY))
g<-
ggplot() +
theme_bw() + theme(
plot.background = element_blank()
,panel.grid.major = element_blank()
,panel.grid.minor = element_blank()
,panel.border = element_rect(color = "black")
,plot.margin = margin(0, 0, 0.5, 0.2, "cm")
) +theme(axis.title.x = element_blank(),legend.position="none", axis.line = element_line(color = 'black'), axis.ticks.length = unit(-0.15, "cm"), axis.text.y = element_text(margin=unit(c(0.35,0.25,0.35,0.25), "cm")), axis.text.x = element_text(margin=unit(c(-0.30,-0.3,-0.3,-0.3), "cm")))+
geom_ribbon(data=prd_dj, aes(x=dj196, ymin=lci, ymax=uci), fill="gray86")+
geom_smooth(data=prd_dj, aes(x=dj196, y= fit), stat = "identity", col= "gray35", linetype="solid")+
geom_ribbon(data=prd_dj2, aes(x=dj196, ymin=lci, ymax=uci), fill="gray86")+
geom_smooth(data=prd_dj2, aes(x=dj196, y= fit), stat = "identity", col= "gray35", linetype="dashed")+
geom_errorbar(data=tabmoy, aes(x=dj196,ymin=sais-sais_sd, ymax=sais+sais_sd), width=range(tabmoy$dj196)[2]/75, col="gray30")+
geom_errorbarh(data=tabmoy, aes(xmin=dj196-dj196_sd, xmax=dj196+dj196_sd, y=sais, height=range(tabmoy$sais)[2]/65), col="gray30")+
geom_point(data = tabmoy, aes(x = dj196, y = sais, fill= Field_site, shape=Field_site), size =5)+scale_shape_manual(values=shape)+ scale_fill_manual(values=pal)+
guides(col = guide_legend(ncol = 2))+theme(axis.title.y = element_text(margin = unit(c(0,0,0,0), "cm")))+
labs(x= "Cumulative thaw degree days on July 15th", y= "Seasonal biomass (mg)")+
theme(
legend.title = element_blank(),
legend.text = element_text( size = 8), legend.key.size = unit(0.5, "cm"))+ annotate("text", x = 415, y=4400, label= "c)", size=4)
}
# print(g2)
return(list(g=g, sum=sum, conf=conf, sl=sl, int=int, pt=pt, models_sais=models_sais))
}