# Conditional IQR simulations for Goldman and Kaplan (2017), "Nonparametric inference on conditional quantile differences and linear combinations, using L-statistics" # Questions? Comments? kaplandm@missouri.edu # Note: for PARALLEL>1, have to "run lines", not "source" SAVE.FLAG <- TRUE # TRUE=save results to file; FALSE=display in console only. # setwd("") tryCatch(source("quantile_inf_np.R"), error=function(w) { cat("First download the file quantile_inf_np.R currently available from http://faculty.missouri.edu/~kaplandm/code/quantile_inf_np.R and save it in the same directory as this file (or elsewhere in R's path)."); stop(w) } ) PARALLEL <- 4 #determine by parallel::detectCores() or prior knowledge. if (!require("np")) stop("Need package np to run Qu and Yoon (2015).") if (!require("fields")) stop("Need package fields to run Qu and Yoon (2015).") if (!require("quantreg")) stop("Need package quantreg to run Qu and Yoon (2015).") if (!require("MatrixModels")) stop("Need package MatrixModels (for rqss) to run Qu and Yoon (2015).") QY.BC <- -1 #0: no bias correction. 1: usual. -1: conservative, like QY15 paper. if (!(QY.BC %in% -1:1)) stop("QY.BC must be -1, 0, or 1") UNIFb <- sqrt(3) #for Findex in 018,019: Uniform(-UNIFb,UNIFb) set.seed(112358) rqssfn <- function(X) (X*(1-X))^(1/2)*sin(2*pi*(1+2^(-7/5))/(X+2^(-7/5))) constfn <- function(X) ifelse(X==X,0,0) if (SAVE.FLAG) { OUTFILE <- paste0(format(Sys.time(),"%Y_%m_%d"),"_quantile_inf_np_sims_iqr",".txt") NREPLIC <- 1000 BREPLIC <- 1 CASES <- list(list(n=400,ALPHA=0.05,Y.fn=rqssfn,Findex=010,XUNIF=T,x0type=2), list(n=400,ALPHA=0.05,Y.fn=rqssfn,Findex=012,XUNIF=T,x0type=2), list(n=400,ALPHA=0.05,Y.fn=rqssfn,Findex=014,XUNIF=T,x0type=2), list(n=400,ALPHA=0.05,Y.fn=rqssfn,Findex=016,XUNIF=T,x0type=2) ) } else { OUTFILE <- "" NREPLIC <- 20 BREPLIC <- 1 CASES <- list(list(n=400,ALPHA=0.05,Y.fn=constfn,Findex=011,XUNIF=T,x0type=2)) } cat(sprintf("NREPLIC=%d,BREPLIC=%d,PARALLEL=%d", NREPLIC,BREPLIC,PARALLEL), file=OUTFILE,sep="\n",append=TRUE) # PARALLEL setup and error checking if (PARALLEL>1) { if (!require("parallel") || !require("foreach")) {warning("Install package foreach in order to run PARALLEL."); PARALLEL <- 1} if (!require("doParallel")) {warning("Install package doParallel in order to run PARALLEL."); PARALLEL <- 1} PARALLEL <- tryCatch({workers <- makeCluster(PARALLEL); registerDoParallel(workers); on.exit(stopCluster(workers),add=TRUE); PARALLEL}, error=function(Err){warning(sprintf("Error creating %d clusters",PARALLEL)); 1}) clusterSetRNGStream(workers,112358) } else if (PARALLEL<0) { warning("PARALLEL must be a non-negative integer: 0 or 1 for not parallel, positive for number of parallel CPUs to use.") } for (icase in 1:length(CASES)) { p <- 0.25 #do (1-p)-quantile minus p-quantile n <- CASES[[icase]]$n ALPHA <- CASES[[icase]]$ALPHA Y.fn <- CASES[[icase]]$Y.fn Findex <- CASES[[icase]]$Findex XUNIF <- CASES[[icase]]$XUNIF x0type <- CASES[[icase]]$x0type if (x0type==1) { NUM.X0 <- 6 #6,11: really hard (multiple points at local max/min); 7 less so TRIM.X0 <- 0.04 #trim off top and bottom TRIM portion of X0 points ONESIDED <- 0 #0:two-sided; 1/-1:upper/lower one-sided x0.ps <- TRIM.X0 + (1-2*TRIM.X0)*c(0:(NUM.X0-1))/(NUM.X0-1) #X points of interest, as quantiles of X dist'n } else if (x0type==2) { if (!XUNIF) stop("Can only use x0type==2 with XUNIF==TRUE") x0.ps <- c(0.050,0.087,0.125,0.181,0.237,0.324,0.411,0.558,0.706,0.853) #local extrema + midpoints NUM.X0 <- length(x0.ps) } else { stop(sprintf("Not supported: x0type=%g",x0type)) } # CUBEADJ <- min(n.c,n.t)^(-1/5) / min(n.c,n.t)^(-7/19) #for local cubic bandwidth adjustment, multiply L-stat bandwidth by this cat(sprintf("p=%g,1-p=%g,n=%d,ALPHA=%g,Y.fn=%s,Findex=%03d,QY.BC=%d,XUNIF=%d,x0type=%d", p,1-p,n,ALPHA,as.character(as.expression(body(Y.fn))),Findex,QY.BC,XUNIF,x0type), file=OUTFILE,sep="\n",append=TRUE) cat(sprintf("Start time for this case is %s",format(Sys.time(), "%X, %A, %d %b %Y")), sep='\n',file=OUTFILE,append=TRUE) if (XUNIF) { Xs <- matrix(runif(NREPLIC*n,0,1),NREPLIC,n) x0s <- qunif(x0.ps,0,1) } else { Xs <- matrix(rnorm(NREPLIC*n,0,1),NREPLIC,n) x0s <- qnorm(x0.ps,0,1) } cat(sprintf("x0s = %s",paste0(sprintf("%5.3f",x0s),collapse=" & ")), file=OUTFILE,sep='\n',append=TRUE) y0s <- Y.fn(x0s) #no add'l effect of heterosk. since p-quantile of U is zero sig0 <- 0.2 #Homoskedastic if Findex is even, heteroskedastic if odd if (((Findex-10*floor(Findex/10)) %% 2) == 0) { sig.fn <- function(X) rep.int(sig0,length(X)) #homoskedastic } else { sig.fn <- function(X) sig0*(1+X) #heteroskedastic } #Different distribution shapes for 010/011 vs. 012/013 vs. ... p0 <- 0.5 #Center w/ median=0 if (floor((Findex-10*floor(Findex/10))/2)==0) { Us <- matrix(rnorm(NREPLIC*n,0,1)-qnorm(p0,0,1), NREPLIC, n) #Gaussian (010/011) IQR0s <- sig.fn(x0s)*(qnorm(1-p)-qnorm(p)) Fstr <- "$N(0,1)$" } else if (floor((Findex-10*floor(Findex/10))/2)==1) { Us <- matrix(rt(NREPLIC*n,3)-qt(p0,3), NREPLIC, n) IQR0s <- sig.fn(x0s)*(qt(1-p,3)-qt(p,3)) Fstr <- "$t_3$" } else if (floor((Findex-10*floor(Findex/10))/2)==2) { Us <- matrix(rt(NREPLIC*n,1)-qt(p0,1), NREPLIC, n) IQR0s <- sig.fn(x0s)*(qt(1-p,1)-qt(p,1)) Fstr <- "Cauchy" } else if (floor((Findex-10*floor(Findex/10))/2)==3) { Us <- matrix(rchisq(NREPLIC*n,3)-qchisq(p0,3), NREPLIC, n) IQR0s <- sig.fn(x0s)*(qchisq(1-p,3)-qchisq(p,3)) Fstr <- "$\\chi^2_3$" } else if (floor((Findex-10*floor(Findex/10))/2)==4) { #018/9 Us <- matrix(runif(NREPLIC*n,-UNIFb,UNIFb)-qunif(p0,-UNIFb,UNIFb), NREPLIC, n) IQR0s <- sig.fn(x0s)*(qunif(1-p,-UNIFb,UNIFb)-qunif(p,-UNIFb,UNIFb)) Fstr <- "Uniform" } else {stop("Not supported.")} Ys <- Y.fn(Xs) + sig.fn(Xs)*Us tmp <- rep.int(NA,NREPLIC*NUM.X0); dim(tmp) <- c(NREPLIC,NUM.X0) CIs <- list(pt.lo=tmp,pt.hi=tmp,jt.lo=tmp,jt.hi=tmp, QYu.pt.lo=tmp,QYu.pt.hi=tmp,QYg.pt.lo=tmp,QYg.pt.hi=tmp) # if (PARALLEL>1) { cat(sprintf("Time elapsed = %g seconds",system.time({ # DONE ABOVE clusterSetRNGStream(workers,112358) #for bootstrap (nothing else randomized inside dopar) writeLines(c(""),"tmp_iqr.txt") #clear file of prev output CIs.raw <- foreach(irep=1:NREPLIC,.combine=rbind,.inorder=FALSE) %dopar% { #,.export=c("Xs"),.packages=c("")) %dopar% { sink("tmp_iqr.txt",append=TRUE) source("quantile_inf_np.R") X <- Xs[irep,]; Y <- Ys[irep,] tmp <- quantile.inf.np(Y=Y, X=X, x0s=x0s, p=c(1-p,p), PSI=c(1,-1), ALPHA=ALPHA, JOINT.FLAG=T, ONESIDED=0, METHOD.TYPE='lincom', NORM.APPROX.CAL=FALSE, NORM.APPROX.Q=F) Lstat.CIs <- c(tmp[[1]]$CI.pointwise.lows,tmp[[1]]$CI.pointwise.highs, tmp[[1]]$CI.joint.lows,tmp[[1]]$CI.joint.highs) Lstat.hs <- tmp[[1]]$bandwidth.pointwise # Qu and Yoon (2015) QYu.CIs <- QYg.CIs <- rep.int(NA,2*NUM.X0) #uniform or Gaussian kernel #kernel-dependent values: intKsq, mu2K for uniform, gaussian intKsq <- c(1,1/(2*sqrt(pi))); mu2K <- c(1/12,1) if (!require("np")) stop("Need package np to run Qu and Yoon (2015).") if (!require("fields")) stop("Need package fields to run Qu and Yoon (2015).") if (!require("quantreg")) stop("Need package quantreg to run Qu and Yoon (2015).") # estimate objects for bandwidth: fX=c(fX.c,fX.t), fYX=c(fYX.c,fYX.t), Qpp=c(Qpp.c,Qpp.t) # "50" means for conditional median, for use in Cor 1 optimal bandwidth # Otherwise, for use in estimating Thm 2 scale factor g50 <- function(lam,y,x) AIC(rqss(y~qss(x, lambda=lam),tau=1/2),k = -1) #k=-1 does BIC (look @ code for AIC.rqss) Q50 <- tryCatch({ tmplamstar <- optimize(g50, interval = c(0.001, 0.5), x=X, y=Y) tmp <- rqss(Y~qss(X, lambda=tmplamstar$min), tau=1/2) predict.rqss(tmp,newdata=data.frame(X=x0s)) }, error=function(err) rep.int(NA,NUM.X0)) #warning=function(err) rep.int(NA,NUM.X0) # g1 <- function(lam,y,x) AIC(rqss(y~qss(x, lambda=lam),tau=p),k = -1) #k=-1 does BIC (look @ code for AIC.rqss) Q1 <- tryCatch({ tmplamstar <- optimize(g1, interval = c(0.001, 0.5), x=X, y=Y) tmp <- rqss(Y~qss(X, lambda=tmplamstar$min), tau=p) predict.rqss(tmp,newdata=data.frame(X=x0s)) }, error=function(err) rep.int(NA,NUM.X0)) #warning=function(err) rep.int(NA,NUM.X0) g2 <- function(lam,y,x) AIC(rqss(y~qss(x, lambda=lam),tau=1-p),k = -1) #k=-1 does BIC (look @ code for AIC.rqss) Q2 <- tryCatch({ tmplamstar <- optimize(g2, interval = c(0.001, 0.5), x=X, y=Y) tmp <- rqss(Y~qss(X, lambda=tmplamstar$min), tau=1-p) predict.rqss(tmp,newdata=data.frame(X=x0s)) }, error=function(err) rep.int(NA,NUM.X0)) #warning=function(err) rep.int(NA,NUM.X0) # fYX50 <- tryCatch(npcdens(npcdensbw(formula=Y~X,bwmethod='normal-reference'),txdat=X,tydat=Y,exdat=x0s,eydat=Q50)$condens, warning=function(w)rep.int(NA,NUM.X0), error=function(err)rep.int(NA,NUM.X0)) fYX1 <- tryCatch(npcdens(npcdensbw(formula=Y~X,bwmethod='normal-reference'),txdat=X,tydat=Y,exdat=x0s,eydat=Q1)$condens, warning=function(w)rep.int(NA,NUM.X0), error=function(err)rep.int(NA,NUM.X0)) fYX2 <- tryCatch(npcdens(npcdensbw(formula=Y~X,bwmethod='normal-reference'),txdat=X,tydat=Y,exdat=x0s,eydat=Q2)$condens, warning=function(w)rep.int(NA,NUM.X0), error=function(err)rep.int(NA,NUM.X0)) # kde.h <- hns(X,deriv.order=0) fX <- kde(X,eval.points=x0s,h=kde.h)$estimate # tmp <- qsreg(X,Y,alpha=1/2); Qpp50 <- predict(tmp,derivative=2,x=x0s) tmp <- qsreg(X,Y,alpha= p ); Qpp1 <- predict(tmp,derivative=2,x=x0s) tmp <- qsreg(X,Y,alpha=1-p); Qpp2 <- predict(tmp,derivative=2,x=x0s) for (ktype in 1:2) { #1=uniform, 2=Gaussian if (ktype==1) Kfn <- function(u)ifelse((-1/2<=u)&(u<=1/2),1,0) else Kfn <- dnorm for (ix in 1:NUM.X0) { # bias from Theorem 1: dtau=c(dtau.c,dtau.t) dtau1 <- (1/2)*Qpp1[ix]*mu2K[ktype] dtau2 <- (1/2)*Qpp2[ix]*mu2K[ktype] # bandwidth from (16) for median and (17) for other quantiles h50 <- n^(-1/5) * ((0.5*(1-0.5)*1*intKsq[ktype])/(fX[ix]*fYX50[ix]^2*Qpp50[ix]^2*mu2K[ktype]^2))^(1/5) hp <- h50 * ((2*p*(1-p))/(pi*dnorm(qnorm(p))^2))^(1/5) #same for p and 1-p x0 <- x0s[ix] tmp <- tryCatch({ kwgts <- Kfn((X-x0)/hp) # use bandwidths to get local linear QR estimators w/ rq() tmp <- rq(Y~X,weights=kwgts,tau=p) ahat1 <- predict(tmp,newdata=data.frame(X=x0)) tmp <- rq(Y~X,weights=kwgts,tau=1-p) ahat2 <- predict(tmp,newdata=data.frame(X=x0)) # use formulas derived from Gaussian process limits and conservative bias correction sig1 <- 1/(sqrt(n*hp*fX[ix])*fYX1[ix]) sig2 <- 1/(sqrt(n*hp*fX[ix])*fYX2[ix]) V <- p*(1-p)*intKsq[ktype]*(sig1^2+sig2^2) - 2*p^2*sig1*sig2*intKsq[ktype] # CI: est +/- bias +/- z*stdev BC <- c(0,0) B <- (dtau2-dtau1)*hp^2 if (QY.BC==1) BC <- c(B,B) else if (QY.BC==-1) BC <- c(max(0,B),min(0,B)) c(ahat2-ahat1 - BC[1] - qnorm(1-ALPHA/2)*sqrt(V), ahat2-ahat1 - BC[2] - qnorm(ALPHA/2)*sqrt(V)) }, error=function(err)c(NA,NA)) CI.lo <- ifelse(is.numeric(tmp[1]),tmp[1],NA) CI.hi <- ifelse(is.numeric(tmp[2]),tmp[2],NA) if (ktype==1) { QYu.CIs[ix] <- CI.lo; QYu.CIs[ix+NUM.X0] <- CI.hi } else if (ktype==2) { QYg.CIs[ix] <- CI.lo; QYg.CIs[ix+NUM.X0] <- CI.hi } else stop("ktype may only be 1 (uniform) or 2 (Gaussian).") } } sink() c(Lstat.CIs, QYu.CIs, QYg.CIs) } })[3]),file=OUTFILE,sep='\n',append=TRUE) CIs$pt.lo <- CIs.raw[,1:NUM.X0] CIs$pt.hi <- CIs.raw[,(NUM.X0+1):(2*NUM.X0)] CIs$jt.lo <- CIs.raw[,(2*NUM.X0+1):(3*NUM.X0)] CIs$jt.hi <- CIs.raw[,(3*NUM.X0+1):(4*NUM.X0)] CIs$QYu.pt.lo <- CIs.raw[,(4*NUM.X0+1):(5*NUM.X0)] CIs$QYu.pt.hi <- CIs.raw[,(5*NUM.X0+1):(6*NUM.X0)] CIs$QYg.pt.lo <- CIs.raw[,(6*NUM.X0+1):(7*NUM.X0)] CIs$QYg.pt.hi <- CIs.raw[,(7*NUM.X0+1):(8*NUM.X0)] } else { cat(sprintf("Time elapsed = %g seconds",system.time({ for (irep in 1:NREPLIC) { tmp <- quantile.inf.np(Y=Y, X=X, x0s=x0s, p=c(1-p,p), PSI=c(1,-1), ALPHA=ALPHA, JOINT.FLAG=T, ONESIDED=0, METHOD.TYPE='lincom', NORM.APPROX.CAL=FALSE, NORM.APPROX.Q=FALSE) CIs$pt.lo[irep,] <- tmp[[1]]$CI.pointwise.lows CIs$pt.hi[irep,] <- tmp[[1]]$CI.pointwise.highs CIs$jt.lo[irep,] <- tmp[[1]]$CI.joint.lows CIs$jt.hi[irep,] <- tmp[[1]]$CI.joint.highs } })[3]),file=OUTFILE,sep='\n',append=TRUE) } IQR0mat <- do.call(rbind, replicate(NREPLIC, IQR0s, simplify=FALSE)) na.mat <- ifelse(is.na(CIs$pt.lo-CIs$pt.hi),NA,TRUE) na.QYg.mat <- ifelse(is.na(CIs$QYg.pt.lo-CIs$QYg.pt.hi),NA,TRUE) keep.QYg.sub.mat <- ifelse(!is.na(na.mat) & !is.na(na.QYg.mat),TRUE,FALSE) CP.pt <- apply(na.mat & CIs$pt.lo<=IQR0mat & IQR0mat<=CIs$pt.hi,2,mean,na.rm=TRUE) CP.jt <- mean(apply(na.mat & CIs$jt.lo<=IQR0mat & IQR0mat<=CIs$jt.hi,1,prod,na.rm=TRUE)) CP.QYu.pt <- apply(na.mat & CIs$QYu.pt.lo<=IQR0mat & IQR0mat<=CIs$QYu.pt.hi,2,mean,na.rm=TRUE) CP.QYg.pt <- apply(na.mat & CIs$QYg.pt.lo<=IQR0mat & IQR0mat<=CIs$QYg.pt.hi,2,mean,na.rm=TRUE) CP.QYg.full.pt <- apply(CIs$QYg.pt.lo<=IQR0mat & IQR0mat<=CIs$QYg.pt.hi,2,mean,na.rm=TRUE) na.mat2 <- ifelse(is.na(na.mat),NA,0) med.len.pt <- apply(na.mat2 + CIs$pt.hi-CIs$pt.lo,2,median,na.rm=TRUE) med.len.QYu.pt <- apply(na.mat2 + CIs$QYu.pt.hi-CIs$QYu.pt.lo,2,median,na.rm=TRUE) med.len.QYg.pt <- apply(na.mat2 + CIs$QYg.pt.hi-CIs$QYg.pt.lo,2,median,na.rm=TRUE) med.len.QYg.full.pt <- apply(CIs$QYg.pt.hi-CIs$QYg.pt.lo,2,median,na.rm=TRUE) cat(sprintf("Keeping (Lstat) {%s} of %d=NREPLIC\n",paste0(sprintf("%d",colSums(na.mat,na.rm=TRUE)),collapse=", "),NREPLIC), file=OUTFILE,sep="",append=TRUE) cat(sprintf("Of those, QYg ok for {%s}\n",paste0(sprintf("%d",colSums(keep.QYg.sub.mat,na.rm=TRUE)),collapse=", "),NREPLIC), file=OUTFILE,sep="",append=TRUE) cat(sprintf("Overall, QYg ok for {%s} of %d=NREPLIC\n",paste0(sprintf("%d",colSums(na.QYg.mat,na.rm=TRUE)),collapse=", "),NREPLIC), file=OUTFILE,sep="",append=TRUE) col2str <- sprintf("%-10s",Fstr) #for LaTeX output below cat(" &&\\multicolumn{10}{c}{\\textit{Coverage Probability}} \\\\\n",file=OUTFILE,sep="",append=T) abbrevs <- data.frame(desc="method abbrevs") for (vn in c("CP.pt","CP.QYg.pt","CP.QYg.full.pt","CP.QYu.pt")) { assign(paste0(vn,".str"), paste0(sprintf("%5.3f",get(vn)),collapse=" & ")) if (vn=="CP.pt") { abbrevs[vn] <- "Lstat" } else if (vn=="CP.QYg.full.pt") { abbrevs[vn] <- "QYg(f)" } else { abbrevs[vn] <- substr(vn,4,nchar(vn)-3) } cat(sprintf("%-7s&%s& %s \\\\\n",abbrevs[vn],col2str,get(paste0(vn,".str"))),file=OUTFILE,sep="",append=TRUE) } # cat(" &&\\multicolumn{10}{c}{\\textit{Median Interval Length}} \\\\\n",file=OUTFILE,sep="",append=T) abbrevs <- data.frame(desc="method abbrevs") for (vn in c("med.len.pt","med.len.QYg.pt","med.len.QYg.full.pt","med.len.QYu.pt")) { assign(paste0(vn,".str"), paste0(sprintf("%5.3f",get(vn)),collapse=" & ")) if (vn=="med.len.pt") { abbrevs[vn] <- "Lstat" } else if (vn=="med.len.QYg.full.pt") { abbrevs[vn] <- "QYg(f)" } else { abbrevs[vn] <- substr(vn,9,nchar(vn)-3) } cat(sprintf("%-7s&%s& %s \\\\\n",abbrevs[vn],col2str,get(paste0(vn,".str"))),file=OUTFILE,sep="",append=TRUE) } cat("CP.pt = ",file=OUTFILE,sep="",append=TRUE) cat(sprintf("%5.3f & ",CP.pt),file=OUTFILE,sep=" ",append=TRUE) cat("\n",file=OUTFILE,append=TRUE) cat("CP.QYg.pt = ",file=OUTFILE,sep="",append=TRUE) cat(sprintf("%5.3f & ",CP.QYg.pt),file=OUTFILE,sep=" ",append=TRUE) cat("\n",file=OUTFILE,append=TRUE) cat("CP.QYg.full.pt = ",file=OUTFILE,sep="",append=TRUE) cat(sprintf("%5.3f & ",CP.QYg.full.pt),file=OUTFILE,sep=" ",append=TRUE) cat("\n",file=OUTFILE,append=TRUE) cat("CP.QYu.pt = ",file=OUTFILE,sep="",append=TRUE) cat(sprintf("%5.3f & ",CP.QYu.pt),file=OUTFILE,sep=" ",append=TRUE) cat("\n",file=OUTFILE,append=TRUE) # cat("med.len.pt = ",file=OUTFILE,sep="=",append=TRUE) cat(sprintf("%5.3f & ",med.len.pt),file=OUTFILE,sep=" ",append=TRUE) cat("\n",file=OUTFILE,append=TRUE) cat("med.len.QYg.pt = ",file=OUTFILE,sep="=",append=TRUE) cat(sprintf("%5.3f & ",med.len.QYg.pt),file=OUTFILE,sep=" ",append=TRUE) cat("\n",file=OUTFILE,append=TRUE) cat("med.len.QYg.full.pt = ",file=OUTFILE,sep="=",append=TRUE) cat(sprintf("%5.3f & ",med.len.QYg.full.pt),file=OUTFILE,sep=" ",append=TRUE) cat("\n",file=OUTFILE,append=TRUE) cat("med.len.QYu.pt = ",file=OUTFILE,sep="=",append=TRUE) cat(sprintf("%5.3f & ",med.len.QYu.pt),file=OUTFILE,sep=" ",append=TRUE) cat("\n",file=OUTFILE,append=TRUE) # cat(sprintf("CP.joint = %5.3f",CP.jt),file=OUTFILE,sep="\n",append=TRUE) } #end of CASES loop cat(sprintf("End time is %s",format(Sys.time(), "%X, %A, %d %b %Y")), sep='\n',file=OUTFILE,append=TRUE) if (PARALLEL>1) stopCluster(workers) #EOF