Revision a8c14dda
Added by Benoit Parmentier about 11 years ago
| climate/research/oregon/interpolation/analyses_papers_methods_comparison_part5.R | ||
|---|---|---|
| 4 | 4 |
#It uses inputs from interpolation objects created at earlier stages... |
| 5 | 5 |
#Note that this is exploratory code i.e. not part of the worklfow. |
| 6 | 6 |
#AUTHOR: Benoit Parmentier # |
| 7 |
#DATE: 09/06/2013 #
|
|
| 7 |
#DATE: 09/13/2013 #
|
|
| 8 | 8 |
#PROJECT: NCEAS INPLANT: Environment and Organisms --TASK#491-- # |
| 9 | 9 |
################################################################################################### |
| 10 | 10 |
|
| ... | ... | |
| 62 | 62 |
################## PARAMETERS ########## |
| 63 | 63 |
|
| 64 | 64 |
|
| 65 |
#"/home/parmentier/Data/IPLANT_project/Oregon_interpolation/Oregon_03142013/output_data_365d_gam_day_lst_comb4_07152013/"
|
|
| 65 |
#path to gam CAI and kriging analyes with hold-out
|
|
| 66 | 66 |
in_dir1 <- "/home/parmentier/Data/IPLANT_project/Oregon_interpolation/Oregon_03142013/output_data_365d_gam_CAI_lst_comb3_08312013/" |
| 67 |
in_dir <- "/home/parmentier/Data/IPLANT_project/Oregon_interpolation/Oregon_03142013/output_data_365d_gam_CAI_lst_comb3_09012013" |
|
| 68 |
in_dir2 <-"/home/parmentier/Data/IPLANT_project/Oregon_interpolation/Oregon_03142013/output_data_365d_gam_CAI_lst_comb3_09032013"
|
|
| 67 |
in_dir2 <- "/home/parmentier/Data/IPLANT_project/Oregon_interpolation/Oregon_03142013/output_data_365d_gam_CAI_lst_comb3_09012013"
|
|
| 68 |
in_dir3 <-"/home/parmentier/Data/IPLANT_project/Oregon_interpolation/Oregon_03142013/output_data_365d_gam_CAI_lst_comb3_09032013"
|
|
| 69 | 69 |
in_dir4 <- "/home/parmentier/Data/IPLANT_project/Oregon_interpolation/Oregon_03142013/output_data_365d_kriging_CAI_lst_comb3_09042013" |
| 70 | 70 |
|
| 71 |
#path to gam fusion and kriging fusion analyes with hold-out |
|
| 72 |
in_dir5 <- "/data/project/layers/commons/Oregon_interpolation/output_data_365d_gam_fus_lst_comb3_09092013" |
|
| 73 |
in_dir6 <- "/data/project/layers/commons/Oregon_interpolation/output_data_365d_gam_fus_lst_comb3_09102013" |
|
| 74 |
in_dir7 <-"/data/project/layers/commons/Oregon_interpolation/output_data_365d_gam_fus_lst_comb3_09112013" |
|
| 75 |
in_dir8 <-"/data/project/layers/commons/Oregon_interpolation/output_data_365d_kriging_fus_lst_comb3_09122013" |
|
| 76 |
in_dir9 <-"/data/project/layers/commons/Oregon_interpolation/output_data_365d_kriging_fus_lst_comb3_09132013" |
|
| 77 |
in_dir10 <-"/data/project/layers/commons/Oregon_interpolation/output_data_365d_kriging_fus_lst_comb3_09142013" |
|
| 78 |
|
|
| 71 | 79 |
raster_prediction_obj1 <-load_obj(file.path(in_dir1,"raster_prediction_obj_gam_CAI_dailyTmax_365d_gam_CAI_lst_comb3_08312013.RData")) |
| 72 |
raster_prediction_obj <-load_obj(file.path(in_dir,"raster_prediction_obj_gam_CAI_dailyTmax_365d_gam_CAI_lst_comb3_09012013.RData"))
|
|
| 73 |
raster_prediction_obj2 <-load_obj(file.path(in_dir2,"raster_prediction_obj_gam_CAI_dailyTmax_365d_gam_CAI_lst_comb3_09032013.RData"))
|
|
| 80 |
raster_prediction_obj2 <-load_obj(file.path(in_dir2,"raster_prediction_obj_gam_CAI_dailyTmax_365d_gam_CAI_lst_comb3_09012013.RData"))
|
|
| 81 |
raster_prediction_obj3 <-load_obj(file.path(in_dir3,"raster_prediction_obj_gam_CAI_dailyTmax_365d_gam_CAI_lst_comb3_09032013.RData"))
|
|
| 74 | 82 |
raster_prediction_obj4 <-load_obj(file.path(in_dir4,"raster_prediction_obj_kriging_CAI_dailyTmax_365d_kriging_CAI_lst_comb3_09042013.RData")) |
| 75 |
|
|
| 83 |
|
|
| 84 |
raster_prediction_obj5 <-load_obj(file.path(in_dir5,"raster_prediction_obj_gam_fusion_dailyTmax_365d_gam_fus_lst_comb3_09092013.RData")) |
|
| 85 |
raster_prediction_obj6 <-load_obj(file.path(in_dir6,"raster_prediction_obj_gam_fusion_dailyTmax_365d_gam_fus_lst_comb3_09102013.RData")) |
|
| 86 |
raster_prediction_obj7 <-load_obj(file.path(in_dir7,"raster_prediction_obj_gam_fusion_dailyTmax_365d_gam_fus_lst_comb3_09112013.RData")) |
|
| 87 |
raster_prediction_obj8 <-load_obj(file.path(in_dir8,"raster_prediction_obj_kriging_fusion_dailyTmax_365d_kriging_fus_lst_comb3_09122013.RData")) |
|
| 88 |
raster_prediction_obj9 <-load_obj(file.path(in_dir9,"raster_prediction_obj_kriging_fusion_dailyTmax_365d_kriging_fus_lst_comb3_09132013.RData")) |
|
| 89 |
raster_prediction_obj10 <-load_obj(file.path(in_dir10,"raster_prediction_obj_kriging_fusion_dailyTmax_365d_kriging_fus_lst_comb3_09142013.RData")) |
|
| 90 |
|
|
| 76 | 91 |
out_dir<-"/home/parmentier/Data/IPLANT_project/paper_multitime_scale__analyses_tables_fig_09032013" |
| 77 | 92 |
setwd(out_dir) |
| 78 | 93 |
y_var_name <- "dailyTmax" |
| 79 | 94 |
y_var_month <- "TMax" |
| 80 | 95 |
#y_var_month <- "LSTD_bias" |
| 81 |
out_suffix <- "_OR_09032013"
|
|
| 96 |
out_suffix <- "_OR_09132013"
|
|
| 82 | 97 |
#script_path<-"/data/project/layers/commons/data_workflow/env_layers_scripts/" |
| 83 | 98 |
#### FUNCTION USED IN SCRIPT |
| 84 | 99 |
|
| 85 |
function_analyses_paper <-"contribution_of_covariates_paper_interpolation_functions_08152013.R"
|
|
| 100 |
function_analyses_paper <-"contribution_of_covariates_paper_interpolation_functions_09092013.R"
|
|
| 86 | 101 |
|
| 87 | 102 |
script_path<-"/home/parmentier/Data/IPLANT_project/env_layers_scripts/" #path to script |
| 88 | 103 |
source(file.path(script_path,function_analyses_paper)) #source all functions used in this script. |
| ... | ... | |
| 90 | 105 |
################################################################################# |
| 91 | 106 |
############ ANALYSES 1: Average accuracy per proportion for monthly hold out in muli-timescale mehtods... ####### |
| 92 | 107 |
|
| 93 |
tb_mv_gam_CAI <-rbind(raster_prediction_obj1$tb_month_diagnostic_v,raster_prediction_obj$tb_month_diagnostic_v,raster_prediction_obj2$tb_month_diagnostic_v) |
|
| 94 |
tb_ms_gam_CAI <-rbind(raster_prediction_obj1$tb_month_diagnostic_s,raster_prediction_obj$tb_month_diagnostic_s,raster_prediction_obj2$tb_month_diagnostic_s) |
|
| 95 |
|
|
| 96 |
tb_v_gam_CAI <-rbind(raster_prediction_obj1$tb_diagnostic_v,raster_prediction_obj$tb_diagnostic_v,raster_prediction_obj2$tb_diagnostic_v) |
|
| 97 |
tb_s_gam_CAI <-rbind(raster_prediction_obj1$tb_diagnostic_s,raster_prediction_obj$tb_diagnostic_s,raster_prediction_obj2$tb_diagnostic_s) |
|
| 108 |
tb_mv_gam_CAI <-rbind(raster_prediction_obj1$tb_month_diagnostic_v,raster_prediction_obj2$tb_month_diagnostic_v,raster_prediction_obj3$tb_month_diagnostic_v) |
|
| 109 |
tb_ms_gam_CAI <-rbind(raster_prediction_obj1$tb_month_diagnostic_s,raster_prediction_obj2$tb_month_diagnostic_s,raster_prediction_obj3$tb_month_diagnostic_s) |
|
| 98 | 110 |
|
| 99 |
prop_obj_gam_CAI_v <- calc_stat_prop_tb_diagnostic(names_mod,names_id,tb_v) |
|
| 111 |
tb_v_gam_CAI <-rbind(raster_prediction_obj1$tb_diagnostic_v,raster_prediction_obj2$tb_diagnostic_v,raster_prediction_obj3$tb_diagnostic_v) |
|
| 112 |
tb_s_gam_CAI <-rbind(raster_prediction_obj1$tb_diagnostic_s,raster_prediction_obj2$tb_diagnostic_s,raster_prediction_obj3$tb_diagnostic_s) |
|
| 113 |
#prop_obj_gam_CAI_v <- calc_stat_prop_tb_diagnostic(names_mod,names_id,tb_v) |
|
| 100 | 114 |
|
| 101 | 115 |
tb_mv_kriging_CAI <- raster_prediction_obj4$tb_month_diagnostic_v |
| 102 | 116 |
tb_ms_kriging_CAI <- raster_prediction_obj4$tb_month_diagnostic_s |
| 103 | 117 |
|
| 104 | 118 |
tb_v_kriging_CAI <- raster_prediction_obj4$tb_diagnostic_v |
| 105 | 119 |
tb_s_kriging_CAI <- raster_prediction_obj4$tb_diagnostic_s |
| 106 |
|
|
| 107 |
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 |
|
| 108 |
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
|
|
| 120 |
|
|
| 121 |
### SAME for gam fusion |
|
| 122 |
|
|
| 123 |
tb_mv_gam_fus <-rbind(raster_prediction_obj5$tb_month_diagnostic_v,raster_prediction_obj6$tb_month_diagnostic_v,raster_prediction_obj7$tb_month_diagnostic_v) |
|
| 124 |
tb_ms_gam_fus <-rbind(raster_prediction_obj5$tb_month_diagnostic_s,raster_prediction_obj6$tb_month_diagnostic_s,raster_prediction_obj7$tb_month_diagnostic_s) |
|
| 125 |
|
|
| 126 |
tb_v_gam_fus <-rbind(raster_prediction_obj5$tb_diagnostic_v,raster_prediction_obj6$tb_diagnostic_v,raster_prediction_obj7$tb_diagnostic_v) |
|
| 127 |
tb_s_gam_fus <-rbind(raster_prediction_obj5$tb_diagnostic_s,raster_prediction_obj6$tb_diagnostic_s,raster_prediction_obj7$tb_diagnostic_s) |
|
| 128 |
|
|
| 129 |
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 |
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 |
|
|
| 132 |
tb_v_kriging_fus <-rbind(raster_prediction_obj8$tb_diagnostic_v,raster_prediction_obj9$tb_diagnostic_v,raster_prediction_obj10$tb_diagnostic_v) |
|
| 133 |
tb_s_kriging_fus <-rbind(raster_prediction_obj8$tb_diagnostic_s,raster_prediction_obj9$tb_diagnostic_s,raster_prediction_obj10$tb_diagnostic_s) |
|
| 134 |
|
|
| 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
|
|
| 137 |
|
|
| 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 |
|
| 109 | 142 |
|
| 110 | 143 |
##### DAILY AVERAGE ACCURACY : PLOT AND DIFFERENCES... |
| 111 | 144 |
|
| 112 | 145 |
for(i in 1:length(list_tb)){
|
| 113 |
i<-i+1
|
|
| 146 |
i <- i+1
|
|
| 114 | 147 |
tb <-list_tb[[i]] |
| 115 | 148 |
plot_name <- names(list_tb)[i] |
| 116 | 149 |
pat_str <- "tb_m" |
| ... | ... | |
| 139 | 172 |
main=paste("rmse ",plot_name,sep=" "),
|
| 140 | 173 |
pch=1:length(avg_tb$pred_mod), |
| 141 | 174 |
par.settings=list(superpose.symbol = list( |
| 142 |
pch=1:length(avg_tb$pred_mod))),
|
|
| 175 |
pch=1:length(avg_tb$pred_mod))), |
|
| 143 | 176 |
auto.key=list(columns=5)) |
| 144 | 177 |
|
| 145 | 178 |
dev.off() |
| ... | ... | |
| 162 | 195 |
return(tb_diff) |
| 163 | 196 |
} |
| 164 | 197 |
|
| 198 |
|
|
| 165 | 199 |
metric_names <- c("mae","rmse","me","r")
|
| 166 | 200 |
diff_kriging_CAI <- diff_df(tb_s_kriging_CAI,tb_v_kriging_CAI,metric_names) |
| 167 | 201 |
diff_gam_CAI <- diff_df(tb_s_gam_CAI,tb_v_gam_CAI,metric_names) |
| 168 | 202 |
|
| 169 |
boxplot(diff_kriging_CAI$rmse,diff_gam_CAI$rmse) |
|
| 203 |
boxplot(diff_kriging_CAI$rmse,diff_gam_CAI$rmse,names=c("kriging_CAI","gam_CAI"),
|
|
| 204 |
main="Difference between training and testing daily rmse") |
|
| 205 |
|
|
| 206 |
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) |
|
| 209 |
|
|
| 210 |
boxplot(diff_kriging_m_CAI$rmse,diff_gam_m_CAI$rmse,names=c("kriging_CAI","gam_CAI"),
|
|
| 211 |
main="Difference between training and monhtly testing rmse") |
|
| 212 |
|
|
| 213 |
### For fusion |
|
| 214 |
|
|
| 215 |
metric_names <- c("mae","rmse","me","r")
|
|
| 216 |
diff_kriging_fus <- diff_df(tb_s_kriging_fus,tb_v_kriging_fus,metric_names) |
|
| 217 |
diff_gam_fus <- diff_df(tb_s_gam_fus,tb_v_gam_fus,metric_names) |
|
| 218 |
|
|
| 219 |
boxplot(diff_kriging_fus$rmse,diff_gam_fus$rmse,names=c("kriging_fus","gam_fus"),
|
|
| 220 |
main="Difference between training and testing daily rmse") |
|
| 221 |
|
|
| 222 |
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) |
|
| 225 |
|
|
| 226 |
boxplot(diff_kriging_m_fus$rmse,diff_gam_m_fus$rmse,names=c("kriging_fus","gam_fus"),
|
|
| 227 |
main="Difference between training and testing FUS rmse") |
|
| 228 |
|
|
| 229 |
### NOW PLOT OF COMPARISON BETWEEN Kriging and GAM |
|
| 230 |
|
|
| 231 |
xyplot(as.formula(plot_formula),group=pred_mod,type="b", |
|
| 232 |
data=avg_tb, |
|
| 233 |
main=paste("rmse ",plot_name,sep=" "),
|
|
| 234 |
pch=1:length(avg_tb$pred_mod), |
|
| 235 |
par.settings=list(superpose.symbol = list( |
|
| 236 |
pch=1:length(avg_tb$pred_mod))), |
|
| 237 |
auto.key=list(columns=5)) |
|
| 170 | 238 |
|
Also available in: Unified diff
analyses paper methods part 5, monthly hold out (0 to 70%) for gam_fus and kriging_fus