1print("#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++") 2print(" ================== Welcome to 1dGC.R ================== ") 3print("AFNI Vector (or Multivariate) Auto-Regressive (VAR or MAR) Modeling Package!") 4print("#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++") 5print("Version 1.2.2, Jan. 19, 2012") 6print("Author: Gang Chen (gangchen@mail.nih.gov)") 7# [PT: Aug 12, 2020] commented out, bc link no longer works 8#print("Website: https://afni.nimh.nih.gov/sscc/gangc/VAR.html") 9print(" SSCC/NIMH, National Institutes of Health, Bethesda MD 20892") 10print("#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++") 11 12first.in.path <- function(file) { 13 ff <- paste(strsplit(Sys.getenv('PATH'),':')[[1]],'/', file, sep='') 14 ff<-ff[lapply(ff,file.exists)==TRUE]; 15 #cat('Using ', ff[1],'\n'); 16 return(gsub('//','/',ff[1], fixed=TRUE)) 17} 18source(first.in.path('AFNIio.R')) 19 20 21# header assumed for multi-column files, but not for one-column ones 22# readMultiFiles <- function(nFiles, dim, inData) { 23# inFile <- vector('list', nFiles) # list of file names with path attached 24# fn <- vector('list', nFiles) # ist of file names 25# for (ii in 1:nFiles) { 26# inFile[[ii]] <- tclvalue( tkgetOpenFile( filetypes = 27# "{{Files} {.1D}} {{All files} *}", 28# title = paste('Choose number', ii, 'input file'))) 29# fn[[ii]] <- strsplit(inFile[[ii]], "/")[[1]][length(strsplit(inFile[[ii]], "/")[[1]])] 30# print(sprintf("No. %i file just read in: %s", ii, fn[[ii]])) 31# if (dim==1) inData[,ii] <- read.table(inFile[[ii]], header=FALSE) 32# if (dim==2) inData[[ii]] <- read.table(inFile[[ii]], header=TRUE) 33# } 34# return(inData) 35# } 36 37#readMultiFiles <- function(nFiles, dim, inData, type) { 38# inFile <- vector('list', nFiles) # list of file names with path attached 39# fn <- vector('list', nFiles) # list of file names 40# if (dim==1) inData <- matrix(data = NA, nrow = dim(read.table(inFile[[ii]], header=FALSE))[1], ncol = nFiles) 41# if (dim==2) inData <- vector('list', nFiles) 42# for (ii in 1:nFiles) { 43# inFile[[ii]] <- readline(sprintf("No. %i %s file name: ", ii, type)) 44# fn[[ii]] <- strsplit(inFile[[ii]], "/")[[1]][length(strsplit(inFile[[ii]], "/")[[1]])] 45# print(sprintf("No. %i file just read in: %s", ii, fn[[ii]])) 46# if (dim==1) inData[,ii] <- read.table(inFile[[ii]], header=FALSE) 47# if (dim==2) inData[[ii]] <- read.table(inFile[[ii]], header=TRUE) 48# } 49# return(inData) 50#} 51 52readMultiFiles <- function(nFiles, dim, type) { 53 inFile <- vector('list', nFiles) # list of file names with path attached 54 fn <- vector('list', nFiles) # list of file names 55 if (dim==1) { 56 inFile[[1]] <- readline(sprintf("No. 1 %s file name: ", type)) 57 inData <- matrix(data = NA, nrow = dim(read.table(inFile[[1]], header=FALSE))[1], ncol = nFiles) 58 } 59 if (dim==2) inData <- vector('list', nFiles) 60 for(ii in 1:nFiles) { 61 if(((dim==1)&(ii>1))|(dim==2)) inFile[[ii]] <- readline(sprintf("No. %i %s file name: ", ii, type)) 62 fn[[ii]] <- strsplit(inFile[[ii]], "/")[[1]][length(strsplit(inFile[[ii]], "/")[[1]])] 63 print(sprintf("No. %i file just read in: %s", ii, fn[[ii]])) 64 #browser() 65 if (dim==1) inData[,ii] <- read.table(inFile[[ii]], header=FALSE)[,1] 66 if (dim==2) inData[[ii]] <- read.table(inFile[[ii]], header=TRUE) 67 } 68 return(inData) 69} 70 71 72plotTS <- function(dataFrame, nCurves, msg) { 73 if (nCurves <= 5) { 74 dev.new(); par(mfrow=c(nCurves, 1)) 75 for (ii in 1:nCurves) { 76 plot(dataFrame[,ii], ann=FALSE, axes=TRUE) 77 if (ii==1) title(msg) 78 lines(dataFrame[,ii]) 79 mtext(sprintf("%s", names(dataFrame)[ii]), side=2, line=2) 80 } 81 mtext("time", side=1, line=2.5) 82 } else for (ii in 1:nCurves) { 83 dev.new() 84 plot(dataFrame[,ii], ann=FALSE, axes=TRUE) 85 title(msg) 86 lines(dataFrame[,ii]) 87 mtext(sprintf("%s", names(dataFrame)[ii]), side=2, line=2) 88 mtext("time", side=1, line=2.5) 89 } 90} 91 92# plot network 93plotNet <- function(net, selfLoop, edgeWd, arrScl, edgeCol, msg) { 94 netData <- network(net, loops=selfLoop, directed=TRUE) 95 dev.new() 96 97 if(sum(abs(net)>10e-10)==1) # only ONE path in the network, it seems the matrix collapses into one number 98 plot.network(netData, displaylabels=selfLoop, mode="circle", edge.lwd=edgeWd, 99 arrowhead.cex=arrScl, edge.col=sum((abs(net)>10e-10)*(edgeCol)), loop.cex=5, 100 boxed.label=FALSE, label.pos=0, vertex.col=3) else 101 plot.network(netData, displaylabels=selfLoop, mode="circle", edge.lwd=edgeWd, 102 arrowhead.cex=arrScl, edge.col=edgeCol, loop.cex=5, 103 boxed.label=FALSE, label.pos=0, vertex.col=3) 104 105 title(msg) 106} 107 108 109### rma version 0.562 by Wolfgang Viechtbauer, Department of Methodology and Statistics, University of Maastricht 110 111rma <- function(yi, vi, mods=NULL, method="REML", addint=TRUE, ci=95, digits=4, btt=NULL, tau2=NULL, knha=FALSE, subset=NULL, ll=FALSE, control = list()) { 112 113 if (length(yi) != length(vi)) { 114 stop("Length of yi and vi is not the same.") 115 } 116 117 if (is.vector(mods)) { 118 mods <- cbind(mods) 119 } 120 121 if (is.data.frame(mods)) { 122 mods <- as.matrix(mods) 123 } 124 125 k <- length(yi) 126 ids <- 1:k 127 128 if (!is.null(subset)) { 129 yi <- yi[subset] 130 vi <- vi[subset] 131 mods <- mods[subset, , drop=FALSE] 132 ids <- ids[subset] 133 k <- length(yi) 134 } 135 136 if (sum(vi <= 0) > 0) { 137 null.ids <- which(vi<=0) 138 yi <- yi[-c(null.ids)] 139 vi <- vi[-c(null.ids)] 140 mods <- mods[-c(null.ids), , drop=FALSE] 141 ids <- ids[-c(null.ids)] 142 k <- length(yi) 143 warning("Outcomes with non-positive sampling variances have been excluded from the analysis.") 144 } 145 146 if (k <= 1) { 147 stop("Processing terminated since k <= 1.") 148 } 149 150 if (is.null(mods) && addint == FALSE) { 151 warning("Must either include an intercept (addint=TRUE) and/or moderators in model. Coerced intercept into the model.") 152 addint <- TRUE 153 } 154 155 Y <- as.matrix(yi) 156 157 if (addint == TRUE) { 158 X <- cbind(intrcpt=rep(1,k), mods) 159 } else { 160 X <- mods 161 } 162 163 p <- dim(X)[2] 164 165 if (method == "FE" || method == "FU") { 166 if (p > k) { 167 stop("The number of parameters to be estimated is larger than the number of observed outcomes.") 168 } 169 } else { 170 if (is.null(tau2)) { 171 if (p > k-1) { 172 stop("The number of parameters to be estimated is larger than the number of observed outcomes.") 173 } 174 } else { 175 if (p > k) { 176 stop("The number of parameters to be estimated is larger than the number of observed outcomes.") 177 } 178 } 179 } 180 181 if ( (p == 1) && (sum(X == 1) == k) ) { 182 int.only <- TRUE 183 } else { 184 int.only <- FALSE 185 } 186 187 tr <- function(X) { 188 sum(diag(X)) 189 } 190 191 con <- list(tau2.init=NULL, threshold=10^-5, maxiter=50, maxtau2=50, verbose=FALSE) 192 con[names(control)] <- control 193 194 se.tau2 <- NA 195 I2 <- NA 196 s2w <- NA 197 198 if (is.null(tau2) == TRUE) { 199 200 if (!is.element(method, c("FE", "FU", "HS", "HE", "DL", "SJ", "ML", "REML", "EB"))) { 201 stop("Unknown 'method' specified.") 202 } 203 204 if (method == "HS") { 205 wi <- 1/vi 206 W <- diag(wi) 207 P <- W - W %*% X %*% solve(t(X) %*% W %*% X) %*% t(X) %*% W 208 RSS <- t(Y) %*% P %*% Y 209 tau2 <- RSS/sum(wi) - k/sum(wi) 210 } 211 212 if (method == "HE") { 213 P <- diag(k) - X %*% solve(t(X) %*% X) %*% t(X) 214 RSS <- t(Y) %*% P %*% Y 215 tau2 <- ( RSS - tr( P %*% diag(vi) ) ) / (k-p) 216 } 217 218 if (method == "DL") { 219 wi <- 1/vi 220 W <- diag(wi) 221 P <- W - W %*% X %*% solve(t(X) %*% W %*% X) %*% t(X) %*% W 222 RSS <- t(Y) %*% P %*% Y 223 tau2 <- ( RSS - (k-p) ) / tr(P) 224 } 225 226 if (method == "SJ") { 227 tau2.0 <- var(yi) * (k-1)/k 228 wi <- 1/(vi + tau2.0) 229 W <- diag(wi) 230 P <- W - W %*% X %*% solve(t(X) %*% W %*% X) %*% t(X) %*% W 231 tau2 <- (tau2.0/(k-p)) * t(Y) %*% P %*% Y 232 } 233 234 if (is.element(method, c("ML", "REML", "EB"))) { 235 236 conv <- 1 237 change <- 1000 238 iter <- 0 239 240 if (is.null(con$tau2.init)) { 241 tau2 <- max(0, var(yi) - mean(vi)) 242 } else { 243 tau2 <- con$tau2.init 244 } 245 246 while (change > con$threshold) { 247 if (con$verbose == TRUE) cat("Iteration:", iter, " Estimate of (Residual) Heterogeneity:", tau2, "\n") 248 iter <- iter + 1 249 tau2.old <- tau2 250 wi <- 1/(vi + tau2) 251 W <- diag(wi) 252 P <- W - W %*% X %*% solve(t(X) %*% W %*% X) %*% t(X) %*% W 253 if (method == "REML") { 254 adj <- solve( tr(P%*%P) ) %*% ( t(Y)%*%P%*%P%*%Y - tr(P) ) 255 } 256 if (method == "ML") { 257 adj <- solve( tr(W%*%W) ) %*% ( t(Y)%*%P%*%P%*%Y - tr(W) ) 258 } 259 if (method == "EB") { 260 adj <- solve( tr(W) ) %*% ( (k/(k-p)) %*% t(Y)%*%P%*%Y - k ) 261 } 262 while (tau2 + adj < 0) { 263 adj <- adj / 2 264 } 265 tau2 <- tau2 + adj 266 change <- abs(tau2.old - tau2) 267 if (iter > con$maxiter) { 268 conv <- 0 269 break 270 } 271 } 272 273 if (conv == 0) { 274 stop("Fisher scoring algorithm did not converge. Try increasing the number of iterations (maxiter), lowering the threshold (threshold), or use a different estimation method.") 275 } 276 277 if (method == "ML") { 278 se.tau2 <- sqrt( 2/sum(wi^2) ) 279 } 280 if (method == "REML") { 281 se.tau2 <- sqrt( 2/tr(P%*%P) ) 282 } 283 284 } 285 286 tau2 <- max(0, c(tau2)) 287 288 } 289 290 if (method == "FE" || method == "FU") { 291 tau2 <- 0 292 } 293 294 wi <- 1/vi 295 W <- diag(wi) 296 P <- W - W %*% X %*% solve(t(X) %*% W %*% X) %*% t(X) %*% W 297 QE <- t(Y) %*% P %*% Y 298 QEp <- 1 - pchisq(QE, df=k-p) 299 300 if (int.only == TRUE) { 301 s2 <- (k-1)*sum(wi) / ( sum(wi)^2 - sum(wi^2) ) 302 I2 <- 100 * tau2 / (tau2 + s2) 303 } 304 305 if (method == "FU") { 306 b <- solve(t(X) %*% X) %*% t(X) %*% Y 307 vb <- solve(t(X) %*% X) %*% t(X) %*% diag(vi) %*% X %*% solve(t(X) %*% X) 308 } else { 309 wi <- 1/(vi + tau2) 310 W <- diag(wi) 311 b <- solve(t(X) %*% W %*% X) %*% t(X) %*% W %*% Y 312 vb <- solve(t(X) %*% W %*% X) 313 } 314 315 if (knha == TRUE) { 316 P <- W - W %*% X %*% solve(t(X) %*% W %*% X) %*% t(X) %*% W 317 s2w <- c( t(Y) %*% P %*% Y ) / (k-p) 318 vb <- s2w * vb 319 if (method == "FE" || method == "FU") { 320 warning("The Knapp & Hartung method is not meant to be used in the context of fixed-effects models.") 321 } 322 } 323 324 if (is.null(btt)) { 325 if (p > 1) { 326 if (addint == TRUE) { 327 btt <- 2:p 328 } else { 329 btt <- 1:p 330 } 331 } else { 332 btt <- 1 333 } 334 } 335 336 m <- length(btt) 337 338 QME <- t(b)[btt] %*% solve(vb[btt,btt]) %*% b[btt] 339 340 if (knha == FALSE) { 341 QMEp <- 1 - pchisq(QME, df=m) 342 } else { 343 QMEp <- 1 - pf( QME/m, df1=m, df2=k-p ) 344 } 345 346 alpha <- (100-ci)/100 347 se <- sqrt(diag(vb)) 348 z <- b / se 349 350 if (knha == FALSE) { 351 zp <- 2*(1-pnorm(abs(z))) 352 crit <- qnorm(1-alpha/2) 353 } else { 354 zp <- 2*(1-pt(abs(z), df=k-p)) 355 crit <- qt(1-alpha/2, df=k-p) 356 } 357 358 ci.lb <- b - crit * se 359 ci.ub <- b + crit * se 360 361 if (method == "FU") { 362 P <- W - W %*% X %*% solve(t(X) %*% X) %*% t(X) - X %*% solve(t(X) %*% X) %*% t(X) %*% W + X %*% solve(t(X) %*% X) %*% t(X) %*% W %*% X %*% solve(t(X) %*% X) %*% t(X) 363 } else { 364 P <- W - W %*% X %*% solve(t(X) %*% W %*% X) %*% t(X) %*% W 365 } 366 ll.ML <- -1/2 *(k) * log(2*pi) - 1/2 * sum( log(vi + tau2) ) - 1/2 * t(Y) %*% P %*% Y 367 ll.REML <- -1/2 *(k-p) * log(2*pi) + 0/2 * log( det(t(X)%*%X) ) - 1/2 * sum( log(vi + tau2) ) - 1/2 * log( det(t(X)%*%W%*%X) ) - 1/2 * t(Y) %*% P %*% Y 368 dev.ML <- -2 * ll.ML 369 dev.REML <- -2 * ll.REML 370 AIC.ML <- -2 * ll.ML + 2*(p + (if (method == "FE" || method == "FU") 0 else 1) ) 371 BIC.ML <- -2 * ll.ML + (p + (if (method == "FE" || method == "FU") 0 else 1) ) * log(k) 372 AIC.REML <- -2 * ll.REML + 2*(p + (if (method == "FE" || method == "FU") 0 else 1) ) 373 BIC.REML <- -2 * ll.REML + (p + (if (method == "FE" || method == "FU") 0 else 1) ) * log(k-p) 374 fit.stats <- c(ll.ML=ll.ML, dev.ML=dev.ML, AIC.ML=AIC.ML, BIC.ML=BIC.ML, ll.REML=ll.REML, dev.REML=dev.REML, AIC.REML=AIC.REML, BIC.REML=BIC.REML) 375 376 res <- list(b, se, z, zp, ci.lb, ci.ub, vb, tau2, se.tau2, k, p, m, fit.stats, QE, QEp, QME, QMEp, s2w, I2, int.only, yi, vi, X, ids, method, knha, btt, addint, digits, ci, ll, control) 377 names(res) <- c("b", "se", "z", "zp", "ci.lb", "ci.ub", "vb", "tau2", "se.tau2", "k", "p", "m", "fit.stats", "QE", "QEp", "QME", "QMEp", "s2w", "I2", "int.only", "yi", "vi", "X", "ids", "method", "knha", "btt", "addint", "digits", "ci", "ll", "control") 378 class(res) <- c("rma") 379 res 380 381} 382 383 384pkgLoad('network') 385 386# [PT: Aug 12, 2020] commented out, bc link no longer works 387#print("~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~") 388#print("Visit https://afni.nimh.nih.gov/sscc/gangc/VAR.html and make sure") 389#print("you've acquired the data for the analysis in desirable data format.") 390#print("~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~") 391 392print("################################################################") 393print("Please consider citing the following if this program is useful for you:") 394cat("\n\tChen et al., Vector autoregression, structural equation modeling, and\n") 395cat("\ttheir synthesis in neuroimaging data analysis, Comput. Biol. Med. (2011),\n") 396cat("\tVolume 41, Issue 12, 1142-1155.\n\n") 397print("################################################################") 398 399print("Use CNTL-C on Unix or ESC on GUI version of R to stop at any moment.") 400 401anaType <- as.integer(readline("Analysis type (0: quit; 1: individual; 2: group)? ")) 402 403if (anaType==1) { 404#libLoad("gsl") # Legendre polynomials 405#libLoad("vars") # VAR modeling 406pkgLoad(c('gsl', 'vars')) 407 408anotherAna <- 1 409while (anotherAna==1) { 410 411 412nROIs <- as.integer(readline("Number of regions/nodes (e.g., 6)? ")) # number of ROIs 413print("Header with one line of labels is optional in multi-column files, but NOT allowed in one-column files.") 414yForm <- as.integer(readline("ROI time series data type (0: MULTIPLE one-column files; 1: ONE multi-column file)? ")) # input format 415 416if (yForm) { # read ROI file (in dataframe), and take label from the header 417# fn <- tclvalue( tkgetOpenFile( filetypes = 418# "{{ROI files in multi-column 1D format} {.1D}} {{All files} *}", 419# title = paste('Choose ROIs time series file'))) 420 fn <- readline("ROIs time series file name: ") 421 print(sprintf("File just read in: %s", strsplit(fn, "/")[[1]][length(strsplit(fn, "/")[[1]])])) 422 yHeader <- as.integer(readline("Does this multi-column file have a header (0: no; 1: yes)? ")) 423 if (yHeader == 1) myData <- read.table(fn, header=TRUE) else { 424 myData <- read.table(fn, header=FALSE) 425 for (ii in 1:nROIs) names(myData)[ii] <- readline(sprintf("Name for region/node number %i? ", ii)) 426 } 427# nROIs <- ncol(myData) 428 nTotal <- nrow(myData) 429} else { 430 431 fn <- vector('list', nROIs) 432# fn[[1]] <- tclvalue(tkgetOpenFile(filetypes = 433# "{{ROI files in one-column 1D format} {.1D}} {{All files} *}", 434# title = paste('Choose number', 1, 'ROI time series file'))) 435 fn[[1]] <- readline(sprintf("No. 1 ROI time series file name: ")) 436 print(sprintf("No. %i file just read in: %s", 1, strsplit(fn[[1]], "/")[[1]][length(strsplit(fn[[1]], "/")[[1]])])) 437 myData <- data.frame(matrix(data=NA, nrow=dim(read.table(fn[[1]], header=FALSE))[1], ncol=nROIs, dimnames = NULL)) 438 myData[,1] <- read.table(fn[[1]], header=FALSE) 439 for (ii in 2:nROIs) { # read ROI 1D file (in 1 column) 440# fn[[ii]] <- tclvalue( tkgetOpenFile( filetypes = 441# "{{ROI files in one-column format} {.1D}} {{All files} *}", 442# title = paste('Choose number', ii, 'ROI time series file'))) 443 fn[[ii]] <- readline(sprintf("No. %i ROI time series file name: ", ii)) 444 print(sprintf("No. %i file just read in: %s", ii, strsplit(fn[[ii]], "/")[[1]][length(strsplit(fn[[ii]], "/")[[1]])])) 445 myData[,ii] <- read.table(fn[[ii]], header=FALSE) 446 } #for (ii in 2:nROIs) 447 for (ii in 1:nROIs) names(myData)[ii] <- readline(sprintf("Name for region/node number %i? ", ii)) 448# take labels from the filenames 449# ROIlab <- vector('list', nROIs) 450# for (ii in 1:nROIs) ROIlab[[ii]] <- strsplit(strsplit(fn[[ii]], "/")[[1]][length(strsplit(fn[[ii]], "/")[[1]])], ".1D")[[1]] 451# names(myData) <- ROIlab 452 nTotal <- dim(myData)[1] 453} 454print("#++++++++++++++++++++++++++++++++++++++++++++") 455print("If there are n consecutive chunks of data, enter n-1 here.") 456print("If all the time series are consecutive, enter 0 breaks.") 457nChunks <- as.integer(readline("Number of breaks in the time series? "))+1 # number of runs 458 459nPts <- array(data = NA, dim=nChunks) 460if (nChunks == 1) nPts[1] <- nTotal else 461for (ii in 1:nChunks) nPts[ii] <- as.integer(readline(paste("Length of number", ii, "run/segment? "))) 462 463print("#++++++++++++++++++++++++++++++++++++++++++++") 464 465print("Prior detrending is NOT recommended due to potential complications. ") 466print("Trend can be modeled through specifying the order of a polynomial here, ") 467print("or can be included as part of covariates later on. If you plan to model the") 468print("trend with your own regressors or don't need to model it, choose -1 here.") 469print("If trend has already been removed (not recommended), choose 0 here.") 470nPoly <- as.integer(readline("Order of polynomials for drifting effects (note: -1 means no trend removal!)? ")) # Legendre 471 472print("#++++++++++++++++++++++++++++++++++++++++++++") 473 474COV <- as.integer(readline("Any other covariates or confounding effects than drifting (0: no; 1: yes)? ")) 475if (as.logical(COV)) { 476 nCOVs <- as.integer(readline("Number of covariates (e.g., 6)? ")) # number of regions: 6 477 print("Header with one line of labels is optional in multi-column files, but NOT allowed in one-column files.") 478 covForm <- as.integer(readline("Covariates data type (0: MULTIPLE one-column files; 1: ONE multi-column file)? ")) # covariates format 479 if (covForm) { 480# fncov <- tclvalue( tkgetOpenFile( filetypes = 481# "{{Covariate File} {.1D}} {{All files} *}", 482# title = paste('Choose covariates file in multi-column format'))) 483 fncov <- readline("Covariates file name: ") 484 covHeader <- as.integer(readline("Does this multi-column file have a header (0: no; 1: yes)? ")) 485 if (covHeader == 1) exData <- read.table(fncov, header=TRUE) else { 486 exData <- read.table(fncov, header=FALSE) 487 for (ii in 1:nCOVs) names(exData)[ii] <- readline(sprintf("Name for covariate number %i? ", ii)) 488 } 489 } else { 490# covn <- vector('list', nCOVs) 491 exData <- data.frame(matrix(data=NA, nrow=nTotal, ncol=nCOVs, dimnames = NULL)) 492 exData <- readMultiFiles(nCOVs, 1, "covariate") # 1: assuming no header 493# covFN <- exTmp[[1]]; exData <- exTmp[[2]] 494 495# for (ii in 1:nCOVs) { 496# covn[[ii]] <- tclvalue( tkgetOpenFile( filetypes = 497# "{{Covariate Files} {.1D}} {{All files} *}", 498# title = paste('Choose number', ii, 'covariate time series file'))) 499# print(sprintf("No. %i file just read in: %s", ii, strsplit(covn[[ii]], "/")[[1]][length(strsplit(covn[[ii]], "/")[[1]])])) 500# exData[,ii] <- read.table(covn[[ii]], header=FALSE) 501# } 502 for (ii in 1:nCOVs) names(exData)[ii] <- readline(sprintf("Name for covariate number %i? ", ii)) 503# COVlab <- vector('list', nCOVs) 504# for (ii in 1:nCOVs) COVlab[[ii]] <- strsplit(covFN[[ii]], ".1D")[[1]] 505# names(exData) <- COVlab 506 } 507} else {exData <- NULL; nCOVs <- 0} 508print("#++++++++++++++++++++++++++++++++++++++++++++") 509 510# plot out the time series, and let the user make sure they look OK 511plotROIs <- as.integer(readline("Plot out the original ROI time series (0: no; 1: yes)? ")) 512#print("Check the fitting for each ROI time series") 513if (plotROIs) plotTS(myData, nROIs, "Original time series") 514 515 516if (as.logical(COV)) { 517 plotCov <- as.integer(readline("Plot out covariates time series provided by you (0: no; 1: yes)? ")) 518 if (plotCov) plotTS(exData, nCOVs, "Covariate time series") 519} # if (as.logical(COV)) 520print("-----------------") 521 522# scaling the data per run per ROI 523newData <- myData 524print(sprintf("If normalization was NOT performed during pre-processing, you can scale the data now.")) 525scaleTS <- as.integer(readline("Scale the ROI time series (0: no; 1: yes)? ")) 526if (scaleTS) { 527 jumpPts <- 0 528 for (ii in 1:nChunks) { 529 for (jj in 1:nROIs) newData[(jumpPts+1):(jumpPts+nPts[ii]),jj] <- 530 myData[(jumpPts+1):(jumpPts+nPts[ii]),jj]/mean(myData[(jumpPts+1):(jumpPts+nPts[ii]),jj]) 531 if (ii < nChunks) jumpPts <- jumpPts+nPts[ii] 532 } 533} 534 535# create exogenous variables with Legendre polynomials from gsl 536if (nPoly > -1) { 537 trendMat <- as.data.frame(array(0, dim = c(nTotal, (nPoly+1)*nChunks))) 538 jumpPts <- 0 539 for (ii in 1:nChunks) { 540 trendMat[(jumpPts+1):(jumpPts+nPts[ii]),(1+(nPoly+1)*(ii-1)):((nPoly+1)*ii)] <- 541 t(legendre_Pl_array(nPoly, seq(from=-1,to=1,len=nPts[ii]))) 542 names(trendMat)[(1+(nPoly+1)*(ii-1)):((nPoly+1)*ii)] <- sprintf("Run%iTrend%i", ii, seq(nPoly+1)-1) 543 if (ii < nChunks) jumpPts <- jumpPts+nPts[ii] 544 } 545 if (is.null(exData)) exMat <- trendMat else exMat <- cbind(trendMat, exData) 546} else exMat <- exData # if no baseline and trend, do nothing 547# plot out those polynomials here??? 548 549maxLags <- as.integer(readline("Highest order (or number of lags) for the VAR model (1,2,...)? ")) 550 551critSel <- VARselect(newData, lag.max = maxLags, type = "none", exogen=exMat) 552print(sprintf("Suggested orders for VAR:")) 553print(critSel$selection) 554print("AIC: Akaike Information Criterion;") 555print("HQ: Hannan-Quinn criterion;") 556print("SC: Schwartz Criterion;") 557print("FPE: Final Prediction Error criterion.") 558print(" ") 559print("Usually consistency exists between AIC and FPE, and between") 560print("HQ and SC. AIC and FPE tend to overestimate the 'true order'.") 561print("Since there is no such a universally best criterion, it may be") 562print("better to try different analysis with various orders within") 563print("the range covered by the 4 criteria.") 564 565anotherLag <- TRUE 566while (anotherLag) { 567 568#print(sprintf("Select the order of VAR model based on above criteria:")) 569nLags <- as.integer(readline("Select order of VAR model based on above criteria (e.g., 3)? ")) 570 571# generate intervention dummy variables for across-run/block breaks: nLags dummies per run 572if (nChunks > 1) { 573 breakMat <- as.data.frame(array(0, dim = c(nTotal, (nChunks-1)*nLags))) 574 jumpPts <- 0 575 for (ii in 1:(nChunks-1)) { 576 jumpPts <- jumpPts+nPts[ii] 577 for (jj in 1:(nLags)) { 578 breakMat[,(ii-1)*nLags+jj] <- c(rep(0, jumpPts+jj-1), 1, rep(0, nTotal-jumpPts-jj)) 579 names(breakMat)[(ii-1)*nLags+jj] <- sprintf("Run%iLag%i", ii, jj) 580 } 581 } 582 if (is.null(exMat)) exMatMod <- breakMat else exMatMod <- cbind(breakMat, exMat) 583 print(sprintf("Suggested orders for VAR updated:")) 584 print(VARselect(newData, lag.max = maxLags, type = "none", exogen=exMatMod)$selection) 585} else exMatMod <- exMat 586print("#++++++++++++++++++++++++++++++++++++++++++++") 587 588fm <- VAR(newData, p=nLags, type="none", exogen=exMatMod) 589 590qualityCheck <- as.integer(readline("Want to run a few quality check tests (0: no; 1: yes)? ")) 591if (qualityCheck) { 592# the modulus of the eigenvalues (presumably less than 1 as stable condition) in the reverse characteristic polynomial; stable process is stationary, but the converse is not true 593#print("Quality check of the model:") 594if (prod(roots(fm)<1)) print("Eigenvalues of the companion coefficient matrix indciate that the VAR(p) process is stable and thus stationary") else print("The VAR(p) process seems unstable and thus is not stationary") 595print("-----------------") 596print("Normality testing of the residuals") 597print(normality.test(fm)) 598print("-----------------") 599print("Serial correlation test:") 600print(serial.test(fm, lags.pt=16, lags.bg=5, type=c("PT.asymptotic"))) 601print(serial.test(fm, lags.pt=16, lags.bg=5, type=c("PT.adjusted"))) 602print(serial.test(fm, lags.pt=16, lags.bg=5, type=c("BG"))) 603print(serial.test(fm, lags.pt=16, lags.bg=5, type=c("ES"))) 604print("-----------------") 605print("Autoregressive conditional heteroskedasticity (ARCH) test") 606print(arch.test(fm)) 607archPlot <- as.integer(readline("Plot out ARCH test result (0: no; 1: yes)? ")) 608if (archPlot) {dev.new(); plot(arch.test(fm))} 609print("-----------------") 610statPlot <- as.integer(readline("Plot out stability test (0: no; 1: yes)? ")) 611if (statPlot) { 612 print("Available empirical fluctuation process (efp) types are:") 613 print("1. OLS-CUSUM: CUmulative SUMs of Ordinary Least Squares residuals;") 614 print("2. Rec-CUSUM: CUmulative SUMs of Recursive residuals;") 615 print("3. OLS-MOSUM: MOving SUMs of Ordinary Least Squares residuals;") 616 print("4. Rec-MOSUM: MOving SUMs of Recursive residuals;") 617 print("5. RE: Recursive Estimates of the path coefficients;") 618 print("6. ME: Moving Estimates of the path coefficients;") 619 print("7. Score-CUSUM: CUmulative SUMs of the ML Scores;") 620 print("8. Score-MOSUM: MOving SUMs of the ML scores.") 621 anotherType <- TRUE 622 while (anotherType) { 623 procNo <- as.integer(readline("Select process type (1-8)? ")) 624 if (procNo==1) procType <- "OLS-CUSUM" 625 if (procNo==2) procType <- "Rec-CUSUM" 626 if (procNo==3) procType <- "OLS-MOSUM" 627 if (procNo==4) procType <- "Rec-MOSUM" 628 if (procNo==5) procType <- "RE" 629 if (procNo==6) procType <- "ME" 630 if (procNo==7) procType <- "Score-CUSUM" 631 if (procNo==8) procType <- "Score-CUSUM" 632 dev.new(); plot(stability(fm, type = procType, h = 0.15, dynamic = FALSE, rescale = TRUE)) 633 anotherType <- as.integer(readline("Want to plot stability with another type (0: no; 1: yes)? ")) 634 } # while (anotherType) 635} # if (statPlot) 636 637print("-----------------") 638checkCov <- as.integer(readline("Check significance of covariates (0: no; 1: yes)? ")) 639if (checkCov) { 640 anotherCovPth <- TRUE 641 totCOVs <- (nPoly+1)*nChunks+nCOVs # total covariates 642 while (anotherCovPth) { 643 pCovThresh <- as.numeric(readline("p-threshold for covariates (e.g., 0.05)? ")) 644 #Info about all the covariates: 645 #lapply(coef(fm), function(x) x[nROIs*nLags+(1:totCOVs),]) 646 covPList <- lapply(coef(fm), function(x) x[nROIs*nLags+(1:totCOVs),4]<=pCovThresh) 647 #covSigList <- vector(mode="logical", totCOVs) 648 #for (ii in 1:nROIs) covSigList <- covSigList+covPList[[ii]] 649 covSigList <- apply(do.call(cbind, covPList), 1, sum) 650 # detailed info: apply(do.call(cbind, covPList), c(1,2), sum) 651 if (length(covSigList[covSigList==0])) { 652 print(sprintf("With a threshold of %f the following covariates don't show significance in the model:", pCovThresh)) 653 print(names(covSigList[covSigList==0])) 654 print("You may consider removing them from the model. However, when polynomial terms") 655 print("show up in the above list, only if the highest order of the polynomials") 656 print("indicates insignificant would you try decreasing the order.") 657 } else print(sprintf("All covarates show significance in the model with a threshold of %f.", pCovThresh)) 658 anotherCovPth <- as.integer(readline("Want to try another p-threshold for covariates (0: no; 1: yes)? ")) 659 } # 660} 661} # model quality check 662 663print("#++++++++++++++++++++++++++++++++++++++++++++") 664 665# spill out the original path matrix with direction going from rows to columns 666 667netMatR <- array(data=NA, dim=c(nLags, nROIs, nROIs)) # original path coefficient matrix 668netMatT <- array(data=NA, dim=c(nLags, nROIs, nROIs)) # t values matrix 669for (ii in 1:nROIs) for (jj in 1:nROIs) for (kk in 1:nLags) { 670 netMatR[kk,jj,ii] <- coef(fm)[[ii]][jj+nROIs*(kk-1), 1] # path coefficients 671 netMatT[kk,jj,ii] <- coef(fm)[[ii]][jj+nROIs*(kk-1), 3] # t values 672} 673 674for (ii in 1:nLags) { 675 print(sprintf("Path coefficient matrix with a lag of %i (direction goes from row to column):", ii)) 676 print(matrix(netMatR[ii,,], nrow = nROIs, ncol = nROIs, dimnames = list(names(myData), names(myData)))) 677 saveMat <- as.integer(readline("Save above path matrix for group analysis (0: no; 1: yes)? ")) 678 if (saveMat) { 679 matName <- as.character(readline("File name prefix (e.g., PathLag1Subj1)? ")) 680 write.table(netMatR[ii,,], file=sprintf("%s.1D", matName), append=FALSE, row.names=names(myData), col.names=names(myData)) 681 } 682 print("-----------------") 683 print(sprintf("Matrix of t values with a lag of %i (direction goes from row to column):", ii)) 684 print(matrix(netMatT[ii,,], nrow = nROIs, ncol = nROIs, dimnames = list(names(myData), names(myData)))) 685 print(sprintf("DFs = %i for null hypothesis H_0: a path coefficient = 0.", summary(fm)$varresult[[1]]$df[2])) 686 saveMatT <- as.integer(readline("Save matrix of t values for group analysis (0: no; 1: yes)? ")) 687 if (saveMatT) { 688 matName <- as.character(readline("File name prefix (e.g., TLag1Subj1)? ")) 689 write.table(netMatT[ii,,], file=sprintf("%s.1D", matName), append=FALSE, row.names=names(myData), col.names=names(myData)) 690 } 691 print("-----------------") 692} 693 694if (nLags>1) { # overall network with all lags collapsed 695 #libLoad("car") # for linear.hypothesis 696 pkgLoad('car') 697 # initialization for overall network across lags 698 netCMatF <- matrix(data=NA, nrow=nROIs, ncol=nROIs, dimnames = list(names(myData), names(myData))) 699 netCMatP <- matrix(data=NA, nrow=nROIs, ncol=nROIs, dimnames = list(names(myData), names(myData))) 700 fTest <- vector("list", nROIs) # initialization for F test name list 701 for (ii in 1:nROIs) for (jj in 1:nLags) 702 fTest[[ii]] <- c(fTest[[ii]], sprintf("%s.l%d",names(myData)[ii],jj)) 703 for (ii in 1:nROIs) for (jj in 1:nROIs) { 704 ltTmp <- linear.hypothesis(fm$varresult[[ii]], fTest[[jj]]) 705 netCMatF[jj, ii] <- ltTmp$F[2]; netCMatP[jj, ii] <- ltTmp$Pr[2] 706 } 707 print(sprintf("Overall F matrix for ALL %i lags (direction goes from row to column):", nLags)) 708 print(netCMatF) 709 print(sprintf("Numerator DFs = %i, and Denominator DFs =%i", -ltTmp$Df[2], ltTmp$Res.Df[1])) 710 print(sprintf("Null hypothesis H_0: path coef from region i to j is 0 for all %i lags.", nLags)) 711 saveCMatF <- as.integer(readline("Save above F matrix (0: no; 1: yes)? ")) 712 if (saveCMatF) { 713 matCFName <- as.character(readline("File name prefix (e.g., FMatSubj1)? ")) 714 write.table(netCMatF, file=sprintf("%s.1D", matCFName), append=FALSE) 715 } 716 print("-----------------") 717 print(sprintf("Overall p-value matrix for ALL %i lags (direction goes from row to column):", nLags)) 718 print(netCMatP) 719 print(sprintf("Numerator DFs = %i, and Denominator DFs =%i", -ltTmp$Df[2], ltTmp$Res.Df[1])) 720 print(sprintf("Null hypothesis H_0: path coef from region i to j is 0 for all %i lags.", nLags)) 721 saveCMatP <- as.integer(readline("Save above p-value matrix (0: no; 1: yes)? ")) 722 if (saveCMatP) { 723 matCPName <- as.character(readline("File name prefix (e.g., PMatSubj1)? ")) 724 write.table(netCMatP, file=sprintf("%s.1D", matCPName), append=FALSE) 725 } 726 print("-----------------") 727# totCOVs <- (nPoly+1)*nChunks+nCOVs 728# fCovTest <- vector("list", totCOVs) # initialization for covariate F test name list 729# for (ii in 1:totCOVs) 730# fCovTest[[ii]] <- names(myData)[ii],jj) 731 732 733} # if (nLags>1) 734 735print("Covariance matrix of residuals:") 736print(summary(fm)$covres) 737saveCovMat <- as.integer(readline("Save above contemporaneous covariance matrix (0: no; 1: yes)? ")) 738if (saveCovMat) { 739 matCovName <- as.character(readline("File name prefix (e.g., CovSubj1)? ")) 740 write.table(summary(fm)$covres, file=sprintf("%s.1D", matCovName), append=FALSE, row.names=TRUE, col.names=TRUE) 741} 742print("-----------------") 743 744print("Correlation matrix of residuals:") 745print(summary(fm)$corres) 746saveCorMat <- as.integer(readline("Save above contemporaneous correlation matrix (0: no; 1: yes)? ")) 747if (saveCorMat) { 748 matCorName <- as.character(readline("File name prefix (e.g., CorSubj1)? ")) 749 write.table(summary(fm)$corres, file=sprintf("%s.1D", matCorName), append=FALSE, row.names=TRUE, col.names=TRUE) 750} 751print("-----------------") 752 753seeModel <- as.integer(readline("Want to see the model details and statistics (0: no; 1: yes)? ")) 754if (seeModel) print(summary(fm)) 755 756print("#++++++++++++++++++++++++++++++++++++++++++++") 757plotFit <- as.integer(readline("Plot out the model fit, residuals, residual ACF and PACF (0: no; 1: yes)? ")) 758print("-----------------") 759 760if (plotFit) { dev.new(); plot(fm) } 761 762anotherPth <- TRUE 763while (anotherPth) { 764 765print(sprintf("There are totally %i paths in the model. Select a low p value ", nLags*nROIs^2)) 766print("if you're concerned about multiple comparisons issue:") 767pThresh <- as.numeric(readline("p-threshold for causal effects (e.g., 0.05)? ")) 768 769# connection goes from row to column, which is the default in network package 770netMat <- array(data=NA, dim=c(nLags+1, nROIs, nROIs)) # thresholded network 771for (ii in 1:nROIs) 772 for (jj in 1:nROIs) { 773 for (kk in 1:nLags) 774 netMat[kk,jj,ii] <- sum(sum(coef(fm)[[ii]][jj+nROIs*(kk-1), 4] 775 <= pThresh)>0) 776 if (nLags>1) { 777 ltTmp <- linear.hypothesis(fm$varresult[[ii]], fTest[[jj]]) 778 netCMatF[jj, ii] <- ltTmp$F[2]; netCMatP[jj, ii] <- ltTmp$Pr[2] 779 } 780 } 781 782print("#++++++++++++++++++++++++++++++++++++++++++++") 783# spill the matrices with those insignificant ones being masked: direction goes from rows to cols 784for (ii in 1:nLags) { 785 print(sprintf("Path coefficient matrix with a lag of %i with insignificant ones masked with NAs:", ii)) 786# print(matrix(netMatR[ii,,]*netMat[ii,,], nrow = nROIs, ncol = nROIs, dimnames = list(names(myData), names(myData)))) 787 print(matrix(mapply(function(x,y) ifelse(x, y, NA), netMat[ii,,], netMatR[ii,,]), nrow = nROIs, ncol = nROIs, 788 dimnames = list(names(myData), names(myData)))) 789 print("-----------------") 790 print(sprintf("Matrix of t values with a lag of %i with insignificant ones masked with NAs:", ii)) 791# print(matrix(netMatT[ii,,]*netMat[ii,,], nrow = nROIs, ncol = nROIs, dimnames = list(names(myData), names(myData)))) 792 print(matrix(mapply(function(x,y) ifelse(x, y, NA), netMat[ii,,], netMatT[ii,,]), nrow = nROIs, ncol = nROIs, 793 dimnames = list(names(myData), names(myData)))) 794 print(sprintf("DFs = %i for null hypothesis H_0: a path coefficient = 0.", summary(fm)$varresult[[1]]$df[2])) 795 print("-----------------") 796} 797 798if (nLags>1) { 799 print(sprintf("Overall F matrix for all %i lags with insignificant ones masked with NAs:", nLags)) 800# print(netCMatF*(netCMatP<=pThresh)) 801 print(matrix(mapply(function(x,y) ifelse(x, y, NA), netCMatP<=pThresh, netCMatF), nrow = nROIs, ncol = nROIs, 802 dimnames = list(names(myData), names(myData)))) 803 print("-----------------") 804 print(sprintf("Overall p-value matrix for all %i lags with insignificant ones masked with NAs:", nLags)) 805# print(netCMatP*(netCMatP<=pThresh)) 806 print(matrix(mapply(function(x,y) ifelse(x, y, NA), netCMatP<=pThresh, netCMatP), nrow = nROIs, ncol = nROIs, 807 dimnames = list(names(myData), names(myData)))) 808} 809 810plotNetwork <- as.integer(readline("Plot out the identified network (0: no; 1: yes)? ")) 811if (plotNetwork) { 812 813set.seed(1702) 814net <- network.initialize(nROIs) 815attr(net, "vertex.names") <- names(myData) 816edgeScale <- as.numeric(readline("Scale factor for path thickness (e.g., 20)? ")) 817arrowScale <- as.numeric(readline("Scale factor for arrows (e.g., 2)? ")) 818print("A self-loop indicates a region/node has effect on itself from t-1 to t.") 819selfLoop <- as.integer(readline("Show self-loops in the network (0: no; 1: yes)? ")) 820 821for (ii in 1:nLags) { 822 netData <- netMat[ii,,] 823 rownames(netData) <- names(myData); colnames(netData) <- names(myData) 824 # if path +, col<-2 (red); if path -, col<-4 (blue) 825 plotNet(netData, selfLoop, netMatR[ii,,]*edgeScale, arrowScale, 3-sign(netMat[ii,,]*netMatR[ii,,]), sprintf("Network with lag = %s", ii)) 826} # for (ii in 1:nLags) 827 828if (nLags>1) { # overall network 829 logF <- log(netCMatF)*(netCMatP<=pThresh) 830 # if path +, col<-2 (red); if path -, col<-4 (blue) 831 plotNet(logF, selfLoop, logF*edgeScale/10, arrowScale, 3-sign(logF), sprintf("Overall network for all %s lags (scaled by ln(F))", nLags)) 832} # if (nLags>1)) 833 834} # if (plotNetwork) 835 836 837print("#++++++++++++++++++++++++++++++++++++++++++++") 838anotherPth <- as.integer(readline("Want to try another p-threshold/plotting set-up for network (0: no; 1: yes)? ")) 839} 840print("#++++++++++++++++++++++++++++++++++++++++++++") 841anotherLag <- as.integer(readline("Want to try another number of lags for VAR (0: no; 1: yes)? ")) 842} 843print("#++++++++++++++++++++++++++++++++++++++++++++") 844 845#causality(fm, cause=names(myData)[1]) 846 847anotherAna <- as.integer(readline("Next (0: quit; 1: another individual analysis; 2: group analysis)? ")) 848} # while (anotherAna==1) 849 850} # if (anaType==1) 851 852# Group analysis below 853 854if (anaType==2 || anotherAna==2) { 855 856print("It's suggested to run a separate group analysis for each lag.") 857print("-----------------") 858print("The following from subjects can be taken for group analysis:") 859print("(1) Both path coefficients and their t-statistics (highly recommended); or") 860print("(2) Path coefficients only (OK but NOT as reliable as (1); or") 861print("(3) t-statistics only (NOT recommended); or") 862print("(4) F-statistics only (NOT recommended).") 863print("-----------------") 864 865grpType <- as.integer(readline("Which type of input files (1: path and t; 2: path only; 3: t only; 4: F only)? ")) 866 867#fishConv <- TRUE # Fisher transformation 868#if(fishConv) libLoad("psych") 869 870doneGrp <- 1 871while (doneGrp) { 872 873nGrp <- as.integer(readline("Number of groups (1 or 2)? ")) 874if (nGrp==1) { #one-sample t 875 876nSubjs <- as.integer(readline("Number of subjects (e.g., 12)? ")) # number of subjects 877#gfn <- vector('list', nSubjs) 878pathList <- vector('list', nSubjs) 879 880if (grpType==1) { # Fisher tranformation not done since it doesn't seem to matter much in testing 881 print("Provide path coef files first:") 882 print("-----------------") 883 pathList <- readMultiFiles(nSubjs, 2, "subject path coef matrix") #2: with header 884 print("-----------------") 885 print("Now provide t-statistic files:") 886 pathTList <- readMultiFiles(nSubjs, 2, "subject t-statistic matrix") # get t values 887 pathSEList <- mapply("/", pathList, pathTList, SIMPLIFY = FALSE) # calculate stand. error 888 zMat <- do.call(rbind, lapply(lapply(pathList, as.matrix), c)) # convert to matrix 889 zSEMat <- do.call(rbind, lapply(lapply(pathSEList, as.matrix), c)) # convert se to matrix 890 runMeta <- TRUE 891} 892 893if (grpType==2) { # Fisher tranformation not done since it doesn't seem to matter much in testing 894 pathList <- readMultiFiles(nSubjs, 2, "subject path coef matrix") #2: with header 895 zMat <- do.call(rbind, lapply(lapply(pathList, as.matrix), c)) # convert to matrix 896 runMeta <- FALSE 897} 898 899if (grpType==3) { # added here simply because other people still practice this approach 900 pathList <- readMultiFiles(nSubjs, 2, "subject t-statistic matrix") #2: with header 901 pathList <- lapply(pathList, function(x) log(x^2)) # convert to F, then take log transform 902 zMat <- do.call(rbind, lapply(lapply(pathList, as.matrix), c)) # convert to matrix 903 runMeta <- FALSE 904} 905 906if (grpType==4) { # added here simply because other people still practice this approach 907 pathList <- readMultiFiles(nSubjs, 2, "subject F-statistic matrix") #2: with header 908 pathList <- lapply(pathList, log) # take log transform 909 zMat <- do.call(rbind, lapply(lapply(pathList, as.matrix), c)) # convert to matrix 910 runMeta <- FALSE 911} 912 913#for (ii in 1:nSubjs) { # read in path matrices 914# fn[[ii]] <- tclvalue( tkgetOpenFile( filetypes = 915# "{{Path coef matrix file} {.1D}} {{All files} *}", 916# title = paste('Choose subject No.', ii, 'subject path coef matrix file'))) 917# print(sprintf("No. %i file just read in: %s", ii, strsplit(fn[[ii]], "/")[[1]][length(strsplit(fn[[ii]], "/")[[1]])])) 918# pathList[[ii]] <- read.table(fn[[ii]], header=TRUE) 919#} 920 921nROIsG <- dim(pathList[[1]])[1] 922roiNames <- names(pathList[[1]]) 923 924if (runMeta) { 925# libLoad("meta") 926# if(fishConv) {zMat <- fisherz(zMat); zSEMat <- fisherz(zSEMat)} 927 resList <- vector('list', nROIsG^2) # memory allocation 928# for (ii in 1:nROIsG^2) resList[[ii]] <- metagen(zMat[,ii], zSEMat[,ii]) # meta analysis 929# for (ii in 1:nROIsG^2) resList[[ii]] <- rma(zMat[,ii], zSEMat[,ii]^2) # meta analysis 930 for (ii in 1:nROIsG^2) resList[[ii]] <- rma(zMat[,ii], zSEMat[,ii]^2, method="REML") # meta analysis 931# zOutList <- lapply(lapply(resList, summary), function(x) x$random$TE) # extract group z-score 932 zOutList <- lapply(resList, function(x) as.numeric(x$b)) 933# pList <- lapply(lapply(resList, summary), function(x) x$random$p) # extract p 934 pList <- lapply(resList, function(x) as.numeric(x$zp)) # extract p 935# qList <- lapply(resList, function(x) x$QE) # Q for homogeneity test 936 qList <- lapply(resList, function(x) x$QEp) # p-value of Q for homogeneity tes 937 grpR <- matrix(unlist(zOutList), nrow=nROIsG, ncol=nROIsG, dimnames = list(roiNames, roiNames)) 938 grpP <- matrix(unlist(pList), nrow=nROIsG, ncol=nROIsG, dimnames = list(roiNames, roiNames)) 939 grpQ <- matrix(unlist(qList), nrow=nROIsG, ncol=nROIsG, dimnames = list(roiNames, roiNames)) 940 941} else { # no meta analysis: one-sample t-test 942# resList <- apply(do.call(rbind, lapply(lapply(zList, as.matrix), c)), 2, t.test) 943# if(fishConv) zMat <- fisherz(zMat) 944# for (ii in 1:nROIsG^2) resList[[ii]] <- t.test(zMat[,ii]) 945 resList <- apply(zMat, 2, t.test) 946 zOutList <- lapply(resList, function(x) as.numeric(x$estimate)) 947 pList <- lapply(resList, function(x) as.numeric(x$p.value)) 948 grpR <- matrix(unlist(zOutList), nrow=nROIsG, ncol=nROIsG, dimnames = list(roiNames, roiNames)) 949# if(fishConv) grpR <- fisherz2r(grpZ) # convert z to r 950 grpP <- matrix(unlist(pList), nrow=nROIsG, ncol=nROIsG, dimnames = list(roiNames, roiNames)) 951} # if (runMeta) 952 print("-----------------") 953 print("Group path matrix (direction goes from row to column):") 954 print(grpR) 955 print("-----------------") 956 saveRMat <- as.integer(readline("Save the above path matrix (0: no; 1: yes)? ")) 957 if (saveRMat) { 958 matRName <- as.character(readline("File name prefix for path matrix? ")) 959 write.table(grpR, file=sprintf("%s.1D", matRName), append=FALSE, row.names=TRUE, col.names=TRUE) 960 } 961 print("-----------------") 962 print("Group p matrix (direction goes from row to column):") 963 print(grpP) 964 print("-----------------") 965 savePMat <- as.integer(readline("Save the above p matrix (0: no; 1: yes)? ")) 966 if (savePMat) { 967 matPName <- as.character(readline("File name prefix for p matrix? ")) 968 write.table(grpP, file=sprintf("%s.1D", matPName), append=FALSE, row.names=TRUE, col.names=TRUE) 969 } 970 print("-----------------") 971 if (runMeta) { 972 print("p-value matrix for cross-subject heterogeneity/variability (direction goes from row to column):") 973 print(grpQ) 974 } 975 print("-----------------") 976 anotherPthG <- TRUE 977 while (anotherPthG) { 978 979 pThreshG <- as.numeric(readline("p-threshold for group analysis (e.g., 0.05)? ")) 980 surviveR <- as.numeric(grpP<=pThreshG)*grpR # for network plotting 981 showSigR <- matrix(mapply(function(x,y) ifelse(x, y, NA), grpP<=pThreshG, grpR), nrow=nROIsG, ncol=nROIsG, 982 dimnames = list(roiNames, roiNames)) 983#surviveP <- as.numeric(grpP<=pThreshG)*grpP 984 showSigP <- matrix(mapply(function(x,y) ifelse(x, y, NA), grpP<=pThreshG, grpP), nrow=nROIsG, ncol=nROIsG, 985 dimnames = list(roiNames, roiNames)) 986 if (runMeta) showSigQ <- matrix(mapply(function(x,y) ifelse(x, y, NA), grpQ<=pThreshG, grpQ), nrow=nROIsG, 987 ncol=nROIsG, dimnames = list(roiNames, roiNames)) 988 print("Group path matrix with insignificant ones masked with NAs:") 989 print(showSigR) 990 print("-----------------") 991 print("Group p matrix with insignificant ones masked with NAs:") 992 print(showSigP) 993 print("-----------------") 994 if (runMeta) { 995 print("p-value matrix for cross-subject hetetrogeneity/variability with insignificant ones masked with NAs:") 996 print(showSigQ) 997 } 998 print("-----------------") 999 1000 plotNetG <- as.integer(readline("Plot out the identified network (0: no; 1: yes)? ")) 1001 if (plotNetG) { 1002 print("The thickness of a path indicates its relative strength (or effect).") 1003 print("A path in red means positive strength (or effect) while blue is the opposite.") 1004 set.seed(1702) 1005 net <- network.initialize(nROIsG) 1006 attr(net, "vertex.names") <- names(pathList[[1]]) 1007 edgeScaleG <- as.numeric(readline("Scale factor for path thickness (e.g., 1)? ")) 1008 arrowScaleG <- as.numeric(readline("Scale factor for arrow size (e.g., 2)? ")) 1009 selfLoop <- as.integer(readline("Show self-loops in the network (0: no; 1: yes)? ")) 1010 plotNet(surviveR, selfLoop, surviveR*edgeScaleG, arrowScaleG, 3-sign(surviveR), sprintf("Group network with %s subjects", nSubjs)) 1011 } # if (plotNetG) 1012 1013 anotherPthG <- as.integer(readline("Want to try another p-threshold/plotting set-up for group network (0: no; 1: yes)? ")) 1014} 1015} else { # more than one group 1016 1017nSubjs <- vector('integer', nGrp) 1018pathList <- vector('list', nGrp) 1019pathTList <- vector('list', nGrp) 1020pathSEList <- vector('list', nGrp) 1021nROIsG <- vector('integer', nGrp) 1022roiNames <- vector('list', nGrp) 1023zList <- vector('list', nGrp) 1024zMat <- vector('list', nGrp) 1025zSEMat <- vector('list', nGrp) 1026DF <- vector('list', nGrp) 1027 1028for (ii in 1:nGrp) { # Fisher tranformation not done since it doesn't seem to matter much in testing 1029 nSubjs[ii] <- as.integer(readline(sprintf("Number of subjects in group %s (e.g., 12)? ", ii))) 1030 pathList[[ii]] <- vector('list', nSubjs[ii]) 1031 if (grpType==1 ) { 1032 print("Provide path coef files first:") 1033 print("-----------------") 1034 pathList[[ii]] <- readMultiFiles(nSubjs[ii], 2, "subject path coef matrix") 1035 print("-----------------") 1036 print("Now provide t-statistic files in the SAME subject order as path coef's:") 1037 print("-----------------") 1038 pathTList[[ii]] <- readMultiFiles(nSubjs[ii], 2, "subject t-statistic matrix") 1039 pathSEList[[ii]] <- mapply("/", pathList[[ii]], pathTList[[ii]], SIMPLIFY = FALSE) # calculate stand. error 1040 zMat[[ii]] <- do.call(rbind, lapply(lapply(pathList[[ii]], as.matrix), c)) # convert to matrix 1041 zSEMat[[ii]] <- do.call(rbind, lapply(lapply(pathSEList[[ii]], as.matrix), c)) # convert se to matrix 1042 runMeta <- TRUE 1043# print("-----------------") 1044# DF[[ii]] <- as.integer(readline(sprintf("Degrees of freedom for the t-statistics in group %s? ", ii))) 1045 print("-----------------") 1046 } 1047 if (grpType==2) { # Fisher tranformation not done here since it doesn't seem to matter much in testing 1048 pathList[[ii]] <- readMultiFiles(nSubjs[ii], 2, "subject path coef matrix") #2: with header 1049 zMat[[ii]] <- do.call(rbind, lapply(lapply(pathList[[ii]], as.matrix), c)) # convert to matrix 1050 runMeta <- FALSE 1051 } 1052 1053 if (grpType==3) { # added here simply because other people still practice this approach 1054 pathList[[ii]] <- readMultiFiles(nSubjs, 2, "subject t-statistic matrix") #2: with header 1055 pathList[[ii]] <- lapply(pathList[[ii]], function(x) log(x^2)) # convert to F, then take log transform 1056 zMat[[ii]] <- do.call(rbind, lapply(lapply(pathList[[ii]], as.matrix), c)) # convert to matrix 1057 runMeta <- FALSE 1058 } 1059 1060 if (grpType==4) { # added here simply because other people still practice this approach 1061 pathList[[ii]] <- readMultiFiles(nSubjs, 2, "subject F-statistic matrix") #2: with header 1062 pathList[[ii]] <- lapply(pathList[[ii]], log) # take log transform 1063 zMat[[ii]] <- do.call(rbind, lapply(lapply(pathList[[ii]], as.matrix), c)) # convert to matrix 1064 runMeta <- FALSE 1065} 1066 1067 nROIsG[ii] <- dim(pathList[[ii]][[1]])[1] 1068 if (ii>1) for (jj in 1:ii) if (nROIsG[jj]!=nROIsG[1]) {sprintf("Group %i has %i ROIs while group 1 has %i", jj, nROIsG[jj], nROIsG[1]); break} 1069 roiNames[[ii]] <- names(pathList[[ii]][[1]]) 1070 if (ii>1) for (jj in 1:ii) if (!identical(roiNames[[ii]], roiNames[[1]])) {sprintf("Group %i has different ROI namess with group 1", jj, roiNames[[ii]]); break} 1071# if (pthType==0) zList[[ii]] <- lapply(pathList[[ii]], fisherz) # sapply(pathList[[ii]], fisherz, simplify = FALSE) 1072# if (pthType==1) zList[[ii]] <- lapply(pathList[[ii]], function(x) log(x^2)) 1073# if (pthType==2) zList[[ii]] <- lapply(pathList[[ii]], log) 1074# zMat[[ii]] <- do.call(rbind, lapply(lapply(zList[[ii]], as.matrix), c)) # convert from list to matrix for easier handling when running 2-sample t-test 1075} 1076 1077# Instead of looping, we can also use the following aggregation approach: No, it doesn't work!!! 1078# resList <- apply(do.call(rbind, lapply(lapply(zList[[1]], as.matrix), c)), 2, t.test, do.call(rbind, lapply(lapply(zList[[2]], as.matrix), c))) 1079 1080 1081if (runMeta) { 1082# libLoad("meta") 1083# if(fishConv) {zMat <- lapply(zMat, fisherz); zSEMat <- lapply(zSEMat, fisherz)} 1084 resList <- vector('list', nROIsG[1]^2) # memory allocation 1085# for (ii in 1:nROIsG[1]^2) resList[[ii]] <- metacont(rep.int(DF[[1]]+1, 19), zMat[[1]][,ii], zSEMat[[1]][,ii], rep.int(DF[[2]]+1, 18), zMat[[2]][,ii], zSEMat[[2]][,ii]) # meta analysis 1086 for (ii in 1:nROIsG[1]^2) resList[[ii]] <- rma(c(zMat[[1]][,ii], zMat[[2]][,ii]), c(zSEMat[[1]][,ii]^2, zSEMat[[2]][,ii]^2), 1087 mods=c(rep(0, nSubjs[1]), rep(1, nSubjs[2])), method="REML") # meta analysis: REML is selected here! Group2-Group1 1088 1089 zOutList <- lapply(resList, function(x) as.numeric(x$b[2])) # extract effect of group2-group1 1090 z1OutList <- lapply(resList, function(x) as.numeric(x$b[1])) # extract effect of group1, intercept 1091 z2OutList <- mapply("+", zOutList, z1OutList, SIMPLIFY = FALSE) # effect of group2 = sum of two above 1092 pList <- lapply(resList, function(x) as.numeric(x$zp[2])) # extract p for effect of group2-group1 1093 p1List <- lapply(resList, function(x) as.numeric(x$zp[1])) # extract p for group1 effect 1094 1095 varList <- mapply("^", lapply(resList, function(x) as.numeric(x$se[2])), 2, SIMPLIFY = FALSE) # variance of group2-group1 1096 var1List <- mapply("^", lapply(resList, function(x) as.numeric(x$se[1])), 2, SIMPLIFY = FALSE) # variance of group1 1097 # se(grp2)=sqrt(se(grp2-grp1)^2-se(grp1)^2) 1098 zval2List <- mapply("/", z2OutList, lapply(mapply("-", varList, var1List, SIMPLIFY = FALSE), sqrt), SIMPLIFY = FALSE) # z-value of group2 1099 p2List <- lapply(zval2List, function(x) 2*(1-pnorm(abs(x)))) # p-value of group2 1100 qList <- lapply(resList, function(x) x$QEp) # p-value of Q for homogeneity test 1101 grpR <- matrix(unlist(zOutList), nrow=nROIsG[1], ncol=nROIsG[1], dimnames = list(roiNames[[1]], roiNames[[1]])) 1102 grp1R <- matrix(unlist(z1OutList), nrow=nROIsG[1], ncol=nROIsG[1], dimnames = list(roiNames[[1]], roiNames[[1]])) 1103 grp2R <- matrix(unlist(z2OutList), nrow=nROIsG[1], ncol=nROIsG[1], dimnames = list(roiNames[[1]], roiNames[[1]])) 1104 grpP <- matrix(unlist(pList), nrow=nROIsG[1], ncol=nROIsG[1], dimnames = list(roiNames[[1]], roiNames[[1]])) 1105 grp1P <- matrix(unlist(p1List), nrow=nROIsG[1], ncol=nROIsG[1], dimnames = list(roiNames[[1]], roiNames[[1]])) 1106 grp2P <- matrix(unlist(p2List), nrow=nROIsG[1], ncol=nROIsG[1], dimnames = list(roiNames[[1]], roiNames[[1]])) 1107 grpQ <- matrix(unlist(qList), nrow=nROIsG[1], ncol=nROIsG[1], dimnames = list(roiNames[[1]], roiNames[[1]])) 1108 1109} else { # no meta analysis 1110 # two-sample t-test: group2 - group1 1111 resList <- vector('list', nROIsG[1]^2) 1112# want equal variance assumption, otherwise DF would be different from one test to another! 1113 for (ii in 1:nROIsG[1]^2) resList[[ii]] <- t.test(zMat[[2]][,ii], zMat[[1]][,ii], paired=FALSE, var.equal=TRUE) # 2nd-1st 1114 zOutList <- lapply(resList, function(x) as.numeric(x$estimate)) 1115 pList <- lapply(resList, function(x) as.numeric(x$p.value)) 1116 grpR <- matrix(unlist(zOutList), nrow=nROIsG[1], ncol=nROIsG[1], dimnames = list(roiNames[[1]], roiNames[[1]])) 1117 grpP <- matrix(unlist(pList), nrow=nROIsG[1], ncol=nROIsG[1], dimnames = list(roiNames[[1]], roiNames[[1]])) 1118} # if (runMeta) 1119 1120 print("NOTE!!! The results here are about the contrast of the two groups: Group2 - Group1.") 1121 print("That is, they should be interpreted as the amount Group2 is greater than Group1.") 1122 print("They make sense only when comparing to the two individual networks from the two groups.") 1123 1124 print("-----------------") 1125 if (runMeta) { 1126 print("Group1 path matrix (direction goes from row to column):") 1127 print(grp1R) 1128 print("-----------------") 1129 saveRMat <- as.integer(readline("Save Group1 path matrix (0: no; 1: yes)? ")) 1130 if (saveRMat) { 1131 matRName <- as.character(readline("File name prefix for path matrix? ")) 1132 write.table(grp1R, file=sprintf("%s.1D", matRName), append=FALSE, row.names=TRUE, col.names=TRUE) 1133 } 1134 print("-----------------") 1135 print("Group2 path matrix (direction goes from row to column):") 1136 print(grp2R) 1137 print("-----------------") 1138 saveRMat <- as.integer(readline("Save Group2 path matrix (0: no; 1: yes)? ")) 1139 if (saveRMat) { 1140 matRName <- as.character(readline("File name prefix for path matrix? ")) 1141 write.table(grp2R, file=sprintf("%s.1D", matRName), append=FALSE, row.names=TRUE, col.names=TRUE) 1142 } 1143 print("-----------------") 1144 } 1145 print("Group2-Group1 path matrix (direction goes from row to column):") 1146 print(grpR) 1147 print("-----------------") 1148 saveRMat <- as.integer(readline("Save Group2-Group1 path matrix (0: no; 1: yes)? ")) 1149 if (saveRMat) { 1150 matRName <- as.character(readline("File name prefix for path matrix? ")) 1151 write.table(grpR, file=sprintf("%s.1D", matRName), append=FALSE, row.names=TRUE, col.names=TRUE) 1152 } 1153 print("-----------------") 1154 if (runMeta) { 1155 print("Group1 p matrix (direction goes from row to column):") 1156 print(grp1P) 1157 print("-----------------") 1158 savePMat <- as.integer(readline("Save Group1 p matrix (0: no; 1: yes)? ")) 1159 if (savePMat) { 1160 matPName <- as.character(readline("File name prefix for p matrix? ")) 1161 write.table(grp1P, file=sprintf("%s.1D", matPName), append=FALSE, row.names=TRUE, col.names=TRUE) 1162 } 1163 print("-----------------") 1164 print("Group2 p matrix (direction goes from row to column):") 1165 print(grp2P) 1166 print("-----------------") 1167 savePMat <- as.integer(readline("Save Group2 p matrix (0: no; 1: yes)? ")) 1168 if (savePMat) { 1169 matPName <- as.character(readline("File name prefix for p matrix? ")) 1170 write.table(grp2P, file=sprintf("%s.1D", matPName), append=FALSE, row.names=TRUE, col.names=TRUE) 1171 } 1172 print("-----------------") 1173 } 1174 print("Group2-Group1 p matrix (direction goes from row to column):") 1175 print(grpP) 1176 print("-----------------") 1177 savePMat <- as.integer(readline("Save Group2-Group1 p matrix (0: no; 1: yes)? ")) 1178 if (savePMat) { 1179 matPName <- as.character(readline("File name prefix for p matrix? ")) 1180 write.table(grpP, file=sprintf("%s.1D", matPName), append=FALSE, row.names=TRUE, col.names=TRUE) 1181 } 1182 print("-----------------") 1183 if (runMeta) { 1184 print("p-value matrix for cross-subject hetetrogeneity/variability (direction goes from row to column):") 1185 print(grpQ) 1186 } 1187 print("-----------------") 1188 1189 anotherPthG <- TRUE 1190 while (anotherPthG) { 1191 1192 pThreshG <- as.numeric(readline("Two-tailed p-threshold for group analysis (e.g., 0.05)? ")) 1193 surviveR <- as.numeric(grpP<=pThreshG)*grpR # for plotting 1194 showSigR <- matrix(mapply(function(x,y) ifelse(x, y, NA), grpP<=pThreshG, grpR), nrow=nROIsG[1], 1195 ncol=nROIsG[1], dimnames = list(roiNames[[1]], roiNames[[1]])) 1196 showSigP <- matrix(mapply(function(x,y) ifelse(x, y, NA), grpP<=pThreshG, grpP), nrow=nROIsG[1], 1197 ncol=nROIsG[1], dimnames = list(roiNames[[1]], roiNames[[1]])) 1198 if (runMeta) { 1199 survive1R <- as.numeric(grp1P<=pThreshG)*grp1R 1200 survive2R <- as.numeric(grp2P<=pThreshG)*grp2R 1201 showSig1R <- matrix(mapply(function(x,y) ifelse(x, y, NA), grp1P<=pThreshG, grp1R), nrow=nROIsG[1], 1202 ncol=nROIsG[1], dimnames = list(roiNames[[1]], roiNames[[1]])) 1203 showSig2R <- matrix(mapply(function(x,y) ifelse(x, y, NA), grp2P<=pThreshG, grp2R), nrow=nROIsG[1], 1204 ncol=nROIsG[1], dimnames = list(roiNames[[1]], roiNames[[1]])) 1205 showSig1P <- matrix(mapply(function(x,y) ifelse(x, y, NA), grp1P<=pThreshG, grp1P), nrow=nROIsG[1], 1206 ncol=nROIsG[1], dimnames = list(roiNames[[1]], roiNames[[1]])) 1207 showSig2P <- matrix(mapply(function(x,y) ifelse(x, y, NA), grp2P<=pThreshG, grp2P), nrow=nROIsG[1], 1208 ncol=nROIsG[1], dimnames = list(roiNames[[1]], roiNames[[1]])) 1209 showSigQ <- matrix(mapply(function(x,y) ifelse(x, y, NA), grpQ<=pThreshG, grpQ), nrow=nROIsG[1], 1210 ncol=nROIsG[1], dimnames = list(roiNames[[1]], roiNames[[1]])) 1211 print("Group1 path matrix with insignificant path masked with NAs:") 1212 print(showSig1R) 1213 print("-----------------") 1214 print("Group2 path matrix with insignificant path masked with NAs:") 1215 print(showSig2R) 1216 print("-----------------") 1217 } 1218 print("Group2-Group1 path matrix with insignificant path differences masked with NAs:") 1219 print(showSigR) 1220 if (!runMeta) print(sprintf("DFs = %i", round(as.numeric(resList[[1]]$parameter)))) 1221 print("-----------------") 1222 if (runMeta) { 1223 print("Group1 p matrix with insignificant path masked with NAs:") 1224 print(showSig1P) 1225 print("-----------------") 1226 print("Group2 p matrix with insignificant path masked with NAs:") 1227 print(showSig2P) 1228 print("-----------------") 1229 } 1230 print("Group2-Group1 p matrix with insignificant path differences masked with NAs:") 1231 print(showSigP) 1232 print("-----------------") 1233 if (runMeta) { 1234 print("p-value matrix for cross-subject hetetrogeneity/variability with insignificant ones masked with NAs:") 1235 print(showSigQ) 1236 } 1237 print("-----------------") 1238 1239 plotNetG <- as.integer(readline("Plot out the identified network (0: no; 1: yes)? ")) 1240 if (plotNetG) { 1241 print("The thickness of a path indicates the difference magnitude.") 1242 print("A path in red means positive difference while blue is the opposite.") 1243 set.seed(1702) 1244 net <- network.initialize(nROIsG[1]) 1245 attr(net, "vertex.names") <- names(pathList[[1]][[1]]) 1246 edgeScaleG <- as.numeric(readline("Scale factor for path thickness (e.g., 1)? ")) 1247 arrowScaleG <- as.numeric(readline("Scale factor for arrows (e.g., 2)? ")) 1248 selfLoop <- as.integer(readline("Show self-loops in the network (0: no; 1: yes)? ")) 1249 # if path +, col<-2 (red); if path -, col<-4 (blue) 1250 if (runMeta) { 1251 plotNet(survive1R, selfLoop, survive1R*edgeScaleG, arrowScaleG, 3-sign(survive1R), "Network of Group1") 1252 plotNet(survive2R, selfLoop, survive2R*edgeScaleG, arrowScaleG, 3-sign(survive2R), "Network of Group2") 1253 } 1254 plotNet(surviveR, selfLoop, surviveR*edgeScaleG, arrowScaleG, 3-sign(surviveR), "Network of Group2-Group1") 1255 } # if (plotNetG) 1256 1257anotherPthG <- as.integer(readline("Want to try another p-threshold/plotting set-up for group network (0: no; 1: yes)? ")) 1258} # while (anotherPthG) 1259 1260} # more than one group 1261 1262doneGrp <- as.integer(readline("Next (0: done; 1: another group analysis)? ")) 1263} # while (doneGrp) 1264 1265} # if (anaType==2 || anotherAna==2) 1266 1267if (anaType==0 || doneGrp==0 || anotherAna==0) { 1268 print("Make sure to save all desirable figures before leaving!") 1269 quitMe <- as.integer(readline("Quit R (0: no; 1: yes)? ")) 1270 print("***********Thanks for using the program!***********") 1271 print("Any feedback will be welcome - Gang Chen (gangchen@mail.nih.gov)") 1272 if (quitMe) {dev.off(); q()} 1273} 1274print("#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++") 1275