00001 /*---------------------------------------------------------------------------*/
00002 /*
00003   This file contains routines for performing deconvolution analysis.
00004 
00005   File:    Deconvolve.c
00006   Author:  B. Douglas Ward
00007   Date:    31 August 1998
00008 
00009   Mod:     Restructured matrix calculations to improve execution speed.
00010   Date:    16 December 1998
00011 
00012   Mod:     Report mean square error from full model.
00013   Date:    04 January 1999
00014 
00015   Mod:     Earlier termination if unable to invert X matrix.
00016            (Avoids redundant error messages.)
00017   Date:    06 January 1999
00018 
00019   Mod:     Modifications for matrix calculation of general linear tests.
00020   Date:    02 July 1999
00021 
00022   Mod:     Additional statistical output (partial R^2 statistics).
00023   Date:    07 September 1999
00024 
00025   Mod:     Allow reading of multiple input stimulus functions from a single
00026            file by selection of individual columns.
00027   Date:    09 November 1999
00028 
00029   Mod:     Added options for writing the fitted full model time series (-fitts)
00030            and the residual error time series (-errts) to 3d+time datasets.
00031   Date:    22 November 1999
00032 
00033   Mod:     Added -censor option to allow operator to eliminate individual
00034            time points from the analysis.
00035   Date:    21 June 2000
00036 
00037   Mod:     Added -censor option to allow operator to eliminate individual
00038            time points from the analysis.
00039   Date:    21 June 2000
00040 
00041   Mod:     Added screen display of p-values.
00042   Date:    22 June 2000
00043 
00044   Mod:     Added -glt_label option for labeling the GLT matrix.
00045   Date:    23 June 2000
00046 
00047   Mod:     Added -concat option for analysis of a concatenated 3d+time dataset.
00048   Date:    26 June 2000
00049 
00050   Mod:     Increased size of screen output buffer.
00051   Date:    27 July 2000
00052 
00053   Mod:     Additional output with -nodata option (norm.std.dev.'s for
00054            GLT linear constraints).
00055   Date:    11 August 2000
00056 
00057   Mod:     Added -stim_nptr option, to allow input stim functions to be sampled
00058            at a multiple of the 1/TR rate.
00059   Date:    02 January 2001
00060 
00061   Mod:     Changes to eliminate constraints on number of stimulus time series,
00062            number of regressors, number of GLT's, and number of GLT linear
00063            constraints.
00064   Date:    10 May 2001
00065 
00066   Mod:     Made fstat_t2p() a static function to avoid conflicts on CYGWIN.
00067   Date:    08 Jan 2002 - RWCox
00068 
00069   Mod:     Added static function student_t2p() for display of p-values
00070            corresponding to the t-statistics.
00071   Date:    28 January 2002
00072 
00073   Mod:     Modifications to glt_analysis and report_results for enhanced screen
00074            output:  Display of p-values for individual stim function regression
00075            coefficients.  Display of t-stats and p-values for individual linear
00076            constraints within a GLT.
00077   Date:    29 January 2002
00078 
00079   Mod:     Allow user to specify which input stimulus functions are part of
00080            the baseline model.
00081   Date:    02 May 2002
00082 
00083   Mod:     Increased size of screen output buffer (again).
00084   Date:    02 December 2002
00085 
00086   Mod:     global variable legendre_polort replaces x^m with P_m(x)
00087   Date     15 Jul 2004 - RWCox
00088 
00089 */
00090 
00091 /*---------------------------------------------------------------------------*/
00092 /*
00093   Include linear regression analysis software.
00094 */
00095 
00096 #include "RegAna.c"
00097 
00098 
00099 /*---------------------------------------------------------------------------*/
00100 static int legendre_polort = 1 ;  /* 15 Jul 2004 */
00101 static int demean_base     = 1 ;  /* 12 Aug 2004 */
00102 
00103 #define SPOL ((legendre_polort) ? "P_" : "t^")  /* string label for polynomials */
00104 
00105 double legendre( double x , int m )   /* Legendre polynomials over [-1,1] */
00106 {
00107    if( m < 0 ) return 1.0 ;    /* bad input */
00108 
00109    switch( m ){                /*** P_m(x) for m=0..20 ***/
00110     case 0: return 1.0 ;
00111     case 1: return x ;
00112     case 2: return (3.0*x*x-1.0)/2.0 ;
00113     case 3: return (5.0*x*x-3.0)*x/2.0 ;
00114     case 4: return ((35.0*x*x-30.0)*x*x+3.0)/8.0 ;
00115     case 5: return ((63.0*x*x-70.0)*x*x+15.0)*x/8.0 ;
00116     case 6: return (((231.0*x*x-315.0)*x*x+105.0)*x*x-5.0)/16.0 ;
00117     case 7: return (((429.0*x*x-693.0)*x*x+315.0)*x*x-35.0)*x/16.0 ;
00118     case 8: return ((((6435.0*x*x-12012.0)*x*x+6930.0)*x*x-1260.0)*x*x+35.0)/128.0;
00119 
00120            /** 07 Feb 2005: this part generated by Maple, then hand massaged **/
00121 
00122     case 9:
00123       return (0.24609375e1 + (-0.3609375e2 + (0.140765625e3 + (-0.20109375e3
00124               + 0.949609375e2 * x * x) * x * x) * x * x) * x * x) * x;
00125 
00126     case 10:
00127       return -0.24609375e0 + (0.1353515625e2 + (-0.1173046875e3 +
00128               (0.3519140625e3 + (-0.42732421875e3 + 0.18042578125e3 * x * x)
00129              * x * x) * x * x) * x * x) * x * x;
00130 
00131     case 11:
00132       return (-0.270703125e1 + (0.5865234375e2 + (-0.3519140625e3 +
00133              (0.8546484375e3 + (-0.90212890625e3 + 0.34444921875e3 * x * x)
00134              * x * x) * x * x) * x * x) * x * x) * x;
00135 
00136     case 12:
00137       return 0.2255859375e0 + (-0.17595703125e2 + (0.2199462890625e3 +
00138              (-0.99708984375e3 + (0.20297900390625e4 + (-0.1894470703125e4
00139              + 0.6601943359375e3 * x * x) * x * x) * x * x) * x * x) * x * x)
00140              * x * x;
00141 
00142     case 13:
00143       return (0.29326171875e1 + (-0.87978515625e2 + (0.7478173828125e3 +
00144              (-0.270638671875e4 + (0.47361767578125e4 + (-0.3961166015625e4
00145              + 0.12696044921875e4 * x * x) * x * x) * x * x) * x * x) * x * x)
00146             * x * x) * x;
00147 
00148     case 14:
00149       return -0.20947265625e0 + (0.2199462890625e2 + (-0.37390869140625e3 +
00150              (0.236808837890625e4 + (-0.710426513671875e4 +
00151              (0.1089320654296875e5 + (-0.825242919921875e4 +
00152             0.244852294921875e4 * x * x) * x * x) * x * x) * x * x) * x * x)
00153            * x * x) * x * x;
00154 
00155     case 15:
00156       return (-0.314208984375e1 + (0.12463623046875e3 + (-0.142085302734375e4
00157             + (0.710426513671875e4 + (-0.1815534423828125e5 +
00158               (0.2475728759765625e5 + (-0.1713966064453125e5 +
00159                0.473381103515625e4 * x * x) * x * x) * x * x) * x * x)
00160              * x * x) * x * x) * x * x) * x;
00161 
00162     case 16:
00163       return 0.196380615234375e0 + (-0.26707763671875e2 + (0.5920220947265625e3
00164             + (-0.4972985595703125e4 + (0.2042476226806641e5 +
00165               (-0.4538836059570312e5 + (0.5570389709472656e5 +
00166                (-0.3550358276367188e5 + 0.9171758880615234e4 * x * x) * x * x)
00167             * x * x) * x * x) * x * x) * x * x) * x * x) * x * x;
00168 
00169     case 17:
00170       return (0.3338470458984375e1 + (-0.169149169921875e3 +
00171              (0.2486492797851562e4 + (-0.1633980981445312e5 +
00172              (0.5673545074462891e5 + (-0.1114077941894531e6 +
00173              (0.1242625396728516e6 + (-0.7337407104492188e5 +
00174               0.1780400253295898e5 * x * x) * x * x) * x * x) * x * x)
00175            * x * x) * x * x) * x * x) * x * x) * x;
00176 
00177     case 18:
00178       return -0.1854705810546875e0 + (0.3171546936035156e2 +
00179             (-0.8880331420898438e3 + (0.9531555725097656e4 +
00180             (-0.5106190567016602e5 + (0.153185717010498e6 +
00181             (-0.2692355026245117e6 + (0.275152766418457e6 +
00182             (-0.1513340215301514e6 + 0.3461889381408691e5 * x * x) * x * x)
00183            * x * x) * x * x) * x * x) * x * x) * x * x) * x * x) * x * x;
00184 
00185     case 19:
00186       return (-0.3523941040039062e1 + (0.2220082855224609e3 +
00187              (-0.4084952453613281e4 + (0.3404127044677734e5 +
00188              (-0.153185717010498e6 + (0.4038532539367676e6 +
00189              (-0.6420231216430664e6 + (0.6053360861206055e6 +
00190              (-0.3115700443267822e6 + 0.6741574058532715e5 * x * x) * x * x)
00191           * x * x) * x * x) * x * x) * x * x) * x * x) * x * x) * x * x) * x;
00192 
00193     case 20:
00194       return 0.1761970520019531e0 + (-0.3700138092041016e2 +
00195             (0.127654764175415e4 + (-0.1702063522338867e5 +
00196             (0.1148892877578735e6 + (-0.4442385793304443e6 +
00197             (0.1043287572669983e7 + (-0.1513340215301514e7 +
00198             (0.1324172688388824e7 + (-0.6404495355606079e6 +
00199              0.1314606941413879e6 * x * x) * x * x) * x * x) * x * x) * x * x)
00200             * x * x) * x * x) * x * x) * x * x) * x * x;
00201    }
00202 
00203 #if 0
00204    /* order out of range: return Chebyshev instead (it's easy) */
00205 
00206         if(  x >=  1.0 ) x = 0.0 ;
00207    else if ( x <= -1.0 ) x = 3.14159265358979323846 ;
00208    else                  x = acos(x) ;
00209    return cos(m*x) ;
00210 #else
00211    /** if here, m > 20 ==> use recurrence relation **/
00212 
00213    { double pk, pkm1, pkm2 ; int k ;
00214      pkm2 = legendre( x , 19 ) ;
00215      pkm1 = legendre( x , 20 ) ;
00216      for( k=21 ; k <= m ; k++ , pkm2=pkm1 , pkm1=pk )
00217        pk = ((2.0*k-1.0)*x*pkm1 - (k-1.0)*pkm2)/k ;
00218      return pk ;
00219    }
00220 #endif
00221 }
00222 
00223 /*---------------------------------------------------------------------------*/
00224 
00225 #ifdef USE_BASIS
00226 # define IBOT(ss) ((basis_stim[ss]!=NULL) ? 0                       : min_lag[ss])
00227 # define ITOP(ss) ((basis_stim[ss]!=NULL) ? basis_stim[ss]->nfunc-1 : max_lag[ss])
00228 #else
00229 # define IBOT(ss) min_lag[ss]
00230 # define ITOP(ss) max_lag[ss]
00231 #endif
00232 
00233 /*---------------------------------------------------------------------------*/
00234 /*
00235    Initialize independent variable X matrix
00236 */
00237 
00238 int init_indep_var_matrix
00239 (
00240   int p,                      /* number of parameters in the full model */
00241   int qp,                     /* number of poly. trend baseline parameters */
00242   int polort,                 /* degree of polynomial for baseline model */
00243   int nt,                     /* number of images in input 3d+time dataset */
00244   int N,                      /* total number of images used in the analysis */
00245   int * good_list,            /* list of usable time points */
00246   int * block_list,           /* list of block (run) starting points */
00247   int num_blocks,             /* number of blocks (runs) */
00248   int num_stimts,             /* number of stimulus time series */
00249   float ** stimulus,          /* stimulus time series arrays */
00250   int * stim_length,          /* length of stimulus time series arrays */
00251   int * min_lag,              /* minimum time delay for impulse response */
00252   int * max_lag,              /* maximum time delay for impulse response */
00253   int * nptr,                 /* number of stim fn. time points per TR */
00254   int * stim_base ,           /* flag for being in the baseline [12 Aug 2004] */
00255   matrix * xgood              /* independent variable matrix */
00256 )
00257 
00258 {
00259   int m;                    /* X matrix column index */
00260   int n;                    /* X matrix row index */
00261   int is;                   /* input stimulus index */
00262   int ilag;                 /* time lag index */
00263   int ib;                   /* block (run) index */
00264   int nfirst, nlast;        /* time boundaries of a block (run) */
00265   int mfirst, mlast;        /* column boundaries of baseline parameters
00266                                for a block (run) */
00267 
00268   float * stim_array;       /* stimulus function time series */
00269   matrix x;                 /* X matrix */
00270 
00271   int mold ;                /* 12 Aug 2004 */
00272   int ibot,itop ;
00273 
00274 ENTRY("init_indep_var_matrix") ;
00275 
00276 
00277   /*----- Initialize X matrix -----*/
00278 
00279 STATUS("create x matrix" ) ;
00280   matrix_initialize (&x);
00281   matrix_create (nt, p, &x);
00282 
00283 
00284   /*----- Set up columns of X matrix corresponding to
00285           the baseline (null hypothesis) signal model -----*/
00286 
00287 STATUS("loop over blocks") ;
00288 
00289   for (ib = 0;  ib < num_blocks;  ib++)
00290     {
00291       nfirst = block_list[ib];       /* start time index for this run */
00292       if (ib+1 < num_blocks)
00293              nlast = block_list[ib+1];    /* last+1 time index for this run */
00294       else
00295              nlast = nt;
00296 
00297 if(PRINT_TRACING){
00298   char str[256] ;
00299   sprintf(str," block #%d = %d .. %d",ib,nfirst,nlast-1) ; STATUS(str) ;
00300 }
00301 
00302       for (n = nfirst;  n < nlast;  n++)
00303            {
00304             mfirst =  ib    * (polort+1);   /* first column index */
00305             mlast  = (ib+1) * (polort+1);   /* last+1 column index */
00306 
00307        if( !legendre_polort ){                /* the old way: powers */
00308               for (m = mfirst;  m < mlast;  m++)
00309                 x.elts[n][m] = pow ((double)(n-nfirst), (double)(m-mfirst));
00310 
00311        } else {            /* 15 Jul 2004: the new way: Legendre - RWCox */
00312 
00313          double xx , aa=2.0/(nlast-nfirst-1.0) ; /* map nfirst..nlast-1 */
00314          for( m=mfirst ; m < mlast ; m++ ){      /* to interval [-1,1] */
00315            xx = aa*(n-nfirst) - 1.0 ;
00316            x.elts[n][m] = legendre( xx , m-mfirst ) ;
00317          }
00318        }
00319       }
00320 
00321       if( mfirst+1 < mlast && demean_base ){  /* 12 Aug 2004: remove means? */
00322         float sum ;
00323         for( m=mfirst+1 ; m < mlast ; m++ ){
00324           sum = 0.0f ;
00325           for( n=nfirst ; n < nlast ; n++ ) sum += x.elts[n][m] ;
00326           sum /= (nlast-nfirst) ;
00327           for( n=nfirst ; n < nlast ; n++ ) x.elts[n][m] -= sum ;
00328         }
00329       }
00330     }
00331 
00332 
00333   /*----- Set up columns of X matrix corresponding to
00334           time delayed versions of the input stimulus -----*/
00335 
00336 STATUS("loop over stimulus time series") ;
00337 
00338   m = qp;
00339   for (is = 0;  is < num_stimts;  is++){
00340 #ifdef USE_BASIS
00341     if( basis_vect[is] != NULL ){                 /* 16 Aug 2004 */
00342       float *bv=MRI_FLOAT_PTR(basis_vect[is]) ;
00343       int nf=basis_vect[is]->ny , jj ;
00344 if( PRINT_TRACING ){
00345   char str[256] ;
00346   sprintf(str," stim #%d: expanding into %d basis vectors",is,nf) ;
00347   STATUS(str) ;
00348 }
00349       for( jj=0 ; jj < nf ; jj++ ){
00350         for( n=0 ; n < nt ; n++ ) x.elts[n][m] = bv[n+jj*nt] ;
00351         m++ ;
00352       }
00353     }
00354     else {
00355 #endif
00356       if (stim_length[is] < nt*nptr[is])
00357            {
00358              DC_error ("Input stimulus time series is too short");
00359              RETURN (0);
00360            }
00361       stim_array = stimulus[is]; mold = m ;  /* mold = col index we start at */
00362       ibot = IBOT(is) ; itop = ITOP(is) ;
00363 if( PRINT_TRACING ){
00364   char str[256] ;
00365   sprintf(str," stim #%d: ibot=%d itop=%d",is,ibot,itop) ;
00366   STATUS(str) ;
00367 }
00368       for( ilag=ibot ; ilag <= itop ; ilag++ )
00369            {
00370              for (n = 0;  n < nt;  n++)
00371              {
00372                if (n*nptr[is] < ilag)
00373                         x.elts[n][m] = 0.0;
00374                else
00375                         x.elts[n][m] = stim_array[n*nptr[is]-ilag];
00376              }
00377              m++;
00378            }
00379 
00380       /* 12 Aug 2004: remove mean of baseline columns? */
00381       /* 07 Feb 2005: Oops -- used [m] instead of [mm] in the for(n) loops! */
00382 
00383       if( stim_base != NULL && stim_base[is] && demean_base ){
00384         int mm ; float sum ;
00385 STATUS("  remove baseline mean") ;
00386         for( mm=mold ; mm < m ; mm++ ){
00387           sum = 0.0f ;
00388           for( n=0 ; n < nt ; n++ ) sum += x.elts[n][mm] ;
00389           sum /= nt ;
00390           for( n=0 ; n < nt ; n++ ) x.elts[n][mm] -= sum ;
00391         }
00392       }
00393 #ifdef USE_BASIS
00394     }
00395 #endif
00396   }
00397 
00398 
00399   /*----- Keep only those rows of the X matrix which correspond to
00400           usable time points -----*/
00401 
00402 STATUS("extract xgood matrix") ;
00403 
00404   matrix_extract_rows (x, N, good_list, xgood);
00405   matrix_destroy (&x);
00406 
00407 
00408   RETURN (1);
00409 
00410 }
00411 
00412 
00413 /*---------------------------------------------------------------------------*/
00414 /*
00415   Initialization for the regression analysis.
00416 */
00417 
00418 int init_regression_analysis
00419 (
00420   int p,                      /* number of parameters in full model */
00421   int qp,                     /* number of poly. trend baseline parameters */
00422   int num_stimts,             /* number of stimulus time series */
00423   int * baseline,             /* flag for stim function in baseline model */
00424   int * min_lag,              /* minimum time delay for impulse response */
00425   int * max_lag,              /* maximum time delay for impulse response */
00426   matrix xdata,               /* independent variable matrix */
00427   matrix * x_full,            /* extracted X matrix    for full model */
00428   matrix * xtxinv_full,       /* matrix:  1/(X'X)      for full model */
00429   matrix * xtxinvxt_full,     /* matrix:  (1/(X'X))X'  for full model */
00430   matrix * x_base,            /* extracted X matrix    for baseline model */
00431   matrix * xtxinvxt_base,     /* matrix:  (1/(X'X))X'  for baseline model */
00432   matrix * x_rdcd,            /* extracted X matrices  for reduced models */
00433   matrix * xtxinvxt_rdcd      /* matrix:  (1/(X'X))X'  for reduced models */
00434 )
00435 
00436 {
00437   int * plist = NULL;         /* list of model parameters */
00438   int ip, it;                 /* parameter indices */
00439   int is, js;                 /* stimulus indices */
00440   int im, jm;                 /* lag index */
00441   int ok;                     /* flag for successful matrix calculation */
00442   matrix xtxinv_temp;         /* intermediate results */
00443   int ibot,itop ;
00444 
00445 
00446 ENTRY("init_regression_analysis") ;
00447 
00448   /*----- Initialize matrix -----*/
00449   matrix_initialize (&xtxinv_temp);
00450 
00451 
00452   /*----- Initialize matrices for the baseline model -----*/
00453   plist = (int *) malloc (sizeof(int) * p);   MTEST(plist);
00454   for (ip = 0;  ip < qp;  ip++)
00455     plist[ip] = ip;
00456   it = ip = qp;
00457   for (is = 0;  is < num_stimts;  is++)
00458     {
00459       ibot = IBOT(is) ; itop = ITOP(is) ;
00460       for (im = ibot;  im <= itop;  im++)
00461       {
00462         if (baseline[is])
00463         {
00464           plist[ip] = it;
00465           ip++;
00466         }
00467              it++;
00468            }
00469     }
00470   ok = calc_matrices (xdata, ip, plist, x_base, &xtxinv_temp, xtxinvxt_base);
00471   if (!ok)  { matrix_destroy (&xtxinv_temp);  RETURN (0); };
00472 
00473 
00474   /*----- Initialize matrices for stimulus functions -----*/
00475   for (is = 0;  is < num_stimts;  is++)
00476     {
00477       for (ip = 0;  ip < qp;  ip++)
00478         {
00479           plist[ip] = ip;
00480         }
00481 
00482       it = ip = qp;
00483 
00484       for (js = 0;  js < num_stimts;  js++)
00485         {
00486           ibot = IBOT(js) ; itop = ITOP(js) ;
00487           for (jm = ibot;  jm <= itop;  jm++)
00488             {
00489               if (is != js){ plist[ip] = it; ip++; }
00490               it++;
00491             }
00492         }
00493 
00494       ibot = IBOT(is) ; itop = ITOP(is) ;
00495       ok = calc_matrices (xdata, p-(itop-ibot+1),
00496                      plist, &(x_rdcd[is]), &xtxinv_temp, &(xtxinvxt_rdcd[is]));
00497       if (!ok)  { matrix_destroy (&xtxinv_temp);  RETURN (0); };
00498     }
00499 
00500 
00501   /*----- Initialize matrices for full model -----*/
00502   for (ip = 0;  ip < p;  ip++)
00503     plist[ip] = ip;
00504   ok = calc_matrices (xdata, p, plist, x_full, xtxinv_full, xtxinvxt_full);
00505   if (!ok)  { matrix_destroy (&xtxinv_temp);  RETURN (0); };
00506 
00507 
00508   /*----- Destroy matrix -----*/
00509   matrix_destroy (&xtxinv_temp);
00510 
00511   if (plist != NULL) { free(plist);  plist = NULL; }
00512 
00513   RETURN (1);
00514 }
00515 
00516 
00517 /*---------------------------------------------------------------------------*/
00518 /*
00519   Initialization for the general linear test analysis.
00520 */
00521 
00522 int init_glt_analysis
00523 (
00524   matrix xtxinv,              /* matrix:  1/(X'X)  for full model */
00525   int glt_num,                /* number of general linear tests */
00526   matrix * glt_cmat,          /* general linear test matrices */
00527   matrix * glt_amat,          /* constant GLT matrices for later use */
00528   matrix * cxtxinvct          /* matrices: C(1/(X'X))C' for GLT */
00529 
00530 )
00531 
00532 {
00533   int iglt;                   /* index for general linear test */
00534   int ok;                     /* flag for successful matrix inversion */
00535 
00536 
00537 ENTRY("init_glt_analysis") ;
00538 
00539   for (iglt = 0;  iglt < glt_num;  iglt++)
00540     {
00541       ok = calc_glt_matrix (xtxinv, glt_cmat[iglt], &(glt_amat[iglt]),
00542                             &(cxtxinvct[iglt]));
00543       if (! ok)  RETURN (0);
00544     }
00545 
00546   RETURN (1);
00547 }
00548 
00549 
00550 /*---------------------------------------------------------------------------*/
00551 /*
00552    Calculate results for this voxel.
00553 */
00554 
00555 void regression_analysis
00556 (
00557   int N,                    /* number of usable data points */
00558   int p,                    /* number of parameters in full model */
00559   int q,                    /* number of parameters in baseline model */
00560   int num_stimts,           /* number of stimulus time series */
00561   int * min_lag,            /* minimum time delay for impulse response */
00562   int * max_lag,            /* maximum time delay for impulse response */
00563   matrix x_full,            /* extracted X matrix    for full model */
00564   matrix xtxinv_full,       /* matrix:  1/(X'X)      for full model */
00565   matrix xtxinvxt_full,     /* matrix:  (1/(X'X))X'  for full model */
00566   matrix x_base,            /* extracted X matrix    for baseline model */
00567   matrix xtxinvxt_base,     /* matrix:  (1/(X'X))X'  for baseline model */
00568   matrix * x_rdcd,          /* extracted X matrices  for reduced models */
00569   matrix * xtxinvxt_rdcd,   /* matrix:  (1/(X'X))X'  for reduced models */
00570   vector y,                 /* vector of measured data */
00571   float rms_min,            /* minimum variation in data to fit full model */
00572   float * mse,              /* mean square error from full model */
00573   vector * coef_full,       /* regression parameters */
00574   vector * scoef_full,      /* std. devs. for regression parameters */
00575   vector * tcoef_full,      /* t-statistics for regression parameters */
00576   float * fpart,            /* partial F-statistics for the stimuli */
00577   float * rpart,            /* partial R^2 stats. for the stimuli */
00578   float * ffull,            /* full model F-statistics */
00579   float * rfull,            /* full model R^2 stats. */
00580   int * novar,              /* flag for insufficient variation in data */
00581   float * fitts,            /* full model fitted time series */
00582   float * errts             /* full model residual error time series */
00583 )
00584 
00585 {
00586   int is;                   /* input stimulus index */
00587   float sse_base;           /* error sum of squares, baseline model */
00588   float sse_rdcd;           /* error sum of squares, reduced model */
00589   float sse_full;           /* error sum of squares, full model */
00590   vector coef_temp;         /* intermediate results */
00591   int ibot,itop ;
00592 
00593 
00594 ENTRY("regression_analysis") ;
00595 
00596   /*----- Initialization -----*/
00597   vector_initialize (&coef_temp);
00598 
00599 
00600   /*----- Calculate regression coefficients for baseline model -----*/
00601   calc_coef (xtxinvxt_base, y, &coef_temp);
00602 
00603 
00604   /*----- Calculate the error sum of squares for the baseline model -----*/
00605   sse_base = calc_sse (x_base, coef_temp, y);
00606 
00607 
00608   /*----- Stop here if variation about baseline is sufficiently low -----*/
00609   if (sqrt(sse_base/N) < rms_min)
00610     {
00611       int it;
00612 
00613       *novar = 1;
00614       vector_create (p, coef_full);
00615       vector_create (p, scoef_full);
00616       vector_create (p, tcoef_full);
00617       for (is = 0;  is < num_stimts;  is++)
00618         {
00619           fpart[is] = 0.0;
00620           rpart[is] = 0.0;
00621         }
00622       for (it = 0;  it < N;  it++)
00623         {
00624           fitts[it] = 0.0;
00625           errts[it] = 0.0;
00626         }
00627       *mse = 0.0;
00628       *rfull = 0.0;
00629       *ffull = 0.0;
00630       vector_destroy (&coef_temp);
00631       EXRETURN;
00632     }
00633   else
00634     *novar = 0;
00635 
00636 
00637   /*----- Calculate regression coefficients for the full model  -----*/
00638   calc_coef (xtxinvxt_full, y, coef_full);
00639 
00640 
00641   /*----- Calculate the error sum of squares for the full model -----*/
00642   sse_full = calc_sse_fit (x_full, *coef_full, y, fitts, errts);
00643   *mse = sse_full / (N-p);
00644 
00645 
00646   /*----- Calculate t-statistics for the regression coefficients -----*/
00647   calc_tcoef (N, p, sse_full, xtxinv_full,
00648               *coef_full, scoef_full, tcoef_full);
00649 
00650 
00651   /*----- Determine significance of the individual stimuli -----*/
00652   for (is = 0;  is < num_stimts;  is++)
00653     {
00654 
00655       /*----- Calculate regression coefficients for reduced model -----*/
00656       calc_coef (xtxinvxt_rdcd[is], y, &coef_temp);
00657 
00658 
00659       /*----- Calculate the error sum of squares for the reduced model -----*/
00660       sse_rdcd = calc_sse (x_rdcd[is], coef_temp, y);
00661 
00662 
00663       /*----- Calculate partial F-stat for significance of the stimulus -----*/
00664       ibot = IBOT(is) ; itop = ITOP(is) ;
00665       fpart[is] = calc_freg (N, p, p-(itop-ibot+1), sse_full, sse_rdcd);
00666 
00667 
00668       /*----- Calculate partial R^2 for this stimulus -----*/
00669       rpart[is] = calc_rsqr (sse_full, sse_rdcd);
00670 
00671     }
00672 
00673 
00674   /*----- Calculate coefficient of multiple determination R^2 -----*/
00675   *rfull = calc_rsqr (sse_full, sse_base);
00676 
00677 
00678   /*----- Calculate the total regression F-statistic -----*/
00679   *ffull = calc_freg (N, p, q, sse_full, sse_base);
00680 
00681 
00682   /*----- Dispose of vector -----*/
00683   vector_destroy (&coef_temp);
00684 
00685   EXRETURN ;
00686 }
00687 
00688 
00689 /*---------------------------------------------------------------------------*/
00690 /*
00691   Perform the general linear test analysis for this voxel.
00692 */
00693 
00694 void glt_analysis
00695 (
00696   int N,                      /* number of data points */
00697   int p,                      /* number of parameters in the full model */
00698   matrix x,                   /* X matrix for full model */
00699   vector y,                   /* vector of measured data */
00700   float ssef,                 /* error sum of squares from full model */
00701   vector coef,                /* regression parameters for full model */
00702   int novar,                  /* flag for insufficient variation in data */
00703   matrix * cxtxinvct,         /* matrices: C(1/(X'X))C' for GLT */
00704   int glt_num,                /* number of general linear tests */
00705   int * glt_rows,             /* number of linear constraints in glt */
00706   matrix * glt_cmat,          /* general linear test matrices */
00707   matrix * glt_amat,          /* constant matrices */
00708   vector * glt_coef,          /* linear combinations from GLT matrices */
00709   vector * glt_tcoef,         /* t-statistics for GLT linear combinations */
00710   float * fglt,               /* F-statistics for the general linear tests */
00711   float * rglt                /* R^2 statistics for the general linear tests */
00712 )
00713 
00714 {
00715   int iglt;                   /* index for general linear test */
00716   int q;                      /* number of parameters in the rdcd model */
00717   float sser;                 /* error sum of squares, reduced model */
00718   vector rcoef;               /* regression parameters for reduced model */
00719   vector scoef;               /* std. devs. for regression parameters  */
00720 
00721 ENTRY("glt_analysis") ;
00722 
00723   /*----- Initialization -----*/
00724   vector_initialize (&rcoef);
00725   vector_initialize (&scoef);
00726 
00727 
00728   /*----- Loop over multiple general linear tests -----*/
00729   for (iglt = 0;  iglt < glt_num;  iglt++)
00730     {
00731       /*----- Test for insufficient variation in data -----*/
00732       if (novar)
00733         {
00734           vector_create (glt_rows[iglt], &glt_coef[iglt]);
00735           vector_create (glt_rows[iglt], &glt_tcoef[iglt]);
00736           fglt[iglt] = 0.0;
00737           rglt[iglt] = 0.0;
00738         }
00739       else
00740         {
00741           /*----- Calculate the GLT linear combinations -----*/
00742           calc_lcoef (glt_cmat[iglt], coef, &glt_coef[iglt]);
00743         
00744           /*----- Calculate t-statistics for GLT linear combinations -----*/
00745           calc_tcoef (N, p, ssef, cxtxinvct[iglt],
00746                       glt_coef[iglt], &scoef, &glt_tcoef[iglt]);
00747 
00748           /*----- Calculate regression parameters for the reduced model -----*/
00749           /*-----   (that is, the model in the column space of X but )  -----*/
00750           /*-----   (orthogonal to the restricted column space of XC')  -----*/
00751           calc_rcoef (glt_amat[iglt], coef, &rcoef);
00752 
00753           /*----- Calculate error sum of squares for the reduced model -----*/
00754           sser = calc_sse (x, rcoef, y);
00755 
00756           /*----- Calculate the F-statistic for this GLT -----*/
00757           q = p - glt_rows[iglt];
00758           fglt[iglt] = calc_freg (N, p, q, ssef, sser);
00759 
00760           /*----- Calculate the R^2 statistic for this GLT -----*/
00761           rglt[iglt] = calc_rsqr (ssef, sser);
00762 
00763         }
00764     }
00765 
00766 
00767   /*----- Dispose of vectors -----*/
00768   vector_destroy (&rcoef);
00769   vector_destroy (&scoef);
00770 
00771   EXRETURN ;
00772 }
00773 
00774 
00775 /*---------------------------------------------------------------------------*/
00776 /*
00777   Convert t-values and F-values to p-value.
00778   These routines were copied and modified from: mri_stats.c
00779   16 May 2005: Change names from '_t2p' to '_t2pp' to avoid library
00780                name conflict on Mac OS X 10.4 (stupid system).
00781 */
00782 
00783 
00784 static double student_t2pp( double tt , double dof )
00785 {
00786    double bb , xx , pp ;
00787 
00788    tt = fabs(tt);
00789 
00790    if( dof < 1.0 ) return 1.0 ;
00791 
00792    if (tt >= 1000.0)  return 0.0;
00793 
00794    bb = lnbeta( 0.5*dof , 0.5 ) ;
00795    xx = dof/(dof + tt*tt) ;
00796    pp = incbeta( xx , 0.5*dof , 0.5 , bb ) ;
00797    return pp ;
00798 }
00799 
00800 
00801 static double fstat_t2pp( double ff , double dofnum , double dofden )
00802 {
00803    int which , status ;
00804    double p , q , f , dfn , dfd , bound ;
00805 
00806    if (ff >= 1000.0)  return 0.0;
00807 
00808    which  = 1 ;
00809    p      = 0.0 ;
00810    q      = 0.0 ;
00811    f      = ff ;
00812    dfn    = dofnum ;
00813    dfd    = dofden ;
00814 
00815    cdff( &which , &p , &q , &f , &dfn , &dfd , &status , &bound ) ;
00816 
00817    if( status == 0 ) return q ;
00818    else              return 1.0 ;
00819 }
00820 
00821 /*---------------------------------------------------------------------------*/
00822 
00823 static char lbuf[65536];  /* character string containing statistical summary */
00824 static char sbuf[512];
00825 
00826 
00827 void report_results
00828 (
00829   int N,                      /* number of usable data points */
00830   int qp,                     /* number of poly. trend baseline parameters */
00831   int q,                      /* number of baseline model parameters */
00832   int p,                      /* number of full model parameters */
00833   int polort,                 /* degree of polynomial for baseline model */
00834   int * block_list,           /* list of block (run) starting points */
00835   int num_blocks,             /* number of blocks (runs) */
00836   int num_stimts,             /* number of stimulus time series */
00837   char ** stim_label,         /* label for each stimulus */
00838   int * baseline,             /* flag for stim function in baseline model */
00839   int * min_lag,              /* minimum time delay for impulse response */
00840   int * max_lag,              /* maximum time delay for impulse response */
00841   vector coef,                /* regression parameters */
00842   vector tcoef,               /* t-statistics for regression parameters */
00843   float * fpart,              /* partial F-statistics for the stimuli */
00844   float * rpart,              /* partial R^2 stats. for the stimuli */
00845   float ffull,                /* full model F-statistic */
00846   float rfull,                /* full model R^2 stat. */
00847   float mse,                  /* mean square error from full model */
00848   int glt_num,                /* number of general linear tests */
00849   char ** glt_label,          /* label for general linear tests */
00850   int * glt_rows,             /* number of linear constraints in glt */
00851   vector *  glt_coef,         /* linear combinations from GLT matrices */
00852   vector *  glt_tcoef,        /* t-statistics for GLT linear combinations */
00853   float * fglt,               /* F-statistics for the general linear tests */
00854   float * rglt,               /* R^2 statistics for the general linear tests */
00855   char ** label               /* statistical summary for ouput display */
00856 )
00857 
00858 {
00859   const int MAXBUF = 65000;   /* maximum buffer string length */
00860   int m;                      /* coefficient index */
00861   int is;                     /* stimulus index */
00862   int ilag;                   /* time lag index */
00863 
00864   int iglt;                   /* general linear test index */
00865   int ilc;                    /* linear combination index */
00866 
00867   double pvalue;              /* p-value corresponding to F-value */
00868   int r;                      /* number of parameters in the reduced model */
00869 
00870   int ib;                   /* block (run) index */
00871   int mfirst, mlast;        /* column boundaries of baseline parameters
00872                                for a block (run) */
00873   int ibot,itop ;
00874 
00875 
00876   lbuf[0] = '\0' ;   /* make this a 0 length string to start */
00877 
00878   /** for each reference, make a string into sbuf **/
00879 
00880 
00881   /*----- Statistical results for baseline fit -----*/
00882   if (num_blocks == 1)
00883     {
00884       sprintf (sbuf, "\nBaseline: \n");
00885       if (strlen(lbuf) < MAXBUF)  strcat(lbuf,sbuf); else goto finisher ;
00886       for (m=0;  m < qp;  m++)
00887         {
00888           sprintf (sbuf, "%s%d   coef = %10.4f    ", SPOL,m, coef.elts[m]);
00889           if (strlen(lbuf) < MAXBUF)  strcat(lbuf,sbuf) ; else goto finisher ;
00890           sprintf (sbuf, "%s%d   t-st = %10.4f    ", SPOL,m, tcoef.elts[m]);
00891           if (strlen(lbuf) < MAXBUF)  strcat(lbuf,sbuf) ; else goto finisher ;
00892           pvalue = student_t2pp ((double)tcoef.elts[m], (double)(N-p));
00893           sprintf (sbuf, "p-value  = %12.4e \n", pvalue);
00894           if (strlen(lbuf) < MAXBUF)  strcat (lbuf, sbuf); else goto finisher ;
00895         }
00896     }
00897   else
00898     {
00899       for (ib = 0;  ib < num_blocks;  ib++)
00900         {
00901           sprintf (sbuf, "\nBaseline for Run #%d: \n", ib+1);
00902           if (strlen(lbuf) < MAXBUF)  strcat(lbuf,sbuf); else goto finisher ;
00903         
00904           mfirst =  ib    * (polort+1);
00905           mlast  = (ib+1) * (polort+1);
00906           for (m = mfirst;  m < mlast;  m++)
00907             {
00908               sprintf (sbuf, "%s%d   coef = %10.4f    ",
00909                        SPOL,m - mfirst, coef.elts[m]);
00910               if (strlen(lbuf) < MAXBUF)  strcat(lbuf,sbuf) ; else goto finisher ;
00911               sprintf (sbuf, "%s%d   t-st = %10.4f    ",
00912                        SPOL,m - mfirst, tcoef.elts[m]);
00913               if (strlen(lbuf) < MAXBUF)  strcat(lbuf,sbuf) ; else goto finisher ;
00914               pvalue = student_t2pp ((double)tcoef.elts[m], (double)(N-p));
00915               sprintf (sbuf, "p-value  = %12.4e \n", pvalue);
00916               if (strlen(lbuf) < MAXBUF)  strcat (lbuf, sbuf); else goto finisher ;
00917             }
00918         }
00919     }
00920 
00921 
00922   /*----- Statistical results for stimulus response -----*/
00923   m = qp;
00924   for (is = 0;  is < num_stimts;  is++)
00925     {
00926       if (baseline[is])
00927         sprintf (sbuf, "\nBaseline: %s \n", stim_label[is]);    
00928       else
00929         sprintf (sbuf, "\nStimulus: %s \n", stim_label[is]);
00930       if (strlen(lbuf) < MAXBUF)  strcat(lbuf,sbuf); else goto finisher ;
00931       ibot = IBOT(is) ; itop = ITOP(is) ;
00932       for (ilag = ibot;  ilag <= itop;  ilag++)
00933         {
00934           sprintf (sbuf,"h[%2d] coef = %10.4f    ", ilag, coef.elts[m]);
00935           if (strlen(lbuf) < MAXBUF)  strcat(lbuf,sbuf) ; else goto finisher ;
00936           sprintf  (sbuf,"h[%2d] t-st = %10.4f    ", ilag, tcoef.elts[m]);
00937           if (strlen(lbuf) < MAXBUF)  strcat (lbuf, sbuf); else goto finisher ;
00938           pvalue = student_t2pp ((double)tcoef.elts[m], (double)(N-p));
00939           sprintf (sbuf, "p-value  = %12.4e \n", pvalue);
00940           if (strlen(lbuf) < MAXBUF)  strcat (lbuf, sbuf); else goto finisher ;
00941           m++;
00942         }
00943 
00944       sprintf (sbuf, "       R^2 = %10.4f    ", rpart[is]);
00945       if (strlen(lbuf) < MAXBUF)  strcat (lbuf, sbuf); else goto finisher ;
00946       r = p - (itop-ibot+1);
00947       sprintf (sbuf, "F[%2d,%3d]  = %10.4f    ", p-r, N-p, fpart[is]);
00948       if (strlen(lbuf) < MAXBUF)  strcat (lbuf, sbuf); else goto finisher ;
00949       pvalue = fstat_t2pp ((double)fpart[is], (double)(p-r), (double)(N-p));
00950       sprintf (sbuf, "p-value  = %12.4e \n", pvalue);
00951       if (strlen(lbuf) < MAXBUF)  strcat (lbuf, sbuf); else goto finisher ;
00952     }
00953 
00954 
00955   /*----- Statistical results for full model -----*/
00956   sprintf (sbuf, "\nFull Model: \n");
00957   if (strlen(lbuf) < MAXBUF)  strcat (lbuf, sbuf); else goto finisher ;
00958 
00959   sprintf (sbuf, "       MSE = %10.4f \n", mse);
00960   if (strlen(lbuf) < MAXBUF)  strcat (lbuf, sbuf); else goto finisher ;
00961 
00962   sprintf (sbuf, "       R^2 = %10.4f    ", rfull);
00963   if (strlen(lbuf) < MAXBUF)  strcat (lbuf, sbuf); else goto finisher ;
00964 
00965   sprintf (sbuf, "F[%2d,%3d]  = %10.4f    ", p-q, N-p, ffull);
00966   if (strlen(lbuf) < MAXBUF)  strcat (lbuf, sbuf); else goto finisher ;
00967   pvalue = fstat_t2pp ((double)ffull, (double)(p-q), (double)(N-p));
00968   sprintf (sbuf, "p-value  = %12.4e  \n", pvalue);
00969   if (strlen(lbuf) < MAXBUF)  strcat (lbuf, sbuf); else goto finisher ;
00970 
00971 
00972   /*----- Statistical results for general linear test -----*/
00973   if (glt_num > 0)
00974     {
00975       for (iglt = 0;  iglt < glt_num;  iglt++)
00976         {
00977           sprintf (sbuf, "\nGeneral Linear Test: %s \n", glt_label[iglt]);
00978           if (strlen(lbuf) < MAXBUF)  strcat (lbuf,sbuf); else goto finisher ;
00979           for (ilc = 0;  ilc < glt_rows[iglt];  ilc++)
00980             {
00981               sprintf (sbuf, "LC[%d] coef = %10.4f    ",
00982                        ilc, glt_coef[iglt].elts[ilc]);
00983               if (strlen(lbuf) < MAXBUF)  strcat (lbuf,sbuf); else goto finisher ;
00984               sprintf (sbuf, "LC[%d] t-st = %10.4f    ",
00985                        ilc, glt_tcoef[iglt].elts[ilc]);
00986               if (strlen(lbuf) < MAXBUF)  strcat (lbuf,sbuf); else goto finisher ;
00987               pvalue = student_t2pp ((double)glt_tcoef[iglt].elts[ilc],
00988                                     (double)(N-p));
00989               sprintf (sbuf, "p-value  = %12.4e \n", pvalue);
00990               if (strlen(lbuf) < MAXBUF)  strcat (lbuf, sbuf); else goto finisher ;
00991             }
00992 
00993           sprintf (sbuf, "       R^2 = %10.4f    ", rglt[iglt]);
00994           if (strlen(lbuf) < MAXBUF)  strcat (lbuf,sbuf); else goto finisher ;
00995 
00996           r = p - glt_rows[iglt];
00997           sprintf (sbuf, "F[%2d,%3d]  = %10.4f    ", p-r, N-p, fglt[iglt]);
00998           if (strlen(lbuf) < MAXBUF)  strcat (lbuf, sbuf); else goto finisher ;
00999           pvalue = fstat_t2pp ((double)fglt[iglt],
01000                               (double)(p-r), (double)(N-p));
01001           sprintf (sbuf, "p-value  = %12.4e  \n", pvalue);
01002           if (strlen(lbuf) < MAXBUF)  strcat (lbuf, sbuf); else goto finisher ;
01003         }
01004     }
01005 
01006 finisher:
01007   if (strlen(lbuf) >= MAXBUF)
01008     strcat (lbuf, "\n\nWarning:  Screen output buffer is full. \n");
01009 
01010   *label = lbuf ;  /* send address of lbuf back in what label points to */
01011 
01012 }