#!/usr/bin/env AFNI_Batch_R #iplots are quite cool, # but when it comes to plotting multiple time series (say 576 x 42) # it gets reaaaaal slow. JAVA SUCKS. #tkrplot has the stupid feature of not allowing autoresize of the plot. # you can do it manually by setting hscale and vscale, but I'd rather # play with a copperhead instead #Clean up rm(list = ls()) first.in.path <- function(file) { ff <- paste(strsplit(Sys.getenv('PATH'),':')[[1]],'/', file, sep='') ff<-ff[lapply(ff,file.exists)==TRUE]; #cat('Using ', ff[1],'\n'); return(gsub('//','/',ff[1], fixed=TRUE)) } source(first.in.path('AFNIio.R')) source(first.in.path('AFNIplot.R')) ExecName <- 'PTA' greeting.PTA <- function () return( "#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ================== Welcome to PTA.R ================== Profile Tracking Analysis #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++" ) reference.PTA <- function () return( " Gang Chen (SSCC/NIMH/NIH)" ) #The help function for PTA batch (command line mode) help.PTA.opts <- function (params, alpha = TRUE, itspace=' ', adieu=FALSE) { intro <- ' ================== Welcome to PTA ================== Program for Profile Tracking Analysis (PTA) #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ Version 0.0.5, Oct 11, 2023 Author: Gang Chen (gangchen@mail.nih.gov) Website - https://afni.nimh.nih.gov/gangchen_homepage SSCC/NIMH, National Institutes of Health, Bethesda MD 20892, USA #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ Introduction ------ Profile Tracking Analysis (PTA) estimates nonlinear trajectories or profiles through smoothing splines. Currently the program PTA only works through a command-line scripting mode. Check the examples below: find one close to your specific scenario and use it as a template. The underlying theory is covered in the following paper: Chen, G., Nash, T.A., Cole, K.M., Kohn, P.D., Wei, S.-M., Gregory, M.D., Eisenberg, D.P., Cox, R.W., Berman, K.F., Shane Kippenhan, J., 2021. Beyond linearity in neuroimaging: Capturing nonlinear relationships with application to longitudinal studies. NeuroImage 233, 117891. https://doi.org/10.1016/j.neuroimage.2021.117891 To be able to run PTA, one needs to have the R packages "mgcv" installed with the following command at the terminal: rPkgsInstall -pkgs "mgcv" Alternatively you may install them in R: install.packages("mgcv") When a factor (e.g, groups, conditions) is involved, numerical coding is required in formulating the data information. See Examples 3 and 4. The following website provides some explanations regarding factor coding that might be useful for modeling formulation: https://stats.idre.ucla.edu/r/library/r-library-contrast-coding-systems-for-categorical-variables/ There are two output files generated by PTA: one (with the affix -stat.txt) contains the information about the statistical evidence for various effects while the other (with the affix -prediction.txt) tabulates the predicted values and their standard errors which can be utilized to illustrate the inferred trajectories or trends (e.g., using graphical tools such as ggplot2 in R). ' ex1 <- "Example 1 --- simplest case: one group of subjects with a between-subject quantitative variable that does not vary within subject. Analysis is set up to model the trajectory or trend along age: PTA -prefix age \\ -input data.txt \\ -model 's(age)' \\ -Y height \\ -prediction pred.txt The function 's(age)' indicates that 'age' is modeled via a smooth curve. No empty space is allowed in the model formulation. The file pred.txt lists all the explanatory variables (excluding lower-level variables such as subject) for prediction. The file should be in a data.frame format as below: age 10 12 14 20 22 24 ... The age step in the above example is 2 years. To obtain smoother graphical appearance in plotted profiles, one can set the age values in pred.txt with a small grid sizer of, for example, 0.5. The file data.txt stores the information for all the variables and input data in a data.frame format as below: Subj age height S1 24 175 S2 14 163 ... The subject labels in the above table can be characters or mixtures of characters and numbers, but they cannot be pure numbers. There will be two output files, one age-stat.txt and the other age-prediction.txt: the former shows the statistical evidence; the latter contains a predicted value for each age plus the associated uncertainty (standard error), which can be plotted using tools such as ggplot2. \n" ex2 <- "Example 2 --- Largely same as Example 1, but with 'age' as a within-subject quantitative variable (varying within each subject). The model is now specified by replacing the line of -model in Example 1 with the following two lines: -model 's(age)+s(Subj,bs=\"re\")' \\ -vt Subj 's(Subj)' \\ The second term 's(Subj,bs=\"re\")' in the model specification means that each subject is allowed to have a varying intercept or random effect ('re'). To estimate the smooth trajectory through the option -prediction, the option -vt has to be included in this case to indicate the varying term (usually subjects). That is, if prediction is desirable, one has to explicitly declare the variable (e.g., Subj) that is associated with the varying term (e.g., s(Subj)). No empty space is allowed in the model formulation and the the varying term. The full script version is PTA -prefix age2 \\ -input data.txt \\ -model 's(age)+s(Subj,bs=\"re\")' \\ -vt Subj 's(Subj)' \\ -prediction pred.txt All the rest remains the same as Example 1. \n" ex3 <- "Example 3 --- two groups and one quantitative variable (age). The analysis is set up to compare the trajectory or trend along age between the two groups, which are quantitatively coded as -1 and 1. For example, if the two groups are females and males, you can code females as -1 and males as 1. The following script applies to the situation when the quantitative variable age does not vary within subject, PTA -prefix age3a \\ -input data.txt \\ -model 's(age)+s(age,by=MvF)' \\ -prediction pred.txt The prediction table in the file data.txt contains the following structure: Subj age grp MvsF S1 27 M 1 S2 21 M 1 S3 28 F -1 S4 18 F -1 ... The column grp above is not necessary for modeling, but it is included to be more indicative for the prediction values in the output file age3a-prediction.txt Similarly, the prediction file pred.txt looks like (set the age values with a small grid so that the graphical illustration would be smooth): age grp MvsF 10 M 1 12 M 1 ... 28 M 1 30 M 1 10 F -1 12 F -1 ... 28 F -1 30 F -1 Note that the age values for prediction have a gap of 2 years: The a smaller the gap, the smoother the plotted predictions. On the other hand, go with the script below when the quantitative variable age varies within subject, PTA -prefix age3b \\ -input data.txt \\ -model 's(age)+s(age,by=grp)+s(Subj,bs=\"re\")' \\ -vt Subj 's(Subj)' \\ -prediction pred.txt \n" ex4 <- "Example 4 --- This example demonstrates the situations where more than two levels are involved in a between-individual factor. Suppose that three groups and one quantitative variable (age). The analysis is set up to compare the trajectory or trend along age between the three groups, A, B and C that are quantitatively represented using dummy coding. PTA -prefix age4a \\ -input data.txt \\ -model 's(age)+s(age,by=AvC)+s(age,by=BvC)' \\ -prediction pred.txt The input table in the file data.txt contains the following structure: Subj age grp AvsC BvC S1 27 A 1 0 S2 21 A 1 0 S3 17 B 0 1 S4 24 B 0 1 S5 28 C 0 0 S6 18 C 0 0 ... The column grp above is not necessary for modeling, but it is included to be more indicative for the prediction values in the output file age4a-prediction.txt On the other hand, go with the script below when the quantitative variable age varies within subject, PTA -prefix age4b \\ -input data.txt \\ -model 's(age)+s(age,by=AvC)+s(age,by=BvC)+s(Subj,bs=\"re\")' \\ -vt Subj 's(Subj)' \\ -prediction pred.txt \n" parnames <- names(params) ss <- vector('character') if (alpha) { parnames <- sort(parnames) ss <- paste('Options in alphabetical order:\n', '------------------------------\n', sep='') } else { ss <- paste('Options:\n', '--------\n', sep='') } for (ii in 1:length(parnames)) { op <- params[parnames[ii]][[1]] if (!is.null(op$help)) { ss <- c(ss , paste(itspace, op$help, sep='')); } else { ss <- c(ss, paste(itspace, parnames[ii], '(no help available)\n', sep='')); } } ex5 <- "Example 5 --- Suppose tht we compare the profiles between two conditions across space or time that is expreessed as a variable x. In this case profile estimation and statistical inference are separated into two steps. First, estimate the profile for each condition using Example 1 or Example 2 as a template. Then, make inference about the contrast between the two conditions. Obtain the contrast at each value of x for each individual, and use the difference values as input. Specify the model as below if there are multiple individuals: -model 's(x)+s(id,bs=\"re\")' \\ -vt id 's(id)' \\ For one individual, change the model to -model 's(x)' \\ \n" ss <- paste(ss, sep='\n'); cat(intro, ex1, ex2, ex3, ex4, ex5, ss, reference.PTA(), sep='\n'); if (adieu) exit.AFNI(); } init.PTA.lop <- function () { lop <- list() lop$input = NULL; lop$interactive = FALSE; lop$pprefix=NULL; lop$cprefix=NULL; lop$select = 'ALL_TASKS'; lop$verb=0 return(lop) } set.both.prefix <- function (lop, pp) { pp <- parse.name(pp) lop$pprefix <- paste(pp$name_noext, "-stat.txt",pp$ext, sep='') lop$cprefix <- paste(pp$name_noext, "-prediction.txt",pp$ext, sep='') return(lop) } #Change command line arguments into an options list read.PTA.opts.batch <- function (args=NULL, verb = 0) { params <- list ( '-input' = apl(n = c(1), d = NA, h = paste( "-input file: input file in a table format (sames as the data frame structure of long", "format in R. Use the first row to specify the column names. The subject column, if", "applicable, should not be purely numbers. On the other hand, factors (groups, tasks)", "should be numerically coded using convenient coding methods such as deviation or", "dummy coding. \n" ) ), '-model' = apl(n = 1, d = 1, h = paste( "-model FORMULA: Specify the model formulation through multilevel smoothing splines", " expression FORMULA with more than one variable has to be surrounded within", " (single or double) quotes. Variable names in the formula should be", " consistent with the ones used in the header of the input file.", " The nonlinear trajectory is specified through the expression of s(x,k=?)", " where s() indicates a smooth function, x is a quantitative variable with", " which one would like to trace the trajectory and k is the number of smooth", " splines (knots). The default (when k is missing) for k is 10, which is good", " enough most of the time when there are more than 10 data points of x. When", " there are less than 10 data points of x, choose a value of k slightly less", " than the number of data points.\n", sep = '\n' ) ), '-Y' = apl(n = 1, d = NA, h = paste( "-Y var_name: var_name is used to specify the column name that is designated as", " as the response/outcome variable. The default (when this option is not", " invoked) is 'Y'.\n", sep = '\n' ) ), '-vt' = apl(n=2, h = paste( "-vt var formulation: This option is for specifying varying smoothing terms. Two components", " are required: the first one 'var' indicates the variable (e.g., subject) around", " which the smoothing will vary while the second component specifies the smoothing", " formulation (e.g., s(age,subject)). When there is no varying smoothing terms (e.g.,", " no within-subject variables), do not use this option.\n", sep='\n')), '-prediction' = apl(n=c(1, 1000000), d=NA, h = paste( "-prediction TABLE: Provide a data table so that predicted values could be generated for", "graphical illustration. Usually the table should contain similar structure as the input", "file except that columns for those varying smoothing terms (e.g., subject) and response", "variable (i.e., Y) should not be included. Try to specify equally-spaced values with a small", "for the quantitative variable of modeled trajectory (e.g., age) so that smooth curves could", "be plotted after the analysis. See Examples in the help for a couple of specific tables used", "for predictions.\n", sep = '\n' ) ), # '-qVars' = apl(n=c(1,100), d=NA, h = paste( # "-qVars variable_list: Identify quantitative variables (or covariates) with", # " this option. The list with more than one variable has to be", # " separated with comma (,) without any other characters such as", # " spaces and should be surrounded within (single or double) quotes.", # " For example, -qVars \"Age,IQ\"", # " WARNINGS:", # " 1) Centering a quantitative variable through -qVarsCenters is", # " very critical when other fixed effects are of interest.", # " 2) Between-subjects covariates are generally acceptable.", # " However EXTREME caution should be taken when the groups", # " differ substantially in the average value of the covariate.", # " \n", # sep = '\n' # ) ), '-prefix' = apl(n = 1, d = NA, h = paste( "-prefix PREFIX: Prefix for output files. \n" ) ), '-interactive' = apl(n = 0, d = NA, h = paste( "-interactive: Currently unavailable.\n" ) ), '-verb' = apl(n=1, d = 0, h = paste( "-verb VERB: VERB is an integer specifying verbosity level.\n", " 0 for quiet (Default). 1 or more: talkative.\n" ) ), '-help' = apl(n=0, h = '-help: this help message\n'), '-h' = apl(n=0, h = '-h: this help message\n'), '-show_allowed_options' = apl(n=0, h= "-show_allowed_options: list of allowed options\n" ), '-dbgArgs' = apl(n=0, h = paste( "-dbgArgs: This option will enable R to save the parameters in a", " file called .PTA.dbg.AFNI.args in the current directory", " so that debugging can be performed.\n", sep='\n')) ); ops <- parse.AFNI.args(args, params, other_ok=FALSE, verb=verb ) if (verb) show.AFNI.args(ops, verb=verb-1, hstr=''); if (is.null(ops)) { errex.AFNI('Error parsing arguments. See PTA -help for details.'); } #Parse dems options #initialize with defaults lop <- init.PTA.lop() lop$model <- NULL #lop$qVars <- NA lop$vt <- NULL #lop$qVarCenters <- NA lop$input <- NULL lop$Y <- 'Y' lop$vt <- NULL lop$prediction <- NULL lop$dbgArgs <- FALSE # for debugging purpose lop$verb <- 0 #Get user's input for (i in 1:length(ops)) { opname <- strsplit(names(ops)[i],'^-')[[1]]; opname <- opname[length(opname)]; switch(opname, input = lop$input <- ops[[i]], prefix = lop <- set.both.prefix(lop,ops[[i]]), pprefix = lop$pprefix <- ops[[i]], cprefix = lop$cprefix <- ops[[i]], interactive = lop$interactive <- TRUE, verb = lop$verb <- ops[[i]], help = help.PTA.opts(params, adieu=TRUE), h = help.PTA.opts(params, adieu=TRUE), show_allowed_options = show.AFNI.args(ops, verb=0, hstr="PTA's",adieu=TRUE), Y = lop$Y <- ops[[i]], vt = lop$vt <- ops[[i]], model = lop$model <- ops[[i]], #qVars = lop$qVars <- ops[[i]], vt = lop$vt <- ops[[i]], #qVarCenters = lop$qVarCenters <- ops[[i]], input = lop$input <- dataTable.AFNI.parse(ops[[i]]), #prediction = lop$prediction <- dataTable.AFNI.parse(ops[[i]]), prediction = lop$prediction <- read.table(ops[[i]], header=TRUE), help = help.MSS.opts(params, adieu=TRUE), dbgArgs = lop$dbgArgs <- TRUE ) } return(lop) }# end of read.PTA.opts.batch #Change options list to PTA variable list process.PTA.opts <- function (lop, verb = 0) { if (is.null(lop$pprefix)) lop$interactive <- TRUE else if(file.exists(lop$pprefix) || file.exists(lop$cprefix)) { errex.AFNI(c("File ", lop$pprefix, " exists! Try a different name.\n")) return(NULL) } return(lop) } imgdev <<- NULL plotdev <<- NULL #The plotting function show_xmat <- function (xmat, isel=1:1:ncol(xmat), descr="", interactive=TRUE) { #cat ('isel in show_xmat') #show(isel) if (length(isel) == 0) { err.AFNI("Nothing to show"); return(0) } #Set colors based on group and offset by 2 so that -1 and 0 get colored colvec <- cg <- attr(xmat,'ColumnGroups') #first color regressors of no interest colvec[cg == -1] = 1 colvec[cg == 0 ] = 2 N_basegroup = 2 #2 types of baseline (RONI) #now color regressors of interest, skip yellow collist = c(3,4,5,6,8) #all but 1st two and yellow N_collist = length(collist) for (i in min(colvec[colvec > 0]):max(colvec[colvec > 0])) { colvec[cg == i] <- collist[(i-N_basegroup) %% N_collist + 1] } if (length(isel) > 10) { text.lym <- colnames(xmat) text.rym <- paste('c', seq(0,ncol(xmat)-1,1), sep='') } else { text.lym <- NULL text.rym <- paste('c', seq(0,ncol(xmat)-1,1), sep='') } descr <- paste( 'Xmat ', attr(xmat,'FileName'), '\n', 'Baseline indices: ', xmat.base.index(xmat)-1, '\n', 'Motion indices: ', xmat.motion.index(xmat)-1, '\n', 'Task indices: ', xmat.alltasks.index(xmat)-1, '\n'); if (length(isel)) { view_cond <- kappa(xmat[,isel], exact=TRUE) }else{ view_cond <- 0 } if (length(xmat.alltasks.index(xmat))) { stim_cond <- kappa(xmat[,xmat.alltasks.index(xmat)], exact=TRUE) } else { stim_cond <- 0 } all_cond <- kappa(xmat, exact=TRUE) if (length(xmat.motion.index(xmat))) { all_cond_no_mot <- kappa(xmat[,-xmat.motion.index(xmat)], exact=TRUE) } else { all_cond_no_mot <- 0 } stit = paste ( 'Rall: ', sprintf('%.2f', all_cond), '', 'Rall-motion: ', sprintf('%.2f', all_cond_no_mot), '', 'Rall-roni: ', sprintf('%.2f', stim_cond),'', 'Rviewed: ', sprintf('%.2f', view_cond),''); maintit <- sprintf('%s\n(%d/%d)', trim.string(attr(xmat,'FileName'), 48), length(isel), ncol(xmat)) maintitcorr <- sprintf('cor(%s)\n(%d/%d)', trim.string(attr(xmat,'FileName'), 48), length(isel), ncol(xmat)) if (!is.null(plotdev) && plotdev != 0) { dev.set(plotdev) } plotdev <<- plot.1D( dmat = xmat, dmat.err=NULL, dmat.colsel = isel, col.ystack = TRUE, ttl.main = maintit, ttl.sub=stit, multi.ncol=1, col.mean.line = FALSE, col.color = colvec, col.plot.char = -1, col.plot.type = 'l', leg.show = FALSE, col.text.lym = text.lym, col.text.lym.at = 'YOFF', col.text.rym = text.rym, col.text.rym.at = 'YOFF', xax.label = 'Time', prefix = lop$pprefix, nodisp = !interactive ) #update report window, if interactive mode if (interactive) { #title (paste('Xmat: ', attr(xmat,'FileName'),'\n')); title (''); ttc <<- condition_report(ttc, xmat, isel, descr=descr) } #Need a device for correlation matrix display imgdev <<- image.corr.1D(dmat = xmat, dmat.colsel = isel, prefix = lop$cprefix, nodisp = !interactive, ttl.main = maintitcorr, dev.this = imgdev, dev.new=TRUE) return(0) } show_xmat_old <- function (xmat, isel=1:1:ncol(xmat), descr="") { #cat ('isel in show_xmat') #show(isel) #Set colors based on group and offset by 2 so that -1 and 0 get colored colvec <- cg <- attr(xmat,'ColumnGroups') #first color regressors of no interest colvec[cg == -1] = 1 colvec[cg == 0 ] = 2 N_basegroup = 2 #2 types of baseline (RONI) #now color regressors of interest, skip yellow collist = c(3,4,5,6,8) #all but 1st two and yellow N_collist = length(collist) for (i in min(colvec[colvec > 0]):max(colvec[colvec > 0])) { colvec[cg == i] <- collist[(i-N_basegroup) %% N_collist + 1] } if (length(isel) > 10) { tp = 'single' #Calculate normalized version for full display ffr <- xmat[,isel] ra = matrix(nrow=ncol(ffr), ncol=3) for (i in 1:ncol(ffr)) { ra[i,] <- as.vector(quantile(ffr[,i], c(0,1,0.5))) } offs <-vector(length = ncol(ffr), mode="numeric") for (i in 1:ncol(ffr)) { ffr[,i] <- (ffr[,i]-ra[i,1])/(ra[i,2]-ra[i,1])+1.4*i offs[i] = mean(ffr[,i]) } if (is.null(dev.list())) x11() plot(ffr, plot.type = tp, xy.labels = colnames(xmat[,isel]), xy.lines = TRUE, panel = lines, nc = 3, yax.flip = FALSE, axes = TRUE, col = colvec[isel], main = '') xoff = ((1:length(isel))%%5)*(length(xmat[,1])/frequency(xmat)/5) text (xoff, offs, colnames(xmat[,isel]), col = colvec[isel], adj=c(0,0)) } else { tp = 'multiple' if (is.null(dev.list())) x11() if (length(isel) == 1) { colm <- colvec[isel] } else { colm <- 1 #ARRRRGH only one col is used for multiple line drawings! } plot(xmat[,isel], plot.type = tp, xy.labels = colnames(xmat[,isel]), xy.lines = TRUE, panel = lines, nc = 1, yax.flip = FALSE, axes = TRUE, col = colm, main = '') } cat ( 'Xmat ', attr(xmat,'FileName'), '\n', 'Baseline indices: ', xmat.base.index(xmat)-1, '\n', 'Motion indices: ', xmat.motion.index(xmat)-1, '\n', 'Task indices: ', xmat.alltasks.index(xmat)-1, '\n'); if (length(isel)) { view_cond <- kappa(xmat[,isel], exact=TRUE) }else{ view_cond <- 0 } if (length(xmat.alltasks.index(xmat))) { stim_cond <- kappa(xmat[,xmat.alltasks.index(xmat)], exact=TRUE) } else { stim_cond <- 0 } all_cond <- kappa(xmat, exact=TRUE) if (length(xmat.motion.index(xmat))) { all_cond_no_mot <- kappa(xmat[,-xmat.motion.index(xmat)], exact=TRUE) } else { all_cond_no_mot <- 0 } stit = paste ( 'Rall : ', sprintf('%.2f', all_cond), '\n', 'Rall-motion : ', sprintf('%.2f', all_cond_no_mot), '\n', 'Rall-roni : ', sprintf('%.2f', stim_cond),'\n', 'Rviewed : ', sprintf('%.2f', view_cond),'\n'); #title (paste('Xmat: ', attr(xmat,'FileName'),'\n')); title (''); #update report window ttc <<- condition_report(ttc, xmat, isel, descr=descr) return(0) } text_window_kill <- function() { tkdestroy(ttc) ttc <<-NA } condition_report <- function(tt=NA, xmat, isel=1:1:ncol(xmat), descr="") { if (is.na(tt)[1]) { #cat ('Have no widget\n') tt <- tktoplevel() #Catch the destroy button to reset ttc to NA tcl( "wm", "protocol", tt, "WM_DELETE_WINDOW", function()text_window_kill()) scr <- tkscrollbar(tt, repeatinterval=5, command=function(...)tkyview(txt,...)) #make sure txt is being kept for next call txt <<- tktext( tt,bg="white",font="courier", yscrollcommand=function(...)tkset(scr,...)) tkgrid(txt,scr) tkgrid.configure(scr,sticky="ns") #Numbers that will stay constant if (length(xmat.alltasks.index(xmat))) { stim_cond <- kappa(xmat[,xmat.alltasks.index(xmat)], exact=TRUE) } else { stim_cond <- 0 } all_cond <- kappa(xmat, exact=TRUE) if (length(xmat.motion.index(xmat))) { all_cond_no_mot <- kappa(xmat[,-xmat.motion.index(xmat)], exact=TRUE) } else { all_cond_no_mot <- all_cond } #attr(xmat,'TR') #attr(xmat,'dim') stitcom = paste ( 'Matrix file: ', attr(xmat,'FileName'), '\n', '* TR: ', attr(xmat,'TR'), '\n', '* N: ', length(xmat[,1]), '\n', '* Baseline indices: ', paste(xmat.base.index(xmat)-1, collapse=", "), '\n', '* Motion indices: ', paste(xmat.motion.index(xmat)-1, collapse=", "), '\n', '* Task indices: ', paste(xmat.alltasks.index(xmat)-1, collapse=", "), '\n', '* Task Labels: ', paste( attr(xmat,'TaskNames'), collapse = ", "),'\n', '* Standard condition numbers:', '\n', ' Rall : ', sprintf('%.2f', all_cond), '\n', ' Rall-motion : ', sprintf('%.2f', all_cond_no_mot), '\n', ' Rall-roni : ', sprintf('%.2f', stim_cond),'\n', '=============================================\n'); } else { #cat ('reusing widget\n') stitcom = paste ('User Selection:', paste(descr, collapse=" ") , '\n') } #cat ('Calculating conditions\n') if (length(isel)) { view_cond <- kappa(xmat[,isel], exact=TRUE) } else { view_cond = 0 } stit = paste ( stitcom, ' Rviewed : ', sprintf('%.2f', view_cond),'\n', '----------------------------------------------\n'); #cat ('inserting string\n') tkconfigure(txt, state="normal") #allow editing tkinsert(txt,"end",stit) tkconfigure(txt, state="disabled") #no more editing tkfocus(txt) return(tt) } # model specification mod_gui <- function (tt) { require(tcltk) Name <- tclVar("") entry.Name <-tkentry(tt,width="40",textvariable=Name) tkgrid(tklabel(tt, text=paste( "Enter model formulation\n", "example: s(age)") ) ) tkgrid(entry.Name) OnOK <- function() { #cat ('OnOK\n') NameVal <- gsub("[,;\'\"]", " ",tclvalue(Name)) #show_sel(NameVal,xmat) tkfocus(tt) } OnDone <- function() { tkdestroy(tt) exit_flag <<- TRUE } OnHelp <- function() { PTA.help() } OK.but <-tkbutton(tt,text=" OK ",command=OnOK) Done.but <- tkbutton(tt,text=" Done ",command=OnDone) Help.but <- tkbutton(tt,text=" Help ",command=OnHelp) tkbind(entry.Name, "", OnOK) tkgrid(OK.but, Done.but, Help.but) tkfocus(tt) exit_flag <<- FALSE } sel_gui <- function (tt, xmat) { require(tcltk) Name <- tclVar("") entry.Name <-tkentry(tt,width="20",textvariable=Name) tkgrid(tklabel(tt, text=paste( "Enter indices or label identifiers\n", "example: 0:3, ", attr(xmat,'TaskNames')[1],", 7 9") ) ) tkgrid(entry.Name) OnOK <- function() { #cat ('OnOK\n') NameVal <- gsub("[,;\'\"]", " ",tclvalue(Name)) show_sel(NameVal,xmat) tkfocus(tt) } OnDone <- function() { tkdestroy(tt) exit_flag <<- TRUE } OnHelp <- function() { PTA.help() } OK.but <-tkbutton(tt,text=" OK ",command=OnOK) Done.but <- tkbutton(tt,text=" Done ",command=OnDone) Help.but <- tkbutton(tt,text=" Help ",command=OnHelp) tkbind(entry.Name, "", OnOK) tkgrid(OK.but, Done.but, Help.but) tkfocus(tt) exit_flag <<- FALSE } PTA.help.select <- function() { sh <- paste ( 'To select regressors, you can use regressor indices or regressor \n', 'labels. For example, say your tasks (and therefore regressors) \n', 'are labeled "house", "face", "random"\n', 'Then to view house and face regressors only, select:\n', '"house, face" or "house face" or "h fa" etc.\n', 'You can also use regressor indices (start at 0)\n', 'So you can select "0, 5, 10:12" for regressors 0, 5, 10, 11, and 12\n', 'You can also combine strings and integers in the same selection.\n', 'Commas, semicolons, and quotes are ignored.\n' ) return(sh) } PTA.help <- function() { sh <- paste ( '===============================================================\n', 'PTA Help:\n', '-----------------\n', 'PTA.R is an interactive tool to examine a design matrix.\n', 'The interface consists of 3 windows: Graph, Control, Results.\n', '\n', 'The Graph window displays the entire design matrix or a subset \n', 'designated by the user.\n', 'The Control window allows users to designate which regressors\n', 'to display.\n', 'The Results window displays properties of the design matrix in \n', 'its entirety and the condition number of the matrix formed \n', 'by the designated (or viewed in the graph) regressors only.\n', '\n', PTA.help.select(), '\n', 'When viewing more than 10 regressors, the regressors get automatically\n', 'rescaled and displayed on the same axis. In this case, the amplitude\n', 'of the Y axis is of no relevance. In this display mode, baseline \n', 'regressors are black, motion are red, and each task group gets a distinct\n', 'color.\n', 'When 10 or less regressors are displayed simultaneously, each one\n', 'gets its own subplot and the y axis represents the regressor amplitude\n', 'as it is in the design matrix. Also, in this display mode, all regressors\n', 'are drawn in black.\n', '===============================================================\n', '\n') cat (sh) #tkmessageBox(message=sh) #ugly... #use results output instead tkconfigure(txt, state="normal") #allow editing tkinsert(txt,"0.0",sh) tkconfigure(txt, state="disabled") #no more editing tkfocus(txt) #tkmessageBox(message="See text in Results window") } show_sel <- function(sel, xmat) { ilst <- xmat.select.indices(sel, xmat) if (length(ilst)) { ilst = unique(ilst[ilst >0 & ilst <=ncol(xmat)]) if (length(ilst)) { cat ('selection is: ', ilst-1, '\n') show_xmat(xmat, ilst, descr=sel) } } } main_loop <- function (fn, select=NULL) { #read file if (is.null(lop$xmat <- read.table(fn, header=T))) { err.AFNI(paste("Failed to load input file", fn)); return(1); } #begin the selection toy, this GUI controls the quitting etc. mod_gui(tktoplevel()) #top level widget ttc <<- NA #Show all of lop$xmatr #browser() require(mgcv) require(ggplot2) isel <- NULL if (!is.null(select)) { isel <- xmat.select.indices(select, lop$xmat) if (length(isel)==0) { warn.AFNI(sprintf("Nothing found for selection %s. Showing all.", select)) isel<-NULL; } } if (is.null(isel)) isel <- 1:1:ncol(lop$xmat) show_xmat(lop$xmat, isel) #begin the selection toy, this GUI controls the quitting etc. sel_gui(tktoplevel(), lop$xmat) #dunno what to do for event loop, all I can do is sleep #Without this, the batch mode exits immediately while (!exit_flag) { Sys.sleep(0.5) } return(0) } ################################################################################# ####################### Begin PTA main ################################## ################################################################################# if (!exists('.DBG_args')) { args = (commandArgs(TRUE)) rfile <- first.in.path(sprintf('%s.R',ExecName)) # save only on -dbg_args if ( '-dbgArgs' %in% args ) { try(save(args, rfile, file=".PTA.dbg.AFNI.args", ascii = TRUE), silent=TRUE) } } else { note.AFNI("Using .DBG_args resident in workspace"); args <- .DBG_args } if (!length(args)) { BATCH_MODE <<- 0 err.AFNI("No parameters"); } else { if (!exists('.DBG_args')) { BATCH_MODE <<- 1 } else { BATCH_MODE <<- 0 } if (is.null(lop <- read.PTA.opts.batch(args, verb = 0))) { stop('Error parsing input'); } if (is.null(lop <- process.PTA.opts(lop, verb = lop$verb))) { stop('Error processing input'); } } if (lop$verb) { str(lop); } if (lop$interactive) { # GUI mode require('graphics') require('tcltk') if (is.null(lop$input)) { lop$input <- tclvalue( tkgetOpenFile( filetypes = "{{Input Files .1D} {.1D}} {{Input Files .txt} {.txt}} {{Input Files .tbl} {.tbl}} {{All files} *}", #title = 'Choose input data file (in data-table format)')) title = 'Currently Not Working!')) } if (nchar(lop$input)) { if (main_loop(lop$input, lop$select)) cat ('error') } } else { # command-line mode if (is.null(lop$dat <- read.table(lop$input, header=TRUE))) { errex.AFNI(paste("Failed to load input file", lop$input)); } lop$model <- as.formula(paste0(lop$Y, '~', lop$model)) nr <- nrow(lop$prediction) if(is.null(lop$vt)) lop$Pred <- lop$prediction else { lop$dat[[lop$vt[1]]] <- as.factor(lop$dat[[lop$vt[1]]]) # declare lop$vt[1] as a factor lop$Pred <- lop$prediction[rep(seq_len(nr), times=nlevels(lop$dat[, lop$vt[1]])), ] lop$Pred[,lop$vt[1]] <- as.factor(rep(levels(lop$dat[,lop$vt[1]]), each = nr)) } require(mgcv) try(fm <- gam(lop$model, data=lop$dat, method='REML'), silent=TRUE) cat('###Parametric effects###', file = lop$pprefix, sep = '\n', append=TRUE) cat(capture.output(summary(fm)$p.table), file = lop$pprefix, sep = '\n', append=TRUE) cat('\n###Nonlinear effects###', file = lop$pprefix, sep = '\n', append=TRUE) cat(capture.output(summary(fm)$s.table), file = lop$pprefix, sep = '\n', append=TRUE) try(pred <- predict(fm, lop$Pred, exclude=lop$vt[2], se.fit = T), silent=TRUE) pp <- data.frame(prediction=pred$fit[1:nr], se=pred$se.fit[1:nr]) lop$prediction <- cbind(lop$prediction, pp) row.names(lop$prediction) <- NULL formals(print.data.frame)$row.names <- FALSE cat(capture.output(print(lop$prediction, row.names = FALSE)), file = lop$cprefix, sep = '\n', append=TRUE) print(sprintf("PTA was performed successfully with output files %s and %s", lop$pprefix, lop$cprefix)) }