############################################### ## ## ## Script for OLS regression and Diagnostics ## ## Model Evaluation System 3.0 ## ## Copyright (c) 2019 Luis Orlindo Tedeschi ## ## July 03, 2019 ## ## ## ############################################### # Clear the memory rm(list=ls(all=TRUE)) #-------------------# # Functions library # #-------------------# #------# # Code # #------# # Information for MES. # Use different values to separate customized R files filesel<-"SC6" graphwidth<-480 graphheight<-480 # Assign file names args<-commandArgs() fn1<-substring(args[4],8) fn2<-substring(fn1,1,nchar(fn1)-2) # Remove the .R extension fndata<-paste(fn2,".csv",sep="") fntext1<-paste(fn2,"_",filesel,"_text1_Script_Output.txt",sep="") if (args[6]=="JPG") { graphext<-"jpg" } else {if (args[6]=="PNG") { graphext<-"png" } else { graphext<-"bmp" } } fngraph1<-paste(fn2,"_",filesel,"_graph1_Scatter_Plot.",graphext,sep="") fngraph2<-paste(fn2,"_",filesel,"_graph2_Residual_vs_Predicted_Plot.",graphext,sep="") fngraph3<-paste(fn2,"_",filesel,"_graph3_Residual_vs_MinMaxPredicted_Plots.",graphext,sep="") fngraph4<-paste(fn2,"_",filesel,"_graph4_Residual_vs_StandardizedPredicted_Plots.",graphext,sep="") fngraph5<-paste(fn2,"_",filesel,"_graph5_Residual_vs_CenteredPredicted_Plots.",graphext,sep="") # Add the output text file here sink(fntext1) # Read the database # Add the input datafile here data1<-read.csv(fndata,header=F) # Assign variables to common used variables x<-data1$V1; sdx<-data1$V2; y<-data1$V3; sdy<-data1$V4; rankx<-data1$V5; ranky<-data1$V6; exclude<-data1$V7 # Linear LS regression Y on X lsregout<-lm(y~x) # Plot linear regression and residuals # Use: x11() for screen output # Use: bmp() to save the graph cat(" *********************** * Scatter Plot Y on X * *********************** \n") if (args[6]=="JPG") { jpeg(filename=fngraph1,bg="white",width=2*graphwidth,height=graphheight) } else {if (args[6]=="PNG") { png(filename=fngraph1,bg="white",width=2*graphwidth,height=graphheight) } else { bmp(filename=fngraph1,bg="white",width=2*graphwidth,height=graphheight) } } layout(matrix(1:2,1,2)) #plot(x,y,xlim=range(min(min(x),min(y)):max(max(x),max(y))),ylim=range(min(min(x),min(y)):max(max(x),max(y))),xlab="Model-Predicted",ylab="Observed",type='p',main="Linear Regression: Y vs. X") plot(x,y,xlab="Model-Predicted",ylab="Observed",type='p',main="Linear Regression: Y vs. X") abline(0,1, lty=1,col=1) abline(lsregout,lty=2,col=2) legend("bottomright",legend=c('Y=X','Linear'),lty=1:2,col=1:2) plot(y,x,ylab="Model-Predicted",xlab="Observed",type='p',main="Linear Regression: X vs. Y") abline(0,1, lty=1,col=1) legend("bottomright",legend=c('Y=X'),lty=1,col=1) #text(1,1,as.expression(substitute(italic(R)^2==r,list(r=round(summary(lsregout)$r.squared,3))))) # LS regression cat(" ****************************************** * Ordinary Least-Square (OLS) Regression * ****************************************** \n") lsregout summary(lsregout) # Plot residuals minmaxvalues<-(lsregout$fitted.values-min(lsregout$fitted.values))/(max(lsregout$fitted.values)-min(lsregout$fitted.values)) zvalues<-scale(lsregout$fitted.values) centervalues<-lsregout$fitted.values-mean(lsregout$fitted.values) cat(" ***************** * Residual Plot * ***************** \n") if (args[6]=="JPG") { jpeg(filename=fngraph2,bg="white",width=graphwidth,height=graphheight) } else {if (args[6]=="PNG") { png(filename=fngraph2,bg="white",width=graphwidth,height=graphheight) } else { bmp(filename=fngraph2,bg="white",width=graphwidth,height=graphheight) } } plot(lsregout$fitted.values,lsregout$residuals,xlab="Regression-Predicted",ylab="Residuals",type='p',main="Residuals (Y) vs Regression-Predicted (X)") abline(h=0,lty=1,col=1) abline(lm(lsregout$residuals~lsregout$fitted.values),lty=2,col=2) legend("bottomright",legend=c('Y=0','Linear'),lty=1:2,col=1:2) dev.off() cat(" ************************************* * Min-Max Normalization of Residues * ************************************* \n") if (args[6]=="JPG") { jpeg(filename=fngraph3,bg="white",width=graphwidth,height=graphheight) } else {if (args[6]=="PNG") { png(filename=fngraph3,bg="white",width=graphwidth,height=graphheight) } else { bmp(filename=fngraph3,bg="white",width=graphwidth,height=graphheight) } } plot(minmaxvalues,lsregout$residuals,xlab="Normalized Regression-Predicted",ylab="Residuals",type='p',main="Residuals (Y) vs Normalized Regression-Predicted (X)") abline(h=0,lty=1,col=1) abline(lm(lsregout$residuals~minmaxvalues),lty=2,col=2) legend("bottomright",legend=c('Y=0','Linear'),lty=1:2,col=1:2) dev.off() cat(" ************************* * Standardized Residues * ************************* \n") if (args[6]=="JPG") { jpeg(filename=fngraph4,bg="white",width=graphwidth,height=graphheight) } else {if (args[6]=="PNG") { png(filename=fngraph4,bg="white",width=graphwidth,height=graphheight) } else { bmp(filename=fngraph4,bg="white",width=graphwidth,height=graphheight) } } plot(zvalues,lsregout$residuals,xlab="Standardized Regression-Predicted",ylab="Residuals",type='p',main="Residuals (Y) vs Standardized Regression-Predicted (X)") abline(h=0,lty=1,col=1) abline(lm(lsregout$residuals~zvalues),lty=2,col=2) legend("bottomright",legend=c('Y=0','Linear'),lty=1:2,col=1:2) dev.off() cat(" ********************* * Centered Residues * ********************* \n") if (args[6]=="JPG") { jpeg(filename=fngraph5,bg="white",width=graphwidth,height=graphheight) } else {if (args[6]=="PNG") { png(filename=fngraph5,bg="white",width=graphwidth,height=graphheight) } else { bmp(filename=fngraph5,bg="white",width=graphwidth,height=graphheight) } } plot(centervalues,lsregout$residuals,xlab="Centered Regression-Predicted",ylab="Residuals",type='p',main="Residuals (Y) vs Centered Regression-Predicted (X)") abline(h=0,lty=1,col=1) abline(lm(lsregout$residuals~zvalues),lty=2,col=2) legend("bottomright",legend=c('Y=0','Linear'),lty=1:2,col=1:2) dev.off() cat(" **************** * End of Script* **************** \n") sink()