############################################## ## ## ## Script for Assessment of Model Agreement ## ## Model Evaluation System 3.1.8 ## ## Copyright (c) 2010 Luis Orlindo Tedeschi ## ## July 3, 2019 ## ## ## ############################################## # Reference: Lin, L. I.-K., A. S. Hedayat, B. Sinha, and M. Yang. 2002. # Statistical methods in assessing agreement: Models, issues, # and tools. Journal of the American Statistical Association. # 97(457):257-270. # Clear the memory rm(list=ls(all=TRUE)) #-------------------# # Functions library # #-------------------# packageexist<-function(a) { packs<-.packages(all.available=TRUE) for (i in 1:length(packs)) { if (packs[i]==a) {return(TRUE)} } return(FALSE) } check.necessary.packages <- function(necessary, updatepackages = FALSE, askrepos = FALSE, selectpackages = FALSE, removetempdir = TRUE) { Unix <- FALSE if (.Platform$OS.type == "unix") { Unix <- TRUE askrepos <- FALSE } if (Unix == FALSE) library(tcltk) ok <- TRUE if (file.exists(Sys.getenv("R_LIBS_USER")) == FALSE) dir.create(Sys.getenv("R_LIBS_USER"), recursive = TRUE) # .libPaths(if(.Machine$sizeof.pointer == 8) '~/Rlibs/x64' else '~/Rlibs') .libPaths(c(Sys.getenv("R_LIBS_USER"))) installed <- necessary %in% installed.packages()[, "Package"] destdir<-NA if (length(necessary[!installed]) >= 1) { destdir<-paste("C:\\temp\\",basename(tempdir()),"\\",sep="") if (dir.exists(destdir) == FALSE) dir.create(destdir, recursive = TRUE) if (askrepos == T) { m <- getCRANmirrors(all = TRUE, local.only = FALSE) res <- tkchoose2("CRAN Mirror", "Choose your country to select the CRAN mirrors.\n'NONE' will use pre-selected CRAN mirrors", unique(m$Country)) if (res == -1) { m2 <- m } else { if (res > 0L) { m2 <- subset(m, m$Country == res) } else { m2 <- data.frame(Name = c("0-cloud1", "0-cloud2", "Imperial", "Los Angeles", "Oxford", "Omegahat", "Berkeley", "Melbourne", "Piracicaba"), URL = c("https://cloud.r-project.org/", "http://cloud.r-project.org/", "http://cran.ma.imperial.ac.uk", "http://cran.stat.ucla.edu", "http://www.stats.ox.ac.uk/pub/RWin", "http://www.omegahat.org/R", "http://cran.cnr.Berkeley.edu", "http://cran.ms.unimelb.edu.au", "http://brieger.esalq.usp.br/CRAN"), stringsAsFactors = F) } } } else { m2 <- data.frame(Name = c("0-cloud1", "0-cloud2", "Imperial", "Los Angeles", "Oxford", "Omegahat", "Berkeley", "Melbourne", "Piracicaba"), URL = c("https://cloud.r-project.org/", "http://cloud.r-project.org/", "http://cran.ma.imperial.ac.uk", "http://cran.stat.ucla.edu", "http://www.stats.ox.ac.uk/pub/RWin", "http://www.omegahat.org/R", "http://cran.cnr.Berkeley.edu", "http://cran.ms.unimelb.edu.au", "http://brieger.esalq.usp.br/CRAN"), stringsAsFactors = F) } local({ for (i in 1:nrow(m2)) { if (length(necessary[!installed]) >= 1) { options(repos = NULL) mirror <- m2$Name[i] URL <- m2$URL[i] r <- getOption("repos") r[mirror] <- gsub("/$", "", URL) options(repos = r) contriburl <- contrib.url(getOption("repos"), type = getOption("pkgType")) # con<-url(paste(contriburl,'/_Info.txt',sep='')) con <- url(contriburl) f <- try(open(con, "r"), silent = TRUE) close(con) if (!inherits(f, "try-error")) { cat(paste("Trying '", contriburl, "'\n", sep = "")) install.packages(necessary[!installed],lib = Sys.getenv("R_LIBS_USER"), destdir=destdir, Ncpus = getOption("Ncpus", 2), dependencies = TRUE, INSTALL_opts = c("--no-lock")) installed <- necessary %in% installed.packages()[, "Package"] } } } }) installed <- necessary %in% installed.packages()[, "Package"] if (length(necessary[!installed]) >= 1) { cat(paste("FAILED to install package(s): '", necessary[!installed], "\n", sep = "")) ok <- FALSE } } if (updatepackages == TRUE) { if (is.na(destdir)==TRUE) { destdir<-paste("C:\\temp\\",basename(tempdir()),"\\",sep="") if (dir.exists(destdir) == FALSE) dir.create(destdir, recursive = TRUE) } if (selectpackages == TRUE) { update.packages(lib = Sys.getenv("R_LIBS_USER"), ask = "graphics") } else { #packages.toupdate<-old.packages(lib = Sys.getenv("R_LIBS_USER")) #packages.toupdate.count<-length(packages.toupdate) #if (packages.toupdate.count>=5) #{ # for (i in 1:packages.toupdate.count) # { # install.packages(pkgs=packages.toupdate[i],lib = Sys.getenv("R_LIBS_USER"), # destdir=destdir, Ncpus = getOption("Ncpus", 2), # dependencies = TRUE, INSTALL_opts = c("--no-lock")) # } #} else #{ update.packages(lib = Sys.getenv("R_LIBS_USER"), ask = FALSE, destdir=destdir, dependencies = TRUE, INSTALL_opts = c("--no-lock")) #} } } if (removetempdir == TRUE) unlink(destdir, recursive=TRUE) return(ok) } #------# # Code # #------# # Information for MES. # Use different values to separate customized R files filesel<-"SC4" # 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 (packageexist("svDialogs")==FALSE) { if (file.exists("SciViews_0.9-5.zip")==TRUE) { install.packages("SciViews_0.9-5.zip", repos=NULL) } else { install.packages("svDialogs", repos = "http://cran.r-project.org") if (packageexist("svDialogs")==FALSE) { cat("\n\n ***** Package svDialogs not found. Modeling halted *****\n\n") quit(status=1) } } } check.necessary.packages("svDialogs", FALSE) library(svDialogs) if (packageVersion("svDialogs")<="0.9-43") { # target can be fixed or random target<-guiDlgInput(message = "Are the factors 'fixed' or 'random'?", title = "Factor type (1/6)", default="fixed") # error can be const or prop error<-guiDlgInput(message = "Are the errors 'const' (constant) or 'prop' (proportional)?", title = "Error type (2/6)", default="const") # least acceptable concordance correlation coefficient (CCC) value la_ccc<-as.numeric(guiDlgInput(message = "What is the least acceptable value for Concordance Correlation Coefficient (CCC)?", title = "Least Acceptable CCC (3/6)", default="0.95")) # least acceptable total deviation index (TDI) value la_tdi<-as.numeric(guiDlgInput(message = "What is the least acceptable value for Total Deviation Index (TDI)?", title = "Least Acceptable TDI (4/6)", default="1")) # least acceptable coverage probability (CP) value la_cp<-as.numeric(guiDlgInput(message = "What is the least acceptable value for Coverage Probability (CP)?", title = "Least Acceptable CP (5/6)", default="0.8")) # confidence interval p-value conf<-as.numeric(guiDlgInput(message = "What is the Confidence Interval (CI) P-value?", title = "CI P-value (6/6)", default="0.95")) } else { # target can be fixed or random target<-dlgInput(message = "Are the factors 'fixed' or 'random'?", title = "Factor type (1/6)", default="fixed")$res # error can be const or prop error<-dlgInput(message = "Are the errors 'const' (constant) or 'prop' (proportional)?", title = "Error type (2/6)", default="const")$res # least acceptable concordance correlation coefficient (CCC) value la_ccc<-as.numeric(dlgInput(message = "What is the least acceptable value for Concordance Correlation Coefficient (CCC)?", title = "Least Acceptable CCC (3/6)", default="0.95")$res) # least acceptable total deviation index (TDI) value la_tdi<-as.numeric(dlgInput(message = "What is the least acceptable value for Total Deviation Index (TDI)?", title = "Least Acceptable TDI (4/6)", default="1")$res) # least acceptable coverage probability (CP) value la_cp<-as.numeric(dlgInput(message = "What is the least acceptable value for Coverage Probability (CP)?", title = "Least Acceptable CP (5/6)", default="0.8")$res) # confidence interval p-value conf<-as.numeric(dlgInput(message = "What is the Confidence Interval (CI) P-value?", title = "CI P-value (6/6)", default="0.95")$res) } # 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; y<-data1$V3 if(error=='prop'){ y<-log(y) x<-log(x) } d<-y-x m1<-mean(y) m2<-mean(x) yyy<-mean(d) v1<-var(y) v2<-var(x) yy<-var(d) e2<-sum(d^2) n<-length(d) mu_d<-m1-m2 d2<-mu_d^2 e2<-e2/n sd2<-(e2-d2)*n/(n-3) s12<-v1*(n-1)/n s22<-v2*(n-1)/n u<-mu_d/sqrt(sqrt(s12*s22)) v<-sqrt(s12/s22) ca<-2/(v+1/v+u^2) rc<-1-e2/(d2+s12+s22) r<-(rc/ca) d2_s2<-(d2/sd2) b1<-r*v b0<-m1-b1*m2 se2<-s12*(1-r^2)*n/(n-3) e2<-e2*n/(n-1) kp<-la_tdi if (error=='prop'){ kp<-log(la_TDI/100+1) } # For CP when the target value is fixed if (target=='fixed'){ di<-b0+(b1-1)*x d2i<-di^2 cpi<-pchisq(kp^2/se2,1,d2i/se2) se<-sqrt(se2) kpm<-(kp+di)/se kmm<-(kp-di)/se c0<-dnorm(-kpm)-dnorm(kmm) c1<-c0*x c2<--kpm*dnorm(-kpm)-kmm*dnorm(kmm) cp<-mean(cpi) c0<-mean(c0) c1<-mean(c1) c2<-mean(c2) s_p<-c0^2+(c0*m2-c1)^2/s22+c2^2/2 s_t<-sqrt(s_p/((n-3)*cp^2*(1-cp)^2)) t<-log(cp/(1-cp)) t_l<-t-qnorm(conf)*s_t cp_l<-exp(t_l)/(1+exp(t_l)) } # For CCC, accuracy and precision l_ca<-log(ca/(1-ca)) z_rc<-log((1+rc)/(1-rc))/2 z_r<-log((1+r)/(1-r))/2 # For MSD & approx TDI & CP w<-log(e2) k<-qnorm(1-(1-la_cp)/2)*sqrt(e2) cp_a<-pchisq(kp^2/e2,1) # Standard errors of CCC, accuracy, precision, and MSD if (target=='random'){ sz_rc<-((1-r^2)*rc^2/(1-rc^2)/r^2 + 2*rc^3*(1-rc)*u^2/r/(1-rc^2)^2-rc^4*u^4/r^2/(1-rc^2)^2/2)/(n-2) sz_r<-1/(n-3) sl_ca<-((ca^2*u^2*(v+1/v-2*r) + ca^2*(v^2+1/v^2+2*r^2)/2 + (1+r^2)*(ca*u^2-1))/((n-2)*(1-ca)^2)) s_w<-2*(1-d2^2/e2^2)/(n-2) } if (target=='fixed'){ sz_rc<-(1-r^2)*rc^2/((n-2)*(1-rc^2)^2*r^2)*(v*u^2*rc^2+(1-rc*r*v)^2 + v^2*rc^2*(1-r^2)/2) sz_r<-(1-r^2/2)/(n-3) sl_ca<-(u^2*v*ca^2*(1-r^2) + (1-v*ca)^2*(1-r^4)/2)/((n-2)*(1-ca)^2) s_w<-2*(1-(d2+s22*(1-b1)^2)^2/e2^2)/(n-2) } sz_rc<-sqrt(sz_rc) sz_r<-sqrt(sz_r) sl_ca<-sqrt(sl_ca) s_w<-sqrt(s_w) z_rc_l<-z_rc-qnorm(conf)*sz_rc z_r_l<-z_r-qnorm(conf)*sz_r l_ca_l<-l_ca-qnorm(conf)*sl_ca rc_l<-tanh(z_rc_l) r_l<-tanh(z_r_l) ca_l<-exp(l_ca_l)/(1+exp(l_ca_l)) w_u<-w+qnorm(conf)*s_w e2_u<-exp(w_u) k_u<-qnorm(1-(1-la_cp)/2)*sqrt(e2_u) cp_a_l<-pchisq(kp^2/e2_u,1) if (error=='prop'){ k<-100*(exp(k)-1) k_u<-100*(exp(k_u)-1) } # For CP when the target value is random if (target=='random'){ cp<-pchisq(kp^2/sd2,1,d2_s2) sd<-sqrt(sd2) kpm<-(kp+mu_d)/sd kmm<-(kp-mu_d)/sd s_t<-((dnorm(-kpm)-dnorm(kmm))^2 + (kmm*dnorm(kmm)+kpm*dnorm(-kpm))^2/2)/((n-3)*cp^2*(1-cp)^2) s_t<-sqrt(s_t) t<-log(cp/(1-cp)) t_l<-t-qnorm(conf)*s_t cp_l<-exp(t_l)/(1+exp(t_l)) } cat(" ********************************************* * Assessment of Model Agreement * * * * TDI = Total Deviation Index * * CP = Coverage Probability * * CCC = Concordance Correlation Coefficient * ********************************************* Parameters: ----------- . Factor type : ", target," . Error type : ", error," . Least acceptable CCC : ", la_ccc," . Least acceptable TDI : ", la_tdi," . Leat acceptable CP : ", la_cp," . Confidence interval : ", conf," Concordance correlation coefficient (rc) : ", format(rc, digits=4, width=10, justify='right')," Lower CCC confidence interval : ", format(rc_l,digits=4,width=10,justify='right')," (This value has to be GREATER than ",la_ccc,") Precision (r) : ", format(r, digits=4, width=10, justify='right')," Lower r confidence interval : ", format(r_l,digits=4,width=10,justify='right')," Accuracy (Cb) : ", format(ca, digits=4, width=10, justify='right')," Lower Cb confidence interval : ", format(ca_l,digits=4,width=10,justify='right')," Mean shift (u) : ", format(u, digits=4, width=10, justify='right')," Variance shift (v) : ", format(v, digits=4, width=10, justify='right')," Mean Square Deviance (MSD) : ", format(e2, digits=4, width=10, justify='right')," Upper MSD confidence interval : ", format(e2_u,digits=4,width=10,justify='right')," Total Deviation Index (TDI) : ", format(k, digits=4, width=10, justify='right')," (Assuming Coverage Protection (CP) of ",la_cp,") Upper TDI confidence interval : ", format(k_u,digits=4,width=10,justify='right')," (This value has to be LOWER than ", la_tdi,") Coverage Protection (CP) : ", format(cp, digits=4, width=10, justify='right')," Lower CP confidence interval : ", format(cp_l,digits=4,width=10,justify='right')," (This value has to be GREATER than ", la_cp,") Relative Bias Square (RBS) : ", format(d2_s2, digits=4, width=10, justify='right')," Attention: The RBS must be less than 1 or 8 for a least acceptable coverage probability of 0.9 and 0.8, respectively, in order for the TDI to be valid \n") cat(" Reference: Lin, L. I.-K., A. S. Hedayat, B. Sinha, and M. Yang. 2002. Statistical methods in assessing agreement: Models, issues, and tools. Journal of the American Statistical Association. 97(457):257-270. \n") sink()