1 /*****************************************************************************
2 Major portions of this software are copyrighted by the Medical College
3 of Wisconsin, 1994-2001, and are released under the Gnu General Public
4 License, Version 2. See the file README.Copyright for details.
5 ******************************************************************************/
6
7 /*
8 This program calculates the single-factor analysis of variance (ANOVA)
9 for 3 dimensional AFNI data sets.
10
11 File: 3dANOVA.c
12 Author: B. D. Ward
13 Date: 09 December 1996
14
15 Mod: Incorporated include file 3dANOVA.h.
16 Date: 15 January 1997
17
18 Mod: Added option to check for required disk space.
19 Date: 23 January 1997
20
21 Mod: Added protection against divide by zero.
22 Date: 13 November 1997
23
24 Mod: Extensive changes required to implement the 'bucket' dataset.
25 Date: 30 December 1997
26
27 Mod: Separated 3dANOVA.h and 3dANOVA.lib files.
28 Date: 5 January 1998
29
30 Mod: Continuation of 'bucket' dataset changes.
31 Date: 9 January 1998
32
33 Mod: Added software for statistical tests of individual cell means.
34 Date: 27 October 1998
35
36 Mod: Added changes for incorporating History notes.
37 Date: 09 September 1999
38
39 Mod: Replaced dataset input code with calls to THD_open_dataset,
40 to allow operator selection of individual sub-bricks for input.
41 Date: 03 December 1999
42
43 Mod: Added call to AFNI_logger.
44 Date: 15 August 2001
45
46 Mod: Modified routine write_afni_data of 3dANOVA.lib so that all output
47 subbricks will now have the scaled short integer format.
48 Date: 14 March 2002
49
50 Mod: Set MAX_NAME_LENGTH equal to THD_MAX_NAME.
51 Date: 02 December 2002
52
53 Mod: Setup to use .1D dataset filenames on output if these are input.
54 Date: 14 March 2003 - RWCox
55
56 Mod: -help menu modified.
57 Date: 21 July 2005 - P Christidis
58
59 Mod: small -help correction
60 Date: 25 Nov 2005 [rickr]
61
62 Mod: Modified contrast t-stat computations, and added -old_method,
63 -OK, -assume_sph and -debug options.
64 Date: 08 Dec 2005 [rickr]
65 */
66
67 /*---------------------------------------------------------------------------*/
68
69 #define PROGRAM_NAME "3dANOVA" /* name of this program */
70 #define PROGRAM_AUTHOR "B. Douglas Ward" /* program author */
71 #define PROGRAM_INITIAL "09 Dec 1996" /* date of initial program release */
72 #define PROGRAM_LATEST "21 Jul 2005" /* date of latest program revision */
73
74 /*---------------------------------------------------------------------------*/
75
76 #if 0
77 #define SUFFIX ".3danova" /* suffix for temporary data files */
78 #else
79 char SUFFIX[1024] = ".3danova." ;
80 #endif
81
82 #include "3dANOVA.h"
83 #include "3dANOVA.lib"
84
85
86 /*---------------------------------------------------------------------------*/
87 /*
88 Routine to display 3dANOVA help menu.
89 */
90
display_help_menu(int detail)91 void display_help_menu(int detail)
92 {
93 printf
94 (
95 "This program performs single factor Analysis of Variance (ANOVA)\n"
96 "on 3D datasets\n"
97 "\n"
98 "---------------------------------------------------------------\n"
99 "\n"
100 "Usage:\n"
101 "-----\n"
102 "\n"
103 "3dANOVA\n"
104 " -levels r : r = number of factor levels\n"
105 "\n"
106 " -dset 1 filename : data set for factor level 1\n"
107 " . . .. . .\n"
108 " -dset 1 filename data set for factor level 1\n"
109 " . . .. . .\n"
110 " -dset r filename data set for factor level r\n"
111 " . . .. . .\n"
112 " -dset r filename data set for factor level r\n"
113 "\n"
114 " [-voxel num] : screen output for voxel # num\n"
115 "\n"
116 " [-diskspace] : print out disk space required for\n"
117 " program execution\n"
118 "\n"
119 " [-mask mset] : use sub-brick #0 of dataset 'mset'\n"
120 " to define which voxels to process\n"
121
122 "\n"
123 " [-debug level] : request extra output\n"
124 "\n"
125 "The following commands generate individual AFNI 2-sub-brick datasets:\n"
126 " (In each case, output is written to the file with the specified\n"
127 " prefix file name.)\n"
128 "\n"
129 " [-ftr prefix] : F-statistic for treatment effect\n"
130 "\n"
131 " [-mean i prefix] : estimate of factor level i mean\n"
132 "\n"
133 " [-diff i j prefix] : difference between factor levels\n"
134 "\n"
135 " [-contr c1...cr prefix] : contrast in factor levels\n"
136 "\n"
137 "Modified ANOVA computation options: (December, 2005)\n"
138 "\n"
139 " ** For details, see %s\n"
140 "\n"
141 "[-old_method] request to perform ANOVA using the previous\n"
142 " functionality (requires -OK, also)\n"
143 "\n"
144 "[-OK] confirm you understand that contrasts that\n"
145 " do not sum to zero have inflated t-stats, and\n"
146 " contrasts that do sum to zero assume sphericity\n"
147 " (to be used with -old_method)\n"
148 "\n"
149 "[-assume_sph] assume sphericity (zero-sum contrasts, only)\n"
150 "\n"
151 " This allows use of the old_method for\n"
152 " computing contrasts which sum to zero (this\n"
153 " includes diffs, for instance). Any contrast\n"
154 " that does not sum to zero is invalid, and\n"
155 " cannot be used with this option (such as\n"
156 " ameans).\n"
157 "\n"
158 "The following command generates one AFNI 'bucket' type dataset:\n"
159 "\n"
160 " [-bucket prefix] : create one AFNI 'bucket' dataset whose\n"
161 " sub-bricks are obtained by\n"
162 " concatenating the above output files;\n"
163 " the output 'bucket' is written to file\n"
164 " with prefix file name\n"
165 "\n", ANOVA_MODS_LINK);
166
167 printf
168 (
169 "N.B.: For this program, the user must specify 1 and only 1 sub-brick\n"
170 " with each -dset command. That is, if an input dataset contains\n"
171 " more than 1 sub-brick, a sub-brick selector must be used,\n"
172 " e.g., -dset 2 'fred+orig[3]'\n"
173 );
174
175 printf
176 ("\n"
177 "Example of 3dANOVA:\n"
178 "------------------\n"
179 "\n"
180 " Example is based on a study with one factor (independent variable)\n"
181 " called 'Pictures', with 3 levels:\n"
182 " (1) Faces, (2) Houses, and (3) Donuts\n"
183 "\n"
184 " The ANOVA is being conducted on the data of subjects Fred and Ethel:\n"
185 "\n"
186 " 3dANOVA -levels 3 \\\n"
187 " -dset 1 fred_Faces+tlrc \\\n"
188 " -dset 1 ethel_Faces+tlrc \\\n"
189 " \\\n"
190 " -dset 2 fred_Houses+tlrc \\\n"
191 " -dset 2 ethel_Houses+tlrc \\\n"
192 " \\\n"
193 " -dset 3 fred_Donuts+tlrc \\\n"
194 " -dset 3 ethel_Donuts+tlrc \\\n"
195 " \\\n"
196 " -ftr Pictures \\\n"
197 " -mean 1 Faces \\\n"
198 " -mean 2 Houses \\\n"
199 " -mean 3 Donuts \\\n"
200 " -diff 1 2 FvsH \\\n"
201 " -diff 2 3 HvsD \\\n"
202 " -diff 1 3 FvsD \\\n"
203 " -contr 1 1 -1 FHvsD \\\n"
204 " -contr -1 1 1 FvsHD \\\n"
205 " -contr 1 -1 1 FDvsH \\\n"
206 " -bucket fred_n_ethel_ANOVA\n"
207 );
208
209 printf("\n" MASTER_SHORTHELP_STRING );
210
211 printf("---------------------------------------------------\n"
212 "Also see HowTo#5 - Group Analysis on the AFNI website:\n"
213 "https://afni.nimh.nih.gov/pub/dist/HOWTO/howto/ht05_group/html/index.shtml\n"
214 "\n" );
215
216 printf(ANOVA_FLOAT_HELP) ;
217
218 PRINT_COMPILE_DATE; return;
219 }
220
221 /*---------------------------------------------------------------------------*/
222 /*
223 Routine to get user specified anova options.
224 */
225
get_options(int argc,char ** argv,anova_options * option_data)226 void get_options (int argc, char ** argv, anova_options * option_data)
227 {
228 int nopt = 1; /* input option argument counter */
229 int ival, rv; /* integer input and return val */
230 int i, j, k; /* factor level counter */
231 int nijk; /* number of data files in cell i */
232 float fval; /* float input */
233 THD_3dim_dataset * dset=NULL; /* test whether data set exists */
234 char message[MAX_NAME_LENGTH]; /* error message */
235
236
237
238 /*----- add to program log -----*/
239 AFNI_logger (PROGRAM_NAME,argc,argv);
240
241 /*----- initialize the input options -----*/
242 initialize_options (option_data);
243
244
245 /*----- main loop over input options -----*/
246 while (nopt < argc)
247 {
248 /*----- does user request help menu? -----*/
249 if (strcmp(argv[nopt], "-help") == 0 ||
250 strcmp(argv[nopt], "-h") == 0) {
251 display_help_menu(strlen(argv[nopt])>3?2:1);
252 exit(0);
253 }
254
255 /*----- allocate memory for storing data file names -----*/
256 if ((option_data->xname == NULL) && (option_data->a > 0))
257 {
258 option_data->xname =
259 (char *****) malloc (sizeof(char ****) * option_data->a);
260 for (i = 0; i < option_data->a; i++)
261 {
262 option_data->xname[i] =
263 (char ****) malloc (sizeof(char ***) * 1);
264 for (j = 0; j < 1; j++)
265 {
266 option_data->xname[i][j] =
267 (char ***) malloc (sizeof(char **) * 1);
268 for (k = 0; k < 1; k++)
269 {
270 option_data->xname[i][j][k] =
271 (char **) malloc (sizeof(char *) * MAX_OBSERVATIONS);
272 }
273 }
274 }
275 }
276
277
278 /*----- -diskspace -----*/
279 if( strncmp(argv[nopt],"-diskspace",5) == 0 )
280 {
281 option_data->diskspace = 1;
282 nopt++ ; continue ; /* go to next arg */
283 }
284
285
286 /*----- -debug level 8 Dec 2005 [rickr] -----*/
287 if( strncmp(argv[nopt],"-debug",4) == 0 ) {
288 if( ++nopt >= argc ) ANOVA_error("need a level after -debug!") ;
289 option_data->debug = atoi(argv[nopt]);
290 nopt++ ; continue ; /* go to next arg */
291 }
292
293
294 /*----- -datum type -----*/
295 if( strncmp(argv[nopt],"-datum",6) == 0 ){
296 if( ++nopt >= argc ) ANOVA_error("need an argument after -datum!") ;
297
298 if( strcmp(argv[nopt],"short") == 0 ){
299 option_data->datum = MRI_short ;
300 } else if( strcmp(argv[nopt],"float") == 0 ){
301 option_data->datum = MRI_float ;
302 } else {
303 char buf[256] ;
304 sprintf(buf,"-datum of type '%s' is not supported in 3dANOVA!",
305 argv[nopt] ) ;
306 ANOVA_error(buf) ;
307 }
308 nopt++ ; continue ; /* go to next arg */
309 }
310
311 /*------------------------------------------------------------*/
312 /*----- Using the old_method: 08 Dec 2005 [rickr] -----*/
313 /* if -old_method
314 if -OK, all contrasts are okay
315 else if -assume_sph, contrasts adding to 0 are okay
316 else complain and fail
317
318 bits: -old_method = 001, -OK = 010, -assume_sph = 100
319
320 valid bit patterns:
321 000 - use the new method
322 011 - use the old method
323 101 - use the old method (only allows zero-sum contrasts)
324 ------------------------------------------------------------*/
325
326 /*----- -old_method 23 Nov 2005 [rickr] -----*/
327 if (strncmp(argv[nopt], "-old_method", 6) == 0)
328 {
329 option_data->old_method |= 1;
330 nopt++;
331 continue;
332 }
333
334 /*----- -OK: denote both OK and old_method by old_method = 3 -----*/
335 if (strncmp(argv[nopt], "-OK", 3) == 0)
336 {
337 option_data->old_method |= 2;
338 nopt++;
339 continue;
340 }
341
342 /*----- -assume_sph: denote assume_sphericity by old_method = 4 ----*/
343 if (strncmp(argv[nopt], "-assume_sph", 11) == 0)
344 {
345 option_data->old_method |= 5; /* also set -old_method bit */
346 nopt++;
347 continue;
348 }
349
350 /*------- end old_method checks ------------------------------*/
351 /*------------------------------------------------------------*/
352
353
354 /*----- -session dirname -----*/
355 if( strncmp(argv[nopt],"-session",6) == 0 ){
356 nopt++ ;
357 if( nopt >= argc ) ANOVA_error("need argument after -session!") ;
358 strcpy(option_data->session , argv[nopt++]) ;
359 continue ;
360 }
361
362
363 /*----- -voxel num -----*/
364 if (strncmp(argv[nopt], "-voxel", 6) == 0)
365 {
366 nopt++;
367 if (nopt >= argc) ANOVA_error ("need argument after -voxel ");
368 rv = sscanf (argv[nopt], "%d", &ival);
369 if (rv < 1 || ival <= 0)
370 ANOVA_error ("illegal argument after -voxel ");
371 option_data->nvoxel = ival;
372 nopt++;
373 continue;
374 }
375
376
377 /*----- -levels r -----*/
378 if (strncmp(argv[nopt], "-levels", 7) == 0)
379 {
380 nopt++;
381 if (nopt >= argc) ANOVA_error ("need argument after -levels ");
382 rv = sscanf (argv[nopt], "%d", &ival);
383 if (rv < 1 || (ival <= 0) || (ival > MAX_LEVELS))
384 ANOVA_error ("illegal argument after -levels ");
385 option_data->a = ival;
386 nopt++;
387 continue;
388 }
389
390
391 /*----- -dset level filename -----*/
392 if (strncmp(argv[nopt], "-dset", 5) == 0)
393 {
394 nopt++;
395 if (nopt+1 >= argc) ANOVA_error ("need 2 arguments after -dset ");
396 rv = sscanf (argv[nopt], "%d", &ival);
397 if ((rv < 1) || (ival <= 0) || (ival > option_data->a))
398 ANOVA_error ("illegal argument after -dset ");
399
400 option_data->na[ival-1] += 1;
401 nijk = option_data->na[ival-1];
402 if (nijk > MAX_OBSERVATIONS)
403 ANOVA_error ("too many data files");
404
405 /*--- check whether input files exist ---*/
406 nopt++;
407 dset = THD_open_dataset( argv[nopt] ); CHECK_OPEN_ERROR(dset,argv[nopt]);
408
409 /*--- check number of selected sub-bricks ---*/
410 if (DSET_NVALS(dset) != 1)
411 {
412 sprintf(message,"Must specify exactly 1 sub-brick for file %s\n",
413 argv[nopt]);
414 ANOVA_error (message);
415 }
416
417 THD_delete_3dim_dataset( dset , False ) ; dset = NULL ;
418
419 option_data->xname[ival-1][0][0][nijk-1]
420 = malloc (sizeof(char) * MAX_NAME_LENGTH);
421 strcpy (option_data->xname[ival-1][0][0][nijk-1],
422 argv[nopt]);
423
424 nopt++;
425 continue;
426 }
427
428
429 /*----- -ftr filename -----*/
430 if (strncmp(argv[nopt], "-ftr", 4) == 0)
431 {
432 nopt++;
433 if (nopt >= argc) ANOVA_error ("need argument after -ftr ");
434 option_data->nftr = 1;
435 option_data->ftrname = malloc (sizeof(char) * MAX_NAME_LENGTH);
436 strcpy (option_data->ftrname, argv[nopt]);
437 nopt++;
438 continue;
439 }
440
441
442 /*----- -mean level filename -----*/
443 if (strncmp(argv[nopt], "-mean", 5) == 0)
444 {
445 nopt++;
446 if (nopt+1 >= argc) ANOVA_error ("need 2 arguments after -mean ");
447
448 option_data->num_ameans++;
449 if (option_data->num_ameans > MAX_MEANS)
450 ANOVA_error ("too many factor level mean estimates");
451
452 rv = sscanf (argv[nopt], "%d", &ival);
453 if ((rv < 1) || (ival <= 0) || (ival > option_data->a)) {
454 fprintf(stderr,"** bad opt #%d: %s\n", nopt, argv[nopt]);
455 ANOVA_error ("illegal argument after -mean ");
456 }
457 option_data->ameans[option_data->num_ameans-1] = ival - 1;
458 nopt++;
459
460 option_data->amname[option_data->num_ameans-1]
461 = malloc (sizeof(char) * MAX_NAME_LENGTH);
462 strcpy (option_data->amname[option_data->num_ameans-1], argv[nopt]);
463 nopt++;
464 continue;
465 }
466
467
468 /*----- -diff level1 level2 filename -----*/
469 if (strncmp(argv[nopt], "-diff", 5) == 0)
470 {
471 nopt++;
472 if (nopt+2 >= argc) ANOVA_error ("need 3 arguments after -diff ");
473
474 option_data->num_adiffs++;
475 if (option_data->num_adiffs > MAX_DIFFS)
476 ANOVA_error ("too many factor level differences");
477
478 rv = sscanf (argv[nopt], "%d", &ival);
479 if ((rv == 0) || (ival <= 0) || (ival > option_data->a)) {
480 fprintf(stderr,"** bad opt #%d: %s\n", nopt, argv[nopt]);
481 ANOVA_error ("illegal argument after -diff ");
482 }
483 option_data->adiffs[option_data->num_adiffs-1][0] = ival - 1;
484 nopt++;
485
486 rv = sscanf (argv[nopt], "%d", &ival);
487 if ((rv == 0) || (ival <= 0) || (ival > option_data->a)) {
488 fprintf(stderr,"** bad opt #%d: %s\n", nopt, argv[nopt]);
489 ANOVA_error ("illegal argument after -diff ");
490 }
491 option_data->adiffs[option_data->num_adiffs-1][1] = ival - 1;
492 nopt++;
493
494 option_data->adname[option_data->num_adiffs-1]
495 = malloc (sizeof(char) * MAX_NAME_LENGTH);
496 strcpy (option_data->adname[option_data->num_adiffs-1], argv[nopt]);
497 nopt++;
498 continue;
499 }
500
501
502 /*----- -contr c1 ... cr filename -----*/
503 if (strncmp(argv[nopt], "-contr", 6) == 0)
504 {
505 nopt++;
506 if (nopt + option_data->a >= argc)
507 ANOVA_error ("need r+1 arguments after -contr ");
508
509 option_data->num_acontr++;
510 if (option_data->num_acontr > MAX_CONTR)
511 ANOVA_error ("too many factor level contrasts");
512
513 for (i = 0; i < option_data->a; i++)
514 {
515 sscanf (argv[nopt], "%f", &fval);
516 option_data->acontr[option_data->num_acontr - 1][i] = fval ;
517 nopt++;
518 }
519
520 option_data->acname[option_data->num_acontr-1]
521 = malloc (sizeof(char) * MAX_NAME_LENGTH);
522 strcpy (option_data->acname[option_data->num_acontr-1], argv[nopt]);
523 nopt++;
524 continue;
525 }
526
527
528 /*----- -bucket filename -----*/
529 if (strncmp(argv[nopt], "-bucket", 4) == 0)
530 {
531 nopt++;
532 if (nopt >= argc) ANOVA_error ("need argument after -bucket ");
533 option_data->bucket_filename = malloc (sizeof(char)*MAX_NAME_LENGTH);
534 strcpy (option_data->bucket_filename, argv[nopt]);
535 nopt++;
536 continue;
537 }
538
539 /*----- -mask filename [11 Mar 2009: RWCox] -----*/
540 if( strncmp(argv[nopt],"-mask",5) == 0 )
541 {
542 THD_3dim_dataset *mset ; byte *amask ;
543 nopt++ ;
544 if( option_data->mask != NULL ) ANOVA_error("Can't have 2 -mask options");
545 if( nopt >= argc ) ANOVA_error("need argument after -mask" );
546 mset = THD_open_dataset( argv[nopt] ) ;
547 CHECK_OPEN_ERROR(mset,argv[nopt]) ;
548 DSET_load(mset) ; CHECK_LOAD_ERROR(mset) ;
549 amask = THD_makemask( mset , 0 , 1.0f , -1.0f ) ;
550 if( amask == NULL ){
551 WARNING_message("Can't create mask from dataset '%s'",argv[nopt]) ;
552 } else {
553 int nmvox = THD_countmask(DSET_NVOX(mset),amask) ;
554 if( nmvox < 1 ){
555 WARNING_message("Mask from dataset '%s' is empty",argv[nopt]) ;
556 free(amask) ;
557 } else {
558 INFO_message("Mask from dataset '%s' has %d voxels",argv[nopt],nmvox);
559 option_data->mask = amask ;
560 option_data->nmask = DSET_NVOX(mset) ;
561 }
562 }
563 DSET_delete(mset) ; nopt++ ; continue ;
564 }
565
566 /*----- unknown command -----*/
567 {
568 int bbot=nopt-3, btop=nopt+3, ib;
569 if( bbot < 0 ) bbot = 0;
570 if( btop >= argc ) btop = argc-1;
571 ERROR_message ("Unrecognized command line option #%d: %s\n",
572 nopt, argv[nopt]);
573 fprintf(stderr,"** bad opt part of:");
574 for( ib=bbot; ib<= btop; ib++ )
575 fprintf(stderr," %s", argv[ib]);
576 fputc('\n', stderr);
577 suggest_best_prog_option(argv[0], argv[nopt]);
578 exit(1);
579 }
580 }
581
582 if (argc < 2) {
583 ERROR_message("Too few options");
584 display_help_menu(0);
585 exit(1);
586 }
587
588 if (option_data->old_method == 1)
589 ANOVA_error("-old_method is insufficient by itself");
590 }
591
592 /*---------------------------------------------------------------------------*/
593 /*
594 Routine to check whether temporary files already exist.
595 */
596
check_temporary_files(anova_options * option_data)597 void check_temporary_files (anova_options * option_data)
598 {
599 char filename[MAX_NAME_LENGTH]; /* temporary file name */
600
601 int i; /* file counter */
602
603 for (i = 0; i < option_data->a; i++)
604 {
605 /*----- temporary file name -----*/
606 sprintf (filename, "y.%d", i);
607
608 /*----- see if file already exists -----*/
609 check_one_temporary_file (filename);
610 }
611
612
613 check_one_temporary_file ("ysum");
614 check_one_temporary_file ("ss");
615 check_one_temporary_file ("ssto");
616 check_one_temporary_file ("sstr");
617 check_one_temporary_file ("sse");
618
619 }
620
621
622 /*---------------------------------------------------------------------------*/
623 /*
624 Routine to check whether output files already exist.
625 */
626
check_output_files(anova_options * option_data)627 void check_output_files (anova_options * option_data)
628 {
629 int i; /* index */
630
631 if (option_data->nftr > 0)
632 check_one_output_file (option_data, option_data->ftrname);
633
634 if (option_data->num_ameans > 0)
635 for (i = 0; i < option_data->num_ameans; i++)
636 check_one_output_file (option_data, option_data->amname[i]);
637
638 if (option_data->num_adiffs > 0)
639 for (i = 0; i < option_data->num_adiffs; i++)
640 check_one_output_file (option_data, option_data->adname[i]);
641
642 if (option_data->num_acontr > 0)
643 for (i = 0; i < option_data->num_acontr; i++)
644 check_one_output_file (option_data, option_data->acname[i]);
645
646 if (option_data->bucket_filename != NULL)
647 check_one_output_file (option_data, option_data->bucket_filename);
648 }
649
650
651 /*---------------------------------------------------------------------------*/
652 /*
653 Routine to check for valid inputs.
654 */
655
check_for_valid_inputs(anova_options * option_data)656 void check_for_valid_inputs (anova_options * option_data)
657 {
658 int i; /* factor level index */
659 char message[MAX_NAME_LENGTH]; /* error message */
660
661 if (option_data->a < 2)
662 ANOVA_error ("must specify number of factor levels (r>1) ");
663
664 if (option_data->nt < option_data->a + 1)
665 ANOVA_error ("too few data sets for ANOVA");
666
667 for (i = 0; i < option_data->a; i++)
668 if (option_data->na[i] == 0)
669 {
670 sprintf (message, "level %d has too few data sets", i+1);
671 ANOVA_error (message);
672 }
673
674 /* check contrasts (show error, and specify 3dANOVA) */
675 if ( !contrasts_are_valid(option_data, 1, 1) )
676 ANOVA_error("invalid contrast(s)\n");
677 }
678 /*---------------------------------------------------------------------------*/
679 /*
680 Routine to calculate the number of data files that have to be stored.
681 Modified to account for 'bucket' dataset output.
682 */
683
required_data_files(anova_options * option_data)684 int required_data_files (anova_options * option_data)
685 {
686 int r; /* number of factor levels */
687 int nftr; /* number of F-treatment output files
688 (0 or 1) */
689 int nmeans; /* number of estimated mean output files */
690 int ndiffs; /* number of difference output files */
691 int ncontr; /* number of contrast output files */
692 int nout; /* number of output files */
693 int nmax; /* maximum number of disk files */
694
695
696 /*----- initialize local variables -----*/
697 r = option_data->a;
698 nftr = option_data->nftr;
699 nmeans = option_data->num_ameans;
700 ndiffs = option_data->num_adiffs;
701 ncontr = option_data->num_acontr;
702 nout = nftr + nmeans + ndiffs + ncontr;
703
704 nmax = r + 2 + nout;
705 if (option_data->bucket_filename != NULL)
706 nmax = max (nmax, 2*nout);
707
708 return (nmax);
709 }
710
711
712 /*---------------------------------------------------------------------------*/
713 /*
714 Routine to perform all ANOVA initialization.
715 */
716
initialize(int argc,char ** argv,anova_options ** option_data)717 void initialize (int argc, char ** argv, anova_options ** option_data)
718 {
719 int i; /* factor level index */
720
721
722
723 /*----- save command line for history notes -----*/
724 commandline = tross_commandline( PROGRAM_NAME , argc,argv ) ;
725
726
727 /*----- allocate memory space for input data -----*/
728 *option_data = (anova_options *) malloc(sizeof(anova_options));
729 if (*option_data == NULL)
730 ANOVA_error ("memory allocation error");
731
732 /*----- get command line inputs -----*/
733 get_options(argc, argv, *option_data);
734
735 /*----- use first data set to get data set dimensions -----*/
736 (*option_data)->first_dataset = (*option_data)->xname[0][0][0][0];
737 get_dimensions (*option_data);
738 printf ("Data set dimensions: nx = %d ny = %d nz = %d nxyz = %d \n",
739 (*option_data)->nx, (*option_data)->ny,
740 (*option_data)->nz, (*option_data)->nxyz);
741 if ((*option_data)->nvoxel > (*option_data)->nxyz)
742 ANOVA_error ("argument of -voxel is too large");
743
744 /*----- count number of observations per treatment level -----*/
745 (*option_data)->nt = 0;
746 for (i = 0; i < (*option_data)->a; i++)
747 (*option_data)->nt += (*option_data)->na[i];
748
749 /*----- check for valid inputs -----*/
750 check_for_valid_inputs (*option_data);
751
752 /*----- check whether temporary files already exist -----*/
753 check_temporary_files (*option_data);
754
755 /*----- check whether output files already exist -----*/
756 if( THD_deathcon() ) check_output_files (*option_data);
757
758 /*----- check whether there is sufficient disk space -----*/
759 if ((*option_data)->diskspace) check_disk_space (*option_data);
760 }
761
762 /*---------------------------------------------------------------------------*/
763 /* 8 Dec 2005 [rickr]
764 routine to compute the degrees of freedom from a contrast
765
766 df = sum_i_in_contr( n_i - 1 )
767
768 */
df_for_contr(anova_options * option_data,float * contr)769 int df_for_contr(anova_options * option_data, float * contr)
770 {
771 int c, df = 0;
772
773 for (c = 0; c < option_data->a; c++)
774 if (contr[c] != 0.0)
775 df += (option_data->na[c] - 1);
776
777 return df;
778 }
779
780
781 /*---------------------------------------------------------------------------*/
782 /* 7 Dec 2005 [rickr,gangc]
783 routine to compute the mean and tstat of an ANOVA contrast
784
785 t = L / sqrt(1/df * SL2 * sum_c2)
786 L = sum_i(contr[i] * Y_i...)
787 SL = sum_i[ step(|contr[i]|) * ( sum_j(Y_ij^2) - n_i*Y_i.^2 ) ]
788 sum_c2 = sum_i[ contr[i]^2/n_i ]
789
790 note: Y_i. is considered the mean Y_i here
791
792 For accuracy, sums are computed using doubles, then copied to float.
793 The contrast is passed in to allow for more uses of this function.
794 */
calc_contr(anova_options * option_data,float * contr,float * mean,float * tmean)795 void calc_contr(anova_options * option_data, float * contr,
796 float * mean, /* save memory: use to read files */
797 float * tmean /* save memory: use to sum over obs. */
798 )
799 {
800 const float EPSILON = 1.0e-15; /* protect against divide by zero */
801 double * L; /* cumulative contrast mean */
802 double * S1; /* sum(n_i*Y_i^2), for mean Y_i */
803 double * S2; /* sum(Y_ij^2) {SL = S2-S1} */
804 double sum_c2; /* sum_i(c_i^2/n_i) */
805 double denom; /* for computations */
806 int i, j; /* indices for levels of factors A and C */
807 int a; /* number of levels */
808 int n; /* number of observations per cell */
809 int df; /* degrees of freedom */
810 int ixyz, nxyz; /* voxel counters */
811 int nvoxel; /* output voxel # */
812
813 /*----- initialize local variables -----*/
814 a = option_data->a;
815 nxyz = option_data->nxyz;
816 nvoxel = option_data->nvoxel;
817
818 /*----- allocate memory space for calculations -----*/
819 L = (double *) malloc(sizeof(double)*nxyz);
820 S1 = (double *) malloc(sizeof(double)*nxyz);
821 S2 = (double *) malloc(sizeof(double)*nxyz);
822 if (!L || !S1 || !S2)
823 ANOVA_error ("calc_contr: unable to allocate sufficient memory");
824
825 for (ixyz = 0; ixyz < nxyz; ixyz++) /* init cumulative sums */
826 L[ixyz] = S1[ixyz] = S2[ixyz] = 0.0;
827 sum_c2 = 0.0;
828
829 if (option_data->debug > 1){
830 fprintf(stderr,"-d contr (%d levels) = ",a);
831 for ( i = 0; i < a; i++ ) fprintf(stderr,"%f ",contr[i]);
832 fprintf(stderr,"\n-d observations = ");
833 for ( i = 0; i < a; i++ ) fprintf(stderr,"%d ",option_data->na[i]);
834 fputc('\n',stderr);
835 }
836
837 /*----- loop over factor levels -----*/
838 for ( i = 0; i < a; i++ )
839 {
840 if (contr[i] == 0.0 ) continue; /* then skip this index */
841
842 n = option_data->na[i]; /* num observations */
843
844 for (ixyz = 0; ixyz < nxyz; ixyz++) /* clear sum over observations */
845 tmean[ixyz] = 0.0;
846
847 /*----- process observations within treatment level -----*/
848 for (j = 0; j < n; j++)
849 {
850 read_afni_data (option_data, option_data->xname[i][0][0][j], mean);
851
852 for (ixyz = 0; ixyz < nxyz; ixyz++)
853 {
854 tmean[ixyz] += mean[ixyz]; /* sum over obs. */
855 S2[ixyz] += (double)mean[ixyz]*mean[ixyz]; /* sum of squares */
856 }
857 }
858
859 /*----- tally sum of contrast and squares for the current k -----*/
860 for (ixyz = 0; ixyz < nxyz; ixyz++)
861 {
862 S1[ixyz] += (double)tmean[ixyz]*tmean[ixyz]/n;/* is n_i*mean_Y_i^2 */
863 L[ixyz] += (double)contr[i]*tmean[ixyz]/n; /* cum. contrast mean */
864 }
865 if (option_data->debug > 1 && nvoxel > 0)
866 fprintf(stderr,"+d mean[%d] = %f, cmean = %f\n",
867 i,tmean[nvoxel-1]/n,L[nvoxel-1]);
868 }
869
870 /* get df and sum_c2 (seems more readable to put in separate loop) */
871 df = 0;
872 for ( i = 0; i < a; i++ )
873 {
874 if (contr[i] == 0.0 ) continue; /* then skip this index */
875 n = option_data->na[i]; /* num observations */
876 sum_c2 += (double)contr[i]*contr[i]/n; /* see above */
877 df += n - 1; /* degrees of freedom */
878 }
879
880 /*----- copy final results to float output -----*/
881 for (ixyz = 0; ixyz < nxyz; ixyz++)
882 {
883 denom = sqrt((S2[ixyz] - S1[ixyz]) * sum_c2 / df); /* t denominator */
884 mean [ixyz] = L[ixyz];
885 tmean[ixyz] = (denom < EPSILON) ? 0.0 : L[ixyz] / denom;
886 }
887 if (option_data->debug > 1 && nvoxel > 0)
888 fprintf(stderr,"+d S1, S2, sum_c2, df = %f, %f, %f, %d\n",
889 S1[nvoxel-1], S2[nvoxel-1], sum_c2, df);
890
891 /*----- fly, be free! {thud} free! {thud} -----*/
892 free(L); free(S1); free (S2);
893 }
894
895
896 /*---------------------------------------------------------------------------*/
897 /*
898 Routine to sum over the observations within one treatment level.
899 Output is stored (temporarily) in data file y"i".3danova, where
900 "i" = treatment level.
901 */
902
calculate_y(anova_options * option_data)903 void calculate_y (anova_options * option_data)
904 {
905 float * x = NULL; /* pointer to original input data */
906 float * y = NULL; /* pointer to sum within treatment data */
907 int i; /* factor level index */
908 int j; /* observation number index */
909 int r; /* number of factor levels (treatments) */
910 int ixyz, nxyz; /* voxel counters */
911 int nvoxel; /* output voxel # */
912 char filename[MAX_NAME_LENGTH]; /* name of file for sum within treatment */
913
914
915 /*----- initialize local variables -----*/
916 r = option_data->a;
917 nxyz = option_data->nxyz;
918 nvoxel = option_data->nvoxel;
919
920 /*----- allocate memory space for calculations -----*/
921 x = (float *) malloc(sizeof(float)*nxyz);
922 y = (float *) malloc(sizeof(float)*nxyz);
923 if ((x == NULL) || (y == NULL))
924 ANOVA_error ("unable to allocate sufficient memory");
925
926
927 /*----- loop over treatment levels -----*/
928 for (i = 0; i < r; i++)
929 {
930 /*----- sum observations within a treatment level -----*/
931 volume_zero (y, nxyz);
932 for (j = 0; j < option_data->na[i]; j++)
933 {
934 read_afni_data (option_data,
935 option_data->xname[i][0][0][j], x);
936 if (nvoxel > 0)
937 printf ("y[%d][%d] = %f \n", i+1, j+1, x[nvoxel-1]);
938 for (ixyz = 0; ixyz < nxyz; ixyz++)
939 y[ixyz] += x[ixyz];
940 }
941
942 /*----- save the sum for this treatment level -----*/
943 if (nvoxel > 0)
944 printf ("y[%d]. = %f \n", i+1, y[nvoxel-1]);
945 sprintf (filename, "y.%d", i);
946 volume_write (filename, y, nxyz);
947 }
948
949 /*----- release memory -----*/
950 free (y); y = NULL;
951 free (x); x = NULL;
952
953 }
954
955
956 /*---------------------------------------------------------------------------*/
957 /*
958 Routine to sum over all observations at all treatment levels.
959 The sum is stored (temporarily) in disk file ysum.3danova.
960 */
961
calculate_ysum(anova_options * option_data)962 void calculate_ysum (anova_options * option_data)
963 {
964 float * y = NULL; /* pointer to sum within treatment data */
965 float * ysum = NULL; /* pointer to overall sum data */
966 int i; /* factor level index */
967 int ixyz, nxyz; /* voxel counters */
968 int nvoxel; /* output voxel # */
969 char filename[MAX_NAME_LENGTH]; /* sum within treatment input file name */
970
971
972 /*----- initialize local variables -----*/
973 nxyz = option_data->nxyz;
974 nvoxel = option_data->nvoxel;
975
976 /*----- allocate memory space for calculations -----*/
977 y = (float *) malloc(sizeof(float)*nxyz);
978 ysum = (float *) malloc(sizeof(float)*nxyz);
979 if ((y == NULL) || (ysum == NULL))
980 ANOVA_error ("unable to allocate sufficient memory");
981
982
983 /*----- sum over all observations -----*/
984 volume_zero (ysum, nxyz);
985 for (i = 0; i < option_data->a; i++)
986 {
987 sprintf (filename, "y.%d", i);
988 volume_read (filename, y, nxyz);
989 for (ixyz = 0; ixyz < nxyz; ixyz++)
990 ysum[ixyz] += y[ixyz];
991 }
992
993 if (nvoxel > 0)
994 printf ("y.. = %f \n", ysum[nvoxel-1]);
995 volume_write ("ysum",ysum, nxyz);
996
997 /*----- release memory -----*/
998 free (ysum); ysum = NULL;
999 free (y); y = NULL;
1000
1001 }
1002
1003
1004 /*---------------------------------------------------------------------------*/
1005 /*
1006 Routine to calculate the sum of squares of all observations (SS).
1007 The output is stored (temporarily) in disk file ss.3danova.
1008 */
1009
calculate_ss(anova_options * option_data)1010 void calculate_ss (anova_options * option_data)
1011 {
1012 float * x = NULL; /* pointer to original input data */
1013 float * ss = NULL; /* pointer to sum of squares data */
1014 int i; /* factor level index */
1015 int j; /* observation data index */
1016 int r; /* number of factor levels */
1017 int ixyz, nxyz; /* voxel counters */
1018 int nvoxel; /* output voxel # */
1019 float xval; /* float input data */
1020
1021
1022 /*----- initialize local variables -----*/
1023 r = option_data->a;
1024 nxyz = option_data->nxyz;
1025 nvoxel = option_data->nvoxel;
1026
1027 /*----- allocate memory space for calculations -----*/
1028 x = (float *) malloc(sizeof(float)*nxyz);
1029 ss = (float *) malloc(sizeof(float)*nxyz);
1030 if ((x == NULL) || (ss == NULL))
1031 ANOVA_error ("unable to allocate sufficient memory");
1032
1033 /*----- sum the squares of all observations -----*/
1034 volume_zero (ss, nxyz);
1035 for (i = 0; i < r; i++)
1036 for (j = 0; j < option_data->na[i]; j++)
1037 {
1038 read_afni_data (option_data,
1039 option_data->xname[i][0][0][j], x);
1040 for (ixyz = 0; ixyz < nxyz; ixyz++)
1041 {
1042 xval = x[ixyz];
1043 ss[ixyz] += xval * xval;
1044 }
1045 }
1046
1047 if (nvoxel > 0)
1048 printf ("SS = %f \n", ss[nvoxel-1]);
1049 volume_write ("ss", ss, nxyz);
1050
1051 /*----- release memory -----*/
1052 free (ss); ss = NULL;
1053 free (x); x = NULL;
1054
1055 }
1056
1057 /*---------------------------------------------------------------------------*/
1058 /*
1059 Routine to calculate total (corrected for the mean) sum of squares (SSTO).
1060 The output is stored (temporarily) in disk file ssto.3danova.
1061 */
1062
calculate_ssto(anova_options * option_data)1063 void calculate_ssto (anova_options * option_data)
1064 {
1065 float * ysum = NULL; /* ptr to sum over all observations */
1066 float * ssto = NULL; /* pointer to ssto data */
1067 int ixyz, nxyz; /* voxel counters */
1068 int nvoxel; /* output voxel # */
1069 int nt; /* total number of observations */
1070 float yval; /* sum over all observations */
1071
1072
1073 /*----- initialize local variables -----*/
1074 nt = option_data->nt;
1075 nxyz = option_data->nxyz;
1076 nvoxel = option_data->nvoxel;
1077
1078 /*----- allocate memory space for calculations -----*/
1079 ysum = (float *) malloc(sizeof(float)*nxyz);
1080 ssto = (float *) malloc(sizeof(float)*nxyz);
1081 if ((ysum == NULL) || (ssto == NULL))
1082 ANOVA_error ("unable to allocate sufficient memory");
1083
1084 /*----- SSTO = SS - ((YSUM)^2)/nt -----*/
1085 volume_read ("ss", ssto, nxyz);
1086 volume_read ("ysum", ysum, nxyz);
1087 for (ixyz = 0; ixyz < nxyz; ixyz++)
1088 {
1089 yval = ysum[ixyz];
1090 ssto[ixyz] -= yval * yval / nt;
1091 }
1092
1093 /*----- ss data file is no longer needed -----*/
1094 volume_delete ("ss");
1095
1096 /*----- save total (corrected) sum of squares -----*/
1097 if (nvoxel > 0)
1098 printf ("SSTO = %f \n", ssto[nvoxel-1]);
1099 volume_write ("ssto", ssto, nxyz);
1100
1101 /*----- release memory -----*/
1102 free (ssto); ssto = NULL;
1103 free (ysum); ysum = NULL;
1104
1105 }
1106
1107 /*---------------------------------------------------------------------------*/
1108 /*
1109 Routine to calculate the sum of squares due to the treatment (SSTR).
1110 The output is stored (temporarily) in file sstr.3danova.
1111 */
1112
calculate_sstr(anova_options * option_data)1113 void calculate_sstr (anova_options * option_data)
1114 {
1115 float * y = NULL; /* input data pointer */
1116 float * sstr = NULL; /* sstr data pointer */
1117 int i; /* factor level index */
1118 int ixyz, nxyz; /* voxel counters */
1119 int nvoxel; /* output voxel # */
1120 int ni; /* number of observations at level i */
1121 int nt; /* total number of observations */
1122 float yval; /* temporary float value */
1123 char filename[MAX_NAME_LENGTH]; /* input data file name */
1124
1125
1126 /*----- assign local variables -----*/
1127 nt = option_data->nt;
1128 nxyz = option_data->nxyz;
1129 nvoxel = option_data->nvoxel;
1130
1131 /*----- allocate memory space for calculations -----*/
1132 y = (float *) malloc(sizeof(float)*nxyz);
1133 sstr = (float *) malloc(sizeof(float)*nxyz);
1134 if ((y == NULL) || (sstr == NULL))
1135 ANOVA_error ("unable to allocate sufficient memory");
1136
1137 volume_zero (sstr, nxyz);
1138 for (i = 0; i < option_data->a; i++)
1139 {
1140 ni = option_data->na[i];
1141 sprintf (filename, "y.%d", i);
1142 volume_read (filename, y, nxyz);
1143 for (ixyz = 0; ixyz < nxyz; ixyz++)
1144 {
1145 yval = y[ixyz];
1146 sstr[ixyz] += yval * yval / ni;
1147 }
1148 }
1149
1150 volume_read ("ysum", y, nxyz);
1151 for (ixyz = 0; ixyz < nxyz; ixyz++)
1152 {
1153 yval = y[ixyz];
1154 sstr[ixyz] -= yval * yval / nt;
1155
1156 /*----- protection against round-off error -----*/
1157 if (sstr[ixyz] < 0.0) sstr[ixyz] = 0.0;
1158 }
1159
1160 /*----- ysum data file is no longer needed -----*/
1161 volume_delete ("ysum");
1162
1163 /*----- save treatment sum of squares -----*/
1164 if (nvoxel > 0)
1165 printf ("SSTR = %f \n", sstr[nvoxel-1]);
1166 volume_write ("sstr", sstr, nxyz);
1167
1168 /*----- release memory -----*/
1169 free (sstr); sstr = NULL;
1170 free (y); y = NULL;
1171
1172 }
1173
1174
1175 /*---------------------------------------------------------------------------*/
1176 /*
1177 Routine to calculate the error sum of squares, SSE = SSTO - SSTR .
1178 The result is stored (temporarily) in disk file sse.3danova.
1179 */
1180
calculate_sse(anova_options * option_data)1181 void calculate_sse (anova_options * option_data)
1182 {
1183 float * sstr = NULL; /* pointer to sstr data */
1184 float * sse = NULL; /* pointer to sse data */
1185 int ixyz, nxyz; /* voxel counters */
1186 int nvoxel; /* output voxel # */
1187
1188
1189 /*----- assign local variables -----*/
1190 nxyz = option_data->nxyz;
1191 nvoxel = option_data->nvoxel;
1192
1193 /*----- allocate memory space for calculations -----*/
1194 sstr = (float *) malloc(sizeof(float)*nxyz);
1195 sse = (float *) malloc(sizeof(float)*nxyz);
1196 if ((sstr == NULL) || (sse == NULL))
1197 ANOVA_error ("unable to allocate sufficient memory");
1198
1199 /*----- read SSTO (total sum of squares) -----*/
1200 volume_read ("ssto", sse, nxyz);
1201
1202 /*----- read SSTR (treatment sum of squares) -----*/
1203 volume_read ("sstr", sstr, nxyz);
1204 for (ixyz = 0; ixyz < nxyz; ixyz++)
1205 {
1206 /*----- SSE = SSTO - SSTR -----*/
1207 sse[ixyz] -= sstr[ixyz];
1208
1209 /*----- protection against round-off error -----*/
1210 if (sse[ixyz] < 0.0) sse[ixyz] = 0.0;
1211 }
1212
1213 /*----- ssto data file is no longer needed -----*/
1214 volume_delete ("ssto");
1215
1216 /*----- save error sum of squares -----*/
1217 if (nvoxel > 0)
1218 printf ("SSE = %f \n", sse[nvoxel-1]);
1219 volume_write ("sse", sse, nxyz);
1220
1221 /*----- release memory -----*/
1222 free (sse); sse = NULL;
1223 free (sstr); sstr = NULL;
1224
1225 }
1226
1227
1228 /*---------------------------------------------------------------------------*/
1229 /*
1230 Routine to calculate the F-statistic for treatment, F = MSTR / MSE,
1231 where MSTR = SSTR / (r-1),
1232 and MSE = SSE / (n-r).
1233
1234 The output is stored as a 2 sub-brick AFNI data set. The first sub-brick
1235 contains the square root of MSTR (mean sum of squares due to treatment),
1236 and the second sub-brick contains the corresponding F-statistic.
1237 */
1238
calculate_f_statistic(anova_options * option_data)1239 void calculate_f_statistic (anova_options * option_data)
1240 {
1241 const float EPSILON = 1.0e-10; /* protect against divide by zero */
1242 float * fin = NULL; /* pointer to input float data */
1243 float * mstr = NULL; /* pointer to MSTR data */
1244 float * ftr = NULL; /* pointer to F due-to-treatment */
1245 int r; /* number of factor levels */
1246 int nt; /* total number of observations */
1247 int ixyz, nxyz; /* voxel counters */
1248 int nvoxel; /* output voxel # */
1249 float fval; /* float value used in calculations */
1250 float mse; /* mean square error */
1251
1252
1253 /*----- initialize local variables -----*/
1254 r = option_data->a;
1255 nt = option_data->nt;
1256 nxyz = option_data->nxyz;
1257 nvoxel = option_data->nvoxel;
1258
1259 /*----- allocate memory space for calculations -----*/
1260 /*----- Note: the order in which memory allocations take place
1261 is important! -----*/
1262 ftr = (float *) malloc(sizeof(float)*nxyz);
1263 mstr = (float *) malloc(sizeof(float)*nxyz);
1264 fin = (float *) malloc(sizeof(float)*nxyz);
1265 if ((fin == NULL) || (ftr == NULL) || (mstr == NULL))
1266 ANOVA_error ("unable to allocate sufficient memory");
1267
1268
1269 /*----- calculate mean SS due to treatments -----*/
1270 volume_read ("sstr", fin, nxyz);
1271 fval = 1.0 / (r - 1.0);
1272 for (ixyz = 0; ixyz < nxyz; ixyz++)
1273 mstr[ixyz] = fin[ixyz] * fval; /*--- MSTR = SSTR / (r-1) ---*/
1274 if (nvoxel > 0)
1275 printf ("MSTR = %f \n", mstr[nvoxel-1]);
1276
1277 /*----- calculate F-statistic FTR = MSTR / MSE -----*/
1278 volume_read ("sse", fin, nxyz);
1279 fval = 1.0 / (nt - r);
1280 for (ixyz = 0; ixyz < nxyz; ixyz++)
1281 {
1282 mse = fin[ixyz] * fval; /*--- MSE = SSE / (nt-r) ---*/
1283 if (mse < EPSILON)
1284 ftr[ixyz] = 0.0;
1285 else
1286 ftr[ixyz] = mstr[ixyz] / mse; /*--- FTR = MSTR / MSE ---*/
1287 }
1288 if (nvoxel > 0)
1289 printf ("FTR = %f \n", ftr[nvoxel-1]);
1290
1291 /*----- release memory -----*/
1292 free (fin); fin = NULL;
1293
1294 /*----- write out afni data file -----*/
1295 for (ixyz = 0; ixyz < nxyz; ixyz++)
1296 mstr[ixyz] = sqrt(mstr[ixyz]); /*-- mstr now holds square root --*/
1297 write_afni_data (option_data, option_data->ftrname, mstr, ftr, r-1, nt-r);
1298
1299 /*----- this data file is no longer needed -----*/
1300 volume_delete ("sstr");
1301
1302 /*----- release memory -----*/
1303 free (mstr); mstr = NULL;
1304 free (ftr); ftr = NULL;
1305
1306 }
1307
1308
1309 /*---------------------------------------------------------------------------*/
1310 /*
1311 Routine to calculate the mean treatment effect at the user specified
1312 treatment level. The output is stored as a 2 sub-brick AFNI data set.
1313 The first sub-brick contains the estimated mean at this treatment level,
1314 and the second sub-brick contains the corresponding t-statistic.
1315 */
1316
calculate_means(anova_options * option_data)1317 void calculate_means (anova_options * option_data)
1318 {
1319 const float EPSILON = 1.0e-10; /* protect against divide by zero */
1320 float * fin = NULL; /* pointer to float input data */
1321 float * mean = NULL; /* pointer to treatment mean data */
1322 float * tmean = NULL; /* pointer to t-statistic data */
1323 float * contr = NULL; /* for new_method contrast */
1324 int imean; /* output mean option index */
1325 int level; /* factor level index */
1326 int ni; /* number obs. at factor level i */
1327 int ixyz, nxyz; /* voxel counters */
1328 int nvoxel; /* output voxel # */
1329 int r; /* number of factor levels */
1330 int nt; /* total number of observations */
1331 int df; /* degrees of freedom */
1332 int num_means; /* number of user requested means */
1333 float fval; /* for calculating std. dev. */
1334 float stddev; /* est. std. dev. of factor mean */
1335 char filename[MAX_NAME_LENGTH]; /* input data file name */
1336
1337
1338 /*----- initialize local variables -----*/
1339 r = option_data->a;
1340 nt = option_data->nt;
1341 num_means = option_data->num_ameans;
1342 nxyz = option_data->nxyz;
1343 nvoxel = option_data->nvoxel;
1344
1345 /*----- allocate memory space for calculations -----*/
1346 mean = (float *) malloc(sizeof(float)*nxyz);
1347 tmean = (float *) malloc(sizeof(float)*nxyz);
1348 contr = (float *) malloc(sizeof(float)*r);
1349 if ((mean == NULL) || (tmean == NULL) || (contr == NULL))
1350 ANOVA_error ("unable to allocate sufficient memory");
1351
1352 /*----- loop over user specified treatment means -----*/
1353 for (imean = 0; imean < num_means; imean++)
1354 {
1355 level = option_data->ameans[imean];
1356 ni = option_data->na[level];
1357
1358 if (option_data->old_method) /* old method 8 Dec 2005 [rickr] */
1359 {
1360 df = nt-r;
1361
1362 /*----- allocate temporary memory space -----*/
1363 fin = (float *) malloc(sizeof(float)*nxyz);
1364 if (fin == NULL)
1365 ANOVA_error ("unable to allocate sufficient memory");
1366
1367 /*----- estimate factor mean for this treatment level -----*/
1368 sprintf (filename, "y.%d", level);
1369 volume_read (filename, fin, nxyz);
1370 for (ixyz = 0; ixyz < nxyz; ixyz++)
1371 mean[ixyz] = fin[ixyz] / ni;
1372
1373 /*----- divide by estimated standard deviation of factor mean -----*/
1374 volume_read ("sse", fin, nxyz);
1375 fval = (1.0 / (nt-r)) * (1.0 / ni);
1376 for (ixyz = 0; ixyz < nxyz; ixyz++)
1377 {
1378 stddev = sqrt(fin[ixyz] * fval);
1379 if (stddev < EPSILON)
1380 tmean[ixyz] = 0.0;
1381 else
1382 tmean[ixyz] = mean[ixyz] / stddev;
1383 }
1384
1385 /*----- release memory -----*/
1386 free (fin); fin = NULL;
1387 } else { /* new method using contrast (for all cases) */
1388 for(ixyz = 0; ixyz < r; ixyz ++) contr[ixyz] = 0.0; /* clear contr */
1389 contr[level] = 1.0;
1390 calc_contr(option_data, contr, mean, tmean);
1391 df = df_for_contr(option_data, contr);
1392 }
1393
1394 if (nvoxel > 0){
1395 printf ("Mean [%d] = %f \n", level+1, mean[nvoxel-1]);
1396 printf ("t-Mean [%d] = %f, df = %d\n", level+1, tmean[nvoxel-1], df);
1397 }
1398
1399 /*----- write out afni data file -----*/
1400 write_afni_data (option_data, option_data->amname[imean],
1401 mean, tmean, df, 0);
1402 }
1403
1404 /*----- release memory -----*/
1405 free (tmean); free (mean); free (contr);
1406
1407 }
1408
1409
1410 /*---------------------------------------------------------------------------*/
1411 /*
1412 Routine to estimate the difference in the means between two user specified
1413 treatment levels. The output is a 2 sub-brick AFNI data set. The first
1414 sub-brick contains the estimated difference in the means. The second
1415 sub-brick contains the corresponding t-statistic.
1416 */
1417
calculate_differences(anova_options * option_data)1418 void calculate_differences (anova_options * option_data)
1419 {
1420 const float EPSILON = 1.0e-10; /* protect against divide by zero */
1421 float * fin = NULL; /* pointer to float input data */
1422 float * diff = NULL; /* pointer to est. diff. in means */
1423 float * tdiff = NULL; /* pointer to t-statistic data */
1424 float * contr; /* for contrast in new method */
1425 int r; /* number of factor levels */
1426 int nt; /* total number of observations */
1427 int ixyz, nxyz; /* voxel counters */
1428 int nvoxel; /* output voxel # */
1429 int num_diffs; /* number of user requested diffs. */
1430 int idiff; /* index for requested differences */
1431 int i=-1, j=-1; /* factor level indices */
1432 int ni, nj; /* number of obs. at levels i and j */
1433 int df; /* degrees of freedom */
1434 float fval; /* for calculating std. dev. */
1435 float stddev; /* est. std. dev. of difference */
1436 char filename[MAX_NAME_LENGTH]; /* input file name */
1437
1438
1439 /*----- initialize local variables -----*/
1440 r = option_data->a;
1441 nt = option_data->nt;
1442 num_diffs = option_data->num_adiffs;
1443 nxyz = option_data->nxyz;
1444 nvoxel = option_data->nvoxel;
1445
1446 /*----- allocate memory space for calculations -----*/
1447 diff = (float *) malloc(sizeof(float)*nxyz);
1448 tdiff = (float *) malloc(sizeof(float)*nxyz);
1449 contr = (float *) malloc(sizeof(float)*r);
1450 if ((diff == NULL) || (tdiff == NULL) || (contr == NULL))
1451 ANOVA_error ("unable to allocate sufficient memory");
1452
1453 /*----- loop over user specified treatment differences -----*/
1454 for (idiff = 0; idiff < num_diffs; idiff++)
1455 {
1456
1457 if (option_data->old_method) /* old method 8 Dec 2005 [rickr] */
1458 {
1459 df = nt - r;
1460
1461 /*----- allocate temporary memory space -----*/
1462 fin = (float *) malloc(sizeof(float)*nxyz);
1463 if (fin == NULL) ANOVA_error ("unable to allocate sufficient memory");
1464
1465 /*----- read first treatment level mean -----*/
1466 i = option_data->adiffs[idiff][0];
1467 ni = option_data->na[i];
1468 sprintf (filename, "y.%d", i);
1469 volume_read (filename, fin, nxyz);
1470 for (ixyz = 0; ixyz < nxyz; ixyz++)
1471 diff[ixyz] = fin[ixyz] / ni;
1472
1473 /*----- subtract second treatment level mean -----*/
1474 j = option_data->adiffs[idiff][1];
1475 nj = option_data->na[j];
1476 sprintf (filename, "y.%d", j);
1477 volume_read (filename, fin, nxyz);
1478 for (ixyz = 0; ixyz < nxyz; ixyz++)
1479 diff[ixyz] -= fin[ixyz] / nj;
1480
1481 /*----- divide by estimated standard deviation of difference -----*/
1482 volume_read ("sse", fin, nxyz);
1483 fval = (1.0 / (nt-r)) * ((1.0/ni) + (1.0/nj));
1484 for (ixyz = 0; ixyz < nxyz; ixyz++)
1485 {
1486 stddev = sqrt (fin[ixyz] * fval);
1487 if (stddev < EPSILON)
1488 tdiff[ixyz] = 0.0;
1489 else
1490 tdiff[ixyz] = diff[ixyz] / stddev;
1491 }
1492 } else { /* new method */
1493 for (ixyz = 0; ixyz < r; ixyz++ ) contr[ixyz] = 0.0; /* clear old */
1494 contr[option_data->adiffs[idiff][0]] = 1.0;
1495 contr[option_data->adiffs[idiff][1]] = -1.0; /* set difference */
1496 calc_contr(option_data, contr, diff, tdiff);
1497 df = df_for_contr(option_data, contr);
1498 }
1499
1500 if (nvoxel > 0) {
1501 printf ("Diff [%d] - [%d] = %f \n", i+1, j+1, diff[nvoxel-1]);
1502 printf ("t-Diff [%d] - [%d] = %f, df = %d \n",
1503 i+1, j+1, tdiff[nvoxel-1], df);
1504 }
1505
1506 /*----- release memory -----*/
1507 free (fin); fin = NULL;
1508
1509 /*----- write out afni data file -----*/
1510 write_afni_data (option_data, option_data->adname[idiff],
1511 diff, tdiff, df, 0);
1512
1513 }
1514
1515 /*----- release memory -----*/
1516 free (tdiff); free (diff); free (contr);
1517
1518 }
1519
1520
1521 /*---------------------------------------------------------------------------*/
1522 /*
1523 Routine to estimate a user specified contrast in treatment levels.
1524 The output is stored as a 2 sub-brick AFNI data set. The first
1525 sub-brick contains the estimated contrast. The second sub-brick contains
1526 the corresponding t-statistic.
1527 */
1528
calculate_contrasts(anova_options * option_data)1529 void calculate_contrasts (anova_options * option_data)
1530 {
1531 const float EPSILON = 1.0e-10; /* protect against divide by zero */
1532 float * fin = NULL; /* pointer to float input data */
1533 float * contr = NULL; /* pointer to contrast estimate */
1534 float * tcontr = NULL; /* pointer to t-statistic data */
1535 int r; /* number of factor levels */
1536 int nt; /* total number of observations */
1537 int ixyz, nxyz; /* voxel counters */
1538 int nvoxel; /* output voxel # */
1539 int num_contr; /* number of user requested contrasts */
1540 int icontr; /* index of user requested contrast */
1541 int level; /* factor level index */
1542 int ni; /* number of obs. at factor level i */
1543 int df; /* degrees of freedom */
1544 float fval; /* for calculating std. dev. */
1545 float c; /* contrast coefficient */
1546 float stddev; /* est. std. dev. of contrast */
1547 char filename[MAX_NAME_LENGTH]; /* input data file name */
1548
1549
1550 /*----- initialize local variables -----*/
1551 r = option_data->a;
1552 nt = option_data->nt;
1553 num_contr = option_data->num_acontr;
1554 nxyz = option_data->nxyz;
1555 nvoxel = option_data->nvoxel;
1556
1557 /*----- allocate memory space for calculations -----*/
1558 contr = (float *) malloc(sizeof(float)*nxyz);
1559 tcontr = (float *) malloc(sizeof(float)*nxyz);
1560 fin = (float *) malloc(sizeof(float)*nxyz);
1561 if ((contr == NULL) || (tcontr == NULL) || (fin == NULL))
1562 ANOVA_error ("unable to allocate sufficient memory");
1563
1564
1565 /*----- loop over user specified constrasts -----*/
1566 for (icontr = 0; icontr < num_contr; icontr++)
1567 {
1568 if (option_data->old_method) /* old method 08 Dec 2005 [rickr] */
1569 {
1570 df = nt - r;
1571
1572 volume_zero (contr, nxyz);
1573 fval = 0.0;
1574
1575 for (level = 0; level < r; level++)
1576 {
1577 c = option_data->acontr[icontr][level];
1578 if (c == 0.0) continue;
1579
1580 /*----- add c * treatment level mean to contrast -----*/
1581 sprintf (filename, "y.%d", level);
1582 volume_read (filename, fin, nxyz);
1583 ni = option_data->na[level];
1584 fval += c * c / ni;
1585 for (ixyz = 0; ixyz < nxyz; ixyz++)
1586 contr[ixyz] += c * fin[ixyz] / ni;
1587 }
1588
1589 /*----- divide by estimated standard deviation of the contrast -----*/
1590 volume_read ("sse", fin, nxyz);
1591 for (ixyz = 0; ixyz < nxyz; ixyz++)
1592 {
1593 stddev = sqrt ((fin[ixyz] / (nt-r)) * fval);
1594 if (stddev < EPSILON)
1595 tcontr[ixyz] = 0.0;
1596 else
1597 tcontr[ixyz] = contr[ixyz] / stddev;
1598 }
1599 } else { /* new method */
1600 calc_contr(option_data, option_data->acontr[icontr], contr, tcontr);
1601 df = df_for_contr(option_data, option_data->acontr[icontr]);
1602 }
1603
1604 if (nvoxel > 0){
1605 printf ("Contr no.%d = %f \n", icontr+1, contr[nvoxel-1]);
1606 printf ("t-Contr no.%d = %f, df = %d\n",icontr+1,tcontr[nvoxel-1],df);
1607 }
1608
1609 /*----- write out afni data file -----*/
1610 write_afni_data (option_data, option_data->acname[icontr],
1611 contr, tcontr, df, 0);
1612 }
1613
1614 /*----- release memory -----*/
1615 free (tcontr); tcontr = NULL;
1616 free (contr); contr = NULL;
1617 free (fin); fin = NULL;
1618 }
1619
1620
1621 /*---------------------------------------------------------------------------*/
1622 /*
1623 Routine to perform single factor ANOVA.
1624 */
1625
calculate_anova(anova_options * option_data)1626 void calculate_anova (anova_options * option_data)
1627 {
1628
1629 /*----- sum observations within treatment -----*/
1630 calculate_y (option_data);
1631
1632 /*----- sum all observations -----*/
1633 calculate_ysum (option_data);
1634
1635 /*----- calculate sum of squares -----*/
1636 calculate_ss (option_data);
1637
1638 /*----- calculate total (corrected for the mean) sum of squares -----*/
1639 calculate_ssto (option_data);
1640
1641 /*----- calculate treatment sum of squares -----*/
1642 calculate_sstr (option_data);
1643
1644 /*----- calculate error sum of squares -----*/
1645 calculate_sse (option_data);
1646
1647 }
1648
1649
1650 /*---------------------------------------------------------------------------*/
1651 /*
1652 Routine to analyze the results from a single factor ANOVA.
1653 */
1654
analyze_results(anova_options * option_data)1655 void analyze_results (anova_options * option_data)
1656 {
1657
1658 /*----- calculate F-statistic for treatment effect -----*/
1659 if (option_data->nftr > 0) calculate_f_statistic (option_data);
1660
1661 /*----- estimate factor level means -----*/
1662 if (option_data->num_ameans > 0) calculate_means (option_data);
1663
1664 /*----- estimate differences in factor level means -----*/
1665 if (option_data->num_adiffs > 0) calculate_differences (option_data);
1666
1667 /*----- calculate contrasts in factor levels -----*/
1668 if (option_data->num_acontr > 0) calculate_contrasts (option_data);
1669
1670 }
1671
1672
1673 /*---------------------------------------------------------------------------*/
1674 /*
1675 Routine to create an AFNI "bucket" output dataset.
1676 */
1677
create_bucket(anova_options * option_data)1678 void create_bucket (anova_options * option_data)
1679
1680 {
1681 char bucket_str[20000]; /* command line for program 3dbucket */
1682 char refit_str[20000]; /* command line for program 3drefit */
1683 /* (changed from 10K to 20K for Shruti) */
1684 THD_3dim_dataset * dset=NULL; /* input afni data set pointer */
1685 THD_3dim_dataset * new_dset=NULL; /* output afni data set pointer */
1686 int i; /* file index */
1687 int ibrick; /* sub-brick index number */
1688 char str[100]; /* temporary character string */
1689
1690
1691 /*----- read first dataset -----*/
1692 dset = THD_open_dataset (option_data->first_dataset) ;
1693 CHECK_OPEN_ERROR(dset,option_data->first_dataset) ;
1694
1695 if( DSET_IS_1D(dset) ) USE_1D_filenames(1) ; /* 14 Mar 2003 */
1696 else if( DSET_IS_3D(dset) ) USE_1D_filenames(3) ; /* 21 Mar 2003 */
1697
1698
1699 /*----- make an empty copy of this dataset -----*/
1700 new_dset = EDIT_empty_copy( dset ) ;
1701 THD_delete_3dim_dataset (dset , False); dset = NULL;
1702 EDIT_dset_items (new_dset,
1703 ADN_directory_name, option_data->session,
1704 ADN_none);
1705
1706
1707 /*----- begin command line for program 3dbucket -----*/
1708 strcpy (bucket_str, "3dbucket ");
1709 strcat (bucket_str, " -prefix ");
1710 strcat (bucket_str, option_data->bucket_filename);
1711
1712
1713 /*----- begin command line for program 3drefit -----*/
1714 strcpy (refit_str, "3drefit -addFDR ");
1715 ibrick = -1;
1716
1717
1718 /*----- make F-statistic sub-bricks -----*/
1719 if (option_data->nftr != 0)
1720 {
1721 add_file_name (new_dset, option_data->ftrname, bucket_str);
1722
1723 ibrick++;
1724 sprintf (str, " -sublabel %d %s:Inten ",
1725 ibrick, label_from_filename(option_data->ftrname));
1726 strcat (refit_str, str);
1727
1728 ibrick++;
1729 sprintf (str, " -sublabel %d %s:F-stat ",
1730 ibrick, label_from_filename(option_data->ftrname));
1731 strcat (refit_str, str);
1732 }
1733
1734
1735 /*----- make factor level mean sub-bricks -----*/
1736 if (option_data->num_ameans > 0)
1737 for (i = 0; i < option_data->num_ameans; i++)
1738 {
1739 add_file_name (new_dset, option_data->amname[i], bucket_str);
1740
1741 ibrick++;
1742 sprintf (str, " -sublabel %d %s:Mean ",
1743 ibrick, label_from_filename(option_data->amname[i]));
1744 strcat (refit_str, str);
1745
1746 ibrick++;
1747 sprintf (str, " -sublabel %d %s:t-stat ",
1748 ibrick, label_from_filename(option_data->amname[i]));
1749 strcat (refit_str, str);
1750 }
1751
1752
1753 /*----- make difference in factor level means sub-bricks -----*/
1754 if (option_data->num_adiffs > 0)
1755 for (i = 0; i < option_data->num_adiffs; i++)
1756 {
1757 add_file_name (new_dset, option_data->adname[i], bucket_str);
1758
1759 ibrick++;
1760 sprintf (str, " -sublabel %d %s:Diff ",
1761 ibrick, label_from_filename(option_data->adname[i]));
1762 strcat (refit_str, str);
1763
1764 ibrick++;
1765 sprintf (str, " -sublabel %d %s:t-stat ",
1766 ibrick, label_from_filename(option_data->adname[i]));
1767 strcat (refit_str, str);
1768 }
1769
1770
1771 /*----- make contrast in factor level means sub-bricks -----*/
1772 if (option_data->num_acontr > 0)
1773 for (i = 0; i < option_data->num_acontr; i++)
1774 {
1775 add_file_name (new_dset, option_data->acname[i], bucket_str);
1776
1777 ibrick++;
1778 sprintf (str, " -sublabel %d %s:Contr ",
1779 ibrick, label_from_filename(option_data->acname[i]));
1780 strcat (refit_str, str);
1781
1782 ibrick++;
1783 sprintf (str, " -sublabel %d %s:t-stat ",
1784 ibrick, label_from_filename(option_data->acname[i]));
1785 strcat (refit_str, str);
1786 }
1787
1788
1789 /*----- invoke program 3dbucket to generate bucket type output dataset
1790 by concatenating previous output files -----*/
1791 printf("Writing `bucket' dataset ");
1792 printf("into %s\n", option_data->bucket_filename);
1793 fprintf(stderr,"RUNNING COMMAND: %s\n",bucket_str);
1794 system (bucket_str);
1795
1796
1797 /*----- invoke program 3drefit to label individual sub-bricks -----*/
1798 add_file_name (new_dset, option_data->bucket_filename, refit_str);
1799 fprintf(stderr,"RUNNING COMMAND: %s\n",refit_str);
1800 system (refit_str);
1801
1802
1803 /*----- release memory -----*/
1804 THD_delete_3dim_dataset (new_dset , False); new_dset = NULL;
1805
1806 }
1807
1808
1809 /*---------------------------------------------------------------------------*/
1810 /*
1811 Routine to release memory and remove any remaining temporary data files.
1812 If the '-bucket' option has been used, create the bucket dataset and
1813 dispose of the 2 sub-brick datasets.
1814 */
1815
terminate(anova_options * option_data)1816 void terminate (anova_options * option_data)
1817 {
1818 int i;
1819 THD_3dim_dataset * dset=NULL; /* input afni data set pointer */
1820 THD_3dim_dataset * new_dset=NULL; /* output afni data set pointer */
1821 char filename[MAX_NAME_LENGTH];
1822
1823
1824 /*----- remove temporary data files -----*/
1825 volume_delete ("sstr");
1826 volume_delete ("sse");
1827 for (i = 0; i < option_data->a; i++)
1828 {
1829 sprintf (filename, "y.%d", i);
1830 volume_delete (filename);
1831 }
1832
1833
1834 /*----- create bucket dataset -----*/
1835 if (option_data->bucket_filename != NULL)
1836 create_bucket (option_data);
1837
1838
1839 /*----- if 'bucket' datset was created, remove the individual 2-subbrick
1840 data files -----*/
1841 if (option_data->bucket_filename != NULL)
1842 {
1843
1844 /*----- read first dataset -----*/
1845 dset = THD_open_dataset (option_data->first_dataset) ;
1846 CHECK_OPEN_ERROR(dset,option_data->first_dataset) ;
1847
1848 /*----- make an empty copy of this dataset -----*/
1849 new_dset = EDIT_empty_copy (dset);
1850 THD_delete_3dim_dataset (dset , False); dset = NULL;
1851 EDIT_dset_items (new_dset,
1852 ADN_directory_name, option_data->session,
1853 ADN_none);
1854
1855 /*----- remove F-statistic data file -----*/
1856 if (option_data->nftr != 0)
1857 remove_dataset (new_dset, option_data->ftrname);
1858
1859 /*----- remove mean factor level data files -----*/
1860 if (option_data->num_ameans > 0)
1861 for (i = 0; i < option_data->num_ameans; i++)
1862 remove_dataset (new_dset, option_data->amname[i]);
1863
1864 /*----- remove difference in factor levels data files -----*/
1865 if (option_data->num_adiffs > 0)
1866 for (i = 0; i < option_data->num_adiffs; i++)
1867 remove_dataset (new_dset, option_data->adname[i]);
1868
1869 /*----- remove contrast in factor levels data files -----*/
1870 if (option_data->num_acontr > 0)
1871 for (i = 0; i < option_data->num_acontr; i++)
1872 remove_dataset (new_dset, option_data->acname[i]);
1873
1874 THD_delete_3dim_dataset (new_dset , False); new_dset = NULL;
1875 }
1876
1877
1878 /*----- deallocate memory -----*/
1879 destroy_anova_options (option_data);
1880
1881 }
1882
1883
1884 /*---------------------------------------------------------------------------*/
1885 /*
1886 Single factor analysis of variance (ANOVA).
1887 */
1888
main(int argc,char ** argv)1889 int main (int argc, char ** argv)
1890 {
1891 anova_options * option_data = NULL;
1892
1893 #if 0
1894 /*----- Identify software -----*/
1895 printf ("\n\n");
1896 printf ("Program: %s \n", PROGRAM_NAME);
1897 printf ("Author: %s \n", PROGRAM_AUTHOR);
1898 printf ("Initial Release: %s \n", PROGRAM_INITIAL);
1899 printf ("Latest Revision: %s \n", PROGRAM_LATEST);
1900 printf ("\n");
1901 #endif
1902
1903 UNIQ_idcode_fill(SUFFIX+strlen(SUFFIX)) ; /* 27 Apr 2012 */
1904
1905 /*-- 20 Apr 2001: addto the arglist, if user wants to [RWCox] --*/
1906
1907 mainENTRY("3dANOVA main"); machdep(); PRINT_VERSION("3dANOVA"); AUTHOR(PROGRAM_AUTHOR);
1908 { int new_argc ; char ** new_argv ;
1909 addto_args( argc , argv , &new_argc , &new_argv ) ;
1910 if( new_argv != NULL ){ argc = new_argc ; argv = new_argv ; }
1911 }
1912
1913 /*----- program initialization -----*/
1914 initialize (argc, argv, &option_data);
1915
1916 /*----- warn user (after any help) -----*/
1917 if( !option_data->old_method )
1918 fprintf(stderr,"\n"
1919 "** Changes have been made for 3dANOVA computations.\n"
1920 " For details, please see:\n"
1921 " %s\n\n", ANOVA_MODS_LINK);
1922
1923 /*----- calculate sums of squares -----*/
1924 calculate_anova (option_data);
1925
1926 /*----- generate requested output -----*/
1927 analyze_results (option_data);
1928
1929 /*----- terminate program -----*/
1930 terminate (option_data);
1931 free (option_data); option_data = NULL;
1932
1933 exit(0);
1934 }
1935