Project

General

Profile

« Previous | Next » 

Revision 992b2eeb

Added by Benoit Parmentier about 11 years ago

analyses for multi-time scale paper, monthly hold-out proportions (0-70%) for all methods

View differences:

climate/research/oregon/interpolation/analyses_papers_methods_comparison_part5.R
75 75
in_dir8 <-"/data/project/layers/commons/Oregon_interpolation/output_data_365d_kriging_fus_lst_comb3_09122013"
76 76
in_dir9 <-"/data/project/layers/commons/Oregon_interpolation/output_data_365d_kriging_fus_lst_comb3_09132013"
77 77
in_dir10 <-"/data/project/layers/commons/Oregon_interpolation/output_data_365d_kriging_fus_lst_comb3_09142013"
78
in_dir11 <-"/data/project/layers/commons/Oregon_interpolation/output_data_365d_gwr_CAI_lst_comb3_09162013"
79
in_dir12 <-"/data/project/layers/commons/Oregon_interpolation/output_data_365d_gwr_CAI_lst_comb3_09172013"
80
in_dir13 <-"/data/project/layers/commons/Oregon_interpolation/output_data_365d_gwr_cai_lst_comb3_09282013"
81
in_dir14 <-"/data/project/layers/commons/Oregon_interpolation/output_data_365d_gwr_fus_lst_comb3_09232013"
82
in_dir15 <-"/data/project/layers/commons/Oregon_interpolation/output_data_365d_gwr_fus_lst_comb3_09262013"
78 83

  
84
#better as list and load one by one specific element from the object
79 85
raster_prediction_obj1 <-load_obj(file.path(in_dir1,"raster_prediction_obj_gam_CAI_dailyTmax_365d_gam_CAI_lst_comb3_08312013.RData"))
80 86
raster_prediction_obj2 <-load_obj(file.path(in_dir2,"raster_prediction_obj_gam_CAI_dailyTmax_365d_gam_CAI_lst_comb3_09012013.RData"))
81 87
raster_prediction_obj3 <-load_obj(file.path(in_dir3,"raster_prediction_obj_gam_CAI_dailyTmax_365d_gam_CAI_lst_comb3_09032013.RData"))
......
87 93
raster_prediction_obj8 <-load_obj(file.path(in_dir8,"raster_prediction_obj_kriging_fusion_dailyTmax_365d_kriging_fus_lst_comb3_09122013.RData"))
88 94
raster_prediction_obj9 <-load_obj(file.path(in_dir9,"raster_prediction_obj_kriging_fusion_dailyTmax_365d_kriging_fus_lst_comb3_09132013.RData"))
89 95
raster_prediction_obj10 <-load_obj(file.path(in_dir10,"raster_prediction_obj_kriging_fusion_dailyTmax_365d_kriging_fus_lst_comb3_09142013.RData"))
96
raster_prediction_obj11 <-load_obj(file.path(in_dir11,"raster_prediction_obj_gwr_CAI_dailyTmax_365d_gwr_CAI_lst_comb3_09162013.RData"))
97
raster_prediction_obj12 <-load_obj(file.path(in_dir12,"raster_prediction_obj_gwr_CAI_dailyTmax_365d_gwr_CAI_lst_comb3_09172013.RData"))
98
raster_prediction_obj13 <-load_obj(file.path(in_dir13,"raster_prediction_obj_gwr_CAI_dailyTmax_365d_gwr_cai_lst_comb3_09282013.RData"))
99
raster_prediction_obj14 <-load_obj(file.path(in_dir14,"raster_prediction_obj_gwr_fusion_dailyTmax_365d_gwr_fus_lst_comb3_09232013.RData"))
100
raster_prediction_obj15 <-load_obj(file.path(in_dir15,"raster_prediction_obj_gwr_fusion_dailyTmax_365d_gwr_fus_lst_comb3_09262013.RData"))
101

  
102
#raster_prediction_obj_gwr_CAI_dailyTmax_365d_gwr_cai_lst_comb3_09282013.RData
103
#raster_prediction_obj_gwr_fusion_dailyTmax_365d_gwr_fus_lst_comb3_09262013.RData
90 104

  
91 105
out_dir<-"/home/parmentier/Data/IPLANT_project/paper_multitime_scale__analyses_tables_fig_09032013"
92 106
setwd(out_dir)
93 107
y_var_name <- "dailyTmax"
94 108
y_var_month <- "TMax"
95 109
#y_var_month <- "LSTD_bias"
96
out_suffix <- "_OR_09132013"
110
out_suffix <- "_OR_09292013"
97 111
#script_path<-"/data/project/layers/commons/data_workflow/env_layers_scripts/"
98 112
#### FUNCTION USED IN SCRIPT
99 113

  
100
function_analyses_paper <-"contribution_of_covariates_paper_interpolation_functions_09092013.R"
114
function_analyses_paper <-"contribution_of_covariates_paper_interpolation_functions_09232013.R"
101 115

  
102 116
script_path<-"/home/parmentier/Data/IPLANT_project/env_layers_scripts/" #path to script
103 117
source(file.path(script_path,function_analyses_paper)) #source all functions used in this script.
......
112 126
tb_s_gam_CAI <-rbind(raster_prediction_obj1$tb_diagnostic_s,raster_prediction_obj2$tb_diagnostic_s,raster_prediction_obj3$tb_diagnostic_s)
113 127
#prop_obj_gam_CAI_v <- calc_stat_prop_tb_diagnostic(names_mod,names_id,tb_v)
114 128

  
129
tb_mv_gwr_CAI <-rbind(raster_prediction_obj11$tb_month_diagnostic_v,raster_prediction_obj12$tb_month_diagnostic_v,raster_prediction_obj13$tb_month_diagnostic_v)
130
tb_ms_gwr_CAI <-rbind(raster_prediction_obj11$tb_month_diagnostic_s,raster_prediction_obj12$tb_month_diagnostic_s,raster_prediction_obj13$tb_month_diagnostic_s)
131

  
132
tb_v_gwr_CAI <-rbind(raster_prediction_obj11$tb_diagnostic_v,raster_prediction_obj12$tb_diagnostic_v,raster_prediction_obj13$tb_diagnostic_v)
133
tb_s_gwr_CAI <-rbind(raster_prediction_obj11$tb_diagnostic_s,raster_prediction_obj12$tb_diagnostic_s,raster_prediction_obj13$tb_diagnostic_s)
134

  
115 135
tb_mv_kriging_CAI <- raster_prediction_obj4$tb_month_diagnostic_v
116 136
tb_ms_kriging_CAI <- raster_prediction_obj4$tb_month_diagnostic_s
117 137

  
......
126 146
tb_v_gam_fus <-rbind(raster_prediction_obj5$tb_diagnostic_v,raster_prediction_obj6$tb_diagnostic_v,raster_prediction_obj7$tb_diagnostic_v)
127 147
tb_s_gam_fus <-rbind(raster_prediction_obj5$tb_diagnostic_s,raster_prediction_obj6$tb_diagnostic_s,raster_prediction_obj7$tb_diagnostic_s)
128 148

  
149
tb_mv_gwr_fus <-rbind(raster_prediction_obj14$tb_month_diagnostic_v,raster_prediction_obj15$tb_month_diagnostic_v)
150
tb_ms_gwr_fus <-rbind(raster_prediction_obj14$tb_month_diagnostic_s,raster_prediction_obj15$tb_month_diagnostic_s)
151

  
152
tb_v_gwr_fus <-rbind(raster_prediction_obj14$tb_diagnostic_v,raster_prediction_obj15$tb_diagnostic_v)
153
tb_s_gwr_fus <-rbind(raster_prediction_obj14$tb_diagnostic_s,raster_prediction_obj15$tb_diagnostic_s)
154

  
129 155
tb_mv_kriging_fus <-rbind(raster_prediction_obj8$tb_month_diagnostic_v,raster_prediction_obj9$tb_month_diagnostic_v,raster_prediction_obj10$tb_month_diagnostic_v)
130 156
tb_ms_kriging_fus <-rbind(raster_prediction_obj8$tb_month_diagnostic_s,raster_prediction_obj9$tb_month_diagnostic_s,raster_prediction_obj10$tb_month_diagnostic_s)
131 157

  
132 158
tb_v_kriging_fus <-rbind(raster_prediction_obj8$tb_diagnostic_v,raster_prediction_obj9$tb_diagnostic_v,raster_prediction_obj10$tb_diagnostic_v)
133 159
tb_s_kriging_fus <-rbind(raster_prediction_obj8$tb_diagnostic_s,raster_prediction_obj9$tb_diagnostic_s,raster_prediction_obj10$tb_diagnostic_s)
134 160

  
135
#list_tb <-list(tb_v_gam_CAI,tb_v_kriging_CAI,tb_s_gam_CAI,tb_s_kriging_CAI,tb_mv_gam_CAI,tb_mv_kriging_CAI,tb_ms_gam_CAI,tb_ms_kriging_CAI) #Add fusion here
136
#names(list_tb) <- c("tb_v_gam_CAI","tb_v_kriging_CAI","tb_s_gam_CAI","tb_s_kriging_CAI","tb_mv_gam_CAI","tb_mv_kriging_CAI","tb_ms_gam_CAI","tb_ms_kriging_CAI") #Add fusion here
161
list_tb <- list(tb_v_gam_CAI,tb_v_kriging_CAI,tb_v_gwr_CAI,tb_s_gam_CAI,tb_s_kriging_CAI,tb_s_gwr_CAI,
162
           tb_mv_gam_CAI,tb_mv_kriging_CAI,tb_mv_gwr_CAI,tb_ms_gam_CAI,tb_ms_kriging_CAI,tb_ms_gwr_CAI,
163
           tb_v_gam_fus,tb_v_kriging_fus,tb_v_gwr_fus,tb_s_gam_fus,tb_s_kriging_fus,tb_s_gwr_fus,
164
           tb_mv_gam_fus,tb_mv_kriging_fus,tb_mv_gwr_fus,tb_ms_gam_fus,tb_ms_kriging_fus,tb_ms_gwr_fus)
165

  
166
names(list_tb) <- c("tb_v_gam_CAI","tb_v_kriging_CAI","tb_v_gwr_CAI","tb_s_gam_CAI","tb_s_kriging_CAI","tb_s_gwr_CAI",
167
        "tb_mv_gam_CAI","tb_mv_kriging_CAI","tb_mv_gwr_CAI","tb_ms_gam_CAI","tb_ms_kriging_CAI","tb_ms_gwr_CAI",
168
        "tb_v_gam_fus","tb_v_kriging_fus","tb_v_gwr_fus","tb_s_gam_fus","tb_s_kriging_fus","tb_s_gwr_fus",
169
        "tb_mv_gam_fus","tb_mv_kriging_fus","tb_mv_gwr_fus","tb_ms_gam_fus","tb_ms_kriging_fus","tb_ms_gwr_fus")
137 170

  
138
list_tb <-list(tb_v_gam_CAI,tb_v_kriging_CAI,tb_s_gam_CAI,tb_s_kriging_CAI,tb_mv_gam_CAI,tb_mv_kriging_CAI,tb_ms_gam_CAI,tb_ms_kriging_CAI, #Add fusion here
139
               tb_v_gam_fus,tb_v_kriging_fus,tb_s_gam_fus,tb_s_kriging_fus,tb_mv_gam_fus,tb_mv_kriging_fus,tb_ms_gam_fus,tb_ms_kriging_fus) #Add fusion here
140
names(list_tb) <- c("tb_v_gam_CAI","tb_v_kriging_CAI","tb_s_gam_CAI","tb_s_kriging_CAI","tb_mv_gam_CAI","tb_mv_kriging_CAI","tb_ms_gam_CAI","tb_ms_kriging_CAI", #Add fusion here
141
                   "tb_v_gam_fus","tb_v_kriging_fus","tb_s_gam_fus","tb_s_kriging_fus","tb_mv_gam_fus","tb_mv_kriging_fus","tb_ms_gam_fus","tb_ms_kriging_fus") #Add fusion here
171
#list_tb <-list(tb_v_gam_CAI,tb_v_kriging_CAI,tb_v_gwr_CAI,tb_s_gam_CAI,tb_s_kriging_CAI,tb_s_gwr_CAI,tb_mv_gam_CAI,tb_mv_kriging_CAI,tb_ms_gam_CAI,tb_ms_kriging_CAI,tb_ms_gwr_CAI,tb_ms_gwr_CAI #Add fusion here
172
#               tb_v_gam_fus,tb_v_kriging_fus,tb_s_gam_fus,tb_s_kriging_fus,tb_mv_gam_fus,tb_mv_kriging_fus,tb_ms_gam_fus,tb_ms_kriging_fus) #Add fusion here
173
#names(list_tb) <- c("tb_v_gam_CAI","tb_v_kriging_CAI","tb_v_gwr_CAI","tb_s_gam_CAI","tb_s_kriging_CAI","tb_s_gwr_CAI","tb_mv_gam_CAI","tb_mv_kriging_CAI","tb_ms_gam_CAI","tb_ms_kriging_CAI","tb_ms_gwr_CAI","tb_ms_gwr_CAI" #Add fusion here
174
#                   "tb_v_gam_fus","tb_v_kriging_fus","tb_s_gam_fus","tb_s_kriging_fus","tb_mv_gam_fus","tb_mv_kriging_fus","tb_ms_gam_fus","tb_ms_kriging_fus") #Add fusion here
142 175

  
143
##### DAILY AVERAGE ACCURACY : PLOT AND DIFFERENCES...
176
##### DAILY AVERAGE ACCURACY : PLOT AND DIFFERENCES...Cd
144 177

  
145 178
for(i in 1:length(list_tb)){
146
  i <- i+1
179
  #i <- i+1
147 180
  tb <-list_tb[[i]]
148 181
  plot_name <- names(list_tb)[i]
149 182
  pat_str <- "tb_m"
......
167 200
  png(paste("Figure__accuracy_rmse_prop_month_",plot_name,out_suffix,".png", sep=""),
168 201
      height=480*layout_m[1],width=480*layout_m[2])
169 202
    
170
  xyplot(as.formula(plot_formula),group=pred_mod,type="b",
203
  p<- xyplot(as.formula(plot_formula),group=pred_mod,type="b",
171 204
          data=avg_tb,
172 205
          main=paste("rmse ",plot_name,sep=" "),
173 206
          pch=1:length(avg_tb$pred_mod),
174 207
          par.settings=list(superpose.symbol = list(
175 208
          pch=1:length(avg_tb$pred_mod))),
176 209
          auto.key=list(columns=5))
210
  print(p)
177 211
  
178 212
  dev.off()
179 213
  
......
199 233
metric_names <- c("mae","rmse","me","r")
200 234
diff_kriging_CAI <- diff_df(tb_s_kriging_CAI,tb_v_kriging_CAI,metric_names)
201 235
diff_gam_CAI <- diff_df(tb_s_gam_CAI,tb_v_gam_CAI,metric_names)
236
diff_gwr_CAI <- diff_df(tb_s_gwr_CAI,tb_v_gwr_CAI,metric_names)
237

  
238
layout_m<-c(1,1) #one row two columns
239
par(mfrow=layout_m)
202 240

  
203
boxplot(diff_kriging_CAI$rmse,diff_gam_CAI$rmse,names=c("kriging_CAI","gam_CAI"),
241
png(paste("Figure__accuracy_rmse_prop_month_",plot_name,out_suffix,".png", sep=""),
242
    height=480*layout_m[1],width=480*layout_m[2])
243
boxplot(diff_kriging_CAI$rmse,diff_gam_CAI$rmse,diff_gwr_CAI$rmse,names=c("kriging_CAI","gam_CAI","gwr_CAI"),
204 244
        main="Difference between training and testing daily rmse")
245
dev.off()
205 246

  
247
#remove prop 0,
248
diff_kriging_CAI <- diff_df(tb_s_kriging_CAI[tb_s_kriging_CAI$prop_month!=0,],tb_v_kriging_CAI[tb_v_kriging_CAI$prop_month!=0,],metric_names)
249
diff_gam_CAI <- diff_df(tb_s_gam_CAI[tb_s_gam_CAI$prop_month!=0,],tb_v_gam_CAI[tb_v_gam_CAI$prop_month!=0,],metric_names)
250
diff_gwr_CAI <- diff_df(tb_s_gwr_CAI[tb_s_gwr_CAI$prop_month!=0,],tb_v_gwr_CAI[tb_v_gwr_CAI$prop_month!=0,],metric_names)
251
boxplot(diff_kriging_CAI$rmse,diff_gam_CAI$rmse,diff_gwr_CAI$rmse,names=c("kriging_CAI","gam_CAI","gwr_CAI"),
252
        main="Difference between training and testing daily rmse")
253

  
254
#now monthly accuracy
206 255
metric_names <- c("mae","rmse","me","r")
207
diff_kriging_m_CAI <- diff_df(tb_ms_kriging_CAI,tb_mv_kriging_CAI,metric_names)
208
diff_gam_m_CAI <- diff_df(tb_ms_gam_CAI,tb_mv_gam_CAI,metric_names)
256
diff_kriging_m_CAI <- diff_df(tb_ms_kriging_CAI[tb_ms_kriging_CAI$prop!=0,],tb_mv_kriging_CAI,metric_names)
257
diff_gam_m_CAI <- diff_df(tb_ms_gam_CAI[tb_ms_gam_CAI$prop!=0,],tb_mv_gam_CAI,metric_names)
258
diff_gwr_m_CAI <- diff_df(tb_ms_gwr_CAI[tb_ms_gwr_CAI$prop!=0,],tb_mv_gwr_CAI,metric_names)
259

  
260
layout_m<-c(1,1) #one row two columns
261
par(mfrow=layout_m)
209 262

  
210
boxplot(diff_kriging_m_CAI$rmse,diff_gam_m_CAI$rmse,names=c("kriging_CAI","gam_CAI"),
263
png(paste("Figure__accuracy_rmse_prop_month_",plot_name,out_suffix,".png", sep=""),
264
    height=480*layout_m[1],width=480*layout_m[2])
265
boxplot(diff_kriging_m_CAI$rmse,diff_gam_m_CAI$rmse,diff_gwr_m_CAI$rmse,names=c("kriging_CAI","gam_CAI","gwr_CAI"),
211 266
        main="Difference between training and monhtly testing rmse")
267
dev.off()
268

  
269
#boxplot(diff_kriging_m_CAI$rmse,diff_gam_m_CAI$rmse,diff_gwr_CAI,names=c("kriging_CAI","gam_CAI","gwr_CAI"),
270
#        main="Difference between training and monhtly testing rmse")
212 271

  
213 272
### For fusion
214 273

  
215 274
metric_names <- c("mae","rmse","me","r")
216 275
diff_kriging_fus <- diff_df(tb_s_kriging_fus,tb_v_kriging_fus,metric_names)
217 276
diff_gam_fus <- diff_df(tb_s_gam_fus,tb_v_gam_fus,metric_names)
277
diff_gwr_fus <- diff_df(tb_s_gwr_fus,tb_v_gwr_fus,metric_names)
278

  
279
layout_m<-c(1,1) #one row two columns
280
par(mfrow=layout_m)
218 281

  
219
boxplot(diff_kriging_fus$rmse,diff_gam_fus$rmse,names=c("kriging_fus","gam_fus"),
282
png(paste("Figure__accuracy_rmse_prop_month_",plot_name,out_suffix,".png", sep=""),
283
    height=480*layout_m[1],width=480*layout_m[2])
284
boxplot(diff_kriging_fus$rmse,diff_gam_fus$rmse,diff_gwr_fus$rmse,names=c("kriging_fus","gam_fus","gwr_fus"),
220 285
        main="Difference between training and testing daily rmse")
286
dev.off()
221 287

  
222 288
metric_names <- c("mae","rmse","me","r")
223
diff_kriging_m_fus <- diff_df(tb_ms_kriging_fus,tb_mv_kriging_fus,metric_names)
224
diff_gam_m_fus <- diff_df(tb_ms_gam_fus,tb_mv_gam_fus,metric_names)
289
diff_kriging_m_fus <- diff_df(tb_ms_kriging_fus[tb_ms_kriging_fus$prop!=0,],tb_mv_kriging_fus[tb_mv_kriging_fus$prop!=0,],metric_names)
290
diff_gam_m_fus <- diff_df(tb_ms_gam_fus[tb_ms_gam_fus$prop!=0,],tb_mv_gam_fus[tb_mv_gam_fus$prop!=0,],metric_names)
291
diff_gwr_m_fus <- diff_df(tb_ms_gwr_fus[tb_ms_gwr_fus$prop!=0,],tb_mv_gwr_fus[tb_mv_gwr_fus$prop!=0,],metric_names)
292

  
293
layout_m<-c(1,1) #one row two columns
294
par(mfrow=layout_m)
225 295

  
226
boxplot(diff_kriging_m_fus$rmse,diff_gam_m_fus$rmse,names=c("kriging_fus","gam_fus"),
296
png(paste("Figure__accuracy_rmse_prop_month_",plot_name,out_suffix,".png", sep=""),
297
    height=480*layout_m[1],width=480*layout_m[2])
298
boxplot(diff_kriging_m_fus$rmse,diff_gam_m_fus$rmse,diff_gwr_m_fus$rmse, names=c("kriging_fus","gam_fus","gwr_fus"),
227 299
        main="Difference between training and testing FUS rmse")
300
dev.off()
228 301

  
229 302
### NOW PLOT OF COMPARISON BETWEEN Kriging and GAM
230 303

  
......
235 308
tb_v_kriging_CAI
236 309
tb_v_kriging_fus
237 310

  
238
methods_names <- c("tb_v_gam_CAI","tb_v_gam_fus","tb_v_kriging_CAI","tb_v_kriging_fus")
311
methods_names <- c("tb_v_gam_CAI","tb_v_gam_fus","tb_v_kriging_CAI","tb_v_kriging_fus","tb_v_gwr_CAI","tb_v_gwr_fus")
239 312
list_prop_obj <- vector("list",length=length(methods_names))
240 313
for(i in 1:length(methods_names)){
241 314
  tb <- list_tb[[methods_names[i]]]
......
248 321

  
249 322
ac_prop_tb_list <- extract_list_from_list_obj(list_prop_obj,"avg_tb")
250 323
nb_rows <- sapply(ac_prop_tb_list,FUN=nrow)
251
method_interp_names<-c("gam_CAI","gam_fus","kriging_CAI","kriging_fus")
324
method_interp_names<-c("gam_CAI","gam_fus","kriging_CAI","kriging_fus","gwr_CAI","gwr_fus")
252 325
for(i in 1:length(methods_names)){
253 326
  avg_tb<-ac_prop_tb_list[[i]]
254 327
  avg_tb$method_interp <-rep(x=method_interp_names[i],times=nb_rows[i])
......
258 331

  
259 332
#names(ac_prop_tb_list) <- names(list_prop_obj)
260 333

  
261
t44 <- do.call(rbind,ac_prop_tb_list)
334
t44 <- do.call(rbind,ac_prop_tb_list) #contains all accuracy by method, proportion, model, sample etc.
262 335
View(t44)
263 336

  
264
t44[which.min(t44$rmse),]
337
t44[which.min(t44$rmse),] #Find the mimum rmse across all models and methods...
265 338

  
266 339
test <- t44[order(t44$rmse),]
340
test[1:24,]
267 341

  
268 342
test2<-test[test$method_interp%in% c("gam_fus","gam_CAI"),]
269 343
test2[1:24,]

Also available in: Unified diff