n_va$Pred <- predict(Max_like, newdata=cbind(n_va,g="a"))
# Calcul of RMSE, Bias, MAPE% each time:
Bias[i] = 100*mean((n_va$D - n_va$Pred)/n_va$D)
RMSE[i] = 100*sqrt(mean(((n_va$D - n_va$Pred)/n_va$D)^2)) MAPE[i] = 100*mean(abs(n_va$D - n_va$Pred)/n_va$D)
}
i
# Errors: mean(Bias) mean(RMSE) mean(MAPE)
# The end
Có thể bạn quan tâm!
- Thiết lập và thẩm định chéo hệ thống mô hình ước tính sinh khối trên mặt đất cây rừng khộp ở Việt Nam - 19
- Thiết lập và thẩm định chéo hệ thống mô hình ước tính sinh khối trên mặt đất cây rừng khộp ở Việt Nam - 20
- Codes Thực Hiện Trong R Để Thẩm Định Chéo Theo Phương Pháp Loocv Mô Hình: Agb = A×D 2 Hwd B
- Thiết lập và thẩm định chéo hệ thống mô hình ước tính sinh khối trên mặt đất cây rừng khộp ở Việt Nam - 23
- Thiết lập và thẩm định chéo hệ thống mô hình ước tính sinh khối trên mặt đất cây rừng khộp ở Việt Nam - 24
- Thiết lập và thẩm định chéo hệ thống mô hình ước tính sinh khối trên mặt đất cây rừng khộp ở Việt Nam - 25
Xem toàn bộ 207 trang tài liệu này.
Phụ lục 4: Chương trình codes để thiết lập và thẩm định chéo K-Fold đồng thời hệ thống mô hình sinh khối cây rừng theo phương pháp SUR chạy trong phần mềm SAS
/* Seemingly Unrelated Regression for Model System of AGB = Bst + Bbr + Ble
+ Bba */ Using K-Fold cross validation
dm'log;clear;output;clear';
options pageno=1 nodate nocenter;
/* Create a path to the folder to get input file and write the output to */
%let infile_path = C:UsersbaohuOneDrive1 - Article Dip ForestSAS Scripts for SUR;
%let outfile_path = C:UsersbaohuOneDrive1 - Article Dip ForestSAS Scripts for SUR;
/* Ensure no dataset exists globaly with name we will append TBwts results to */ proc datasets library = work nodetails nolist;
delete B_AGB; run;
quit;
/* Code to load 10 random datasets for DF and EBLF -- datasets generated from script Generate_Random_200_datasets */
data B;
infile "&infile_path.Mix2.csv" dlm=',' firstobs=2 dsd truncover;
input Selected Plot_ID$ ID$ Tree_ID$ FT$ D H WD Bst Bbr Ble Bba AGB D2H D2HWD simNum;
keep Selected D H WD Bst Bbr Ble Bba AGB D2H D2HWD simNum; run;
/* Repeat model fitting and validation statistic calculation 200 times */
%macro runit_var1(eq_num=);
/* Step 1: Subset testing (validation) and training data (model fitting) */ data B_train; set B;
if Selected = 1;
run;
data B_test; set B; if Selected = 0; run;
/* Step 2: Fit sur weighted nls models:
/* Should change Weight variable) */
/* Monitor in case it can not run 200 times by simNum in Log window > Run agian with different Weight. Finally record the simNum */
/* Parms based on AGB, BGB model fitting results outside of SUR (coded in R)*/
proc model data=B_train NOPRINT; by simNum;
parms a1=0.032 b11=2.305 b12=0.506 b13=0.796 a2=0.007 b21=2.846
a3=0.019 b31=1.884 a4=0.039 b41=2.216;
Bst = a1*(D**b11)*(H**b12)*(WD**b13); h.Bst = 1/D;
Bbr = a2*(D**b21); h.Bbr = 1/D**-1; Ble = a3*(D**b31); h.Ble = 1/D**2; Bba = a4*(D**b41); h.Bba = 1/D**-2;
AGB = a1*(D**b11)*(H**b12)*(WD**b13) + a2*(D**b21) + a3*(D**b31) + a4*(D**b41);
h.AGB = 1/D**-1;
fit Bst Bbr Ble Bba AGB / sur outest = B_est
/* look at help for PROC MODEL to see how to fetch R2 and AIC from this line
*/;
run;
/* Step 3.1: Calculate summary statistics for each for each estimate */ data B_Pred; merge B_test B_est;
by simNum;
pBst = a1*(D**b11)*(H**b12)*(WD**b13); pBbr = a2*(D**b21);
pBle = a3*(D**b31); pBba = a4*(D**b41);
pAGB = a1*(D**b11)*(H**b12)*(WD**b13) + a2*(D**b21)+ a3*(D**b31) + a4*(D**b41);
RelDif_Bst = (Bst - pBst)/Bst; RelDif_Bbr = (Bbr - pBbr)/Bbr; RelDif_Ble = (Ble - pBle)/Ble; RelDif_Bba = (Bba - pBba)/Bba; RelDif_AGB = (AGB - pAGB)/AGB;
AbsRelDif_Bst = abs(Bst - pBst)/Bst; AbsRelDif_Bbr = abs(Bbr - pBbr)/Bbr; AbsRelDif_Ble = abs(Ble - pBle)/Ble; AbsRelDif_Bba = abs(Bba - pBba)/Bba; AbsRelDif_AGB = abs(AGB - pAGB)/AGB;
RelDif2_Bst = RelDif_Bst**2; RelDif2_Bbr = RelDif_Bbr**2; RelDif2_Ble = RelDif_Ble**2; RelDif2_Bba = RelDif_Bba**2; RelDif2_AGB = RelDif_AGB**2;
run;
/* Step 3.2: aggregate summary statistics to a per simulation level */ proc summary data = B_Pred mean NOPRINT;
by simNum;
var RelDif_Bst RelDif_Bbr RelDif_Ble RelDif_Bba RelDif_AGB AbsRelDif_Bst AbsRelDif_Bbr AbsRelDif_Ble AbsRelDif_Bba
AbsRelDif_AGB
RelDif2_Bst RelDif2_Bbr RelDif2_Ble RelDif2_Bba RelDif2_AGB; output out = B_Results
run;
mean = pbias_Bst pbias_Bbr pbias_Ble pbias_Bba pbias_AGB mape_Bst mape_Bbr mape_Ble mape_Bba mape_AGB mspe_Bst mspe_Bbr mspe_Ble mspe_Bba mspe_AGB;
/* Step 3.3: get rmspe -- couldn't figure out how to get it in the proc summary/ proc means step */
data B_Results; set B_Results; by simNum;
rmspe_Bst = sqrt(mspe_Bst); rmspe_Bbr = sqrt(mspe_Bbr); rmspe_Ble = sqrt(mspe_Ble); rmspe_Bba = sqrt(mspe_Bba); rmspe_AGB = sqrt(mspe_AGB); run;
/* Step 3.4: take mean of Pbias/MAPE/RMSPE across 10 simulations */ proc summary data = B_Results mean NOPRINT;
var pbias_Bst pbias_Bbr pbias_Ble pbias_Bba pbias_AGB mape_Bst mape_Bbr mape_Ble mape_Bba mape_AGB rmspe_Bst rmspe_Bbr rmspe_Ble rmspe_Bba rmspe_AGB;
output out = B_Results
mean = pbias_Bst pbias_Bbr pbias_Ble pbias_Bba pbias_AGB mape_Bst mape_Bbr mape_Ble mape_Bba mape_AGB rmspe_Bst rmspe_Bbr rmspe_Ble rmspe_Bba rmspe_AGB
nmiss(pbias_AGB) = N_convg;
run;
/* Step 3.5: multiply results by 100 for percentage */ data B_Results; set B_Results;
pbias_Bst = pbias_Bst *100; mape_Bst = mape_Bst *100; rmspe_Bst = rmspe_Bst *100;
pbias_Bbr = pbias_Bbr *100; mape_Bbr = mape_Bbr *100; rmspe_Bbr = rmspe_Bbr *100;
pbias_Ble = pbias_Ble *100; mape_Ble = mape_Ble *100; rmspe_Ble = rmspe_Ble *100;
pbias_Bba = pbias_Bba *100; mape_Bba = mape_Bba *100; rmspe_Bba = rmspe_Bba *100;
pbias_AGB = pbias_AGB *100; mape_AGB = mape_AGB *100; rmspe_AGB = rmspe_AGB *100;
run;
/* Step 4: create name column in dataset */ data B_Results; set B_Results;
fname ="&eq_num.";
run;
/* Step 5: Export results as .csv file */
proc append base= B_AGB data= B_Results FORCE; run;
%mend runit_var1;
%runit_var1(eq_num=e1);
/* AGB: Export results as .csv file*/ proc export data = B_AGB
outfile ="&outfile_path.B_AGB.csv";
run;
************************************/
/* SUR: Estimate all parameters using the entire dataset*/ dm'log;clear;output;clear';
options pageno=1 nodate nocenter;
/* Create a path to the folder where *.csv data is stored */
%let infile_path = C:UsersbaohuOneDrive1 - Article Dip ForestSAS Scripts for SUR;
data AllData;
infile "&infile_path.Mixed.csv" dlm=',' firstobs=2 dsd truncover;
input Plot_ID$ ID$ Tree_ID$ FT$ D H WD Bst Bbr Ble Bba AGB D2H D2HWD;
run;
proc model data=AllData;
parms a1=0.032 b11=2.305 b12=0.506 b13=0.796 a2=0.007 b21=2.846
a3=0.019 b31=1.884 a4=0.039 b41=2.216;
Bst = a1*(D**b11)*(H**b12)*(WD**b13); h.Bst = 1/D;
Bbr = a2*(D**b21); h.Bbr = 1/D**-1; Ble = a3*(D**b31); h.Ble = 1/D**2; Bba = a4*(D**b41); h.Bba = 1/D**-2;
AGB = a1*(D**b11)*(H**b12)*(WD**b13) + a2*(D**b21) + a3*(D**b31) + a4*(D**b41);
h.AGB = 1/D**-1;
fit Bst Bbr Ble Bba AGB / sur outest = AllData_est; run;
quit;
proc print data = AllData_est; run;
*/ The end */
Phụ lục 5: Codes để thiết lập đồ thị Bland- Altman so sánh sai lệch ước tính AGB của hai mô hình thiết lập theo SUR và độc lập với độ tin cậy 95%
# Start rm(list=ls())
# Clean dev.off()
# Install of Paskages library(ggplot2) library(cowplot)
# Dỉrectory path
setwd("C:/Users/baohu/OneDrive/PhD Master/PhD Tinh/1. Huong dan Tinh sau cung/Data")
# Enter data t
t <- read.table("tAll.txt", header=T,sep="t",stringsAsFactors = FALSE)
# Prediction of AGB from 2 models SUR and Non-SUR:
t$AGBnonsur = 0.0801995*t$D^2.4584050*t$H^0.1576263*t$WD^0.6434186 t$AGBsur = 0.02055*t$D^2.35241*t$H^0.59142 *t$WD^0.69609 + 0.00669
*t$D^2.85742 + 0.03701 * t$D^1.68095 + 0.01541 * t$D^2.43959
# Bland Altman Plot: library(BlandAltmanLeh)
p1 <- bland.altman.plot(t$AGBsur,t$AGBnonsur, graph.sys ="ggplot2") p1<- p1+ geom_point(cex=3)
p1 <- p1 + xlab("Trung bình ước tính AGB theo SUR và độc lập (kg)") + ylab("Sai lệch (kg)") + theme_bw()
p1 = p1 + theme(axis.title.y = element_text(size = rel(1.7))) p1 = p1 + theme(axis.title.x = element_text(size = rel(1.7))) p1 = p1 + theme(axis.text.x = element_text(size=15))
p1 = p1 + theme(axis.text.y = element_text(size=15)) p1 = p1 + xlim(0,1500)
p1 = p1 + ylim(-200,200) p1
# The end