#!/usr/bin/env AFNI_Batch_R ##!/usr/bin/env afni_run_R # Command line to run this script: 3dMVM.R dataStr.txt diary.txt & # (Output is a file in which the running progress including # error messages will be stored) 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')) ExecName <- '3dMVM' # Global variables iterPar <- 'matrPar' respVar <- c('InputFile', 'Inputfile', 'Ausgang_val', 'ausgang_val') ################################################################################# ##################### Begin MVM Input functions ################################ ################################################################################# #The help function for 3dMVM batch (AFNI-style script mode) help.MVM.opts <- function (params, alpha = TRUE, itspace=' ', adieu=FALSE) { intro <- ' ================== Welcome to 3dMVM ================== AFNI Group Analysis Program with Multi-Variate Modeling Approach #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ Version 3.2.8, Sept 18, 2014 Author: Gang Chen (gangchen@mail.nih.gov) Website - http://afni.nimh.nih.gov/sscc/gangc/MVM.html SSCC/NIMH, National Institutes of Health, Bethesda MD 20892 #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ Usage: ------ 3dMVM is a group-analysis program that performs traditional ANOVA- and ANCOVA- style computations. In addition, it can run multivariate modeling in the sense of multiple simultaneous response variables. For univariate analysis, no bound is imposed on the numbers of explanatory variables, and these variables can be either categorical (factor) or numerical/quantitative (covariate). F-statistics for all main effects and interactions are automatically included in the output. In addition, general linear tests (GLTs) can be requested via symbolic coding. Input files for 3dMVM can be in AFNI, NIfTI, or surface (niml.dset) format. Note that unequal number of subjects across groups are allowed, but scenarios with missing data for a within-subject factor are better modeled with 3dLME. Cases with quantitative variables (covariates) that vary across the levels of a within-subject variable are also better handled with 3dLME. Computational cost with 3dMVM is very high compared to 3dttest++ or 3dANOVAx, but it has the capability to correct for sphericity violations when within-subject variables with more than two levels are involved. If you want to cite the analysis approach, use the following at this moment: Chen, G., Adleman, N.E., Saad, Z.S., Leibenluft, E., Cox, R.W. (in press). Applications of Multivariate Modeling to Neuroimaging Group Analysis: A Comprehensive Alternative to Univariate General Linear Model. NeuroImage. 10.1016/j.neuroimage.2014.06.027 http://afni.nimh.nih.gov/pub/dist/HBM2014/Chen_in_press.pdf In addition to R installation, the following two R packages need to be acquired in R first before running 3dMVM: install.packages("afex") install.packages("phia") The snow package is also needed if one wants to take advantage of parallel computing: install.packages("snow") Once the 3dMVM command script is constructed, it can be run by copying and pasting to the terminal. Alternatively (and probably better) you save the script as a text file, for example, called MVM.txt, and execute it with the following (assuming on tc shell), tcsh -x MVM.txt & or, tcsh -x MVM.txt > diary.txt & The advantage of the latter command is that the progression is saved into the text file diary.txt and, if anything goes awry, can be examined later. More details about 3dMVM can be found at http://afni.nimh.nih.gov/sscc/gangc/MVM.html Thank the R community and Henrik Singmann in particular for the strong technical support.' ex1 <- "\n-------------------------------- Example 1 --- three between-subjects (genotype, sex, and scanner) and two within-subject (condition and emotion) variables: 3dMVM -prefix Example1 -jobs 4 \\ -bsVars 'genotype*sex+scanner' \\ -wsVars \"condition*emotion\" \\ -wsE2 \\ -num_glt 14 \\ -gltLabel 1 face_pos_vs_neg -gltCode 1 'condition : 1*face emotion : 1*pos -1*neg' \\ -gltLabel 2 face_emot_vs_neu -gltCode 2 'condition : 1*face emotion : 1*pos +1*neg -2*neu' \\ -gltLabel 3 sex_by_condition_interaction -gltCode 3 'sex : 1*male -1*female condition : 1*face -1*house' \\ -gltLabel 4 3way_interaction -gltCode 4 'sex : 1*male -1*female condition : 1*face -1*house emotion : 1*pos -1*neg' \\ ... -dataTable \\ Subj genotype sex scanner condition emotion InputFile \\ s1 TT male scan1 face pos s1+tlrc\'[face_pos_beta]\' \\ s1 TT male scan1 face neg s1+tlrc\'[face_neg_beta]\' \\ s1 TT male scan1 face neu s1+tlrc\'[face_neu_beta]\' \\ s1 TT male scan1 house pos s1+tlrc\'[house_pos_beta]\' \\ ... s68 TN female scan2 house pos s68+tlrc\'[face_pos_beta]\' \\ s68 TN female scan2 house neg s68+tlrc\'[face_neg_beta]\' \\ s68 TN female scan2 house neu s68+tlrc\'[house_pos_beta]\' NOTE: 1) Option -wsE2 is used to combine both the univariate testing and the within-subject multivariate approach. This option only makes sense if a within-subject factor has more than 3 level. 2) The 3rd GLT is for the 2-way 2 x 2 interaction between sex and condition, which is essentially a t-test (or one degree of freedom for the numerator of F-statistic). Multiple degrees of freedom for the numerator of F-statistic is currently unavailable. 3) Similarly, the 4th GLT is a 3-way 2 x 2 x 2 interaction, which is a partial (not full) interaction between the three factors because 'emotion' has three levels. The F-test for the full 2 x 2 x 3 interaction is automatically spilled out by 3dMVM.\n" ex2 <- "-------------------------------- Example 2 --- two between-subjects (genotype and sex), onewithin-subject (emotion) factor, plus two quantitative variables (age and IQ).f 3dMVM -prefix Example2 -jobs 24 \\ -bsVars \"genotype*sex+age+IQ\" \\ -wsVars emotion \\ -wsE2 \\ -qVars \"age,IQ\" \\ -qVarCenters '25,105' \\ -num_glt 10 \\ -gltLabel 1 pos_F_vs_M -gltCode 1 'sex : 1*female -1*male emotion : 1*pos' \\ -gltLabel 2 age_pos_vs_neg -gltCode 2 'emotion : 1*pos -1*neg age :' \\ -gltLabel 3 genotype_by_sex -gltCode 3 'genotype : 1*TT -1*NN sex : 1*male -1*female' \\ -gltLabel 4 genotype_by_sex_emotion -gltCode 4 'genotype : 1*TT -1*NN sex : 1*male -1*female emotion : 1*pos -1*neg' \\ ... -dataTable \\ Subj genotype sex age IQ emotion InputFile \\ s1 TT male 24 107 pos s1+tlrc\'[pos_beta]\' \\ s1 TT male 24 107 neg s1+tlrc\'[neg_beta]\' \\ s1 TT male 24 107 neu s1+tlrc\'[neu_beta]\' \\ ... s63 NN female 29 110 pos s63+tlrc\'[pos_beta]\' \\ s63 NN female 29 110 neg s63+tlrc\'[neg_beta]\' \\ s63 NN female 29 110 neu s63+tlrc\'[neu_beta]\' NOTE: 1) Option -wsE2 is used to combine both the univariate testing and the within-subject multivariate approach. This option only makes sense if a within-subject factor has more than 3 level. 2) The 3rd GLT is for the 2-way 2 x 2 interaction between genotype and sex, which is essentially a t-test (or one degree of freedom for the numerator of F-statistic). Multiple degrees of freedom for the numerator of F-statistic is currently unavailable. 3) Similarly, the 4th GLT is a 3-way 2 x 2 x 2 interaction, which is a partial (not full) interaction between the three factors because 'emotion' has three levels. The F-test for the full 2 x 2 x 3 interaction is automatically spilled out by 3dMVM.\n" ex3 <- "--------------------------------- Example 3 --- BOLD response was modeled with multiple basis functions at individual subject level. In addition, there are one between-subjects (Group) and one within- subject (Condition) variable. Furthermore, the variable corresponding to the number of basis functions, Time, is also a within-subject variable. In the end, the F- statistics for the interactions of Group:Condition:Time, Group:Time, and Condition:Time are of specific interest. And these interactions can be further explored with GLTs in 3dMVM. 3dMVM -prefix Example3 -jobs 12 \\ -bsVars Group \\ -wsVars 'Condition*Time' \\ -num_glt 32 \\ -gltLabel 1 old_t0 -gltCode 1 'Group : 1*old Time : 1*t0' \\ -gltLabel 2 old_t1 -gltCode 2 'Group : 1*old Time : 1*t1' \\ -gltLabel 3 old_t2 -gltCode 3 'Group : 1*old Time : 1*t2' \\ -gltLabel 4 old_t3 -gltCode 4 'Group : 1*old Time : 1*t3' \\ -gltLabel 5 yng_t0 -gltCode 5 'Group : 1*yng Time : 1*t0' \\ -gltLabel 6 yng_t1 -gltCode 6 'Group : 1*yng Time : 1*t1' \\ -gltLabel 7 yng_t2 -gltCode 7 'Group : 1*yng Time : 1*t2' \\ -gltLabel 8 yng_t3 -gltCode 8 'Group : 1*yng Time : 1*t3' \\ ... -gltLabel 17 old_face_t0 -gltCode 17 'Group : 1*old Condition : 1*face Time : 1*t0' \\ -gltLabel 18 old_face_t1 -gltCode 18 'Group : 1*old Condition : 1*face Time : 1*t1' \\ -gltLabel 19 old_face_t2 -gltCode 19 'Group : 1*old Condition : 1*face Time : 1*t2' \\ -gltLabel 20 old_face_t3 -gltCode 20 'Group : 1*old Condition : 1*face Time : 1*t3' \\ ... -dataTable \\ Subj Group Condition Time InputFile \\ s1 old face t0 s1+tlrc\'[face#0_beta]\' \\ s1 old face t1 s1+tlrc\'[face#1_beta]\' \\ s1 old face t2 s1+tlrc\'[face#2_beta]\' \\ s1 old face t3 s1+tlrc\'[face#3_beta]\' \\ ... s40 yng house t0 s40+tlrc\'[house#0_beta]\' \\ s40 yng house t1 s40+tlrc\'[house#1_beta]\' \\ s40 yng house t2 s40+tlrc\'[house#2_beta]\' \\ s40 yng house t3 s40+tlrc\'[house#3_beta]\' NOTE: The model for the analysis can also be set up as and is equivalent to 'Group*Condition*Time'.\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='')) } ss <- paste(ss, sep='\n') cat(intro, ex1, ex2, ex3, ss, sep='\n') if (adieu) exit.AFNI(); } #Change command line arguments into an options list read.MVM.opts.batch <- function (args=NULL, verb = 0) { params <- list ( '-prefix' = apl(n = 1, d = NA, h = paste( "-prefix PREFIX: Output file name. For AFNI format, provide prefix only,", " no view+suffix needed. Filename for NIfTI format should have", " .nii attached, while file name for surface data is expected", " to end with .niml.dset. The sub-brick labeled with the '(Intercept)',", " if present, should be interpreted as the overall average", " across factor levels.\n", sep = '\n' ) ), '-mask' = apl(n=1, d = NA, h = paste( "-mask MASK: Process voxels inside this mask only.\n", " Default is no masking.\n" ) ), '-jobs' = apl(n = 1, d = 1, h = paste( "-jobs NJOBS: On a multi-processor machine, parallel computing will speed ", " up the program significantly.", " Choose 1 for a single-processor computer.\n", sep = '\n' ) ), '-model' = apl(n = 1, d = 1, h = paste( "-model FORMULA: This option will phase out at some point. So use -bsVars", " instead. Specify the fixed effects for between-subjects factors ", " and quantitative variables. When no between-subject factors", " are present, simply put 1 for FORMULA. The 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 -dataTable. A+B represents the", " additive effects of A and B, A:B is the interaction between A", " and B, and A*B = A+B+A:B. The effects of within-subject", " factors, if present under -wsVars are automatically assumed", " to interact with the ones specified here. Subject as a variable", " should not occur in the model specification here.\n", sep = '\n' ) ), '-bsVars' = apl(n = 1, d = 1, h = paste( "-bsVars FORMULA: Specify the fixed effects for between-subjects factors ", " and quantitative variables. When no between-subject factors", " are present, simply put 1 for FORMULA. The expression FORMULA", " with more than one variable has to be surrounded within (single or", " double) quotes. No spaces are allowed in the FORMULA expression.", " Variable names in the formula should be consistent with the ones", " used in the header underneath -dataTable. A+B represents the", " additive effects of A and B, A:B is the interaction between A", " and B, and A*B = A+B+A:B. The effects of within-subject", " factors, if present under -wsVars are automatically assumed", " to interact with the ones specified here. Subject as a variable", " should not occur in the model specification here.\n", sep = '\n' ) ), '-wsVars' = apl(n=c(1,100), d=NA, h = paste( "-wsVars FORMULA: Within-subject factors, if present, have to be listed", " here; otherwise the program will choke. If no within-subject ", " exists, don't include this option in the script. Coding for", " additive effects and interactions is the same as in -bsVars. The", " FORMULA with more than one variable has to be surrounded ", " within (single or double) quotes. Note that the within-subject", " variables are assumed to interact with those between-subjects", " variables specified under -bsVars. The hemodynamic response", " time course are better modeled as simultaneous outcomes through", " option -mVar, and not as the levels of a within-subject factor.", " The varialbes under -wsVars and -mVar are exclusive from each", " other.\n", sep = '\n' ) ), '-wsMVT' = apl(n=0, h = paste( "-wsMVT: If at least one within-subject factor is involved in the model, this", " option provides within-subject multivariate testing for any effect", " associated with a within-subject variable. The testing strategy is", " different from the conventional univariate GLM, see more details in", " Chen et al., ...... If all the within-subject factors", " have two levels, the multivariate testing would render the same", " results as the univariate version. So use the option only if at", " least one within-subject factor has more than two levels.", " The F-statistics from the multivariate testing are labeled", " with -wsMVT- in the sub-brick names. \n", sep='\n')), '-wsE2' = apl(n=0, h = paste( "-wsE2: If at least one within-subject factor is involved in the model, any", " omnibus F-test associated with a within-subject factor is assessed", " with both univariate and within-subject multivariate tests. Use", " the option only if at least one within-subject factor has more", " than two levels.\n", sep='\n')), '-mvE4' = apl(n=0, h = paste( "", sep='\n')), '-mvE4a' = apl(n=0, h = paste( "", sep='\n')), '-parSubset' = apl(n=c(1,100), d=NA, h = paste( "", sep='\n')), '-mVar' = apl(n=c(1,100), d=NA, h = paste( "-mVar variable: With this option, the levels of the within-subject factor", " will be treated as simultaneous variables in a multivariate model.", " For example, when the hemodynamic response time course is modeled", " through multiple basis functions such as TENT, TENTzero, CSPLIN,", " CSPLINzero, SPMG2/3, etc., the effect estimates at the multiple", " time points can be treated as simultaneous response variables in", " a multivariate model. Only one within-subject variable is allowed", " currently under -mVar. In addition, in the presence of -mVar, no", " other within-subject factors should be included. If modeling", " extra within-subject factors with -mVar is desirable, consider", " flattening such factors; that is, perform multiple analyses", " at each level or their contrasts of the factor. The output", " for multivariate testing are labeled with -MV0- in the sub-brick", " names. \n", sep = '\n' ) ), '-SC' = apl(n=0, h = paste( "-SC: If a within-subject factor with more than *two* levels is", " involved in the model, 3dMVM automatically provides the", " F-statistics for main and interaction effects with", " sphericity assumption. If the assumption is violated,", " the F-statistics could be inflated to some extent. This option,", " will enable 3dMVM to additionally output the F-statistics of", " sphericity correction for main and interaction effects, which", " are labeled with -SC- in the sub-brick names.", " NOTE: this option should be used only when at least one", " within-subject factor has more than TWO levesl.\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 significantly in the average value of the covariate.", " 3) Within-subject covariates vary across the levels of a", " within-subject factor, and can be analyzed with 3dLME,", " but not 3dMVM.\n", sep = '\n' ) ), '-vVars' = apl(n=c(1,100), d=NA, h = paste( "-vVars variable_list: Identify voxel-wise covariates with this option.", " Currently one voxel-wise covariate is allowed only, but this", " may change if demand occurs...", # " 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 \"Var1,Var2\". " By default mean centering is performed voxel-wise across all", " subjects. Alternatively centering can be specified through a", " global value under -vVarsCenters. If the voxel-wise covariates", " have already been centered, set the centers at 0 with -vVarsCenters.\n", sep = '\n' ) ), '-qVarCenters' = apl(n=1, d=NA, h = paste( "-qVarCenters VALUES: Specify centering values for quantitative variables", " identified under -qVars. Multiple centers are separated by ", " commas (,) within (single or double) quotes. The order of the", " values should match that of the quantitative variables in -qVars.", " Default (absence of option -qVarsCetners) means centering on the", " average of the variable across ALL subjects regardless their", " grouping. If within-group centering is desirable, center the", " variable YOURSELF first before the values are fed into -dataTable.\n", sep = '\n' ) ), '-vVarCenters' = apl(n=1, d=NA, h = paste( "-vVarCenters VALUES: Specify centering values for voxel-wise covariates", " identified under -vVars. Multiple centers are separated by ", " commas (,) within (single or double) quotes. The order of the", " values should match that of the quantitative variables in -qVars.", " Default (absence of option -vVarsCetners) means centering on the", " average of the variable across ALL subjects regardless their", " grouping. If within-group centering is desirable, center the", " variable YOURSELF first before the files are fed into -dataTable.\n", sep = '\n' ) ), '-num_glt' = apl(n=1, d=0, h = paste( "-num_glt NUMBER: Specify the number of general linear tests (GLTs). A glt", " is a linear combination of a factor levels. See details in ", " -gltLabel.\n", sep = '\n' ) ), '-gltLabel' = apl(n=c(1,1000), d=NA, h = paste( "-gltLabel k label: Specify the label for the k-th general linear test", " (GLT). A symbolic coding for the GLT is assumed to follow with", " each -gltLabel.\n", sep = '\n' ) ), '-gltCode' = apl(n=c(1,1000), d=NA, h = paste( "-gltCode k CODING: Specify the k-th general linear test (GLT) through a", " weighted combination among factor levels. The symbolic coding has", " to be within (single or double) quotes. For example, the following", " 'Condition : 2*House -3*Face Emotion : 1*positive '", " requests for a test of comparing 2 times House condition", " with 3 times Face condition while Emotion is held at positive", " valence.\n", " NOTE:\n", " 1) The weights for a variable do not have to add up to 0.\n", " 2) When a quantitative variable is present, other effects are", " tested at the center value of the covariate.\n", " 3) The effect for a quantitative variable can be specified with,", " for example, 'Group : 1*Old Age : ', or ", " 'Group : 1*Old - 1*Young Age : '\n", " 4) The absence of a categorical variable in a coding means the", " levels of that factor are averaged (or collapsed) for the GLT.\n", " 5) The appearance of a categorical variable has to be followed", " by the linear combination of its levels. Only a quantitative", " is allowed to have a dangling coding as seen in 'Age :'\n", sep = '\n' ) ), '-dataTable' = apl(n=c(1, 1000000), d=NA, h = paste( "-dataTable TABLE: List the data structure with a header as the first line.\n", " NOTE:\n", " 1) This option has to occur last; that is, no other options are", " allowed thereafter. Each line should end with a backslash except for", " the last line.\n", " 2) The first column is fixed and reserved with label 'Subj', and the", " last is reserved for 'InputFile'. Each row should contain only one", " effect estimate in the table of long format (cf. wide format) as", " defined in R. The level labels of a factor should contain at least", " one character. Input files can be in AFNI, NIfTI or surface format.", " AFNI files can be specified with sub-brick selector (square brackets", " [] within quotes) specified with a number or label. Unequal number of", " subjects across groups is allowed, but situations with missing data", " for a within-subject factor are better handled with 3dLME.\n", " 3) It is fine to have variables (or columns) in the table that are", " not modeled in the analysis.\n", " 4) The context of the table can be saved as a separate file, e.g.,", " called table.txt. Do not forget to include a backslash at the end of", " each row. In the script specify the data with '-dataTable @table.txt'.", " This option is useful: (a) when there are many input files so that", " the program complains with an 'Arg list too long' error; (b) when", " you want to try different models with the same dataset (see 3) above).\n", sep = '\n' ) ), '-help' = apl(n=0, h = '-help: this help message\n'), '-show_allowed_options' = apl(n=0, h= "-show_allowed_options: list of allowed options\n" ), '-cio' = apl(n=0, h = paste( "-cio: Use AFNI's C io functions, which is default. Alternatively -Rio", " can be used.\n", sep='\n')), '-Rio' = apl(n=0, h = "-Rio: Use R's io functions. The alternative is -cio.\n") ) ops <- parse.AFNI.args(args, params, other_ok=FALSE) if (verb) show.AFNI.args(ops, verb=0, hstr='') if (is.null(ops)) errex.AFNI('Error parsing arguments. See 3dMVM -help for details.') #Parse dems options #initialize with defaults com_history<-AFNI.command.history(ExecName, args,NULL) lop <- list (com_history = com_history) lop$nNodes <- 1 lop$model <- 1 lop$maskFN <- NA lop$wsVars <- NA lop$mVar <- NA lop$qVars <- NA lop$vVars <- NA lop$vQV <- NA lop$qVarCenters <- NA lop$vVarCenters <- NA lop$num_glt <- 0 lop$gltLabel <- NULL lop$gltCode <- NULL lop$dataTable <- NULL lop$SC <- FALSE lop$wsMVT <- FALSE lop$wsE2 <- FALSE # combining UVT and wsMVT, and then replacing UVT: only applicable for an # effect associated with a within-subject factor with more than 2 levels lop$mvE4 <- FALSE # combining 5 effects: 1st half UVT (AUC), 2nd half UVT (EXC), 1st half # wsMVT (AUC), 2nd half wsMVT (EXC), and MVT: assuming ONLY one within-subject # factor (e.g., effects from basis functions). Those 5 individual results # are not output (cf., lop$mvE4a). lop$mvE4a <- FALSE # combining 5 effects: 1st half UVT (AUC), 2nd half UVT (EXC), 1st half # wsMVT (AUC), 2nd half wsMVT (EXC), and MVT: assuming ONLY one within-subject # factor (e.g., effects from basis functions). Results added (or appended) to # others (cf., lop$mvE4a). lop$parSubset <- NA # parameter subset from option -matrPar for DTI data analysis lop$iometh <- 'clib' 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, prefix = lop$outFN <- pprefix.AFNI.name(ops[[i]]), mask = lop$maskFN <- ops[[i]], jobs = lop$nNodes <- ops[[i]], model = lop$model <- ops[[i]], bsVars = lop$model <- ops[[i]], wsVars = lop$wsVars <- ops[[i]], mVar = lop$mVar <- ops[[i]], qVars = lop$qVars <- ops[[i]], vVars = lop$vVars <- ops[[i]], qVarCenters = lop$qVarCenters <- ops[[i]], vVarCenters = lop$vVarCenters <- ops[[i]], num_glt = lop$num_glt <- ops[[i]], gltLabel = lop$gltLabel <- ops[[i]], gltCode = lop$gltCode <- ops[[i]], dataTable = lop$dataTable <- dataTable.AFNI.parse(ops[[i]]), parSubset = lop$parSubset <- ops[[i]], help = help.MVM.opts(params, adieu=TRUE), SC = lop$SC <- TRUE, wsMVT = lop$wsMVT <- TRUE, wsE2 = lop$wsE2 <- TRUE, mvE4 = lop$mvE4 <- TRUE, mvE4a = lop$mvE4a <- TRUE, cio = lop$iometh <- 'clib', Rio = lop$iometh <- 'Rlib' ) } return(lop) }# end of read.MVM.opts.batch # construct a glt list for testInteraction in phia # NEED to solve the problem when a quantitative variable is tested alone: # with pairwise = NULL!!! gltConstr <- function(cStr, dataStr) { pos <- which(cStr==":") vars <- cStr[pos-1] nvar <- length(vars) pos <- c(pos, length(cStr)+2) # add an artificial one for convenient usage below varsOK <- vars %in% colnames(dataStr) if(all(varsOK)) { gltList <- vector('list', nvar) for(ii in 1:nvar) { #browser() lvl <- levels(dataStr[,vars[ii]]) gltList[[ii]] <- rep(0, length(lvl)) sepTerms <- unlist(lapply(cStr[(pos[ii]+1):(pos[ii+1]-2)], strsplit, '\\*')) lvlInv <- sepTerms[seq(2,length(sepTerms),2)] # levels involved lvlOK <- lvlInv %in% lvl if(all(lvlOK)) { sq <- match(lvlInv, lvl) gltList[[ii]][sq] <- as.numeric(sepTerms[seq(1,length(sepTerms),2)]) } else errex.AFNI(paste("Incorrect level coding in variable", vars[ii], ": ", lvlInv[which(!lvlOK)], " \n ")) } names(gltList) <- vars return(gltList) } else errex.AFNI(paste("Incorrect variable name in GLT coding: ", vars[which(!varsOK)], " \n ")) } # test: pairwise = NULL # gltConstr(clist[[1]], lop$dataStr) # QVpos <- which(lop$gltCode[[n]] %in% lop$QV) # gltConstr(lop$gltCode[[n]][-c(QVpos, QVpos+1)], lop$dataStr) #Change options list to 3dMVM variable list process.MVM.opts <- function (lop, verb = 0) { if(is.null(lop$outFN)) errex.AFNI(c("Output filename not specified! Add filename with -prefix.\n")) an <- parse.AFNI.name(lop$outFN) if(an$type == "NIML") { if(file.exists(lop$outFN)) errex.AFNI(c("File ", lop$outFN, " exists! Try a different name.\n")) } else if(file.exists(paste(lop$outFN,"+tlrc.HEAD", sep="")) || file.exists(paste(lop$outFN,"+tlrc.BRIK", sep="")) || file.exists(paste(lop$outFN,"+orig.HEAD", sep="")) || file.exists(paste(lop$outFN,"+orig.BRIK", sep=""))) { errex.AFNI(c("File ", lop$outFN, " exists! Try a different name.\n")) return(NULL) } #Make sure new io must be used with anything but BRIK format an <- parse.AFNI.name(lop$outFN) if(an$type != 'BRIK' && lop$iometh != 'clib') errex.AFNI(c('Must of use -cio option with any input/output ', 'format other than BRIK')) if(!is.na(lop$qVars[1])) lop$QV <- strsplit(lop$qVars, '\\,')[[1]] if(!is.na(lop$vVars[1])) lop$vQV <- strsplit(lop$vVars, '\\,')[[1]] if(!is.na(lop$parSubset[1])) lop$parSubsetVector <- strsplit(lop$parSubset, '\\,')[[1]] if(!is.null(lop$gltLabel)) { sq <- as.numeric(unlist(lapply(lop$gltLabel, '[', 1))) if(identical(sort(sq), as.numeric(seq(1, lop$num_glt)))) { lop$gltLabel <- unlist(lapply(lop$gltLabel, '[', 2)) lop$gltLabel[sq] <- lop$gltLabel } else errex.AFNI(c('The number of -gltLabel is not consistent with that \n', 'specified by -num_glt ', lop$num_glt)) } if(!is.null(lop$gltCode)) { sq <- as.numeric(unlist(lapply(lop$gltCode, '[', 1))) if(identical(sort(sq), as.numeric(seq(1, lop$num_glt)))) { lop$gltCode <- lapply(lop$gltCode, '[',-1) lop$gltCode[sq] <- lop$gltCode } else errex.AFNI(c('The number of -gltLabel is not consistent with that \n', 'specified by -num_glt ', lop$num_glt)) } len <- length(lop$dataTable) # in case the user misspells #wd <- which(lop$dataTable == "InputFile" | lop$dataTable == "Inputfile") wd <- which(lop$dataTable %in% respVar) hi <- len / wd - 1 #browser() if(len %% wd != 0) errex.AFNI(paste('The content under -dataTable is not rectangular !', len, wd)) else { lop$dataStr <- NULL for(ii in 1:wd) lop$dataStr <- data.frame(cbind(lop$dataStr, lop$dataTable[seq(wd+ii, len, wd)])) names(lop$dataStr) <- lop$dataTable[1:wd] # wow, terrible mistake here with as.numeric(lop$dataStr[,jj]) #if(!is.na(lop$qVars)) for(jj in lop$QV) lop$dataStr[,jj] <- as.numeric(lop$dataStr[,jj]) if(!is.na(lop$qVars[1])) for(jj in lop$QV) lop$dataStr[,jj] <- as.numeric(as.character(lop$dataStr[,jj])) # or if(!is.na(lop$qVars)) for(jj in lop$QV) lop$dataStr[,jj] <- as.numeric(levels(lop$dataStr[,jj]))[as.integer(lop$dataStr[,jj])] if(!is.na(lop$vVars[1])) for(jj in lop$vQV) lop$dataStr[,jj] <- as.character(lop$dataStr[,jj]) } # set the covariate default values at their centers lop$covVal <- NULL if(length(lop$QV)>0) { lop$covVal <- rep(0, length(lop$QV)) names(lop$covVal) <- lop$QV } if (lop$num_glt > 0) { lop$gltList <- vector('list', lop$num_glt) lop$slpList <- vector('list', lop$num_glt) for (n in 1:lop$num_glt) { # assuming each GLT has one slope involved and placed last #if(!is.na(lop$qVars)) { if(any(lop$QV %in% lop$gltCode[[n]])) { if(!is.na(lop$qVars) & any(lop$QV %in% lop$gltCode[[n]])) { QVpos <- which(lop$gltCode[[n]] %in% lop$QV) if(QVpos > 1) lop$gltList[[n]] <- gltConstr(lop$gltCode[[n]][-c(QVpos, QVpos+1)], lop$dataStr) if(QVpos == 1) lop$gltList[[n]] <- NA lop$slpList[[n]] <- lop$gltCode[[n]][QVpos] } else if(!is.na(lop$vVars) & any(lop$vQV %in% lop$gltCode[[n]])) { vQVpos <- which(lop$gltCode[[n]] %in% lop$vQV) if(vQVpos > 1) lop$gltList[[n]] <- gltConstr(lop$gltCode[[n]][-c(vQVpos, vQVpos+1)], lop$dataStr) if(vQVpos == 1) lop$gltList[[n]] <- NULL lop$slpList[[n]] <- lop$gltCode[[n]][vQVpos] } else lop$gltList[[n]] <- gltConstr(lop$gltCode[[n]], lop$dataStr) } } if(lop$iometh == 'Rlib') lop$outFN <- paste(lop$outFN, "+tlrc", sep="") else { an <- parse.AFNI.name(lop$outFN) if(an$type == "BRIK" && an$ext == "" && is.na(an$view)) lop$outFN <- paste(lop$outFN, "+tlrc", sep="") if (exists.AFNI.name(lop$outFN) || exists.AFNI.name(modify.AFNI.name(lop$outFN,"view","+tlrc"))) errex.AFNI(c("File ", lop$outFN, " exists! Try a different name.\n")) } if(lop$nNodes < 1) lop$nNodes <- 1 if(!is.na(lop$maskFN)) { if(verb) cat("Will read ", lop$maskFN,'\n') if(is.null(mm <- read.AFNI(lop$maskFN, verb=lop$verb, meth=lop$iometh))) { warning("Failed to read mask", immediate.=TRUE) return(NULL) } lop$maskData <- mm$brk if(verb) cat("Done read ", lop$maskFN,'\n') } if(!is.na(lop$maskFN)) if(!all(dim(lop$maskData[,,,1])==lop$myDim[1:3])) stop("Mask dimensions don't match the input files!") return(lop) } ################################################################################# ################# MVM Computation functions ############################## ################################################################################# # used only in read.MVM.opts.from.file scanLine <- function(file, lnNo=1, marker="\\:") unlist(strsplit(unlist(scan(file, what= list(""), skip=lnNo-1, strip.white=TRUE, nline=1)), marker))[2] # function to obtain multivariate ANOVA # may totally abandon this and switch to linearHypothesis # Pillai test only now # does this work for multiple between-subjects variables (factors and quantitative variables)????????? #maov <- function(mvfm) # Pillai test with type = 3, need to remove the intercept: -1 below # return(stats:::Pillai(Re(eigen(qr.coef(qr(mvfm$SSPE), mvfm$SSP[[-1]]), symmetric = FALSE)$values), mvfm$df[-1], mvfm$error.df)) # general linear test for between-subjects variables only # takes components from Anova(fm$lm, type=3, test='Pillai') as input maov <- function(SSPE, SSP, DF, error.DF) # Pillai test with type = 3 return(stats:::Pillai(Re(eigen(qr.coef(qr(SSPE), SSP), symmetric = FALSE)$values), DF, error.DF)) # Combines 3 tests, and outputs the minimum p-values among 4 tests: # 2 from UVT, 1 from within-subject MVT, and 1 from MVT. # takes model object from aov.car from afex as input: assuming that # model fitting with afex: ONE within-subject factor ONLY! mvCom4 <- function(fm, nF_mvE4) { uvfm <- tryCatch(univ(fm$Anova), error=function(e) NULL) # univariate model uvP <- rep(1, 2*nF_mvE4) p_wsmvt <- rep(1, nF_mvE4) p_mvt <- rep(1, nF_mvE4) if(!is.null(uvfm)) { #nTerms <- nrow(uvfm$anova) # totaly number of effect estimates #outTerms <- nTerms/2 # half of them # UVT p-values uvP <- uvfm$anova[,'Pr(>F)'] # p-values for UVT # within-subject MVT: one set #p_wsmvt <- rep(1, nTerms) # initiation for within-subject MVT for(ii in 1:nF_mvE4) { jj <- nF_mvE4 + ii wsmvt <- maov(fm$Anova$SSPE[[jj]], fm$Anova$SSP[[jj]], fm$Anova$df[jj], fm$Anova$error.df) #p-value for upper F p_wsmvt[ii] <- pf(wsmvt[2], wsmvt[3], wsmvt[4], lower.tail = FALSE) } } # true MVT #mvfm <- Anova(fm$lm, type=3, test='Pillai') mvfm <- tryCatch(Anova(fm$lm, type = 3, test = "Pillai"), error=function(e) NULL) if(is.null(mvfm)) out_p <- apply(cbind(uvP[1:nF_mvE4], uvP[(nF_mvE4+1):(2*nF_mvE4)], p_wsmvt), 1, min) else { for(kk in 1:nF_mvE4) { mvt <- stats:::Pillai(Re(eigen(qr.coef(qr(mvfm$SSPE), mvfm$SSP[[kk]]), symmetric = FALSE)$values), mvfm$df[[kk]], mvfm$error.df) #p-value for upper F p_mvt[kk] <- pf(mvt[2], mvt[3], mvt[4], lower.tail = FALSE) } out_p <- apply(cbind(uvP[1:nF_mvE4], uvP[(nF_mvE4+1):(2*nF_mvE4)], p_wsmvt, p_mvt), 1, min) } return(out_p) } runAOV <- function(inData, dataframe, ModelForm, pars) { # assign voxel-wise covariate values assVV <- function(DF, vQV, value, c) { # centering if(is.na(c)) cvalue <- scale(value, center=TRUE, scale=F) else cvalue <- scale(value, center=c, scale=F) for(ii in 1:length(unique(DF[,vQV]))) { # ii-th subject's rows fn <- paste(vQV, '_fn', sep='') sr <- which(unique(DF[,fn])[ii] == DF[,fn]) #browser() DF[sr, vQV] <- rep(cvalue[ii], length(sr)) } DF[, vQV] <- as.numeric(DF[, vQV]) return(DF) } #pars[[10]] <- list(lop$vQV, ?) #lop$dataStr <- assVV(lop$dataStr, lop$vQV[1], inData[ii,jj,kk,(NoFile+1):(NoFile+nSubj)]) out <- pars[[1]] options(warn = -1) #browser() if (!all(abs(inData) < 10e-8)) { dataframe$Beta<-inData[1:pars[[10]][[3]]] if(any(!is.na(pars[[10]][[1]]))) { dataframe <- assVV(dataframe, pars[[10]][[1]], inData[(pars[[10]][[3]]+1):(pars[[10]][[3]]+pars[[10]][[4]])], pars[[10]][[2]]) #if(pars[[10]][[3]]) dataframe[,pars[[10]][[1]]] <- scale(dataframe[,pars[[10]][[1]]], center=TRUE, scale=F) else # dataframe[,pars[[10]][[1]]] <- scale(dataframe[,pars[[10]][[1]]], center=pars[[10]][[2]], scale=F) } fm <- NULL try(fm <- aov.car(ModelForm, data=dataframe, factorize=FALSE, return='full'), silent=TRUE) if(!is.null(fm)) { uvfm <- tryCatch(univ(fm$Anova), error=function(e) NULL) # univariate model if(!is.null(uvfm)) { if(pars[[6]][1] & pars[[7]]) { # between-subjects factors/variables only tryCatch(Fvalues <- uvfm[1:pars[[2]][1],3], error=function(e) NULL) if(!is.null(Fvalues)) if(!any(is.nan(Fvalues))) out[1:pars[[2]][2]] <- Fvalues } else if(pars[[6]][6]) { # combine 4 tests: assuming only one within-subject factor tryCatch(out[(1:pars[[2]][1])] <- qchisq(mvCom4(fm, pars[[2]][5]), 1, lower.tail = FALSE), error=function(e) NULL) } else {# contain within-subject variable(s) #tryCatch(Fvalues <- unname(uvfm$anova[-1,5]), error=function(e) NULL) tryCatch(Fvalues <- unname(uvfm$anova[,5]), error=function(e) NULL) if(!is.null(Fvalues)) if(!any(is.nan(Fvalues))) { out[1:pars[[2]][2]] <- Fvalues # univariate Fs: no spherecity correction if(pars[[6]][2]) { # sphericity correction getGG <- uvfm$sphericity.correction[,'HF eps'] < pars[[8]][1] GG <- uvfm$sphericity.correction[,'Pr(>F[GG])'] HF <- uvfm$sphericity.correction[,'Pr(>F[HF])'] #Fsc <- ifelse(getGG, GG, HF) Fsc <- ifelse(uvfm$mauchly[, 'p-value'] < 0.05, ifelse(getGG, GG, HF), uvfm$anova[, 'Pr(>F)']) tryCatch(out[(pars[[2]][2]+1):(pars[[2]][2]+pars[[2]][3])] <- qf(Fsc, pars[[8]][[2]], pars[[8]][[3]], lower.tail = FALSE), error=function(e) NULL) #qf(Fsc, uvfm$anova[dimnames(uvfm$sphericity.correction)[[1]], 'num Df'], # uvfm$anova[dimnames(uvfm$sphericity.correction)[[1]], 'den Df'], # lower.tail = FALSE), error=function(e) NULL) } #if(pars[[6]][2]) if(pars[[6]][3]) { # within-subject MVT is requested #for(ii in 1:length(fm$Anova$SSPE)) { for(ii in 1:length(pars[[9]])) tryCatch(out[pars[[2]][2]+pars[[2]][3]+ii] <- maov(fm$Anova$SSPE[[pars[[9]][ii]]], fm$Anova$SSP[[pars[[9]][ii]]], fm$Anova$df[pars[[9]][ii]], fm$Anova$error.df)[2], error=function(e) NULL) } if(pars[[6]][4]) { # option -wsE2: UVT and wsMVT combined #for(ii in 1:length(fm$Anova$SSPE)) { for(ii in pars[[9]]) { wsmvt <- NULL tryCatch(wsmvt <- maov(fm$Anova$SSPE[[ii]], fm$Anova$SSP[[ii]], fm$Anova$df[ii], fm$Anova$error.df), error=function(e) NULL) if(!is.null(wsmvt)) { p_wsmvt <- pf(wsmvt[2], wsmvt[3], wsmvt[4], lower.tail = FALSE) # replace value at index ii if(p_wsmvt < uvfm$anova[ii,'Pr(>F)']) out[ii] <- qf(p_wsmvt, uvfm$anova[ii,'num Df'], uvfm$anova[ii,'den Df'], lower.tail = FALSE) } } } # if(pars[[6]][4]) } #if(!any(is.nan(Fvalues))) } #if(pars[[6]][1] & pars[[7]]) } #if(!is.null(uvfm)) if(!pars[[7]]) { # reall MVM: currently not allowed to mix within-subject variables tryCatch(mvfm <- Anova(fm$lm, type=3, test='Pillai'), error=function(e) NULL) # need to add options for type and test! #if(pars[[6]][1]) tryCatch(out[(pars[[2]][2]+pars[[2]][3]+1):pars[[2]][1]] <- maov(mvfm)[2], error=function(e) NULL) for(ii in 1:pars[[2]][4]) # pars[[2]][4] equals length(mvfm$terms) #tryCatch(out[pars[[2]][2]+pars[[2]][3]+ii] <- tryCatch(out[pars[[2]][2]+pars[[2]][3]+length(pars[[9]])+ii] <- maov(mvfm$SSPE, mvfm$SSP[[ii]], mvfm$df[ii], mvfm$error.df)[2], error=function(e) NULL) } #if(!pars[[7]]) if(pars[[6]][5]) { # combine 4 tests: assuming only one within-subject factor tryCatch(out[(pars[[2]][2]+pars[[2]][3]+length(pars[[9]])+pars[[2]][4]+1): (pars[[2]][2]+pars[[2]][3]+length(pars[[9]])+pars[[2]][4]+pars[[2]][5])] <- qchisq(mvCom4(fm, pars[[2]][5]), 1, lower.tail = FALSE), error=function(e) NULL) } # redundant computations in mvCom4 for option mvE4a # GLT part below if(pars[[3]]>=1) for(ii in 1:pars[[3]]) { if(all(is.na(pars[[4]][[ii]]))) glt <- tryCatch(testInteractions(fm$lm, pair=NULL, slope=pars[[5]][[ii]], covariates=lop$covVal, adjustment="none", idata = fm[["idata"]]), error=function(e) NULL) else glt <- tryCatch(testInteractions(fm$lm, custom=pars[[4]][[ii]], slope=pars[[5]][[ii]], covariates=lop$covVal, adjustment="none", idata = fm[["idata"]]), error=function(e) NULL) if(!is.null(glt)) { out[pars[[2]][1]+2*ii-1] <- glt[1,1] out[pars[[2]][1]+2*ii] <- sign(glt[1,1]) * sqrt(glt[1,4]) # convert F to t } #if(!is.null(glt)) } #if(pars[[3]]>=1) for(ii in 1:pars[[3]]) } } return(out) } # covariates=pars[[6]][7], adjustment="none", idata = fm[["idata"]]), error=function(e) NULL) else # covariates=pars[[6]][7], adjustment="none", idata = fm[["idata"]]), error=function(e) NULL) ################################################################################# ########################## Read information from a file ######################### ################################################################################# #A function to parse all user input from a file read.MVM.opts.from.file <- function (modFile='model.txt', verb = 0) { lop <- list() #List to hold all options from user input lop$verb <- verb lop$iometh <- 'clib' # Line 1: data type - volume or surface #datatype <- scanLine(modFile, lnNo=1) # Line 2: Output filename lop$outFN <-scanLine(modFile, lnNo=2) lop$outFN <- paste(lop$outFN, "+tlrc", sep="") # check if the filename exists already if(any(file.exists(paste(lop$outFN,'.BRIK', sep=''), paste(lop$outFN,'.HEAD', sep='')))) errex.AFNI("Output filename already exists! Please rename it...") # Line 3: MASK lop$maskFN <- scanLine(modFile, lnNo=3) # Line 4: #jobs lop$nNodes<- as.integer(scanLine(modFile, lnNo=4)) # Line 5: model formula lop$model <- scanLine(modFile, lnNo=5) # Line 6: within-subject variables lop$wsVars <- scanLine(modFile, lnNo=6) # Line 7: covariates lop$QV is NA if no covariates lop$qVars <- unlist(strsplit(scanLine(modFile, lnNo=7), "\\*")) lop$QV <- strsplit(lop$qVars, '\\+')[[1]] # Line 8: covariates center value: NA means mean lop$qVarCenters <- unlist(strsplit(scanLine(modFile, lnNo=8), "\\*")) #lop$qVarCenters <- as.numeric(strsplit(lop$qVarCenters, '\\,')[[1]]) # header position (hp) defined by column name InputFile hp <- grep("InputFile", readLines(modFile)) lop$dataStr <- read.table(modFile, skip=hp[1]-1, header=TRUE) # number of contrasts (pair-wise only right now) lop$num_glt <- (hp-12)%/%2 if (lop$num_glt > 0) { lop$gltLabel <- array(data=NA, dim=lop$num_glt) clist <- vector('list', lop$num_glt) lop$gltList <- vector('list', lop$num_glt) lop$slpList <- vector('list', lop$num_glt) for (n in 1:lop$num_glt) { lop$gltLabel[n] <- paste(unlist(scan(file=modFile, what= list(""), skip=10+2*n-1, strip.white=TRUE, nline=1)), collapse="") clist[[n]] <- scan(file=modFile, what= list(""), skip=10+2*n, nline=1)[[1]] if(any(lop$QV %in% clist[[n]])) { QVpos <- which(clist[[n]] %in% lop$QV) lop$gltList[[n]] <- gltConstr(clist[[n]][-c(QVpos, QVpos+1)], lop$dataStr) lop$slpList[[n]] <- clist[[n]][QVpos] } else lop$gltList[[n]] <- gltConstr(clist[[n]], lop$dataStr) } } #browser() return(lop) } # end of read.MVM.from.file ################################################################################# ########################## Begin MVM main ###################################### ################################################################################# if(!exists('.DBG_args')) { args = (commandArgs(TRUE)) rfile <- first.in.path(sprintf('%s.R',ExecName)) try(save(args, rfile, file=".3dMVM.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 cat(greeting.MVM(), "Use CNTL-C on Unix or ESC on GUI version of R to stop at any moment.\n", sep='\n') #browser() if(length(args)<6) modFile <- "model.txt" else modFile <- args[6] if (is.null(lop <- read.MVM.opts.from.file(modFile, verb=0))) { stop('Error parsing input from file!'); } if(0) str(lop) } else { if(!exists('.DBG_args')) { BATCH_MODE <<- 1 } else { BATCH_MODE <<- 0 } if(is.null(lop <- read.MVM.opts.batch(args, verb = 0))) stop('Error parsing input') #str(lop); if(is.null(lop <- process.MVM.opts(lop, verb = lop$verb))) stop('Error processing input') } #if(lop$verb > 1) { #Too much output, big dump of header structs of input dsets.. # str(lop) #} ######################################################################## if(!is.na(lop$qVarCenters)) lop$qVarCenters <- as.numeric(strsplit(as.character(lop$qVarCenters), '\\,')[[1]]) if(!is.na(lop$vVarCenters)) lop$vVarCenters <- as.numeric(strsplit(as.character(lop$vVarCenters), '\\,')[[1]]) library("afex") library("phia") #comArgs <- commandArgs() #if(length(comArgs)<6) modFile <- "model.txt" else #modFile <- comArgs[6] #paste(commandArgs()) wd <- which(lop$dataTable %in% respVar) # update respVar respVar <- lop$dataTable[wd] # even if lop$wsVars is NA (no within-subject factors), it would be still OK for Error(Subj/NA) if(is.na(lop$mVar)) { if(is.na(lop$wsVars)) ModelForm <- as.formula(paste("Beta ~", lop$model, '+Error(Subj)')) else { ModelForm <- as.formula(paste("Beta ~", lop$model, '+Error(Subj/(', lop$wsVars, '))')) ModelForm2 <- as.formula(paste("Beta ~ (", lop$model, ')*', lop$wsVars, '+Error(Subj/(', lop$wsVars, '))')) } } else if(is.na(lop$wsVars)) ModelForm <- as.formula(paste("Beta ~", lop$model, '+Error(Subj/(', lop$mVar, '))')) else ModelForm <- as.formula(paste("Beta ~", lop$model, '+Error(Subj/(', lop$wsVars, '*', lop$mVar, '))')) # Maybe not list for these two, or yes? lop$dataStr$Subj <- as.factor(lop$dataStr$Subj) lop$dataStr[[respVar]] <- as.character(lop$dataStr[[respVar]]) # center on user-speficied value or mean if(any(!is.na(lop$qVars))) if(all(is.na(lop$qVarCenters))) lop$dataStr[,lop$QV] <- scale(lop$dataStr[,lop$QV], center=TRUE, scale=F) else lop$dataStr[,lop$QV] <- scale(lop$dataStr[,lop$QV], center=lop$qVarCenters, scale=F) cat('\n++++++++++++++++++++++++++++++++++++++++++++++++++++\n') cat('***** Summary information of data structure *****\n') #print(sprintf('%i subjects: ', nlevels(lop$dataStr$Subj))) #for(ii in 1:nlevels(lop$dataStr$Subj)) print(sprintf('%s ', levels(lop$dataStr$Subj)[ii])) nSubj <- nlevels(lop$dataStr$Subj) cat(nSubj, 'subjects : ', levels(lop$dataStr$Subj), '\n') cat(length(lop$dataStr[[respVar]]), 'response values\n') for(ii in 2:(dim(lop$dataStr)[2]-1)) if(class(lop$dataStr[,ii]) == 'factor') cat(nlevels(lop$dataStr[,ii]), 'levels for factor', names(lop$dataStr)[ii], ':', levels(lop$dataStr[,ii]), '\n') else if(class(lop$dataStr[,ii]) == 'matrix' | class(lop$dataStr[,ii]) == 'numeric') # numeric may not work cat(length(lop$dataStr[,ii]), 'centered values for numeric variable', names(lop$dataStr)[ii], ':', lop$dataStr[,ii], '\n') cat(lop$num_glt, 'post hoc tests\n') cat('\nContingency tables of subject distributions among the categorical variables:\n\n') if(is.na(lop$mVar)) if(is.na(lop$wsVars)) showTab <- paste('~', lop$model) else showTab <- paste('~', gsub("\\*", "+", lop$model), '+', gsub("\\*", "+", lop$wsVars)) else if(is.na(lop$wsVars)) showTab <- as.formula(paste('~', gsub("\\*", "+", lop$model), "+", gsub("\\*", "+", lop$mVar))) else showTab <- paste('~', lop$model, "+", gsub("\\*", "+", lop$wsVars), "+", gsub("\\*", "+", lop$mVar)) if(!is.na(lop$qVars)) for(ii in 1:length(lop$QV)) showTab <- gsub(paste('\\*',lop$QV[ii], sep=''), '', gsub(paste('\\+',lop$QV[ii], sep=''), '', showTab)) showTab <- as.formula(gsub("\\*", "+", showTab)) # in case there are still some *'s like between-subjects factors if(!(showTab == '~1')) print(xtabs(showTab, data=lop$dataStr)) cat('\nTabulation of subjects against each of the categorical variables:') all_vars <- names(lop$dataStr) for(var in all_vars[-c(1, length(all_vars))]) if(!(var %in% lop$QV)) { cat('\n~~~~~~~~~~~~~~') cat('\nSubj vs ', var, ':\n', sep='') print(table(lop$dataStr$Subj, lop$dataStr[,var])) } cat('\n***** End of data structure information *****\n') cat('++++++++++++++++++++++++++++++++++++++++++++++++++++\n\n') # Assume the last column is input files #FileCol <- length(colnames(lop$dataStr)) FileCol <- dim(lop$dataStr)[2] #Number of input files NoFile <- dim(lop$dataStr[1])[1] # Repeated-measures (use lme) or not (use lm) #if (length(unique(lop$dataStr$Subj)) != length(lop$dataStr$Subj)) RM <- TRUE else RM <- FALSE cat('Reading input files now...\n\n') if(any(is.na(suppressWarnings(as.numeric(lop$dataStr[, FileCol]))))) { # not elegant because "NAs introduced by coercion" # Read in the 1st input file so that we have the dimension information inData <- read.AFNI(lop$dataStr[1, FileCol], verb=lop$verb, meth=lop$iometh) dimx <- inData$dim[1] dimy <- inData$dim[2] dimz <- inData$dim[3] # for writing output purpose head <- inData # ww <- inData$NI_head #myHist <- inData$header$HISTORY_NOTE; myOrig <- inData$origin; myDelta <- inData$delta # Read in all input files inData <- unlist(lapply(lapply(lop$dataStr[,FileCol], read.AFNI, verb=lop$verb, meth=lop$iometh), '[[', 1)) dim(inData) <- c(dimx, dimy, dimz, NoFile) # voxel-wise covariate files if(any(!is.na(lop$vVars))) { for(ii in lop$vQV) if(length(unique(lop$dataStr[,ii])) != nlevels(lop$dataStr$Subj)) errex.AFNI(c("Error with voxel-wise covariate ", ii, ": Each subject is only\n", "allowed to have one volume; that is, the covariate has to be at the\n", "subject level.")) else { # currently consider one voxel-wise covariate only: may generalize later? vQV <- unlist(lapply(lapply(unique(lop$dataStr[,lop$vQV[1]]), read.AFNI, verb=lop$verb, meth=lop$iometh), '[[', 1)) dim(vQV) <- c(dimx, dimy, dimz, length(unique(lop$dataStr[,lop$vQV[1]]))) inData <- c(inData, vQV) dim(inData) <- c(dimx, dimy, dimz, NoFile+nSubj) } } else vQV <- NULL if (!is.na(lop$maskFN)) { Mask <- read.AFNI(lop$maskFN, verb=lop$verb, meth=lop$iometh)$brk[,,,1] inData <- array(apply(inData, 4, function(x) x*Mask), dim=c(dimx,dimy,dimz,NoFile+(!is.na(lop$vQV[1]))*nSubj)) } # assign voxel-wise covariate values assVV <- function(DF, vQV, value, c) { if(is.na(c)) cvalue <- scale(value, center=TRUE, scale=F) else cvalue <- scale(value, center=c, scale=F) for(ii in 1:length(unique(DF[,vQV]))) { # ii-th subject's rows fn <- paste(lop$vQV[1], '_fn', sep='') sr <- which(unique(DF[,fn])[ii] == DF[,fn]) #browser() DF[sr, vQV] <- rep(cvalue[ii], length(sr)) } DF[, vQV] <- as.numeric(DF[, vQV]) return(DF) } # if(all(is.na(lop$vVarCenters))) # lop$dataStr[,lop$QV] <- scale(lop$dataStr[,lop$QV], center=TRUE, scale=F) else # lop$dataStr[,lop$QV] <- scale(lop$dataStr[,lop$QV], center=lop$vVarCenters, scale=F) # add a new column to store the voxel-wise filenames if(any(!is.na(lop$vVars))) { lop$dataStr <- cbind(lop$dataStr, lop$dataStr[, lop$vQV[1]]) names(lop$dataStr)[length(names(lop$dataStr))] <- paste(lop$vQV[1], '_fn', sep='') } #lop$dataStr <- assVV(lop$dataStr, paste(lop$vQV[1], '_fn', sep=''), vQV[20,20,20,]) # DF: lop$dataStr[,lop$vQV[1]] # try out a few voxels and see if the model is OK, and find out the number of F tests and DF's # for t tests (and catch potential problems as well) #ii<-dimx%/%3; jj<-dimy%/%3; kk<-dimz%/%3 ############################### xinit <- dimx%/%3 if(dimy==1) yinit <- 1 else yinit <- dimy%/%2 if(dimz==1) zinit <- 1 else zinit <- dimz%/%2 ii <- xinit; jj <- yinit; kk <- zinit fm<-NULL gltRes <- vector('list', lop$num_glt) cat('If the program hangs here for more than, for example, half an hour,\n') cat('kill the process because the model specification or the GLT coding\n') cat('is likely inappropriate.\n\n') while(is.null(fm)) { fm<-NULL if (all(abs(inData[ii, jj, kk,]) < 10e-8)) fm<-NULL else { lop$dataStr$Beta<-inData[ii, jj, kk,1:NoFile] if(any(!is.na(lop$vVars))) { #lop$dataStr <- assVV(lop$dataStr, lop$vQV[1], vQV[ii,jj,kk,]) lop$dataStr <- assVV(lop$dataStr, lop$vQV[1], inData[ii,jj,kk,(NoFile+1):(NoFile+nSubj)], lop$vVarCenters[1]) #if(all(is.na(lop$vVarCenters))) # lop$dataStr[,lop$QV] <- scale(lop$dataStr[,lop$QV], center=TRUE, scale=F) else # lop$dataStr[,lop$QV] <- scale(lop$dataStr[,lop$QV], center=lop$vVarCenters, scale=F) } #options(warn=-1) try(fm <- aov.car(ModelForm, data=lop$dataStr, factorize=FALSE, return='full'), silent=TRUE) # if(!is.null(fm)) if(!is.na(lop$mVar)) { # if(is.na(lop$wsVars)) mvfm <- Anova(fm$lm, type=2, test='Pillai') # } else { # uvfm <- univ(fm[[1]]) # univariate modeling # } # need to add options for type and test! if(!is.null(fm)) { uvfm <- univ(fm$Anova) # univariate modeling if(!is.na(lop$mVar)) if(is.na(lop$wsVars)) mvfm <- Anova(fm$lm, type=3, test='Pillai') } if(!is.null(fm)) if (lop$num_glt > 0) { n <- 1 while(!is.null(fm) & (n <= lop$num_glt)) { if(all(is.na(lop$gltList[[n]]))) gltRes[[n]] <- tryCatch(testInteractions(fm$lm, pair=NULL, covariates=lop$covVal, slope=lop$slpList[[n]], adjustment="none", idata = fm[["idata"]]), error=function(e) NA) else gltRes[[n]] <- tryCatch(testInteractions(fm$lm, custom=lop$gltList[[n]], slope=lop$slpList[[n]], covariates=lop$covVal, adjustment="none", idata = fm[["idata"]]), error=function(e) NULL) if(any(is.null(gltRes[[n]]))) fm <- NULL n <- n+1 } } } if(!is.null(fm)) { print(sprintf("Great, test run passed at voxel (%i, %i, %i)!", ii, jj, kk)) } else if(ii 0)) brickNames <- c(brickNames, paste(corTerms, '-SC-', 'F')) # if(lop$wsMVT & (nF_MVT > 0)) brickNames <- c(brickNames, paste(names(fm$Anova$SSPE), '-wsMVT-', 'F')) if(lop$wsMVT & (nF_MVT > 0)) brickNames <- c(brickNames, paste(rownames(uvfm$sphericity.correction), '-wsMVT-', 'F')) } #brickNames <- ifelse(is.na(lop$wsVars) & is.na(lop$mVar), # paste(dimnames(uvfm)[[1]][2:(length(dimnames(uvfm)[[1]])-1)], 'F'), # paste(dimnames(uvfm$anova)[[1]][-1], 'F')) if(!is.na(lop$mVar)) brickNames <- c(brickNames, paste(dimnames(uvfm$anova)[[1]][1:nF_mvE4], '-MV0-', 'F')) if(lop$mvE4a) brickNames <- c(brickNames, paste(dimnames(uvfm$anova)[[1]][1:nF_mvE4], '-mvE4', 'Chisq')) if(lop$mvE4) brickNames <- paste(dimnames(uvfm$anova)[[1]][1:nF_mvE4], '-mvE4', 'Chisq') # no appending for(ii in 1:lop$num_glt) { brickNames <- c(brickNames, lop$gltLabel[ii]) brickNames <- c(brickNames, paste(lop$gltLabel[ii], 't')) } if(lop$SC & (nFsc > 0)) { scTerms <- dimnames(uvfm$anova)[[1]] %in% corTerms numDF <- uvfm$anova[,'num Df'][scTerms] denDF <- uvfm$anova[,'den Df'][scTerms] } else { numDF <- NULL denDF <- NULL } # DFs for t-stat t_DF <- NULL if(lop$num_glt>0) for(ii in 1:lop$num_glt) t_DF <- c(t_DF, ifelse(is.na(lop$wsVars) & is.na(lop$mVar), gltRes[[ii]][2,2], gltRes[[ii]][,6])) # DFs for F-stat F_DF <- vector('list', nF) for(ii in 1:nFu) if(is.na(lop$mVar) & is.na(lop$wsVars)) # between-subjects variables only F_DF[[ii]] <- c(uvfm[ii, 'Df'], uvfm[nF+1, 'Df']) else # having within-subject factor F_DF[[ii]] <- c(unname(uvfm$anova[ii,'num Df']), unname(uvfm$anova[ii,'den Df'])) # skip the intercept: ii+1 if(nFsc > 0) for(ii in 1:nFsc) F_DF[[nFu+ii]] <- c(numDF[ii], denDF[ii]) if(nF_MVT > 0) for(ii in 1:nF_MVT) F_DF[[nFu+nFsc+ii]] <- c(unname(maov(fm$Anova$SSPE[[mvtInd[ii]]], fm$Anova$SSP[[mvtInd[ii]]], fm$Anova$df[mvtInd[ii]], fm$Anova$error.df)[3]), unname(maov(fm$Anova$SSPE[[mvtInd[ii]]], fm$Anova$SSP[[mvtInd[ii]]], fm$Anova$df[mvtInd[ii]], fm$Anova$error.df)[4])) if(nFm > 0) for(ii in 1:nFm) { mvtest <- maov(mvfm$SSPE, mvfm$SSP[[ii]], mvfm$df[ii], mvfm$error.df) F_DF[[nFu+nFsc+nF_MVT+ii]] <- c(mvtest[3], mvtest[4]) } ############################### pars <- vector("list", 10) #pars[[1]] <- NoBrick pars[[1]] <- outInit #pars[[2]] <- nF pars[[2]] <- c(nF, nFu, nFsc, nFm, nF_mvE4, numDF, denDF) pars[[3]] <- lop$num_glt pars[[4]] <- lop$gltList pars[[5]] <- lop$slpList pars[[6]] <- c(is.na(lop$wsVars), lop$SC, lop$wsMVT, lop$wsE2, lop$mvE4a, lop$mvE4) # any within-subject factors? pars[[7]] <- is.na(lop$mVar) # any real multivariate modeling: currently for basis functions pars[[8]] <- list(0.75, numDF, denDF) # switching threshold between GG and HF: 0.6 pars[[9]] <- mvtInd # which indices for wsMVT #pars[[9]] <- !is.null(mvtInd) # which indices for wsMVT pars[[10]] <- list(lop$vQV, all(is.na(lop$vVarCenters)), NoFile, nSubj) # only run wsMVT for those terms associated with a within-subject factor: # which(names(fm$Anova$SSPE) %in% dimnames(uvfm$sphericity.correction)[[1]]) #if(any(!is.na(lop$vVars))) lop$dataStr <- assVV(lop$dataStr, lop$vQV[1], vQV[ii,jj,kk,]) #runAOV(inData[ii, jj, kk,], dataframe=lop$dataStr, ModelForm=ModelForm, pars=pars) print(sprintf("Start to compute %s slices along Z axis. You can monitor the progress", dimz)) print("and estimate the total run time as shown below.") print(format(Sys.time(), "%D %H:%M:%OS3")) ############################### #options(warn = -1) # suppress warnings! #getOption('warn') #for(ii in 1:dimx) for(jj in 1:dimy) # aa<-runAOV(inData[ii,jj,kk,], dataframe=lop$dataStr, ModelForm=ModelForm, pars=pars) if(dimy == 1 & dimz == 1) { nSeg <- 20 # drop the dimensions with a length of 1 inData <- inData[, , ,] # break into 20 segments, leading to 5% increamental in parallel computing dimx_n <- dimx%/%nSeg + 1 # number of datasets need to be filled fill <- nSeg-dimx%%nSeg # pad with extra 0s inData <- rbind(inData, array(0, dim=c(fill, NoFile))) # declare output receiver out <- array(0, dim=c(dimx_n, nSeg, NoBrick)) # break input multiple segments for parrel computation dim(inData) <- c(dimx_n, nSeg, NoFile) if (lop$nNodes==1) for(kk in 1:nSeg) { if(NoBrick > 1) out[,kk,] <- aperm(apply(inData[,kk,], 1, runAOV, dataframe=lop$dataStr, ModelForm=ModelForm, pars=pars), c(2,1)) else out[,kk,] <- aperm(apply(inData[,kk,], 1, runAOV, dataframe=lop$dataStr, ModelForm=ModelForm, pars=pars), dim=c(dimx_n, 1)) cat("Computation done: ", 100*kk/nSeg, "%", format(Sys.time(), "%D %H:%M:%OS3"), "\n", sep='') } if (lop$nNodes>1) { library(snow) cl <- makeCluster(lop$nNodes, type = "SOCK") clusterEvalQ(cl, library(afex)); clusterEvalQ(cl, library(phia)) clusterExport(cl, c("mvCom4", "maov", "lop"), envir=environment()) for(kk in 1:nSeg) { if(NoBrick > 1) out[,kk,] <- aperm(parApply(cl, inData[,kk,], 1, runAOV, dataframe=lop$dataStr, ModelForm=ModelForm, pars=pars), c(2,1)) else out[,kk,] <- aperm(parApply(cl, inData[,kk,], 1, runAOV, dataframe=lop$dataStr, ModelForm=ModelForm, pars=pars), dim=c(dimx_n, 1)) cat("Computation done ", 100*kk/nSeg, "%: ", format(Sys.time(), "%D %H:%M:%OS3"), "\n", sep='') } stopCluster(cl) } # convert to 4D dim(out) <- c(dimx_n*nSeg, 1, 1, NoBrick) # remove the trailers (padded 0s) out <- out[-c((dimx_n*nSeg-fill+1):(dimx_n*nSeg)), 1, 1,,drop=F] } else { # Initialization out <- array(0, dim=c(dimx, dimy, dimz, NoBrick)) # set up a voxel @ ii jj kk to test #runAOV(inData[ii, jj, kk,], dataframe=lop$dataStr, ModelForm=ModelForm, pars=pars) if (lop$nNodes==1) for (kk in 1:dimz) { if(NoBrick > 1) out[,,kk,] <- aperm(apply(inData[,,kk,], c(1,2), runAOV, dataframe=lop$dataStr, ModelForm=ModelForm, pars=pars), c(2,3,1)) else out[,,kk,] <- array(apply(inData[,,kk,], c(1,2), runAOV, dataframe=lop$dataStr, ModelForm=ModelForm, pars=pars), dim=c(dimx, dimy, 1)) cat("Z slice ", kk, "done: ", format(Sys.time(), "%D %H:%M:%OS3"), "\n") } if (lop$nNodes>1) { library(snow) cl <- makeCluster(lop$nNodes, type = "SOCK") clusterEvalQ(cl, library(afex)); clusterEvalQ(cl, library(phia)) clusterExport(cl, c("mvCom4", "maov", "lop"), envir=environment()) #clusterCall(cl, maov) # let all clusters access to function maov() #clusterExport(cl, c("maov"), envir=environment()) # let all clusters access to function maov() for (kk in 1:dimz) { if(NoBrick > 1) out[,,kk,] <- aperm(parApply(cl, inData[,,kk,], c(1,2), runAOV, dataframe=lop$dataStr, ModelForm=ModelForm, pars=pars), c(2,3,1)) else out[,,kk,] <- array(parApply(cl, inData[,,kk,], c(1,2), runAOV, dataframe=lop$dataStr, ModelForm=ModelForm, pars=pars), dim=c(dimx, dimy, 1)) cat("Z slice ", kk, "done: ", format(Sys.time(), "%D %H:%M:%OS3"), "\n") } stopCluster(cl) } } # test on Z slice #if (lop$nNodes>1) { # library(snow) # cl <- makeCluster(lop$nNodes, type = "SOCK") # clusterEvalQ(cl, library(afex)); #clusterEvalQ(cl, library(contrast)) # kk<- 32 # # out[,,kk,] <-aperm(parApply(cl, inData[,,kk,], c(1,2), runAOV, dataframe=lop$dataStr, ModelForm=ModelForm, pars=pars), c(2,3,1)) # cat("Z slice ", kk, "done: ", format(Sys.time(), "%D %H:%M:%OS3"), "\n") # stopCluster(cl) #} # avoid overflow Top <- 100 out[is.nan(out)] <- 0 out[out > Top] <- Top out[out < (-Top)] <- -Top ############################### #if(is.na(lop$mVar)) { #if(is.na(lop$wsVars)) outLabel <- paste(dimnames(uvfm)[[1]][2:(1+nF)], " F") else # if(pars[[6]][2]) { # sphericity correction for within-subject factors # outLabel <- NULL # for(ii in (1:nF)[-pars[[9]]+1]) outLabel <- append(outLabel, sprintf("%s F", dimnames(uvfm$lm)[[1]][ii+1])) # for(ii in (pars[[9]]-1)) outLabel <- append(outLabel, sprintf("%s Z", dimnames(uvfm[[1]])[[1]][ii+1])) # } else # outLabel <- paste(dimnames(uvfm[[1]])[[1]][-1], " F") #} else { #if(is.na(lop$wsVars)) outLabel <- paste(mvfm$terms[-1], " F") else # need improvement here for multivariate case # {} #} #for(ii in 1:lop$num_glt) { # outLabel <- append(outLabel, sprintf("%s", lop$gltLabel[ii])) # outLabel <- append(outLabel, sprintf("%s:t", lop$gltLabel[ii])) #} #statpar <- "3drefit" statsym <- NULL #for(ii in 1:nF) statpar <- paste(statpar, " -substatpar ", ii-1, " fift ", F_DF[[ii]][1], F_DF[[ii]][2]) if(lop$mvE4) for(ii in 1:nF) statsym <- c(statsym, list(list(sb=ii-1, typ="fict", par=1))) else if(lop$mvE4a) { for(ii in 1:(nF-nF_mvE4)) statsym <- c(statsym, list(list(sb=ii-1, typ="fift", par=c(F_DF[[ii]][1], F_DF[[ii]][2])))) for(ii in (nF-nF_mvE4+1):nF) statsym <- c(statsym, list(list(sb=ii, typ="fict", par=1))) } else for(ii in 1:nF) statsym <- c(statsym, list(list(sb=ii-1, typ="fift", par=c(F_DF[[ii]][1], F_DF[[ii]][2])))) if(lop$num_glt>0) for(ii in 1:lop$num_glt) statsym <- c(statsym, list(list(sb=nF+2*ii-1, typ="fitt", par=t_DF[ii]))) # statpar <- paste(statpar, " -substatpar ", nF+2*ii-1, " fitt", t_DF[ii]) #if(is.na(lop$mVar)) { # if(is.na(lop$wsVars)) for(ii in 1:nF) statpar <- paste(statpar, " -substatpar ", ii-1, " fift ", # uvfm[ii+1, 'Df'], " ", uvfm[2+nF, 'Df']) else if(pars[[6]][2]) { # sphericity correction for within-subject factors # for(ii in (1:nF)[-pars[[9]]+1]) statpar <- paste(statpar, " -substatpar ", ii-1, " fift ") # for(ii in (pars[[9]]-1)) statpar <- paste(statpar, " -substatpar ", ii-1, " fizt ") } else # for(ii in 1:nF) statpar <- paste(statpar, " -substatpar ", ii-1, " fift ", # unname(uvfm[[1]][ii+1,'num Df']), " ", unname(uvfm[[1]][ii+1,'den Df'])) #} else { # if(is.na(lop$wsVar)) { # for(ii in 1:nF) statpar <- paste(statpar, " -substatpar ", ii-1, " fift ", maov(mvfm)[3], maov(mvfm)[4]) # } else {} #} #if(lop$num_glt>0) for(ii in 1:lop$num_glt) statpar <- paste(statpar, " -substatpar ", nF+2*ii-1, " fitt", # ifelse(is.na(lop$wsVars) & is.na(lop$mVar), gltRes[[ii]][2,2], gltRes[[ii]][,6])) #statpar <- paste(statpar, " -addFDR -newid ", lop$outFN) #write.AFNI(lop$outFN, out, outLabel, note=myHist, origin=myOrig, delta=myDelta, idcode=newid.AFNI()) #write.AFNI(lop$outFN, out, outLabel, defhead=head, idcode=newid.AFNI(), type='MRI_short') #write.AFNI(lop$outFN, out, brickNames, defhead=head, idcode=newid.AFNI(), type='MRI_short') write.AFNI(lop$outFN, out, brickNames, defhead=head, idcode=newid.AFNI(), com_hist=lop$com_history, statsym=statsym, addFDR=1, type='MRI_short') cat("\nCongratulations! You have got an output ", lop$outFN, ".\n\n", sep='') } else { # if(is.numeric(lop$dataStr[, FileCol])): the last column is values insteadof input file names lop$dataStr$Beta <- as.numeric(lop$dataStr[, FileCol]) # convert characters to values lop$outFN <- paste(strsplit(lop$outFN, '\\+tlrc')[[1]], '.txt', sep='') capture.output(cat(''), file = lop$outFN, append = FALSE) if(iterPar %in% dimnames(lop$dataStr)[[2]]) nPar <- nlevels(as.factor(lop$dataStr[[iterPar]])) else nPar <- 1 for(nn in 1:nPar) { fm <- NULL if(nPar == 1) inData <- lop$dataStr else if((levels(as.factor(lop$dataStr[[iterPar]]))[nn] %in% lop$parSubsetVector) | is.null(lop$parSubsetVector)) inData <- lop$dataStr[lop$dataStr[[iterPar]] == levels(as.factor(lop$dataStr[[iterPar]]))[nn],] else inData <- NULL if(!is.null(inData)) try(fm <- aov.car(ModelForm, data=inData, factorize=FALSE, return='full'), silent=TRUE) else fm <-NULL #fm <- aov.car(ModelForm, data=inData, factorize=FALSE, return='full') if(!is.null(fm)) { #out_p <- apply(cbind(uvP[1:outTerms], uvP[(outTerms+1):length(uvP)], p_wsmvt[1:outTerms], # p_wsmvt[(outTerms+1):length(uvP)], p_mvt), 1, min) options(warn = -1) nF_mvE4 <- tryCatch(nrow(univ(fm$Anova)$anova)/2, error=function(e) NULL) if(is.null(nF_mvE4)) { tryCatch(fm2 <- aov(ModelForm2, data=inData), error=function(e) NULL) if(is.null(fm2)) errex.AFNI('Model failure...') else { nF_aov <- dim(summary(fm2)[[1]][[1]])[1]-1 out_p <- apply(cbind(summary(fm2)[[1]][[1]][1:nF_aov,'Pr(>F)'], summary(fm2)[[2]][[1]][2:(nF_aov+1),'Pr(>F)']), 1, min) names(out_p) <- dimnames(summary(fm2)[[1]][[1]])[[1]][1:nF_aov] } } else out_p <- mvCom4(fm, nF_mvE4) # chisq out_chisq <- qchisq(out_p, 1, lower.tail = FALSE) out <- cbind(out_chisq, 1, out_p) dimnames(out)[[2]] <- c('# Chisq', 'DF', 'Pr(>Chisq)') dimnames(out)[[1]] <- sprintf('# %s', dimnames(out)[[1]]) nC <- max(nchar(row.names(out))) term <- formatC(row.names(out), width=-nC) #or term <- sprintf("%-11s", row.names(out)) if(lop$num_glt>=1) { out_post <- matrix(0, nrow = lop$num_glt, ncol = 4) for(ii in 1:lop$num_glt) { if(all(is.na(lop$gltList[[ii]]))) glt <- tryCatch(testInteractions(fm$lm, pair=NULL, slope=lop$slpList[[ii]], covariates=lop$covVal, adjustment="none", idata = fm[["idata"]]), error=function(e) NULL) else glt <- tryCatch(testInteractions(fm$lm, custom=lop$gltList[[ii]], slope=lop$slpList[[ii]], covariates=lop$covVal, adjustment="none", idata = fm[["idata"]]), error=function(e) NULL) if(!is.null(glt)) { out_post[ii,1] <- glt[1,1] out_post[ii,2] <- sign(glt[1,1]) * sqrt(glt[1,4]) # convert F to t out_post[ii,3] <- glt[1,6] out_post[ii,4] <- glt[1,7] } #if(!is.null(glt)) } # for(ii in 1:lop$num_glt) dimnames(out_post)[[1]] <- sprintf('# %s', lop$gltLabel) dimnames(out_post)[[2]] <- c('# value', 't-stat', 'DF', '2-sided-P') nC2 <- max(nchar(row.names(out_post))) term2 <- formatC(row.names(out_post), width=-nC2) } # if(lop$num_glt>=1) options(width = 800) # include the width so that each line has enough capacity if(nPar==1) cat('# RESULTS: ANOVA table\n') else cat('# RESULTS: ANOVA table -', levels(as.factor(lop$dataStr[[iterPar]]))[nn], '\n') cat('-------------------------------------\n') print(setNames(data.frame(unname(out), term,stringsAsFactors=F), c(colnames(out), formatC("",width=-nC))), row.names=F) if(lop$num_glt>=1) { if(nPar==1) cat('\n# RESULTS: Post hoc tests\n') else cat('\n# RESULTS: Post hoc tests -', levels(as.factor(lop$dataStr[[iterPar]]))[nn], '\n') cat('-------------------------------------\n') print(setNames(data.frame(unname(out_post), term2,stringsAsFactors=F), c(colnames(out_post), formatC("",width=-nC2))), row.names=F) } #cat('-------------------------------------\n\n') cat('#####################################\n\n') if(nPar==1) capture.output(cat('\n# RESULTS: ANOVA table\n'), file = lop$outFN, append = TRUE) else capture.output(cat('\n# RESULTS: ANOVA table -', levels(as.factor(lop$dataStr[[iterPar]]))[nn], '\n'), file = lop$outFN, append = TRUE) capture.output(cat(dim(out)[1], '# Number of effects\n'), file = lop$outFN, append = TRUE) capture.output(print(setNames(data.frame(unname(out), term,stringsAsFactors=F), c(colnames(out), formatC("",width=-nC))), row.names=F), file = lop$outFN, append = TRUE) if(lop$num_glt>=1) { if(nPar==1) capture.output(cat('\n# RESULTS: Post hoc tests\n'), file = lop$outFN, append = TRUE) else capture.output(cat('\n# RESULTS: Post hoc tests -', levels(as.factor(lop$dataStr[[iterPar]]))[nn], '\n'), file = lop$outFN, append = TRUE) capture.output(cat(dim(out_post)[1], '# Number of tests\n'), file = lop$outFN, append = TRUE) capture.output(print(setNames(data.frame(unname(out_post), term2,stringsAsFactors=F), c(colnames(out_post), formatC("",width=-nC2))), row.names=F), file = lop$outFN, append = TRUE) } } # if(!is.null(fm)) } # for(nn in unique(lop$dataStr[[iterPar]])) cat("\nCongratulations! The above results are saved in file ", lop$outFN, ".\n\n", sep='') }