*********************************************************************; * Version 3.2 January, 1998 *; *********************************************************************; START _TSTGLH_(C,U,THETA0,THETARNM,THETACNM, _EPSL_, _ECODE_, _OPT_,_PARM_, _SS_, _TGLHP1_, _VNAME_,_VTYPE_, _THETA_,_MSH_,_MSE_,_MID_, _SDTHTA_,_TTHTA_,_PVTHTA_,_FSTATS_, _HEIVAL_,_CANVEC_,_CANRSQ_,_STMAT1_,_TPARM1_, _NONM_,_FSTRNM_,_STMCNM_,_STMRNM_,_TPCNM1_, _THRNM_,_THCNM_,_CANNM_, _URESUL_,_UCOLNM_,_UROWNM_ ); *___________________________________________________________________* This module estimates secondary parameters of the form: _THETA_=C*_B_*U-THETA0, and computes the associated standard errors, test statistics, p-values, etc. Inputs C,U,THETA0,THETARNM,THETACNM, _EPSL_,_ECODE_,_OPT_,_PARM_,_SS_,_TGLHP1_,_VNAME_,_VTYPE_, Outputs _THETA_,_MSH_,_MSE_,_MID_, _SDTHTA_,_TTHTA_,_PVTHTA_,_FSTATS_, _HEIVAL_,_CANVEC_,_CANRSQ_,_STMAT1_,_TPARM1_, _NONM_,_FSTRNM_,_STMCNM_,_STMRNM_,_TPCNM1_, _THRNM_,_THCNM_, _URESUL_,_UCOLNM_,_UROWNM_ *___________________________________________________________________*; FREE _THETA_ _MSH_ _MSE_ _MID_ _SDTHTA_ _TTHTA_ _PVTHTA_ _FSTATS_ _HEIVAL_ _CANVEC_ _CANRSQ_ _STMAT1_ _TPARM1_ _NONM_ _FSTRNM_ _STMCNM_ _STMRNM_ _TPCNM1_ _THRNM_ _THCNM_ _CANNM_ _URESUL_ _UCOLNM_ _UROWNM_ ; * Some modest input validity checks; IF (NROW(_ECODE_)=0)|(NROW(_OPT_ )=0)| (NROW(_EPSL_ )=0)|(NROW(_TGLHP1_)=0) THEN DO; IF _OPT_(|1,6|)=0 THEN PRINT "One or more system matrix does not exist. " , "_ECODE_, _OPT_, _EPSL_, _TGLHP1_ were checked. " ,, "-------------TESTGLH, LINMOD, and IML stopped.------------" , "-------------TESTGLH, LINMOD, and IML stopped.------------" , "-------------TESTGLH, LINMOD, and IML stopped.------------" ; ABORT; END; IF (NROW(_PARM_ )=0)|(NROW(_SS_ )=0)| (NROW(_VNAME_ )=0)|(NROW(_VTYPE_)=0) THEN DO; IF _OPT_(|1,6|)=0 THEN PRINT "One or more system matrix does not exist. " , "_PARM_, _SS_, _VNAME_, _VTYPE_ were checked. " ,, "-------------TESTGLH, LINMOD, and IML stopped.------------" , "-------------TESTGLH, LINMOD, and IML stopped.------------" , "-------------TESTGLH, LINMOD, and IML stopped.------------" ; ABORT; END; _ECODE_ (| 5 , 1 |) = 0; * If first invocation of TESTGLH then print header ; *; IF (_OPT_(| 1 , 5 |)=0)&(_OPT_(|1,6|)=0) THEN DO; NOTET_0 = { "**************************************" , "** BEGIN TESTGLH **" , "**************************************" } ; PRINT , NOTET_0 , ; END; IF _TGLHP1_ & MAX(_OPT_(| 5, 12:16 |)) & ( _OPT_(| 1,5 |)=0) & (_OPT_(|1,6|)=0) THEN DO; PRINT "*****************************************************", "** The following multivariate test statistics are **", "** provided in TESTGLH, along with approximate **", "** P Values: **", "** LROOT: Roy's Largest Root Criterion **", "** LROOT=ROOT1/(1+ROOT1), where **", "** ROOT1=MAX(EVAL(_H_*INV(_E_))) **", "** LAMBDA: Wilks' LAMBDA, The Likelihood **", "** Ratio Criterion **", "** LAMBDA=|_E_|/|_H_+_E_| **", "** HLTRACE: The Hotelling-Lawley Trace **", "** HLTRACE=TRACE(_H_*INV(_E_)) **", "** PILTRACE: Pillai's Trace **", "** PILTRACE=TRACE(_H_*INV(_H_+_E_)) **", "** See the Linmod Language Reference Manual, **", "** Appendix A for definitions of the associated **", "** parameters, and related approximate tests and **", "** P Values. **", "*****************************************************" ; IF _OPT_(| 5 , 17 |) THEN PRINT "*****************************************************", "** **", "** The 'ASSOCITN' column below contains measures **", "** of multivariate association. For the four test **", "** statistics defined above 'ASSOCITN' is equal to:**", "** (1) LROOT **", "** (2) 1 - (LAMBDA**(1/S)) **", "** (3) (HLTRACE/S)/(1 + HLTRACE/S) **", "** (4) PILTRACE/S **", "** For a discussion of these measures of **", "** association see Muller and Peterson, Paper **", "** presented at SUGI, 1983. **", "*****************************************************" ; END; _TGLHP1_ = 0; * Check for errors from FITMODEL *; IF ( _ECODE_(| 5 , 1 |) > 2000 ) THEN DO; NOTET_E1={"TESTGLH: This secondary parameter not estimated due to" " : previous errors in FITMODEL noted above "}; IF _OPT_(|1,6|)=0 THEN PRINT , NOTET_E1, ; _ECODE_ (| 5 , 1 |) = 3001; RETURN; END; * Check to see if input matrices defined *; IF NROW(C) = 0 THEN DO; NOTET_E2={"TESTGLH: C matrix not defined." , " Module skipped. " } ; IF _OPT_(|1,6|)=0 THEN PRINT , NOTET_E2, ; _ECODE_ (| 5 , 1 |) = 3002; RETURN; END; * Create dummy matrix to wipe out unwanted row and column labels ; _NONM_ = { "" }; * Create _XINDEX_ & _YINDEX_ *; _LOC_ = ( _VTYPE_ = 1 ); _LOC_ = LOC( _LOC_ ); IF NCOL( _LOC_ ) = 0 THEN DO; NOTET_E3={"TESTGLH: No model has been fitted." , " : Module skipped. "}; IF _OPT_(|1,6|)=0 THEN PRINT , NOTET_E3, ; _ECODE_ (| 5 , 1 |) = 3007; RETURN; END; ELSE _XINDEX_ = _LOC_; _LOC_ = ( _VTYPE_ = -1 ); _LOC_ = LOC( _LOC_ ); IF NCOL( _LOC_ ) = 0 THEN DO; NOTET_E4={"TESTGLH: No dependent variables have been specified" , " : secondary parameter cannot be estimated " }; IF _OPT_(|1,6|)=0 THEN PRINT , NOTET_E4, ; _ECODE_ (| 5 , 1 |) = 3008; RETURN; END; ELSE _YINDEX_ = _LOC_; * Check for conformability of C U & THETA0 with BETA **; _NRB1_ = NCOL( _XINDEX_ ); _NCB1_ = NCOL( _YINDEX_ ); IF _NRB1_ ^= NCOL( C ) THEN DO; IF _OPT_(|1,6|)=0 THEN PRINT , "TESTGLH : C and BETA do not conform" , " : NCOL( C ) must = NROW( BETA )" ,, "Current values are :" ; _T1_ = NCOL( C ) || _NRB1_ ; _NM1_ = { "NCOL( C )" "NROW( B )" }; IF _OPT_(|1,6|)=0 THEN PRINT , _T1_ (| ROWNAME = _NONM_ COLNAME = _NM1_ |); _ECODE_ (| 5 , 1 |) = 3003; END; _NROWC_=NROW(C); IF NCOL(U)=0 THEN _NCOLU_=_PARM_ (| 1 , 2 |); ELSE _NCOLU_=NCOL(U); IF (NCOL(U)>0) & (NROW(U)^=_NCB1_) THEN DO; IF _OPT_(|1,6|)=0 THEN PRINT , "TESTGLH : U and BETA do not conform" , " : NROW( U ) must = NCOL( BETA )" ,, "Current values are :"; _NM1_ = { "NCOL( B )" "NROW( U )" }; _T1_ = _NCB1_ || NROW( U ); IF _OPT_(|1,6|)=0 THEN PRINT , _T1_ (| ROWNAME = _NONM_ COLNAME = _NM1_ |); _ECODE_ (| 5 , 1 |) = 3004; END; IF NROW( THETA0 ) > 0 THEN DO; IF ( _NCOLU_ ^= NCOL(THETA0)) | ( _NROWC_ ^= NROW(THETA0)) THEN DO; IF _OPT_(|1,6|)=0 THEN PRINT , "TESTGLH: THETA0 and THETA do not conform" , " : They must have identical dimensions" ,, "CURRENT VALUES ARE: "; _T1_ = ( _NROWC_ || _NCOLU_ ) // ( NROW(THETA0) || NCOL(THETA0) ); _NM1_ = { "# Rows" "# Cols" }; _NM2_ = { "THETA" "THETA0" }; IF _OPT_(|1,6|)=0 THEN PRINT , _T1_ (| ROWNAME = _NM2_ COLNAME = _NM1_ |); _ECODE_ (| 5 , 1 |) = 3005; END; END; * End of THETA0 check *; IF ( NCOL( THETACNM ) > 0 ) THEN DO; _NC1_ = NCOL( THETACNM ); IF ( _NC1_ ^= _NCOLU_ ) THEN DO; IF _OPT_(|1,6|)=0 THEN PRINT , "TESTGLH: THETACNM and THETA do not conform" , "Current values are: "; _NM1_ = { "NCOL(THETACNM)" "NCOL(THETA)" }; _T1_ = _NC1_ || _NCOLU_; IF _OPT_(|1,6|)=0 THEN PRINT , _T1_ (| ROWNAME = _NONM_ COLNAME = _NM1_ |) ; _ECODE_ (| 5 , 1 |) = 3009; END; END; * End of THETACNM check *; IF ( NCOL( THETARNM ) > 0 ) THEN DO; _NR1_ = NCOL( THETARNM ); IF ( _NR1_ ^= _NROWC_ ) THEN DO; IF _OPT_(|1,6|)=0 THEN PRINT , "TESTGLH: THETARNM and THETA do not conform" "Current values are: "; _NM1_ = { "NCOL(THETARNM)" "NROW(THETA)" }; _T1_ = _NR1_ || _NROWC_; IF _OPT_(|1,6|)=0 THEN PRINT , _T1_ (| ROWNAME = _NONM_ COLNAME = _NM1_ |) ; _ECODE_ (| 5 , 1 |) = 3010; END; END; * End of THETARNM check *; * If contrast and name matrices ok, then estimate THETA; IF _ECODE_(| 5 , 1|) > 3000 THEN RETURN; * Create name vectors for labeling matrices *; _XNM_ = _VNAME_(| 1 , _XINDEX_ |); _YNM_ = _VNAME_(| 1 , _YINDEX_ |); IF NROW(THETACNM)=0 THEN _THCNM_=_NMLIST_("Col_U_",_NCOLU_); ELSE _THCNM_=THETACNM; IF NROW(THETARNM)=0 THEN _THRNM_=_NMLIST_("Row_C_",_NROWC_); ELSE _THRNM_=THETARNM; _CANNM_=_NMLIST_("CanVar",NCOL(_THCNM_)); * Print C, U, & THETA0 *; IF _OPT_(| 5 , 1 |) & (_OPT_(|1,6|)=0) THEN PRINT , C (|ROWNAME=_THRNM_ COLNAME=_XNM_ &_FMT4_|); IF _OPT_(| 5 , 2 |)&(_OPT_(|1,6|)=0) THEN DO; IF NCOL(U)>0 THEN PRINT , U (| ROWNAME=_YNM_ COLNAME=_THCNM_ &_FMT4_ |); ELSE PRINT , "U not specified, default is I(ncol(BETA))"; END; IF _OPT_(|5,3|) & (NROW(THETA0)>0)&(_OPT_(|1,6|)=0) THEN PRINT , THETA0 (|ROWNAME=_THRNM_ COLNAME=_THCNM_ &_FMT4_|); ** Calculate secondary parameters **; _B_ = _SS_(| _XINDEX_ , _YINDEX_ |); _XPXINV_ = _SS_(| _XINDEX_ , _XINDEX_ |); _THETA_ = C * _B_; _DF_ = _PARM_(| 1 , 1 |) - _PARM_(| 1 , 3 |); IF NCOL( U ) > 0 THEN _THETA_ = _THETA_ * U; IF NCOL( THETA0 ) > 0 THEN _THETA_ = _THETA_ - THETA0; _MID_ = C * _XPXINV_ * C`; _E_=_SS_(|_YINDEX_,_YINDEX_|); * Error SSCP for original Y's; IF NCOL(U)>0 THEN _E_=U`*_E_*U;* Error SSCP for transformed Y's; _MSE_ = _E_ / _DF_; * _MSE_ is error covariance matrix estimate for transformed Y's; *Avoid zero variances; ZERO_VAR=LOC(VECDIAG(_MSE_=0)); IF NCOL(ZERO_VAR)>0 THEN DO; IF _OPT_(|1,6|)=0 THEN PRINT "TSTGLH: Zero variances in _MSE_=U`*_E_*U. Module stops", _MSE_ (| COLNAME=_THCNM_ ROWNAME=_THCNM_ &_FMT5_ |) ; _ECODE_ (| 5 , 1 |) = 3011; RETURN; END; _TOLR_ = _PARM_(| 1 , 5 |); _DIAG_ = ( VECDIAG( _XPXINV_ ) < _TOLR_ ); _LOC_ = _DIAG_`; _LOC_ = LOC( _LOC_ ); _EST_ = 0; IF NCOL( _LOC_ ) = 0 THEN _EST_ = 1; ELSE DO; _W2_ = _XPXINV_(| , _LOC_ |); _I2_ = I( NCOL( _LOC_ ) ); _K1_ = 1; DO _J1_ = 1 TO NROW( _W2_ ); IF _DIAG_(| _J1_ , 1 |) = 1 THEN DO; _W2_ (| _J1_ , |) =- _I2_(| _K1_ , |); _K1_ = _K1_ + 1; END; END; _ZERO_ = C * _W2_; IF ( ABS( _ZERO_ ) < _EPSL_ ) THEN _EST_ = 1; ELSE DO; IF _OPT_(|1,6|)=0 THEN PRINT , "TESTGLH: _THETA_ is not estimable"; _ECODE_ (| 5 , 1 |) = 3006; IF (_OPT_(|1,5|)=0)&(_OPT_(|1,6|)=0) THEN PRINT / ; RETURN; END; END; * of do group for NCOL(_LOC_) > 0 ; * Check THETA for testability **; _TST_ = 1; _MIDINV_ = SWEEP( _MID_ ); IF ANY( VECDIAG( _MIDINV_ ) < _TOLR_ ) THEN DO; IF _OPT_(|1,6|)=0 THEN PRINT , "TESTGLH: C*(( X` * X ) **-) * C ` is singular" , " : only scalar tests can be performed"; _TST_ = 0; _ECODE_ (| 5 , 1 |) = 2001; END; ELSE DO; _H_=_THETA_`*_MIDINV_*_THETA_; _MSH_=_H_/_NROWC_; END; * Print secondary parameter matrices if requested; IF _OPT_(| 5 , 4 |)&(_OPT_(|1,6|)=0) THEN IF NROW(THETA0)>0 THEN PRINT ,, "THETA is the estimate of CBU-THETA0" , _THETA_ (|ROWNAME=_THRNM_ COLNAME=_THCNM_ &_FMT4_|); ELSE PRINT , "THETA is the estimate of CBU" , _THETA_ (|ROWNAME=_THRNM_ COLNAME=_THCNM_ &_FMT4_|); IF _OPT_(| 5,5|)&(_OPT_(|1,6|)=0) THEN PRINT ,, "_MID_ = C * ((X`*X)**-1) * C`" , _MID_ (|ROWNAME=_THRNM_ COLNAME=_THCNM_ &_FMT3_|); * Compute stats for columns of THETA if they are to be printed ; * or saved in matrix form ; IF ( _OPT_(| 5 , 6 |) = 1 ) | ( _OPT_(| 5 , 7 |) = 1 ) THEN DO; *Here compute the stats one column @ time if they are wanted; DO _I1_ = 1 TO NCOL( _THETA_ ); _SD1_ = SQRT( VECDIAG( _MID_ ) # _MSE_(| _I1_, _I1_ |) ); _T1_ = _THETA_(| , _I1_ |) # ( J(NROW(_SD1_),NCOL(_SD1_),1) / ( _SD1_ ) ); * reciprocal of _SD1_ *; _DF1_ = J( NROW( _THETA_ ), 1, _DF_ ); _PV1_ = (J( NROW( _THETA_ ),1,1)-PROBT(ABS(_T1_),_DF1_))#2; IF _OPT_(| 5 , 17 |) THEN DO; _RSQ1_ = ( _T1_ # _T1_ ) / _DF1_; _RSQ1_ = _RSQ1_ / ( 1 + _RSQ1_ ); END; IF _OPT_(| 5 , 6 |) = 1 THEN DO; * print column stats; _STAT1_=_THETA_(|,_I1_ |) ||_SD1_|| _T1_||_DF1_||_PV1_; IF _OPT_(| 5 , 17 |) THEN _STAT1_ = _STAT1_ || _RSQ1_; _ST1CNM_=_THCNM_(|1,_I1_|) || {"Std Err" "t Value" "DF" "2tail p" "R Sqrd"}; IF _OPT_(|1,6|)=0 THEN PRINT ,, "Column of THETA with Associated Statistics" , _STAT1_ (| COLNAME=_ST1CNM_ ROWNAME=_THRNM_ &_FMT4_ |); END; IF _OPT_(| 5 ,7 |)=1 THEN DO; *save column stats in matrix form; _SDTHTA_=_SDTHTA_ || _SD1_; _TTHTA_=_TTHTA_ || _T1_; _PVTHTA_=_PVTHTA_ || _PV1_; IF _OPT_(| 5 ,17 |) THEN _RSTHTA_=_RSTHTA_ || _RSQ1_; END; * of column stat saving do group *; END; * of iterative do for columns of THETA *; END; * of do group for printing or storing column stats *; IF ( ^ _TST_ ) THEN DO; *Any tests requested?; IF (_OPT_(| 1 ,5 |)=0)&(_OPT_(|1,6|)=0) THEN PRINT /; RETURN; END; * Compute and print univariate tests for the columns of THETA ; IF ( _OPT_(| 5 ,8 |)=1) | (_NCOLU_=1) THEN DO; _VMSH_ = VECDIAG(_MSH_) ; _VMSE_ = VECDIAG(_MSE_) ; _FTHTA_= _VMSH_ / _VMSE_ ; _NDF_ =J(_NCOLU_,1,_NROWC_); _DDF_ =J(_NCOLU_,1,_DF_ ); _FPVAL_=J(_NCOLU_,1,1) - PROBF(_FTHTA_,_NDF_,_DDF_); _FSTATS_= _FTHTA_||_NDF_||_DDF_||_FPVAL_; IF _OPT_(| 5 , 17 |) THEN DO; _RSQ1_ = _VMSH_ / (_VMSH_ + (_VMSE_/_NROWC_)#_DF_ ) ; _FSTATS_=_FSTATS_ || _RSQ1_; END; _FSTRNM_={ "F Value" "Num df" "Den df" "p Value" "R Sqrd" }; IF _OPT_(|1,6|)=0 THEN PRINT ,, "Univariate Tests for columns of THETA=C*BETA*U - THETA0" , _FSTATS_ (| COLNAME=_FSTRNM_ ROWNAME=_THCNM_ &_FMT4_ |); END; * of do group for computing and printing col(THETA) F tests; * Print MSH, MSE, and R if requested ; IF _OPT_(| 5 ,9 |)&(_OPT_(|1,6|)=0) THEN PRINT , "_MSH_=_H_/nrow(C) = Mean Square Hypothesis" _MSH_ (| COLNAME=_THCNM_ ROWNAME=_THCNM_ &_FMT5_ |); IF _OPT_(| 5 ,10 |)&(_OPT_(|1,6|)=0) THEN PRINT , "_MSE_=U`*_SIGMA_*U=_E_/_DF_, Error covariance for transformed Y's" , _MSE_ (| COLNAME=_THCNM_ ROWNAME=_THCNM_ &_FMT5_ |); IF (_OPT_(| 5 ,11 |)=1)&(_NCOLU_>1) THEN DO; _T1_=DIAG({1}/SQRT(VECDIAG(_E_))); _R_=_T1_ * _E_ * _T1_; IF _OPT_(|1,6|)=0 THEN PRINT , "Estimated Correlation Matrix based on E" , _R_ (|ROWNAME=_THCNM_ COLNAME=_THCNM_ &_FMT6_ |); END; * This section computes and prints multivariate statistics & tests ; IF (NCOL(_THETA_)>1) & MAX( _OPT_(| 5,12:16 |))=1 THEN DO; *Insure that _E_ will be perfectly symmetric; _E_=( _E_ + _E_` ) / 2 ; _R1_=SOLVE(HALF(_E_) , I(_NCOLU_) ); _T1_=_R1_` * _H_ * _R1_; *Insure that _T1_ will be perfectly symmetric; _T1_=( _T1_ + _T1_` ) / 2 ; CALL EIGEN( _HEIVAL_ ,_E_VEC_ , _T1_ ); _RS1_=MIN( NCOL( _THETA_ ) || NROW( C ) ); _HEIVAL_=_HEIVAL_(| 1 : _RS1_ , |); _E_VEC_=_E_VEC_(| , 1 : _RS1_ |); IF _OPT_(| 5 , 12 |)&(_OPT_(|1,6|)=0) THEN PRINT , "Eigenvalues of _H_ * INV( _E_ )" , _HEIVAL_ (| ROWNAME=_CANNM_ COLNAME=_NONM_ &_FMT3_ |); IF _OPT_(| 5 , 13 |)=1 THEN DO; IF (_OPT_(|1,5|)=0)&(_OPT_(|1,6|)=0) THEN PRINT " Left eigenvectors of _H_ * INV( _E_ ) " , " which equal left eigenvectors of _H_ * INV( _T_ ) " , " normalized such that transpose ( _A_ ) * _E_ * _A_=I. ) " , "JTH column proportional to JTH canonical variate raw weights"; _CANVEC_=_R1_ * _E_VEC_; IF _OPT_(|1,6|)=0 THEN PRINT _CANVEC_ (| ROWNAME=_THCNM_ COLNAME=_CANNM_ &_FMT5_ |); END; IF _OPT_(| 5 , 14 |)=1 THEN DO; _CANRSQ_=_HEIVAL_ / ( 1 + _HEIVAL_ ); IF _OPT_(|1,6|)=0 THEN PRINT , "Generalized squared canonical correlations:" _CANRSQ_ (|ROWNAME=_CANNM_ COLNAME=_NONM_ &_FMT6_ |); END; IF _OPT_(| 5 , 15 |)&(_OPT_(|1,6|)=0) THEN PRINT , "EVEC2 option has not yet been implemented"; * This section computes multivariate test criteria ; * Set up matrix to hold test statistics; IF _OPT_(| 5 , 17 |) THEN _STMAT1_=J( 4, 6, . ); ELSE _STMAT1_=J( 4, 5, . ); * Check for special cases in which an exact F test can be; * computed rather than approximations ; _SCASE_=0; IF NROW( _THETA_ )=2 THEN _SCASE_=2; IF NCOL( _THETA_ )=2 THEN _SCASE_=4; IF NROW( _THETA_ )=1 THEN _SCASE_=1; * Compute basic parameters & store for printing later *; IF _OPT_(| 5 , 16 |)=1 THEN DO; _RP1_=NCOL( _THETA_ ); _RQ1_=NROW( C ); _RS1_=MIN( _RP1_ || _RQ1_ ); _RM1_=( ABS( _RP1_ - _RQ1_ ) - 1 ) / 2; _RN1_=( _DF_ - _RP1_ - 1 ) / 2; _TPARM1_=_DF_||_RQ1_||_RP1_||_RS1_||_RM1_||_RN1_; _RLROOT_=MAX( _HEIVAL_ ); _RLROOT_=_RLROOT_ # 1.0 / ( 1 + _RLROOT_ ); ****************************************************************; *The next block of code uses an algorithm due to Pillai (1965) *; *to approximate the right tail probability of largest root test*; *Code in appendix in Harris (1975),Primer of Multivariate *; *Statistics. A missing value is returned if P > .10 because *; *the algorithm only works for small P. *; *_RLROOT_ is largest root of H*INV(T) *; *_RS1_ is S=MIN(RK(C),RK(U)) *; *_RM1_ is M *; *_RN1_ is N *; *_P1_ is P-Value of largest root test *; ****************************************************************; IF (_RS1_=1) THEN DO; * Compute P value for special case S=1 ; _DDF_=_DF_ + 1 - _RP1_; _NDF_=_RP1_; IF ( _DDF_ <= 0 ) THEN DO; _P1_=.; _F1_=.; _ECODE_ (| 5 , 1 |)=3011; END; ELSE DO; _F1_=_HEIVAL_ # _DDF_ / _NDF_; _P1_=1 - PROBF( _F1_, _NDF_, _DDF_ ); END; END; ELSE DO; * Compute P-value for S>1 (multivariate) ; _P1_=.; _F1_=.; _DDF_=.; _NDF_=.; IF (_RN1_<=-1) THEN _ECODE_(| 5, 1 |)=3011; ELSE DO; IF (ABS(_RLROOT_)<=.000001) THEN _P1_=1; ELSE IF (ABS(_RLROOT_)<=.000001) THEN _P1_=1; ELSE DO; _LOGR_= LGAMMA( _RM1_ + _RN1_ + _RS1_ + 1 ) - LGAMMA( _RM1_ + _RS1_ / 2 + .5 ) - LGAMMA( _RN1_ + _RS1_ / 2 + .5 ) + LGAMMA( _RM1_ + _RN1_ + _RS1_ + .5 ) - LGAMMA( _RM1_ + _RN1_ + _RS1_ / 2 + 1 ) - LGAMMA( _RS1_ / 2 ) + .5 # 1.1447298858494 ; * ='s ln(pi); _COEFK_=0; _BINCO_=(_RLROOT_ ##(_RM1_+_RS1_)) # ( (1-_RLROOT_) ## (_RN1_+1) ); _ELOGR_=EXP( _LOGR_ ); _SUM_=0; _SGN_=1; _PROD_=1; _COMB_=1; _A1_=_RM1_ + _RS1_ + 1; _A2_=_RM1_ + _RN1_ + _RS1_ + 1; _A3_=_RM1_ + .5 + _RS1_ / 2; _A4_=_RM1_ + .5 + _RN1_ + _RS1_; DO _I1_=1 TO ( _RS1_ - 1 ); _TERML_=_ELOGR_#_PROD_ # _COMB_/( _A2_-_I1_ ); _TERMR_=_COEFK_#( _A1_ - _I1_ )/( _A2_-_I1_ ); _COEFK_=_TERML_ - _TERMR_; _BINCO_=_BINCO_ / _RLROOT_; _SGN_ =- _SGN_; _SUM_=_SUM_ + _COEFK_ # _BINCO_ # _SGN_; _PROD_=_PROD_ # (_A3_-(_I1_#.5)) / (_A4_-(_I1_#.5)) ; _COMB_=_COMB_ # ( _RS1_ - _I1_ ) / _I1_; END; IF (2#INT(_RS1_/2)=_RS1_) THEN _P1_=-_SUM_; ELSE _P1_=1-_SUM_ - PROBBETA(_RLROOT_,_RM1_+1,_RN1_+1); IF _P1_ < 0 THEN _P1_=0; IF _P1_ > .1 THEN _P1_=.; END; * of non-boundary pt, S>1, valid P value ; END; * of valid P value for S>1 ; END; * of S>1 block ; _T1_=_RLROOT_ || _RS1_ || _RM1_ || _RN1_ || _P1_; _FSTATS_=( _RLROOT_ || _F1_ || _NDF_ || _DDF_ || _P1_ ); IF _OPT_(| 5 , 17 |) THEN _FSTATS_=_FSTATS_ || _RLROOT_; _STMAT1_ (| 1 , |)=_FSTATS_; _L1_=1; DO _I1_=1 TO _RS1_; _L1_=_L1_ * ( 1 + _HEIVAL_(| _I1_ , 1|) ); END; _L1_=J( NROW(_L1_), NCOL(_L1_), 1 ) / ( _L1_ ); IF (_SCASE_=0) THEN DO; * Compute Rao's F approximation for the LR criteria ; _V1_=_DF_ - ( _RP1_ - _RQ1_ + 1 ) / 2; _Z1_=_RP1_ # _RP1_ # _RQ1_ # _RQ1_ - 4; _Z1_=_Z1_/( ( _RP1_ # _RP1_ )+( _RQ1_ # _RQ1_ )-5 ); _Z1_=SQRT( _Z1_ ); _B1_=( _RP1_ # _RQ1_ - 2 ) / 4; _F1_=( _V1_ # _Z1_ - 2 # _B1_ ) / ( _RP1_ # _RQ1_ ); _F1_=_F1_ # ( 1 - _L1_ ## ( 1 / _Z1_ ) ) / ( _L1_ ## ( 1 / _Z1_ ) ); _NDF_=_RP1_ # _RQ1_; _DDF_=_V1_ # _Z1_ - 2 # _B1_; IF ( _DDF_ <= 0 ) THEN DO; _FPROB1_=.; _F1_=.; _ECODE_ (| 5 , 1 |)=3011; END; ELSE _FPROB1_=1 - PROBF( _F1_, _NDF_, _DDF_ ); _FSTATS_=_F1_ || _NDF_ || _DDF_ || _FPROB1_; END; ELSE DO; *Compute exact F statistic for special cases; IF _SCASE_=1 THEN DO; _F1_=( _DF_ + _RQ1_ - _RP1_ ) / _RP1_; _F1_=_F1_ # ( ( 1 - _L1_ ) / _L1_ ); _NDF_=_RP1_; _DDF_=_DF_ + _RQ1_ - _RP1_; END; IF _SCASE_=2 THEN DO; _F1_=( _DF_ + _RQ1_ - _RP1_ - 1 ) / _RP1_; _F1_=_F1_ # ( ( 1 - SQRT( _L1_ )) / SQRT( _L1_ )); _NDF_=2 # _RP1_; _DDF_=2 # ( _DF_ + _RQ1_ - _RP1_ - 1 ); END; IF _SCASE_=4 THEN DO; _F1_=( _DF_ - 1 ) / _RQ1_; _F1_=_F1_ # ( ( 1 - SQRT( _L1_ )) / SQRT( _L1_ )); _NDF_=2 # _RQ1_; _DDF_=2 # ( _DF_ - 1 ); END; * Print F stat & P value ; IF ( _DDF_ <= 0 ) THEN DO; _PVAL1_=.; _F1_=.; _ECODE_ (| 5 , 1 |)=3011; END; ELSE _PVAL1_=1 - PROBF( _F1_, _NDF_, _DDF_ ); _FSTATS_=_F1_ || _NDF_ || _DDF_ || _PVAL1_; END; * of do group for special cases ; IF _OPT_(| 5 , 17 |) THEN DO; _RSQ1_=1 - ( _L1_ ## ( 1 / _RS1_ ) ); _FSTATS_=_FSTATS_ || _RSQ1_; END; _STMAT1_ (| 2 , |)=( _L1_ || _FSTATS_ ); * tr(_H_*inv(_E_)), Hotelling-Lawley trace ; _TR1_=SUM( _HEIVAL_ ); IF ( _SCASE_ ^= 1 ) THEN DO; _F1_=2 # ( ( _RS1_ # _RN1_ ) + 1 ) # _TR1_; _F1_=_F1_/( _RS1_#_RS1_#( ( 2 # _RM1_ ) + _RS1_ + 1)); _NDF_=_RS1_ # ( ( 2 # _RM1_ ) + _RS1_ + 1 ); _DDF_=2 # ( ( _RS1_ # _RN1_ ) + 1 ); *Use McKeon (1974, Biometric, 61,p381-383) degrees of freedom; *Use Pillai F approximation if McKeon not defined; MCK1=(_DF_)#(_DF_) - (_DF_)#(2#_RP1_ + 3) + _RP1_#(_RP1_ + 3); MCK2=MCK1/ ((_DF_)#(_RQ1_+_RP1_+1)-(_RQ1_+2#_RP1_+_RP1_##2-1)); MCKDF2= 4 + (_RQ1_#_RP1_ + 2)# MCK2; IF (MCKDF2<=0)&(_OPT_(|1,5|)=0)&(_OPT_(|1,6|)=0) THEN PRINT "Warning: McKeon method not defined due to small N." , "Pillai F approximation used for Hotelling-Lawley. "; ELSE DO; _DDF_=MCKDF2; MCK_GAM=_NDF_#(_DDF_-2)/((_DF_-_RP1_-1)#_DDF_); _F1_=_TR1_/MCK_GAM; END; IF ( _DDF_ <= 0 ) THEN DO; _PVAL1_=.; _F1_=.; _ECODE_ (| 5 , 1 |)=3011; END; ELSE _PVAL1_=1 - PROBF( _F1_, _NDF_, _DDF_ ); _FSTATS_=_F1_ || _NDF_ || _DDF_ || _PVAL1_; IF _OPT_(| 5 , 17 |) THEN DO; _RSQ1_=_TR1_ / ( _RS1_ + _TR1_ ); _FSTATS_=_FSTATS_ || _RSQ1_; END; END; _STMAT1_ (| 3 , |)=( _TR1_ || _FSTATS_ ); * Pillai's Trace, tr( _H_*inv(_H_ + _E_)); _TR2_=SUM( _HEIVAL_ / ( 1 + _HEIVAL_ )); IF ( _SCASE_ ^= 1 ) THEN DO; _S1_=MIN(_RP1_ || _RQ1_); _NDF_=_RQ1_ # _RP1_ ; _DDF_=_S1_ # ( _DF_ + _S1_ - _RP1_ ); *Modify Pillai-Bartlett trace degrees of freedom as per Muller (1998, Journal of Computational & Graphical Statistics); PBTK=( _S1_#(_DF_+_S1_-_RP1_)#(_DF_+_RQ1_+2)#(_DF_+_RQ1_-1) / (_DF_#(_DF_+_RQ1_-_RP1_)) -2 )/ (_S1_#(_DF_+_RQ1_)); _NDF_=PBTK#_NDF_; _DDF_=PBTK#_DDF_; PBTRATIO= _TR2_ / ( _S1_ - _TR2_ ); _F1_=_DDF_#PBTRATIO/_NDF_; IF ( _DDF_ <= 0 ) THEN DO; _PVAL1_=.; _F1_=.; _ECODE_ (| 5 , 1 |)=3011; END; ELSE _PVAL1_=1 - PROBF( _F1_, _NDF_, _DDF_ ); _FSTATS_=_F1_ || _NDF_ || _DDF_ || _PVAL1_; IF _OPT_(| 5 , 17 |) THEN DO; _RSQ1_=_TR2_ / _RS1_; _FSTATS_=_FSTATS_ || _RSQ1_; END; END; _STMAT1_ (| 4 , |)=( _TR2_ || _FSTATS_ ); IF (_OPT_(|1,5|)=0)&(_OPT_(|1,6|)=0) THEN DO; PRINT "Multivariate Test Criteria for Ho: THETA=0"; IF _SCASE_=1 THEN PRINT "Since NROW ( THETA )=1 F Tests are exact"; ELSE IF _SCASE_ ^= 0 THEN PRINT "Since MIN( NROW( THETA ), NCOL( THETA ) )=2," , "F Test for LAMBDA is exact"; END; IF (_ECODE_(|5,1|)=3011)&(_OPT_(|1,5|)=0)&(_OPT_(|1,6|)=0) THEN PRINT "At least one of the four multivariate test statistics" , "had degrees of freedom <= 0. For these statistics " , "the P value has been set to missing. " ; _STMCNM_={"VALUE" "APPROX F" "NUM DF" "DENOM DF" "P VALUE" "ASSOCITN" }; _STMRNM_={ "LROOT" "LAMBDA" "HLTRACE" "PILTRACE" }; IF _OPT_(|5,16|)&(_OPT_(|1,6|)=0) THEN DO; PRINT _STMAT1_ (|COLNAME=_STMCNM_ ROWNAME=_STMRNM_ &_FMT4_|); _TPCNM1_={ "DF" "A" "B" "S" "M" "N" }; PRINT "Parameters for the Criteria" , _TPARM1_ (| ROWNAME=_NONM_ COLNAME=_TPCNM1_ |); IF _STMAT1_(|1,5|)=. THEN PRINT "Missing value given for largest root P value if P > .10 " , "since algorithm usually inaccurate there ( Pillai, 1965 )" ; END; END; * MULTTEST option, multivariate tests ; END; * of section for all multivariate statistics ; * Run UNIREP analysis * ; IF _OPT_(| 6, 1 |)=1 THEN RUN _UNIREP_(U, _ECODE_,_EPSL_,_OPT_,_PARM_, _H_,_E_,_NROWC_,_NCOLU_,_THCNM_, _URESUL_,_UCOLNM_,_UROWNM_ ); IF (_OPT_(| 1 , 5 |)=0)&(_OPT_(|1,6|)=0) THEN PRINT /; FINISH _TSTGLH_ ;