*************************************************************************; * August 26, 2008 *; *************************************************************************; * 23456789A123456789B123456789C123456789D123456789E123456789F123456789G ; *************************************************************************; * *; * This SAS code, used inside PROC IML via an INCLUDE statement, *; * performs a power analysis of multivariate linear hypotheses *; * of the form C*BETA*U - THETA0 = 0. The associated model is *; * Y=X*BETA+E. Users must enter PROC IML to create the following *; * matrices. *; * *; * Note: Vectors may be entered as rows or columns. *; * *; * A) REQUIRED MATRICES *; * *; * (1) ESSENCEX, the essence design matrix. Users unfamiliar with *; * this matrix may simply enter the full design matrix and *; * specify that REPN below be equal to 1. *; * *; * (2) SIGMA, the hypothesized covariance matrix *; * *; * (3) BETA, the matrix of hypothesized regression coefficients *; * *; * (4) C, "between" subject contrast for pre-multiplying BETA *; * *; * WARNING: C must be full rank. *; * By default, ESSENCEX must be full rank, unless option *; * LTFR is on *; * *; * B) OPTIONAL MATRICES *; * *; * (1) U, "within" subject contrast for post-multiplying BETA *; * Default of I(P), where P = NCOL(BETA). *; * *; * WARNING: U must be full rank. *; * U`*SIGMA*U must be symmetric and positive definite. *; * For univariate repeated measures tests, U must be *; * proportional to an orthonormal matrix and orthogonal *; * to a column of 1's, unless option ORTHU is specified. *; * If ORTHU is specified, the program computes a U matrix *; * with the required properties from the user upplied U. *; * *; * (2) THETA0, the matrix of constants to be subtracted from *; * C*BETA*U (CBU). Default matrix of zeroes. *; * *; * (3) REPN, (vector), the # of times each row of ESSENCEX is *; * duplicated. Default of {1}. *; * *; * (4) BETASCAL, (vector), multipliers of BETA. *; * Default of {1}. *; * *; * (5) SIGSCAL, (vector), multipliers of SIGMA. *; * Default of {1}. *; * *; * (6) RHOSCAL, (vector), multipliers of RHO (correlation matrix *; * computed internally from SIGMA). *; * Default of {1}. *; * *; * (7) ALPHA, (vector), the Type I error rate (test size). *; * Default of .05. *; * *; * (8) ROUND, (scalar) # of places to which power values will *; * be rounded. Default of 3. *; * *; * (9) TOLERANCE, (scalar) value not tolerated, numeric zero, *; * used for checking singularity, division problems, etc. *; * singularity. Default of 1E-12. *; * *; * (10) UCDF, a 1x1, 4x1 or 1x4 vector of choices of calculation *; * for UNIREP test statistic CDF. *; * *; * If a 4x1 or 1x4 matrix is specified, the elements correspond *; * to choices for calculation of the CDF of: *; * *; * UCDF[1] --> Uncorrected test *; * UCDF[2] --> Huynh-Feldt test *; * UCDF[3] --> Geisser-Greenhouse test *; * UCDF[4] --> Box conservative test *; *; * Where: *; * *; * UCDF[J] =1 --> Muller and Barton (1989) approximation *; * =2 --> Muller, Edwards and Taylor (2004) approximation *; * =3 --> Muller, Edwards and Taylor (2004) exact, via *; * Davies' algorithm. This algorithm may fail -- if it *; * does, this choice will yield missing values. *; * =4 --> Muller, Edwards and Taylor (2004) exact as in *; * choice 3, except if the Davies' algorithm fails, *; * this replaces the missing value with the value *; * resulting from choice 2 *; * *; * Specifying a 1x1 matrix with choice 1, 2, 3, or 4 as described *; * above gives the same method for all four tests. *; * *; * If UCDF is not specified, the default is {2 2 2 2}`. *; * *; * (11) UMETHOD, a 1x1, 2x1, or 1x2 vector of choices of approximation*; * of mean epsilon hat and mean epsilon tilde for the Huynh-Feldt*; * and Geisser-Greenhouse UNIREP tests. *; * *; * If a 2x1 or 1x2 matrix is specified, the elements correspond *; * to choices for: *; * *; * UMETHOD[1] --> Huynh-Feldt test *; * UMETHOD[2] --> Geisser-Greenhouse test *; * *; * and *; * *; * UMETHOD[J] =1 --> Muller and Barton (1989) approximation *; * =2 --> Method 1, Muller, Edwards, and Taylor (2004)*; * *; * Specifying a 1x1 matrix causes the same approximation to be *; * used for both tests, with choice 1 or 2 as described above. *; * *; * If UMETHOD is not specified, the default is {2, 2} *; * *; * *; * O'Brien & Shieh (1992, OS) recommended multipling noncentrality by *; * N/[N-rankk(X)] if min(rk(C),rk(U))>1 for multivariate tests. *; * *; * (12) MMETHOD default= {4,2,2}; *; * MMETHOD[1] Choices for Hotelling-Lawley Trace *; * = 1 Pillai (1954, 55) 1 moment null approx *; * = 2 McKeon (1974) two moment null approx *; * = 3 Pillai (1959) one moment null approx+ OS noncen mult*; * = 4 McKeon (1974) two moment null approx+ OS noncen mult*; * MMETHOD[2] Choices for Pillai-Bartlett Trace *; * = 1 Pillai (1954, 55) one moment null approx *; * = 2 Muller (1998) two moment null approx *; * = 3 Pillai (1959) one moment null approx + OS noncen mult*; * = 4 Muller (1998) two moment null approx + OS noncen mult*; * MMETHOD[3] Choices for Wilks' Lambda *; * = 1 Rao (1951) two moment null approx *; * = 2 Rao (1951) two moment null approx *; * = 3 Rao (1951) two moment null approx + OS noncen mult *; * = 4 Rao (1951) two moment null approx + OS noncen mult *; * *; * *; * The following 5 parameters control creation of confidence intervals*; * for power values of univariate and MULTIREP tests. Uses methods of *; * Taylor and Muller (1995) for univariate models. **This theory has *; * not yet been extended for the UNIREP tests.** *; * *; * (13) CLTYPE, (scalar) *; * *; * <1 or 0 by 0 matrix (never defined, or FREE-ed) *; * if confidence intervals for power **not** desired *; * =1 if BETA known , SIGMA estimated and CL desired *; * =2 if BETA estimated, SIGMA estimated and CL desired *; * *; * (14) N_EST, (scalar) # of observations in analysis which yielded *; * BETA and SIGMA estimates (required if CLTYPE>=1) *; * *; * (15) RANK_EST, (scalar) rank of design matrix in analysis which *; * yielded BETA and SIGMA estimates (required if CLTYPE>=1) *; * *; * (16) ALPHA_CL, (scalar) lower tail probability for power C.L. *; * default=.025 *; * *; * (17) ALPHA_CU, (scalar) upper tail probability for power C.L. *; * default=.025 *; * *; * *; * (18) DSNAME, (1 row, 2-4 cols), name of output SAS data file. *; * Default WORK.PWRDT###. *; * *; * *; * (19) OPT_ON, OPT_OFF, (vectors), users can specify the options *; * they wish to turn on or off. The possible options are as *; * follows: *; * *; * a) Power calculations Default *; * *; * COLLAPSE, special case** ON *; * HLT, Hotelling-Lawley trace ON *; * PBT, Pillai-Bartlett trace ON *; * WLK, Wilk's Lambda ON *; * UN, uncorrected univariate OFF *; * approach to repeated measures *; * HF, Huynh-Feldt OFF *; * GG, Geisser-Greenhouse ON *; * BOX, Box Conservative OFF *; * *; * **With rank(U)=1 powers of all tests coincide. *; * The COLLAPSE option produces one column of power *; * calculations labeled "POWER" instead of separate *; * columns for WLK, HLT, UNI, etc. *; * *; * With rank(C)=1 and rank(U)>1, powers of all 3 *; * multivariate tests (WLK HLT PBT) coincide. *; * The COLLAPSE option produces one column for all 3. *; * UN, GG, HF, and BOX powers are calculated if requested. *; * *; * LTFR, allow less than full rank *; * ESSENCEX matrix OFF *; * *; * FRACREPN, allow fractional REPN *; * values OFF *; * *; * ORTHU, allow non-orthonormal U *; * from user for UNIREP *; * calculations - program creates *; * orthonormal U for user OFF *; * *; * b) Print/options for output power matrix (SAS dataset) *; * *; * TOTAL_N Total number observations ON *; * SIGSCAL Multiplier for SIGMA ON *; * RHOSCAL Multiplier for RHO OFF *; * BETASCAL Multiplier for BETA ON *; * ALPHA Size of test ON *; * POWERCASE Row # of power matrix OFF *; * MAXRHOSQ Max canonical rho squared OFF *; * FMETHOD Method used for F CDF OFF *; * UCDF Method used for UNIREP test CDF OFF *; * UMETHOD Method used for approx expected value OFF *; * of epilon estimator in GG and HF OFF *; * *; * All of the following are OFF if CLTYPE < 1, otherwise *; * they appear as below: *; * *; * CLTYPE Type of BETA and SIGMA for CL ON *; * N_EST Number of obs in analysis creating *; * BETA/SIGMA estimates OFF *; * RANK_EST Rank of X in analysis creating *; * BETA/SIGMA estimates OFF *; * ALPHA_CL CL Lower tail probability ON *; * ALPHA_CU CL Upper tail probability ON *; * *; * c) Print options - additional matrices *; * *; * ESSENCEX EssenceX matrix ON *; * BETA Beta matrix ON *; * SIGMA Variance-covariance matrix ON *; * RHO Correlation matrix ON *; * C "Between" subjects contrast ON *; * U "Within" subject contrast ON *; * THETA0 Constant to compare to CBU ON, if ^= 0 *; * CBETAU Value of initial C*BETA*U ON *; * RANKX Rank of design matrix X OFF *; * RANKC Rank of C OFF *; * RANKU Rank of U OFF *; * *; * d) Power data options *; * *; * NOPRINT Printed output suppressed OFF *; * DS Request output SAS dataset OFF *; * *; * e) WARN Allows warning messages ON *; * to print. *; * *; * C) SAS DATA FILE *; * *; * If the DS option is selected, the output data file will contain *; * all options requested in a) and b). The user can name the data *; * file by defining the matrix DSNAME= { libref dataset }. For *; * example, if DSNAME = {IN1 MYDATA}, the output data file will be *; * called IN1.MYDATA. IN1 refers to a library defined by a DD *; * statement or a LIBNAME statement. *; * *; * When using a library other than WORK as the default, define *; * DSNAME = {libref dataset defaultlib}; *; * *; * If DSNAME is not defined or if the "dataset" already exists in *; * the library specified by "libref", a default file name is used. *; * The default file names are numbered and of the form PWRDT### *; * (where ### is a number). The program scans the library for the *; * largest numbered data file and assigns the next number to the new *; * data file. The maximum ### is 999.If PWRDT999 exists no more data *; * files can be created. *; * *; * NOTE: The program uses the name _PWRDTMP as an intermediate *; * dataset. If this dataset already exists in the specified library *; * no datasets can be created. *; * *; *************************************************************************; *****MODULE DEFINITIONS*****; *************************************************************************; * *; * The following are natural groupings for the modules required for *; * this software if you wish to split this code into multiple files. *; * *; * 1) POWER, _POWER, *; * 2) _RANKSYMM, _SIZECHK, _TYPEMISSCHK, _SETOPT *; * 3) _HLT, _PBT, _WLK, _SPECIAL, _FIRSTUNI, _HFEXEPS, _GGEXEPS, *; * _LASTUNI, _GLMMPCL *; * 4) _PROBF, _FINV, _FWARN, _SASDS *; * 5) _QPROB, _AS, _IROUND, _COUNTR, _ALOG1, _EXP1, _ORDER, _ERRBD, *; * _CTFF, _TRUNCN, _FINDU, _INTEGR, _CFE *; * 6) NAMELIST, UMEAN, UPOLY1, UPOLY2, UPOLY3 *; * *; * The modules contained in block 5 perform exact computations via *; * Davies' algorithm from Muller Edwards Tayler (2004). _QPROB is the *; * module in block 5 called directly in this power software. The other *; * modules in block 5 are called only by _QPROB or within the *; * modules _QPROB calls. If the exact methods are not needed, then *; * the code for all the block 5 modules can be deleted. *; * This will speed up the time required to read this file because *; * the block 5 modules constitute roughly 20% of the lines of code. *; * *; *************************************************************************; *************************************************************************; * *; * POWER *; * *; * Define POWER command. *; * User will only have to type RUN POWER. The input matrices will not *; * have to be listed in the RUN statement. *; * *; * INPUTS All global, and exactly the same as _POWER( ) *; * documented immediately below *; * *; * OUTPUTS _HOLDPOWER, _HOLDPOWERLBL, _POWERWARN, _POWERWARNLBL *; * which are *not* global. *; * *; *************************************************************************; START POWER; CALL _POWER (_HOLDPOWER, _HOLDPOWERLBL, _POWERWARN, _POWERWARNLBL); FINISH POWER; *************************************************************************; * *; * Primary module: _POWER( ) *; * *; * INPUTS *; * *; * 1) User supplied - required - GLOBAL *; * *; * ESSENCEX, essence design matrix *; * SIGMA, hypothesized variance-covariance matrix *; * BETA, matrix of hypothesized regression coefficients *; * C, "between" subject contrast *; * *; * 2) User supplied - optional (program defaults supplied) - GLOBAL *; * *; * U, "within" subject contrast *; * THETA0, the matrix of constants to be subtracted from CBU *; * REPN, (vector), the # of times each row of ESSENCEX is *; * duplicated *; * BETASCAL, (vector) multipliers for BETA *; * SIGSCAL, (vector) multipliers for SIGMA *; * RHOSCAL, (vector) multipliers for RHO *; * ALPHA, (vector) Type I error rates *; * ROUND, (scalar) # of places to round power calculations *; * TOLERANCE, (scalar) value not tolerated, numeric zero, *; * used for checking singularity. *; * TOLERANC, name for TOLERANCE in previous version of software *; * UCDF, (vector) choices of approximation for UNIREP CDF *; * UMETHOD, (vector) choices of approximation of mean epsilon tilde*; * and mean epsilon hat for the Huynh-Feldt and *; * Geisser-Greenhouse UNIREP tests. *; * CLTYPE, (scalar) *; * <1 or 0 by 0 matrix (never defined, or FREE-ed) *; * if confidence intervals for power **not** desired *; * =1 if BETA known , SIGMA estimated and CL desired *; * =2 if BETA estimated, SIGMA estimated and CL desired *; * N_EST, (scalar) # of observations in analysis which yielded *; * BETA and SIGMA estimates (required if CLTYPE>=1) *; * RANK_EST, (scalar) design matrix rank in analysis which *; * yielded BETA and SIGMA estimates (required if CLTYPE>=1) *; * ALPHA_CL, (scalar) lower tail probability for power C.L. *; * default=.025 *; * ALPHA_CU, (scalar) upper tail probability for power C.L. *; * default=.025 *; * DSNAME, (row) name of output dataset *; * OPT_ON, (column) selected options to turn on *; * OPT_OFF, (column) selected options to turn off *; * *; * OUTPUTS *; * *; * HOLDPOWER, matrix of power calculations *; * HOLDPOWERLBL, column labels for HOLDPOWER *; * POWERWARN, vector of power calculation warning counts *; * POWERWARNLBL, labels for POWERWARN warnings *; * *; *************************************************************************; START _POWER (HOLDPOWER, HOLDPOWERLBL, POWERWARN, POWERWARNLBL) GLOBAL (ESSENCEX, SIGMA, BETA, C, U, THETA0, REPN, BETASCAL, SIGSCAL, RHOSCAL, ALPHA, ROUND, TOLERANCE, TOLERANC, UCDF, UMETHOD, MMETHOD, LIMCASE, CLTYPE, N_EST, RANK_EST, ALPHA_CL, ALPHA_CU, DSNAME, OPT_ON, OPT_OFF, NONCENCL); ***A) Initial Setup; ***A.1) Check for SAS version 7 or higher; If &SYSVER < 7 THEN DO; NOTE = CONCAT("ERROR 1: The power software uses >8 character ", "matrix names, which are allowed only in SAS ", "version 7.0 and higher. This is SAS version ", "&sysver.. The power software will not run."); PRINT NOTE[LABEL=" " ROWNAME=" " COLNAME=" "]; GO TO ERROREND; END; ***A.2) Insure everthing printed by power software goes to output file; RESET NOLOG; ***A.3) Initialize warnings; *** All warnings are stored in vector POWERWARN; *** Labels for the warnings are stored in vector POWERWARNLBL; POWERWARN = J(22,1,0); POWERWARNLBL = CONCAT( "1: A Tiku approximation (see manual) was used in calculating ", "power.") // "2: A Z approximation (see manual) was used in calculating power." // CONCAT( "3: A Z approximation (see manual) was used in calculating power ", "and |Z| > 6, so that the power returned is exactly 0 or 1.") // CONCAT( "4: Power is missing because the FINV function returned a ", "missing value.") // "5: The lower confidence limit on power is conservative." // CONCAT( "6: A Tiku approximation (see manual) was used in calculating ", "the lower confidence limit on power.") // CONCAT( "7: A Z approximation (see manual) was used in calculating the ", "lower confidence limit on power.") // CONCAT( "8: A Z approximation (see manual) was used in calculating the ", "lower confidence limit on power and |Z| > 6, so that the ", "power returned is exactly 0 or 1.") // CONCAT( "9: The lower confidence limit on power is missing because the ", "FINV function returned a missing value.") // "10: The upper confidence limit on power is conservative." // CONCAT("11: A Tiku approxmiation (see manual) was used in calculating ", "the upper confidence limit on power.") // CONCAT("12: A Z approximation (see manual) was used in calculating the ", "upper confidence limit on power.") // CONCAT("13: A Z approximation (see manual) was used in calculating the ", "upper confidence limit on power and |Z| > 6, so that the ", "power returned was exactly 0 or 1.") // CONCAT("14: The upper confidence limit on power is missing because the ", "FINV function returned a missing value.") // CONCAT("15: Power is missing because because the noncentrality could not ", "be computed.") // "16: Confidence limits are missing because power is missing." // CONCAT("17: The approximate expected value of estimated epsilon was ", "truncated up to 1/B.") // CONCAT("18: The approximate expected value of estimated epsilon was ", "truncated down to 1.") // "19: Power missing due to Davies' algorithm fail." // "20: Inputs give off-diagonal correlation = 1 in RHO." // CONCAT("21: (N - R) <= 5, so power approximations may be inaccurate, ", "especially Huynh-Feldt. See the manual.") // CONCAT("22: Power values were rounded to 1 using the value contained ", "in ROUND and should not be be reported as Power = 1. For ", "example, if ROUND = 3 then report Power > 0.999 .") ; ***A.4) Delete previous versions of HOLDPOWER and HOLDPOWERLBL; IF NROW(HOLDPOWER) > 0 THEN FREE HOLDPOWER; IF NROW(HOLDPOWERLBL) > 0 THEN FREE HOLDPOWERLBL; ***A.5) Check for required input matrices, and that they are numeric; IF (NROW(C)=0) | (NROW(BETA)=0) | (NROW(SIGMA)=0) | (NROW(ESSENCEX)=0) THEN DO; IF (NROW(C) = 0) THEN DO; NOTE = "ERROR 2: The matrix C has not been supplied."; PRINT NOTE[LABEL=" " ROWNAME=" " COLNAME=" "]; END; If (NROW(BETA) = 0) THEN DO; NOTE = "ERROR 3: The matrix BETA has not been supplied."; PRINT NOTE[LABEL=" " ROWNAME=" " COLNAME=" "]; END; IF (NROW(SIGMA) = 0) THEN DO; NOTE = "ERROR 4: The matrix SIGMA has not been supplied."; PRINT NOTE[LABEL=" " ROWNAME=" " COLNAME=" "]; END; IF (NROW(ESSENCEX) = 0) THEN DO; NOTE = "ERROR 5: The matrix ESSENCEX has not been supplied."; PRINT NOTE[LABEL=" " ROWNAME=" " COLNAME=" "]; END; GO TO ERROREND; END; IF (_TYPEMISSCHK(C ,"C" ,"N") + _TYPEMISSCHK(BETA ,"BETA" ,"N") + _TYPEMISSCHK(SIGMA ,"SIGMA" ,"N") + _TYPEMISSCHK(ESSENCEX,"ESSENCEX","N")) > 0 THEN GO TO ERROREND; ***B) Define default matrices; IF NROW(U) = 0 THEN DO; U = I(NCOL(BETA)); HITLIST = HITLIST // {U}; END; IF NROW(THETA0) = 0 THEN DO; THETA0 = J(NROW(C),NCOL(U),0); HITLIST = HITLIST // {THETA0}; END; IF NROW(REPN) = 0 THEN DO; REPN = 1; HITLIST = HITLIST // {REPN}; END; IF NCOL(BETASCAL) = 0 THEN DO; BETASCAL = {1}; HITLIST = HITLIST // {BETASCAL}; END; IF NCOL(SIGSCAL) = 0 THEN DO; SIGSCAL = {1}; HITLIST = HITLIST // {SIGSCAL}; END; IF NCOL(RHOSCAL) = 0 THEN DO; RHOSCAL = {1}; HITLIST = HITLIST // {RHOSCAL}; END; IF NROW(ALPHA) = 0 THEN DO; ALPHA = {.05}; HITLIST = HITLIST // {ALPHA}; END; IF NROW(ROUND) = 0 THEN DO; ROUND = 3; HITLIST = HITLIST // {ROUND}; END; IF NROW(TOLERANCE) = 0 THEN DO; TOLERANCE = 1E-12; HITLIST = HITLIST // {TOLERANCE}; END; IF NROW(UCDF) = 0 THEN DO; UCDF = {2, 2, 2, 2}; HITLIST = HITLIST // {UCDF}; END; IF NROW(UMETHOD) = 0 THEN DO; UMETHOD = {2, 2}; HITLIST = HITLIST // {UMETHOD}; END; IF NROW(MMETHOD) = 0 THEN DO; MMETHOD = {4, 2, 2}; HITLIST = HITLIST // {MMETHOD}; END; IF NROW(CLTYPE) = 0 THEN DO; CLTYPE = -1 ; HITLIST = HITLIST // {CLTYPE}; END; IF NROW(ALPHA_CL) = 0 THEN DO; ALPHA_CL =.025; HITLIST = HITLIST // {ALPHA_CL}; END; IF NROW(ALPHA_CU) = 0 THEN DO; ALPHA_CU =.025; HITLIST = HITLIST // {ALPHA_CU}; END; IF NROW(DSNAME) = 0 THEN DO; DSNAME = {WORK DODFAULT WORK}; HITLIST = HITLIST // {DSNAME}; END; ***C) Initial checks on all input matrices ***C.1) Create column vectors from user inputs; *** Check that have only vectors or scalars, not matrices; *** Check type (character or numeric); *** Check that matrices do not contain missing values; CHECKER = 0; CHECKER = CHECKER + _TYPEMISSCHK(U, "U", "N"); CHECKER = CHECKER + _TYPEMISSCHK(THETA0,"THETA0","N"); IF NCOL(REPN) > 1 THEN REPNTEMP = REPN`; ELSE REPNTEMP = REPN; CHECKER = CHECKER + _SIZECHK(REPNTEMP,"REPN") + _TYPEMISSCHK(REPNTEMP,"REPN","N"); IF NCOL(BETASCAL) > 1 THEN BETASCALTEMP = BETASCAL`; ELSE BETASCALTEMP = BETASCAL; CHECKER = CHECKER + _SIZECHK(BETASCALTEMP,"BETASCAL") + _TYPEMISSCHK(BETASCALTEMP,"BETASCAL","N"); IF NCOL(SIGSCAL) > 1 THEN SIGSCALTEMP = SIGSCAL`; ELSE SIGSCALTEMP = SIGSCAL; CHECKER = CHECKER + _SIZECHK(SIGSCALTEMP,"SIGSCAL") + _TYPEMISSCHK(SIGSCALTEMP,"SIGSCAL","N"); IF NCOL(RHOSCAL) > 1 THEN RHOSCALTEMP = RHOSCAL`; ELSE RHOSCALTEMP = RHOSCAL; CHECKER = CHECKER + _SIZECHK(RHOSCALTEMP,"RHOSCAL") + _TYPEMISSCHK(RHOSCALTEMP,"RHOSCAL","N"); IF NCOL(ALPHA) > 1 THEN ALPHATEMP = ALPHA`; ELSE ALPHATEMP = ALPHA; CHECKER = CHECKER + _SIZECHK(ALPHATEMP,"ALPHA") + _TYPEMISSCHK(ALPHATEMP,"ALPHA","N"); CHECKER = CHECKER + _TYPEMISSCHK(ROUND,"ROUND","N"); CHECKER = CHECKER + _TYPEMISSCHK(TOLERANCE,"TOLERANCE","N"); IF NCOL(UCDF) > 1 THEN UCDFTEMP = UCDF`; ELSE UCDFTEMP = UCDF; CHECKER = CHECKER + _SIZECHK(UCDFTEMP,"UCDF") + _TYPEMISSCHK(UCDFTEMP,"UCDF","N"); IF NCOL(UMETHOD) > 1 THEN UMETHODTEMP = UMETHOD`; ELSE UMETHODTEMP = UMETHOD; CHECKER = CHECKER + _SIZECHK(UMETHODTEMP,"UMETHOD") + _TYPEMISSCHK(UMETHOD,"UMETHOD","N"); IF NCOL(MMETHOD) > 1 THEN MMETHODTEMP = MMETHOD`; ELSE MMETHODTEMP = MMETHOD; *Expand 1x1 to 3x1; IF (NCOL(MMETHODTEMP)={1}) & (NROW(MMETHODTEMP)={1}) THEN MMETHODTEMP=J({3},{1},MMETHODTEMP); CHECKER = CHECKER + _SIZECHK(MMETHODTEMP,"MMETHOD") + _TYPEMISSCHK(MMETHODTEMP,"MMETHOD","N"); CHECKER = CHECKER + _TYPEMISSCHK(CLTYPE ,"CLTYPE" ,"N"); CHECKER = CHECKER + _TYPEMISSCHK(ALPHA_CL,"ALPHA_CL","N"); CHECKER = CHECKER + _TYPEMISSCHK(ALPHA_CU,"ALPHA_CU","N"); CHECKER = CHECKER + _TYPEMISSCHK(DSNAME, "DSNAME", "C"); IF CHECKER > 0 THEN GO TO ERROREND; ***D) Define initial parameters for power calculation; XPX = (ESSENCEX` * ESSENCEX + (ESSENCEX` * ESSENCEX)`) / 2; CPC = (C * C` + (C * C`)`) / 2; UPU = (U` * U + (U` * U)`) / 2; R = _RANKSYMM(XPX, "X`X"); ***Compute rank of ESSENCEX; RANKC = _RANKSYMM(CPC, "C`C"); ***Compute rank of C; RANKU = _RANKSYMM(UPU, "U`U"); ***Compute rank of U; Q = NCOL(ESSENCEX); *** Q is the number of columns of X; A = NROW(C); *** A is the number of rows of C; B = NCOL(U); *** B is the number of columns of U; S = MIN(RANKC, RANKU); *** S is the min of Rank(C) and Rank(U); P = NCOL(BETA); *** P is the # of response variables; ***E) Check and set user options; ***E.1) Check OPT_ON and OPT_OFF; CHECKER = 0; IF NROW(OPT_ON) > 0 THEN DO; IF NCOL(OPT_ON) > 1 THEN OPT_ONTEMP = OPT_ON`; ELSE OPT_ONTEMP = OPT_ON ; CHECKER = CHECKER + _SIZECHK(OPT_ONTEMP,"OPT_ON") + _TYPEMISSCHK(OPT_ONTEMP,"OPT_ON","C"); END; IF NROW(OPT_OFF) > 0 THEN DO; IF NCOL(OPT_OFF) > 1 THEN OPT_OFFTEMP = OPT_OFF`; ELSE OPT_OFFTEMP = OPT_OFF ; CHECKER = CHECKER + _SIZECHK(OPT_OFFTEMP,"OPT_OFF") + _TYPEMISSCHK(OPT_OFFTEMP,"OPT_OFF","C"); END; IF CHECKER > 0 THEN GO TO ERROREND; ***E.2) Define default options; ***OPTMATPRINT contains options for printing matrices; IF ANY(THETA0) = 1 THEN OPTMATPRINT = { 1 0 1 1 1 1 0 1 0 1 1 }`; ELSE OPTMATPRINT = { 1 0 1 1 1 1 0 1 0 0 1 }`; OPTMATPRINTLBL = {ESSENCEX RANKX BETA SIGMA RHO C RANKC U RANKU THETA0 CBETAU}`; ***OPTWARNPRINT specifies whether to print warning messages; ***Does NOT control printing error messages - they always print; OPTWARNPRINT = {1}`; OPTWARNPRINTLBL = {WARN}`; ***OPTPOWERCALC1 and OPTPOWERCALC2 contain power calculation options; OPTPOWERCALC1 = { 1 1 1 1 0 0 1 0 }`; OPTPOWERCALC1LBL = {COLLAPSE HLT PBT WLK UN HF GG BOX}`; OPTPOWERCALC2 = { 0 0 0 }`; OPTPOWERCALC2LBL = {LTFR FRACREPN ORTHU}`; ***OPTPOWERMAT1 and OPTPOWERMAT2 specify what includeD in the output matrix; OPTPOWERMAT1A = { 0 1 1 0 1 1 0 }`; OPTPOWERMAT1ALBL = {POWERCASE ALPHA SIGSCAL RHOSCAL BETASCAL TOTAL_N MAXRHOSQ}`; IF CLTYPE < 1 THEN OPTPOWERMAT1B = { 0 0 0 0 0 }`; ELSE OPTPOWERMAT1B = { 1 0 0 1 1 }`; OPTPOWERMAT1BLBL = {CLTYPE N_EST RANK_EST ALPHA_CL ALPHA_CU}`; OPTPOWERMAT1 = OPTPOWERMAT1A // OPTPOWERMAT1B; OPTPOWERMAT1LBL = OPTPOWERMAT1ALBL // OPTPOWERMAT1BLBL; OPTPOWERMAT2 = { 0 0 0 0 }`; OPTPOWERMAT2LBL = {FMETHOD UCDF UMETHOD MMETHOD}`; ***OPTPRINT specifies whether to print matrices specified in OPTMATPRINT ***and whether to print the output power matrix; OPTPRINT = {0}`; OPTPRINTLBL = {NOPRINT}`; ***OPTDATSET specifies whether to create a dataset from the output ***power matrix; OPTDATASET = {0}`; OPTDATASETLBL = {DS}`; ***E.3) Check and set user selected options; IF (( NROW(OPT_ONTEMP) > 0) | (NROW(OPT_OFFTEMP) > 0 )) THEN DO; *Check; OPTALL = OPTMATPRINTLBL // OPTWARNPRINTLBL // OPTPOWERCALC1LBL // OPTPOWERCALC2LBL // OPTPOWERMAT1LBL // OPTPOWERMAT2LBL // OPTPRINTLBL // OPTDATASETLBL; OPTSELECT = OPT_ONTEMP // OPT_OFFTEMP; OPTERR = SETDIF(OPTSELECT, OPTALL); IF NROW(OPTERR) ^= 0 THEN DO; *Option error; OPTOLDLBL = {UNI UNIHF UNIGG UNIBOX}; OPTNEWOLDLBL = {UN HF GG BOX} // OPTOLDLBL; OPTCOMPARELBL = XSECT(OPTERR, OPTOLDLBL); IF NROW(OPTCOMPARELBL) ^= 0 THEN DO; NOTE = CONCAT("ERROR 6: The names for the options ", "that specify which UNIREP test(s) ", "to perform have changed in this new ", "version. You specified the ", "option(s):"); NOTE2 = CONCAT("The new UNIREP options, ", "corresponding to the old ones are:"); PRINT NOTE[LABEL=" " ROWNAME=" " COLNAME=" "], OPTCOMPARELBL [LABEL=" " ROWNAME=" " COLNAME=" "], NOTE2 [LABEL=" " ROWNAME=" " COLNAME=" "], OPTNEWOLDLBL [LABEL = " " ROWNAME = {"NEW: " "OLD: "} COLNAME = " "]; GO TO ERROREND; END; NOTE = CONCAT("ERROR 7: Invalid options requested in ", "OPT_ON or OPT_OFF:"); PRINT NOTE[LABEL=" " ROWNAME=" " COLNAME=" "] , OPTERR [LABEL=" " ROWNAME=" " COLNAME=" "]; GO TO ERROREND; END; *Option error; IF (NROW(OPT_ONTEMP) > 0) & (NROW(OPT_OFFTEMP) > 0) THEN DO; *Compare; OPTNOTBOTH = SETDIF(OPT_ONTEMP, OPT_OFFTEMP); IF (NROW(OPTNOTBOTH) > 0) & (NCOL(OPTNOTBOTH) ^= NROW(OPT_ONTEMP)) THEN DO; OPTBOTH = SETDIF(OPT_ONTEMP, OPTNOTBOTH`); IF NROW(OPTBOTH) > 0 THEN DO; OPTBOTH = OPTBOTH`; NOTE = COMPBL(CONCAT("ERROR 8: The option(s) ", OPTBOTH, " is(are) specified in both ", "OPT_ON and OPT_OFF.")); PRINT NOTE[LABEL=" " ROWNAME=" " COLNAME=" "]; GO TO ERROREND; END; END; IF (NROW(OPTNOTBOTH) = 0) & (NROW(OPT_ONTEMP) >= 1) THEN DO; OPTBOTH = OPT_ONTEMP; NOTE = COMPBL(CONCAT("ERROR 8: The option ", OPTBOTH, " is specified in both OPT_ON ", "and OPT_OFF.")); PRINT NOTE[LABEL=" " ROWNAME=" " COLNAME=" "]; GO TO ERROREND; END; END; *Compare; END; *Check options; *Set options; IF ((NROW(OPT_ONTEMP) > 0) | (NROW(OPT_OFFTEMP) > 0)) THEN DO; CALL _SETOPT(OPTMATPRINT, OPTMATPRINTLBL, OPT_ONTEMP, OPT_OFFTEMP); CALL _SETOPT(OPTWARNPRINT, OPTWARNPRINTLBL, OPT_ONTEMP, OPT_OFFTEMP); CALL _SETOPT(OPTPOWERCALC1, OPTPOWERCALC1LBL, OPT_ONTEMP, OPT_OFFTEMP); CALL _SETOPT(OPTPOWERCALC2, OPTPOWERCALC2LBL, OPT_ONTEMP, OPT_OFFTEMP); CALL _SETOPT(OPTPOWERMAT1, OPTPOWERMAT1LBL, OPT_ONTEMP, OPT_OFFTEMP); CALL _SETOPT(OPTPOWERMAT2, OPTPOWERMAT2LBL, OPT_ONTEMP, OPT_OFFTEMP); CALL _SETOPT(OPTPRINT, OPTPRINTLBL, OPT_ONTEMP, OPT_OFFTEMP); CALL _SETOPT(OPTDATASET, OPTDATASETLBL, OPT_ONTEMP, OPT_OFFTEMP); END; ***E.4) Reset CLTYPE, N_EST, RANK_EST, ALPHA_CL, ALPHA_CU options off; *** when CLTYPE < 1; IF (CLTYPE < 1) & ANY( OPTPOWERMAT1[8:12] = 1) THEN OPTPOWERMAT1[8:12] = {0 0 0 0 0}`; ***E.5) Extract requested variable names for output; KEEPOUT = LOC(OPTPOWERMAT1); IF NROW(KEEPOUT) > 0 THEN KEEPOUTLBL = OPTPOWERMAT1LBL[KEEPOUT, ]`; ELSE IF NROW(KEEPOUTLBL) > 0 THEN FREE KEEPOUTLBL; ***E.6) Insure that at least one test chosen; IF ANY(OPTPOWERCALC1) ^= 1 THEN DO; NOTE = CONCAT("ERROR 9: OPT_OFF combined with defaults implies ", "no power calculation selected."); PRINT NOTE[LABEL=" " ROWNAME=" " COLNAME=" "]; PRINT OPT_OFF; GO TO ERROREND; END; ***E.7) Identify special cases; IF OPTPOWERCALC1[1,1] = 1 THEN DO; ***if COLLAPSE is on; IF S > 1 THEN DO; OPTPOWERCALC1[1,1] = 0; NOTE = CONCAT("WARNING 1: Rank(C*BETA*U) > 1, so ", "COLLAPSE option ignored."); IF OPTWARNPRINT THEN PRINT NOTE[LABEL=" " ROWNAME=" " COLNAME=" "]; END; IF B = 1 THEN DO; *S=1 special case; OPTPOWERCALC1= { 1 0 0 0 0 0 0 0}`; NOTE = CONCAT("WARNING 2: B = 1, so that all tests ", "coincide (univariate case). Since ", "COLLAPSE = 1, power is given as one ", "value with the heading POWER. To ", "print powers with a heading for ", "each test, specify OPT_OFF = {COLLAPSE}."); IF OPTWARNPRINT THEN PRINT NOTE[LABEL=" " ROWNAME=" " COLNAME=" "]; END; IF (B >1) & (A = 1) THEN DO; *other S=1 special case; IF MAX(OPTPOWERCALC1[{2}:{4}]) = 0 THEN DO; *No mult request; OPTPOWERCALC1[1,1] = 0; IF ANY(OPTPOWERCALC1) ^= 1 THEN DO; NOTE = CONCAT("ERROR 9: OPT_OFF combined with ", "defaults implies no power ", "calculation selected."); PRINT NOTE[LABEL=" " ROWNAME=" " COLNAME=" "] , OPT_OFF; GO TO ERROREND; END; END; ELSE DO; *At least one mult request; OPTPOWERCALC1[{2}:{4},1] = 0; NOTE = CONCAT("WARNING 3: A = 1, B > 1, so that all ", "MULTIREP tests coincide. Since ", "COLLAPSE = 1, power for all MULTIREP ", "tests is given as one value with the ", "heading POWER_MULT. To print powers ", "with a heading for each test, specify ", "OPT_OFF = {COLLAPSE};."); IF OPTWARNPRINT THEN PRINT NOTE[LABEL=" " ROWNAME=" " COLNAME=" "]; END; END; *(B >1) & (A = 1); END; *COLLAPSE is "on"; ***COLLAPSE is "off" treated in following; IF (OPTPOWERCALC1[1,1] = 0) & (S = 1) & OPTWARNPRINT THEN DO; IF (B = 1) THEN DO; NOTE = CONCAT("WARNING 4: B = 1, so that all tests ", "coincide (univariate case). To print ", "powers with one overall heading of POWER ", "(thereby reducing the amount of output) ", "remove OPT_OFF = {COLLAPSE}."); PRINT NOTE[LABEL=" " ROWNAME=" " COLNAME=" "]; END; IF (B > 1) THEN DO; NOTE = CONCAT("WARNING 5: A = 1, B > 1, so that all ", "MULTIREP tests coincide. To print powers ", "with one overall heading of POWER_MULT ", "(thereby reducing the amount of output) ", "remove OPT_OFF = {COLLAPSE}."); PRINT NOTE[LABEL=" " ROWNAME=" " COLNAME=" "]; END; END; ***F) Check input matrices; ***F.1) Check REPN; IF MIN(REPNTEMP) <= TOLERANCE THEN DO; NOTE = "ERROR 10: All REPN values must be >= TOLERANCE > 0."; PRINT NOTE[LABEL=" " ROWNAME=" " COLNAME=" "] , REPN; GO TO ERROREND; END; IF (OPTPOWERCALC2[2,1] = 0) & (MAX(INT(REPNTEMP) ^= REPNTEMP)= 1) THEN DO; NOTE = CONCAT("ERROR 11: All REPN values must be positive ", "integers. To allow fractional REPN values, specify ", "OPT_ON = {FRACREPN}."); PRINT NOTE[LABEL=" " ROWNAME=" " COLNAME=" "] , REPN; GO TO ERROREND; END; ***F.2) Check SIGSCAL; IF MIN(SIGSCALTEMP) <= TOLERANCE THEN DO; NOTE = "ERROR 12: All SIGSCAL values must be > TOLERANCE > 0."; PRINT NOTE[LABEL=" " ROWNAME=" " COLNAME=" "] , SIGSCAL , TOLERANCE; GO TO ERROREND; END; ***F.3) Check ALPHA; IF (MIN(ALPHATEMP) <= TOLERANCE) | (MAX(ALPHATEMP) >= 1) THEN DO; NOTE = CONCAT("ERROR 13: All ALPHA values must be > TOLERANCE > 0 ", "and < 1. "); PRINT NOTE[LABEL=" " ROWNAME=" " COLNAME=" "] , ALPHA , TOLERANCE; GO TO ERROREND; END; ***F.4) Check ROUND and set ROUNDOFF; IF MAX(NCOL(ROUND),NROW(ROUND)) > 1 THEN DO; NOTE = "ERROR 14: ROUND must be a scalar (1x1)."; PRINT NOTE[LABEL=" " ROWNAME=" " COLNAME=" "] , ROUND; GO TO ERROREND; END; IF (ROUND < 1) | (ROUND > 15) THEN DO; NOTE = " ERROR 15: User specified ROUND < 1 OR ROUND >15"; PRINT NOTE[LABEL=" " ROWNAME=" " COLNAME=" "] , ROUND; GO TO ERROREND; END; ROUNDOFF = 1 / 10 ** ROUND; ***F.5) Check TOLERANCE; ***Check that TOLERANC doesn't exist (changed the name from TOLERANC to TOLERANCE with this new verion); IF NROW(TOLERANC) > 0 THEN DO; NOTE = CONCAT("ERROR 16: With this version, numeric zero is ", "specified with the scalar TOLERANCE not TOLERANC ", "as in previous versions. You have entered the ", "scalar TOLERANC. Please rename this to TOLERANCE."); PRINT NOTE[LABEL=" " ROWNAME=" " COLNAME=" "] , TOLERANC; GO TO ERROREND; END; IF TOLERANCE <= 0 THEN DO; NOTE = "ERROR 17: User specified TOLERANCE <= zero."; PRINT NOTE[LABEL=" " ROWNAME=" " COLNAME=" "] , TOLERANCE; GO TO ERROREND; END; IF MAX(NCOL(TOLERANCE), NROW(TOLERANCE)) > 1 THEN DO; NOTE = "ERROR 18: TOLERANCE must be a scalar (1x1)."; PRINT NOTE[LABEL=" " ROWNAME=" " COLNAME=" "] , TOLERANCE; GO TO ERROREND; END; IF TOLERANCE >= .01 THEN DO; NOTE = CONCAT("WARNING 6: User specified TOLERANCE > 0.01. This is ", "the value assumed to be numeric zero and affects ", "many calculations. Please check that this value ", "is correct."); IF OPTWARNPRINT THEN DO; PRINT NOTE[LABEL=" " ROWNAME=" " COLNAME=" "] , TOLERANCE; END; END; ***F.6) Check UCDF; IF (NROW(UCDFTEMP) ^= 1) & (NROW(UCDFTEMP) ^= 4) THEN DO; NOTE = "ERROR 19: Matrix UCDF must be 1x1, 4x1, or 1x4"; PRINT NOTE[LABEL=" " ROWNAME=" " COLNAME=" "] , UCDF; GO TO ERROREND; END; IF NROW(SETDIF(UCDFTEMP,{1 2 3 4})) > 0 THEN DO; NOTE = "ERROR 20: Values in matrix UCDF must be 1, 2, 3, or 4"; PRINT NOTE[LABEL=" " ROWNAME=" " COLNAME=" "] , UCDF; GO TO ERROREND; END; IF NROW(UCDFTEMP) = 1 THEN UCDFTEMP = UCDF // UCDF // UCDF // UCDF; IF (ANY(UCDFTEMP = 1) > 0) & OPTWARNPRINT THEN DO; NOTE = CONCAT("WARNING 7: You have chosen the Muller, Barton ", "(1989) approximation for the UNIREP statistic CDF. ", "Muller, Edwards, Taylor (2004) found NO condition ", "where their approximation was not superior to this ", "Muller, Barton approximation. Suggest specifying ", "UCDF={2}; unless you are performing a backwards ", "comparison calculation."); PRINT NOTE[LABEL=" " ROWNAME=" " COLNAME=" "] , UCDF; END; ***F.7) Check UMETHOD; IF (NROW(UMETHODTEMP) ^=1) & (NROW(UMETHODTEMP) ^= 2) THEN DO; NOTE = "ERROR 21: Matrix UMETHOD must be 1x1, 2x1, or 1x2"; PRINT NOTE[LABEL=" " ROWNAME=" " COLNAME=" "] , UMETHOD; GO TO ERROREND; END; IF (NROW(SETDIF(UMETHODTEMP,{1 2})) > 0) THEN DO; NOTE = "ERROR 22: Values in matrix UMETHOD must be 1 or 2"; PRINT NOTE[LABEL=" " ROWNAME=" " COLNAME=" "] , UMETHOD; GO TO ERROREND; END; IF NROW(UMETHODTEMP) = 1 THEN UMETHODTEMP = UMETHOD // UMETHOD; ***F.8) Check MMETHOD; IF (NROW(MMETHODTEMP) ^= 1) & (NROW(MMETHODTEMP) ^= 3) THEN DO; NOTE = "ERROR 23: Matrix UMETHOD must be 1x1, 3x1, or 1x3"; PRINT NOTE[LABEL=" " ROWNAME=" " COLNAME=" "] , MMETHOD; GO TO ERROREND; END; IF (NROW(SETDIF(MMETHODTEMP, {1 2 3 4})) > 0) THEN DO; NOTE = CONCAT("ERROR 24: Values in matrix MMETHOD ", "specifying options for HLT, PBT and WLK ", "tests must be 1, 2, 3, or 4."); PRINT NOTE[LABEL=" " ROWNAME=" " COLNAME=" "] , MMETHOD; GO TO ERROREND; END; ***F.9) Check DSNAME; IF (NROW(DSNAME) > 1) | (NCOL(DSNAME) > 3) THEN DO; NOTE = "ERROR 25: DSNAME must have only 1 row and 2-3 columns."; PRINT NOTE[LABEL=" " ROWNAME=" " COLNAME=" "]; GO TO ERROREND; END; ***F.10) Check CLTYPE, N_EST, RANK_EST, ALPHA_CL, ALPHA_CU; IF NROW(LIMCASE) > 0 THEN DO; NOTE = CONCAT("WARNING 8: In POWERLIM.IML, ", "the choice of type of confidence interval was ", "specified using the vector LIMCASE. In this ", "version, two of the old LIMCASE options have ", "been deleted, the remaining options renumbered, and ", "instead of LIMCASE, the vector CLTYPE now ", "specifies choice of type of confidence interval. ", "The following table displays the renumbering ", "scheme used in CLTYPE: "); TEMP = CONCAT("VALUE ", " LIMCASE ", " CLTYPE") // " " // CONCAT("0 ", "BETA known, SIGMA known ", "NOT VALID") // CONCAT("1 ", "BETA estimated, SIGMA known ", "BETA known, SIGMA estimated") // CONCAT("2 ", "BETA known, SIGMA estimated ", "BETA estimated, SIGMA estimated") // CONCAT("3 ", "BETA estimated, SIGMA estimated ", "NOT VALID"); IF OPTWARNPRINT THEN DO; PRINT NOTE[LABEL=" " ROWNAME=" " COLNAME=" "] , TEMP [LABEL = " " ROWNAME = " " COLNAME = " "]; END; END; IF (NROW(CLTYPE) + NCOL(CLTYPE)) ^= 2 THEN DO; NOTE = "ERROR 26: CLTYPE must be 1 by 1 matrix" ; PRINT NOTE[LABEL=" " ROWNAME=" " COLNAME=" "] , CLTYPE; GO TO ERROREND; END; IF CLTYPE >= 1 THEN DO; IF (1 - ANY(CLTYPE = {1 2})) THEN DO; NOTE = CONCAT("ERROR 27: Invalid CLTYPE. If ", "CLTYPE >= 1, it must be 1 or 2."); PRINT NOTE[LABEL=" " ROWNAME=" " COLNAME=" "]; GO TO ERROREND; END; IF (NROW(N_EST) = 0) | (NROW(RANK_EST) = 0) THEN DO; IF (NROW(N_EST) = 0) THEN DO; NOTE = CONCAT("ERROR 28: N_EST not supplied. ", "This is required whenever ", "CLTYPE >= 1."); PRINT NOTE[LABEL=" " ROWNAME=" " COLNAME=" "]; END; IF (NROW(RANK_EST) = 0) THEN DO; NOTE = CONCAT("ERROR 29: RANK_EST not ", "supplied. This is required whenever ", "CLTYPE >= 1."); PRINT NOTE[LABEL=" " ROWNAME=" " COLNAME=" "]; END; GO TO ERROREND; END; IF (_TYPEMISSCHK(N_EST,"N_EST","N") + _TYPEMISSCHK(RANK_EST,"RANK_EST","N") ) > 0 THEN GO TO ERROREND; IF (NROW(N_EST) + NCOL(N_EST)) ^= 2 THEN DO; NOTE = CONCAT("ERROR 30: A 1x1 matrix N_EST is required if ", "CLTYPE >= 1."); PRINT NOTE[LABEL=" " ROWNAME=" " COLNAME=" "] , N_EST; GO TO ERROREND; END; IF (NROW(RANK_EST) + NCOL(RANK_EST)) ^= 2 THEN DO; NOTE = CONCAT("ERROR 31: A 1x1 matrix RANK_EST is required ", "if CLTYPE >= 1."); PRINT NOTE[LABEL=" " ROWNAME=" " COLNAME=" "] , RANK_EST; GO TO ERROREND; END; DFEE = N_EST - RANK_EST; IF (DFEE <= 0) THEN DO; NOTE = "ERROR 32: (N_EST - RANK_EST) must be > 0"; PRINT NOTE[LABEL=" " ROWNAME=" " COLNAME=" "] , N_EST , RANK_EST; GO TO ERROREND; END; IF (NROW(ALPHA_CL) + NCOL(ALPHA_CL)) > 2 THEN DO; NOTE = "ERROR 33: ALPHA_CL must be a 1x1 matrix" ; PRINT NOTE[LABEL=" " ROWNAME=" " COLNAME=" "]; GO TO ERROREND; END; IF (NROW(ALPHA_CU) + NCOL(ALPHA_CU)) > 2 THEN DO; NOTE = "ERROR 34: ALPHA_CU must be a 1x1 matrix" ; PRINT NOTE[LABEL=" " ROWNAME=" " COLNAME=" "]; GO TO ERROREND; END; IF (ALPHA_CL < 0) | (ALPHA_CU < 0) | (ALPHA_CL >= 1) | (ALPHA_CU >= 1) | ((ALPHA_CL + ALPHA_CU) >= 1) THEN DO; NOTE = CONCAT("ERROR 35: ALPHA_CL and ALPHA_CU must both ", "be >= 0 and <= 1."); PRINT NOTE[LABEL=" " ROWNAME=" " COLNAME=" "] , ALPHA_CL , ALPHA_CU; GO TO ERROREND; END; END; ***G) Checks on matrices to prepare for power calculation; ***G.1) Checks on size and conformity of matrices; CHECKER = 0; IF MIN(VECDIAG(SIGMA)) <= TOLERANCE THEN DO; CHECKER = 1; NOTE = "ERROR 36: At least one variance <= TOLERANCE. Check SIGMA."; PRINT NOTE[LABEL=" " ROWNAME=" " COLNAME=" "] , TOLERANCE , SIGMA ; END; IF (NCOL(SIGMA) = 1) & OPTWARNPRINT THEN DO; NOTE = CONCAT("WARNING 10: SIGMA is a scalar. For this program, ", "a scalar SIGMA must equal the variance, NOT the ", "standard deviation."); PRINT NOTE[LABEL=" " ROWNAME=" " COLNAME=" "]; END; IF (NROW(U) ^= NROW(SIGMA)) THEN DO; CHECKER = 1; NOTE = CONCAT("ERROR 37: The number of rows of U is not equal to ", "number of rows of SIGMA."); PRINT NOTE[LABEL=" " ROWNAME=" " COLNAME=" "] , U , SIGMA; END; IF (NCOL(C) ^= NROW(BETA)) THEN DO; CHECKER = 1; NOTE = CONCAT("ERROR 38: The number of columns of C is not equal ", "to number of rows of BETA."); PRINT NOTE[LABEL=" " ROWNAME=" " COLNAME=" "] , C , BETA; END; IF (NCOL(BETA) ^= NROW(U)) THEN DO; CHECKER = 1; NOTE = CONCAT("ERROR 39: The number of columns of BETA is not ", "equal to number of rows of U."); PRINT NOTE[LABEL=" " ROWNAME=" " COLNAME=" "] , BETA , U; END; IF (Q ^= NROW(BETA)) THEN DO; CHECKER = 1; NOTE = CONCAT("ERROR 40: The number of columns of ESSENCEX is not ", "equal to number of rows of BETA."); PRINT NOTE[LABEL=" " ROWNAME=" " COLNAME=" "] , ESSENCEX , BETA; END; IF (A > NCOL(C)) THEN DO; CHECKER = 1; NOTE = CONCAT("ERROR 41: The number of rows of C is greater than ", "the number of columns of C."); PRINT NOTE[LABEL=" " ROWNAME=" " COLNAME=" "] , C; END; IF (B > NROW(U)) THEN DO; CHECKER = 1; NOTE = CONCAT("ERROR 42: The number of columns of U is greater ", "than the number of rows of U."); PRINT NOTE[LABEL=" " ROWNAME=" " COLNAME=" "] , U; END; IF (NROW(THETA0) ^= A) | (NCOL(THETA0) ^= B) THEN DO; CHECKER = 1; NOTE = "ERROR 43: The THETA0 matrix does not conform to CBETAU."; CBETAU = C * BETA * U; PRINT NOTE[LABEL=" " ROWNAME=" " COLNAME=" "] , THETA0 , CBETAU; END; IF MIN(SIGSCAL) <= TOLERANCE THEN DO; CHECKER = 1; NOTE = "ERROR 44: Smallest value in SIGSCAL <= TOLERANCE (too small)"; PRINT NOTE[LABEL=" " ROWNAME=" " COLNAME=" "] , TOLERANCE , SIGSCAL; END; IF CHECKER = 1 THEN GO TO ERROREND; ***G.2) Define more matrices needed for power calculations and, when needed, their ranks; XPXGINV = GINV(XPX); M = C * XPXGINV * C`; RANKM = _RANKSYMM(M, "C * GINV(X`X) * C`"); ***Compute rank of M; CBETAU = C * BETA * U; SD = DIAG(SQRT(VECDIAG(SIGMA))); RHO = INV(SD) * SIGMA * INV(SD); *** Correlation matrix to be printed; ***G.3) Check and warning for less than full rank ESSENCEX; IF (R ^= Q) & (OPTPOWERCALC2[1,1] ^= 1) THEN DO; NOTE = CONCAT("ERROR 45: You have specified a less than full rank ", "ESSENCEX matrix. Specifying OPT_ON={LTFR}; before ", "RUN POWER; will allow calculations to complete in ", "this case. In doing so, you accept responsibility ", "for the results being the ones desired."); PRINT NOTE[LABEL=" " ROWNAME=" " COLNAME=" "] , R [LABEL = "RANK OF X"] Q [LABEL = "NUMBER OF COLUMNS IN X"] , ESSENCEX; GO TO ERROREND; END; IF (R ^= Q) & (OPTPOWERCALC2[1,1] = 1) THEN DO; NOTE = CONCAT("WARNING 11: You are performing calculations using ", "a less than full rank ESSENCEX matrix. In doing ", "so, you accept responsibility for the results ", "being the ones desired."); IF OPTWARNPRINT THEN PRINT NOTE[LABEL=" " ROWNAME=" " COLNAME=" "] , ESSENCEX; END; ***G.4) Checks for testable GLH and estimable CBETAU; ***Per Lemma 15.1 part b, BIOS 262 notes, if X LTFR then the following are ; ***required for testability: ; *** Full rank C (Rank of C = A) ; *** Full rank M (Rank of M = A) ; *** Full rank U (Rank of U = B) ; *** Check C = C * XPXGINV * XPX ; IF RANKC ^= A THEN DO; NOTE = COMPBL(CONCAT("ERROR 46: C must be full rank to ensure ", "testability of CBETAU = THETA0. Your C matrix has ", CHAR(A,1), " rows, but has rank ", CHAR(RANKC), ".")); PRINT NOTE[LABEL=" " ROWNAME=" " COLNAME=" "] , U ; GO TO ERROREND; END; IF RANKU ^= B THEN DO; NOTE = COMPBL(CONCAT("ERROR 47: U must be full rank to ensure ", "testability of CBETAU = THETA0. Your U matrix has ", CHAR(B,1), " columns, but has rank ", CHAR(RANKU), ".")); PRINT NOTE[LABEL=" " ROWNAME=" " COLNAME=" "] , U ; GO TO ERROREND; END; IF MAX(ABS(C - C*XPXGINV*XPX)) > TOLERANCE THEN DO; NOTE = CONCAT("ERROR 48: A requirement for estimability is that ", "C = C * GINV(X`X) * X`X. Your choices of C and ", "ESSENCEX do not provide this. We suggest using ", "full rank coding for ESSENCEX."); PRINT NOTE[LABEL=" " ROWNAME=" " COLNAME=" "] , C , ESSENCEX; GO TO ERROREND; END; IF RANKM ^= A THEN DO; NOTE = COMPBL(CONCAT("ERROR 49: M = C(GINV(X`X))C` must be full ", "rank to ensure testability of CBETAU = ", "THETA0. Your M matrix is ", CHAR(A,1), "x", CHAR(A,1), "but has rank ", CHAR(RANKM), ".")); PRINT NOTE[LABEL=" " ROWNAME=" " COLNAME=" "] , M ; GO TO ERROREND; END; ***G.5) Create Orthonormal U for UNIREP test power calculations; UTEMP = U; IF (OPTPOWERCALC1[5,1] = 1) | (OPTPOWERCALC1[6,1] = 1) | (OPTPOWERCALC1[7,1] = 1) | (OPTPOWERCALC1[8,1] = 1) THEN DO; IF UPU[1 ,1] ^= 0 THEN UPU = UPU / UPU[1 ,1]; UDIF = ABS( UPU - I(B)); IF ( MAX(UDIF) > SQRT(TOLERANCE) ) | ( (B > 1) * ( MAX( U`*J(P,1,1) ) > SQRT(TOLERANCE) ) > 0 ) THEN DO; IF OPTPOWERCALC2[3,1] = 0 THEN DO; NOTE = CONCAT("ERROR 50: For univariate repeated ", "measures, U must be proportional ", "to an orthonormal matrix [U`U = cI] ", "and orthogonal to a Px1 column ", "of 1's [U`1 = 0]. The U matrix ", "specified does not have these ", "properties. To have this program ", "provide a U matrix with the required", "properties, specify OPT_ON= {ORTHU}; ", " If you do not wish to compute power", " for UNIREP tests, specify OPT_OFF = ", "{GG HF UN BOX};"); PRINT NOTE[LABEL=" " ROWNAME=" " COLNAME=" "] ,U ; GO TO ERROREND; END; IF OPTPOWERCALC2[3,1] = 1 THEN DO; CALL GSORTH(UORTH, T, FLAG, U); UTEMP = UORTH; IF ( (B > 1) * ( MAX( UORTH`*J(P,1,1) ) > SQRT(TOLERANCE) ) > 0 ) THEN DO; NOTE = CONCAT("ERROR 51: You have ", "specified option ORTHU so ", "that the program will ", "provide a U ", "that is proportional to ", "to an orthonormal matrix ", "and orthogonal to a Px1 ", "column of 1's. ", "The original U given cannot ", "be made to be orthogonal ", "to a Px1 column of 1's. ", "Choose a different U ", "matrix."); PRINT NOTE[LABEL=" " ROWNAME=" " COLNAME=" "] , U; GO TO ERROREND; END; CBETAU = C * BETA * UTEMP; NOTE = CONCAT("WARNING 12: For univariate repeated ", "measures, U must be proportional ", "to an orthonormal matrix [U`U = cI] ", "and orthogonal to a Px1 column ", "of 1's [U`1 = 0]. The U matrix ", "specified does not have these ", "properties. A new U matrix with ", "these properties was created from ", "your input and ", "used in UNIREP calculations:"); IF OPTWARNPRINT THEN DO; PRINT NOTE[LABEL=" " ROWNAME=" " COLNAME=" "] , UTEMP[LABEL="U" ROWNAME=" "]; END; END; END; END; ***H) Loop to select ALPHA, BETA, SIGMA, N, and compute power; POWERCASE = 0; ***H.1) Select alpha; DO ALPHACNT = 1 TO NROW(ALPHATEMP); ALPHASCALAR = ALPHATEMP[ALPHACNT,1]; ***H.2) select sigma; DO SIGMACNT = 1 TO NROW(SIGSCALTEMP); SIGMASCALAR = SIGSCALTEMP[SIGMACNT,1]; *** Compute sigscaled covariance matrix; SIGMASCALED1 = SIGMA # SIGMASCALAR; *** Compute rho from sigscaled covariance matrix; VARSCALED = VECDIAG(SIGMASCALED1); IF MIN(VARSCALED) <= TOLERANCE THEN DO; NOTE = COMPBL(CONCAT("ERROR 52: The covariance matrix formed with ", "SIGSCAL element ", CHAR(SIGMACNT,1), " has a variance <= TOLERANCE (too small).")); PRINT NOTE[LABEL=" " ROWNAME=" " COLNAME=" "] , SIGSCAL; GO TO ERROREND; END; STDEVSCALED = SQRT(VARSCALED); RHOSCALED1 = (DIAG(J(P,1,1) / (STDEVSCALED))) * SIGMASCALED1 * (DIAG(J(P,1,1) / (STDEVSCALED))); ***H.3) Create new rhos; DO RHOCNT = 1 TO NROW(RHOSCALTEMP); RHOSCALAR = RHOSCALTEMP[RHOCNT,1]; RHOJUNK = (RHOSCALED1 # RHOSCALAR); *** Diagonal elements not =1; RHOOFFDIAG = RHOJUNK - DIAG(VECDIAG(RHOJUNK)); *** Off-diagonals; RHOOFFDIAG = (RHOOFFDIAG + RHOOFFDIAG`) / 2; ***To insure symmetry; IF MAX(ABS(RHOOFFDIAG)) > 1 THEN DO; NOTE = COMPBL(CONCAT("ERROR 53: SIGSCAL = ", CHAR(SIGMASCALAR,1), " and RHOSCAL = ", CHAR(RHOSCALAR,1), " produce a correlation with absolute ", "value > 1 .")); PRINT NOTE[LABEL=" " ROWNAME=" " COLNAME=" "] , SIGSCAL , RHOSCAL; GO TO ERROREND; END; IF MAX(ABS(RHOOFFDIAG)) = 1 THEN DO; NOTE = COMPBL(CONCAT("WARNING 13: SIGSCAL = ", CHAR(SIGMASCALAR,1), " and RHOSCAL = ", CHAR(RHOSCALAR,1), " produce a correlation with absolute ", "value = 1 .")); POWERWARN[20,1] = POWERWARN[20,1] + 1; IF OPTWARNPRINT THEN PRINT NOTE[LABEL=" " ROWNAME=" " COLNAME=" "] , SIGSCAL , RHOSCAL; GO TO ERROREND; END; RHOSCALED2 = RHOOFFDIAG + I(P); *** Create new sigma from rho; SIGMASCALED2 = DIAG(STDEVSCALED) * RHOSCALED2 * DIAG(STDEVSCALED); SIGMASTAR = UTEMP` * SIGMASCALED2 * UTEMP; SIGMASTAR = (SIGMASTAR + SIGMASTAR`) / 2; ***To insure symmetry; CALL EIGEN(SIGMASTAREVAL, SIGMASTAREVEC, SIGMASTAR); SIGMASTARLBL = {"U`*SIGMA*U"}; RANKSIGMASTAR = _RANKSYMM(SIGMASTAR, SIGMASTARLBL); IF RANKSIGMASTAR ^= B THEN DO; NOTE = COMPBL(CONCAT("ERROR 54: SIGMASTAR = U`*SIGMA*U must be ", "full rank to ensure testability of CBETAU = ", "THETA0. SIGMASTAR is", CHAR(B,1), "x", CHAR(B,1), "but has rank ", CHAR(RANKSIGMASTAR), ".")); PRINT NOTE[LABEL=" " ROWNAME=" " COLNAME=" "] , SIGMASTAR; GO TO ERROREND; END; ***H.4) Select Beta; DO BETACNT = 1 TO NROW(BETASCALTEMP); BETASCALAR = BETASCALTEMP[BETACNT,1]; BETASCALED = BETA # BETASCALAR; THETA = C * BETASCALED * UTEMP; ESSH = (THETA - THETA0)` * SOLVE(M, I(NROW(M))) * (THETA - THETA0); ***H.5) Select N; DO NCNT = 1 TO NROW(REPNTEMP); POWERCASE = POWERCASE + 1; N = REPNTEMP[NCNT,1] # NROW(ESSENCEX); IF (INT(N) ^= N) & OPTWARNPRINT THEN DO; NOTE = CONCAT("WARNING 14: You have specified a design ", "with non-integer N. In doing so, you are ", "responsible for correct interpretation of ", "the results."); PRINT NOTE[LABEL=" " ROWNAME=" " COLNAME=" "]; END; E = SIGMASTAR # (N - R); ***Not used if (N - R) <= 0 ; H = REPNTEMP[NCNT,1] # ESSH; ***H.6) Eigenvalues for H*INV(E); EPS=.; IF (N - R) <= 0 THEN DO; EVAL = J(S,1,.); RHOSQ = . ; D = 1; MTP = B; ***EPS = .; DEIGVAL = .; SLAM1 = .; SLAM2 = .; SLAM3 = .; CANCORRMAX=.; END; ELSE DO; FINV = INV(HALF(E))`; HEIORTH = FINV * H * FINV`; HEIORTH = (HEIORTH + HEIORTH`) / 2; *** Insure symmetry; EVAL = EIGVAL(HEIORTH); EVAL = EVAL[1:S,1]; *** At most S nonzero; EVAL = EVAL # (EVAL > TOLERANCE); *** Set evals= 1 THEN OUTPUT = OUTPUT || CLTYPE || N_EST || RANK_EST || ALPHA_CL || ALPHA_CU; IF NROW(KEEPOUT) = 0 THEN DO; FREE OUTPUT; IF NROW(HOLDPOWERLBL) > 0 THEN FREE HOLDPOWERLBL; END; ELSE DO; OUTPUT = OUTPUT[1, KEEPOUT]; HOLDPOWERLBL = KEEPOUTLBL; END; ***H.8) Compute requested power calculations; ***H.8.a) Start UNIREP power calcs, if any; IF MAX( OPTPOWERCALC1[ {5}:{8} ] ) THEN DO; CALL _FIRSTUNI(D, MTP, EPS, DEIGVAL, SLAM1, SLAM2, SLAM3, SIGMASTAR, B); *** Find first requested UNIREP test so include epsilon; FIRSTUNI= {4} + MIN( LOC( OPTPOWERCALC1[ {5}:{8} ] ) ) ; END; ***H.8.b) Select requested power calculations; *** POWERCALC = 1 --> Univariate test or all mult tests coincide; *** POWERCALC = 2 --> HL Trace; *** POWERCALC = 3 --> PB Trace; *** POWERCALC = 4 --> Wilks LR; *** POWERCALC = 5 --> UNIREP Uncorrected; *** POWERCALC = 6 --> UNIREP Huynh-Feldt; *** POWERCALC = 7 --> UNIREP Geisser-Greenhouse; *** POWERCALC = 8 --> UNIREP Box Conservative; POWERCALCNUM = (LOC(OPTPOWERCALC1))`; DO POWERCNT = 1 TO NROW(POWERCALCNUM); POWERCALC = POWERCALCNUM[POWERCNT,1]; IF POWERCALC=1 THEN CALL _SPECIAL(POWER, LBL, POWERWARN, A, B, S, N, R, EVAL, ALPHASCALAR, OPTPOWERMAT2); IF POWERCALC = 2 THEN CALL _HLT(POWER, LBL, POWERWARN, A, B, S, N, R, EVAL, ALPHASCALAR, MMETHODTEMP, OPTPOWERMAT2); IF POWERCALC = 3 THEN CALL _PBT(POWER, LBL, POWERWARN, A, B, S, N, R, EVAL, ALPHASCALAR, MMETHODTEMP, OPTPOWERMAT2); IF POWERCALC = 4 THEN CALL _WLK(POWER, LBL, POWERWARN, A, B, S, N, R, EVAL, ALPHASCALAR, MMETHODTEMP, OPTPOWERMAT2); IF POWERCALC = 5 THEN DO; EXEPS = 1; CALL _LASTUNI(POWER, LBL, POWERWARN, EXEPS, H, E, A, B, R, N, EPS, POWERCALC, OPTPOWERMAT2, ALPHASCALAR, FIRSTUNI, SIGMASTAREVAL, SIGMASTAREVEC, UCDFTEMP, UMETHODTEMP); END; IF POWERCALC = 6 THEN DO; IF (N - R) <= 0 THEN EXEPS=.; ELSE EXEPS = _HFEXEPS(B, N, R, D, MTP, EPS, DEIGVAL, SLAM1, SLAM2, SLAM3, SIGMASTAREVAL, SIGMASTAREVEC, UMETHODTEMP); CALL _LASTUNI(POWER, LBL, POWERWARN, EXEPS, H, E, A, B, R, N, EPS, POWERCALC, OPTPOWERMAT2, ALPHASCALAR, FIRSTUNI, SIGMASTAREVAL, SIGMASTAREVEC, UCDFTEMP, UMETHODTEMP); END; IF POWERCALC = 7 THEN DO; IF (N - R) <= 0 THEN EXEPS=.; ELSE EXEPS = _GGEXEPS(B, N, R, D, MTP, EPS, DEIGVAL, SLAM1, SLAM2, SLAM3, SIGMASTAREVAL, SIGMASTAREVEC, UMETHODTEMP); CALL _LASTUNI(POWER, LBL, POWERWARN, EXEPS, H, E, A, B, R, N, EPS, POWERCALC, OPTPOWERMAT2, ALPHASCALAR, FIRSTUNI, SIGMASTAREVAL, SIGMASTAREVEC, UCDFTEMP, UMETHODTEMP); END; IF POWERCALC = 8 THEN DO; IF (N - R) <= 0 THEN EXEPS = .; ELSE EXEPS = 1 / B; CALL _LASTUNI(POWER, LBL, POWERWARN, EXEPS, H, E, A, B, R, N, EPS, POWERCALC, OPTPOWERMAT2, ALPHASCALAR, FIRSTUNI, SIGMASTAREVAL, SIGMASTAREVEC, UCDFTEMP, UMETHODTEMP); END; ***H.8.c) Append all requested power calculations in a row; OUTPUT = OUTPUT || POWER; HOLDPOWERLBL = HOLDPOWERLBL || LBL; ***Should only do this once; POWEREQ1 = ((OUTPUT # ( (SUBSTR(HOLDPOWERLBL,1,5) = "POWER") & (SUBSTR(HOLDPOWERLBL,1,6) ^= "POWERC"))) = 1); OUTPUT = ROUND(OUTPUT, ROUNDOFF); POWERROUNDEQ1 = ((OUTPUT # ( (SUBSTR(HOLDPOWERLBL,1,5) = "POWER") & (SUBSTR(HOLDPOWERLBL,1,6) ^= "POWERC"))) = 1); POWERWARN[22,1] = POWERWARN[22,1] + SUM(POWERROUNDEQ1) - SUM(POWEREQ1); FREE POWER LBL; END; ***Powercalc; ***H.8.d) Stack calculations; HOLDPOWER = HOLDPOWER // OUTPUT; ***F.9) Check for sufficient sample size; IF ((N - R) <= 5) THEN DO; IF ((B > 1) & (ANY(OPTPOWERCALC1[{5 8},1]) = 1) & (EPS ^= 1)) | ((S > 1) & (ANY(OPTPOWERCALC1[6:7,1]) = 1)) | ((S > 1) & (ANY(OPTPOWERCALC1[2:4,1]) = 1)) THEN POWERWARN[21,1] = POWERWARN[21,1] + 1; END; END; *** Repn; END; *** Beta; END; *** Rho; END; *** Sigma; END; *** Alpha; ***I) Print warnings from loop; IF (NROW(LOC(POWERWARN > 0)) > 0) & OPTWARNPRINT THEN DO; WARNNONZERO = POWERWARN[LOC(POWERWARN > 0)]; WARNNONZEROLBL = POWERWARNLBL[LOC(POWERWARN > 0)]; WARNINGS = COMPBL(CONCAT("POWER CALCULATION WARNING ", TRIM(WARNNONZEROLBL), " This occurred ", char(warnnonzero,2), " time(s).")); DO I = 1 TO NROW(WARNINGS); PRINT (WARNINGS[I,]) [COLNAME = " " ROWNAME=" " LABEL=" "]; END; END; ***J) Printed Output; ***J.1) Print requested matrices; IF (OPTPRINT = 0) & ANY(OPTMATPRINT) THEN DO; PRINT /; PRINT "****REQUESTED PRINTED MATRICES***" ; IF OPTMATPRINT[ 1,1] = 1 THEN PRINT ESSENCEX; IF OPTMATPRINT[ 2,1] = 1 THEN PRINT R[LABEL = "RANK OF X"]; IF OPTMATPRINT[ 3,1] = 1 THEN PRINT BETA; IF OPTMATPRINT[ 4,1] = 1 THEN PRINT SIGMA; IF OPTMATPRINT[ 5,1] = 1 THEN PRINT RHO; IF OPTMATPRINT[ 6,1] = 1 THEN PRINT C; IF OPTMATPRINT[ 7,1] = 1 THEN PRINT RANKC[LABEL = "RANK OF C"]; IF OPTMATPRINT[ 8,1] = 1 THEN PRINT U[LABEL = "U"]; IF OPTMATPRINT[ 9,1] = 1 THEN PRINT B[LABEL = "RANK OF U"]; IF OPTMATPRINT[10,1] = 1 THEN PRINT THETA0; IF OPTMATPRINT[11,1] = 1 THEN PRINT CBETAU[LABEL="C * BETA * U"]; END; ***J.2) Print power calculations; ***The following 8 statements deal with a user specifying; ***ROUND greater than the default field width; IF ROUND = 8 THEN MATTRIB HOLDPOWER FORMAT=BEST10.; IF ROUND = 9 THEN MATTRIB HOLDPOWER FORMAT=BEST11.; IF ROUND = 10 THEN MATTRIB HOLDPOWER FORMAT=BEST12.; IF ROUND = 11 THEN MATTRIB HOLDPOWER FORMAT=BEST13.; IF ROUND = 12 THEN MATTRIB HOLDPOWER FORMAT=BEST14.; IF ROUND = 13 THEN MATTRIB HOLDPOWER FORMAT=BEST15.; IF ROUND = 14 THEN MATTRIB HOLDPOWER FORMAT=BEST16.; IF ROUND = 15 THEN MATTRIB HOLDPOWER FORMAT=BEST17.; POWERCASETOT = NROW(HOLDPOWER); IF OPTPRINT = 0 THEN DO; IF POWERCASETOT <= 40 THEN PRINT / "****POWER CALCULATIONS****" , HOLDPOWER[ROWNAME=" " COLNAME=HOLDPOWERLBL]; ELSE DO; BRKPT = DO(36, POWERCASETOT, 36); IF MAX(BRKPT) ^= POWERCASETOT THEN BRKPT = BRKPT || POWERCASETOT; START = 1; DO CNT = 1 TO NCOL(BRKPT); STP = BRKPT[1, CNT]; HOLDPOW = HOLDPOWER[START:STP,]; IF CNT=1 THEN PRINT / "****POWER CALCULATIONS****" , HOLDPOW[ROWNAME=" " COLNAME=HOLDPOWERLBL]; ELSE PRINT / HOLDPOW[ROWNAME=" " COLNAME=HOLDPOWERLBL]; START = STP + 1; FREE HOLDPOW; END; END; END; ***J.3) Create SAS dataset if option was requested; IF OPTDATASET THEN CALL _SASDS(HOLDPOWER, HOLDPOWERLBL, DSNAME,OPTWARNPRINT); ***K) Free matrices on the hit list and FINISH; GO TO NOERROREND; ERROREND: ; NOTE = {"PROGRAM TERMINATED WITH ONE OR MORE MAJOR ERRORS." , "PROGRAM TERMINATED WITH ONE OR MORE MAJOR ERRORS." , "PROGRAM TERMINATED WITH ONE OR MORE MAJOR ERRORS."} ; PRINT NOTE[LABEL=" " ROWNAME=" " COLNAME=" "]; NOERROREND: ; DO CNT = 1 TO NROW(HITLIST); IF HITLIST[CNT] = {U} THEN FREE U; ELSE IF HITLIST[CNT] = {THETA0} THEN FREE THETA0; ELSE IF HITLIST[CNT] = {REPN} THEN FREE REPN; ELSE IF HITLIST[CNT] = {BETASCAL} THEN FREE BETASCAL; ELSE IF HITLIST[CNT] = {SIGSCAL} THEN FREE SIGSCAL; ELSE IF HITLIST[CNT] = {RHOSCAL} THEN FREE RHOSCAL; ELSE IF HITLIST[CNT] = {ALPHA} THEN FREE ALPHA; ELSE IF HITLIST[CNT] = {ROUND} THEN FREE ROUND; ELSE IF HITLIST[CNT] = {TOLERANCE} THEN FREE TOLERANCE; ELSE IF HITLIST[CNT] = {UCDF} THEN FREE UCDF; ELSE IF HITLIST[CNT] = {UMETHOD} THEN FREE UMETHOD; ELSE IF HITLIST[CNT] = {MMETHOD} THEN FREE MMETHOD; ELSE IF HITLIST[CNT] = {CLTYPE} THEN FREE CLTYPE; ELSE IF HITLIST[CNT] = {ALPHA_CL} THEN FREE ALPHA_CL; ELSE IF HITLIST[CNT] = {ALPHA_CU} THEN FREE ALPHA_CU; ELSE IF HITLIST[CNT] = {DSNAME} THEN FREE DSNAME; END; FINISH _POWER; *************************************************************************; * *; * _RANKSYMM *; * *; * This function computes the rank of a square symmetric *; * nonnegative definite matrix via eigenvalues. *; * *; * INPUTS *; * *; * 1) Program supplied *; * *; * MATRIX, matrix which will be checked *; * MATRIXLBL, label to identify the matrix *; * TOLERANCE, value not tolerated, numeric zero (global) *; * *; * RETURNS *; * *; * . if MATRIX is not symmetric or positive definite *; * otherwise returns 0<= RANKMATRIX <=dimension of MATRIX *; * *; *************************************************************************; START _RANKSYMM ( MATRIX, MATRIXLBL) GLOBAL (TOLERANCE); ***The following two checks should be performed if this module is used; ***outside of _POWER(), which guarantees the conditions are fulfilled; ***IF NCOL(MATRIX)=0 THEN DO; *** RANKMATRIX = .; *** PRINT CONCAT("ERROR 55: MATRIX ", MATRIXLBL, "does not exist."); *** RETURN(RANKMATRIX); *** END; ***IF NCOL(MATRIX)^=NROW(MATRIX) THEN DO; *** RANKMATRIX = .; *** PRINT CONCAT("ERROR 56: MATRIX ", MATRIXLBL, "is not square."); *** RETURN(RANKMATRIX); *** END; IF ALL(MATRIX = .) THEN DO; RANKMATRIX = .; NOTE = COMPBL(CONCAT("ERROR 57: Matrix ", MATRIXLBL, " is all ", "missing values")); PRINT NOTE[LABEL=" " ROWNAME=" " COLNAME=" "] , MATRIX [LABEL = MATRIXLBL]; RETURN(RANKMATRIX); END; MAXABSVAL = MAX(ABS(MATRIX)); IF MAXABSVAL <= TOLERANCE THEN DO; RANKMATRIX = .; NOTE = COMPBL(CONCAT("ERROR 58: Matrix ", MATRIXLBL, " has MAX(ABS(all elements)) <= TOLERANCE.")); PRINT NOTE[LABEL=" " ROWNAME=" " COLNAME=" "] , MATRIX [LABEL=MATRIXLBL]; RETURN(RANKMATRIX); END; NMATRIX = MATRIX / MAXABSVAL; ***sqrt in next line due to numerical inaccuracy in sscp calculation with some data via sweep; IF MAX(ABS(NMATRIX - NMATRIX`)) >= SQRT(TOLERANCE) THEN DO; RANKMATRIX = .; NOTE = COMPBL(CONCAT("ERROR 59: Matrix ", MATRIXLBL, " is not symmetric within SQRT(TOLERANCE).")); PRINT NOTE[LABEL=" " ROWNAME=" " COLNAME=" "] , MATRIX [LABEL=MATRIXLBL]; RETURN(RANKMATRIX); END; EVALS = EIGVAL(NMATRIX); IF MIN(EVALS) < (-SQRT(TOLERANCE)) THEN DO; RANKMATRIX = .; NOTE = COMPBL(CONCAT("ERROR 60: Matrix ", MATRIXLBL, " is *NOT* nonnegative definite (and has at ", "least one eigenvalue strictly less than ", "zero). This may happen due to programming ", "error or rounding error of a nearly LTFR ", "matrix. This may be able to be fixed using ", "usual scaling/centering techniques. The ", "Eigenvalues/MAX(ABS(original matrix)) are: ")); PRINT NOTE[LABEL=" " ROWNAME=" " COLNAME=" "] , EVALS [LABEL=" " ROWNAME=" "]; NOTE = CONCAT(" The MAX(ABS(original matrix)) is ", CHAR(MAXABSVAL,1), ".") ; PRINT NOTE[LABEL=" " ROWNAME=" " COLNAME=" "]; RETURN(RANKMATRIX); END; RANKMATRIX = SUM(EVALS > SQRT(TOLERANCE)); RETURN(RANKMATRIX); FINISH _RANKSYMM; *************************************************************************; * *; * _SIZECHK *; * *; * This function checks if a matrix has more than one column. *; * *; * INPUTS *; * *; * 1) Program supplied *; * *; * MATRIX, matrix which will be checked *; * MATRIXLBL, label to identify the matrix *; * *; * RETURNS *; * *; * 1, if MATRIX has more than one column *; * 0, otherwise *; * *; *************************************************************************; START _SIZECHK ( MATRIX, MATRIXLBL); CHECKER = 0; IF NCOL(MATRIX) > 1 THEN DO; CHECKER = 1; NOTE = COMPBL(CONCAT("ERROR 61: Input ", MATRIXLBL, " has more than one row and more than one ", "column. We require that ", MATRIXLBL, " be a vector not a matrix.")); MATRIXPRINT = MATRIX`; PRINT NOTE[LABEL=" " ROWNAME=" " COLNAME=" "] , MATRIXPRINT[LABEL=MATRIXLBL]; RETURN(CHECKER); END; RETURN(CHECKER); FINISH _SIZECHK; *************************************************************************; * *; * _TYPEMISSCHK *; * *; * This function verifies that a matrix is of the required type, *; * either character or numeric. Also checks for that a matrix does not *; * contain missing values. *; * *; * INPUTS *; * *; * 1) Program supplied *; * *; * MATRIX, matrix which will be checked *; * MATRIXLBL, label to identify the matrix *; * TARGET, one character, "N" (numeric) or "C" (character) *; * *; * RETURNS *; * *; * 1, if MATRIX is the required type *; * 0, otherwise *; * *; *************************************************************************; START _TYPEMISSCHK ( MATRIX, MATRIXLBL, TARGET); CHECKER = 0; IF TARGET = "N" THEN DO; IF TYPE(MATRIX) ^= TARGET THEN DO; CHECKER = 1; NOTE = COMPBL(CONCAT("ERROR 62: Matrix ", MATRIXLBL, " must be numeric.")); PRINT NOTE[LABEL=" " ROWNAME=" " COLNAME=" "] , MATRIX [LABEL = MATRIXLBL]; RETURN(CHECKER); END; IF NROW(LOC(MATRIX = .)) > 0 THEN DO; CHECKER = 1; NOTE = COMPBL(CONCAT("ERROR 63: Matrix ", MATRIXLBL, " cannot contain missing values.")); PRINT NOTE[LABEL=" " ROWNAME=" " COLNAME=" "] , MATRIX [LABEL = MATRIXLBL]; RETURN(CHECKER); END; END; IF TARGET = "C" THEN DO; IF TYPE(MATRIX) ^= TARGET THEN DO; CHECKER = 1; NOTE = COMPBL(CONCAT("ERROR 64: Matrix ", MATRIXLBL, " must be character.")); PRINT NOTE[LABEL=" " ROWNAME=" " COLNAME=" "] , MATRIX [LABEL = MATRIXLBL]; RETURN(CHECKER); END; IF NROW(LOC(MATRIX = ".")) > 0 THEN DO; CHECKER = 1; NOTE = COMPBL(CONCAT("ERROR 65: Matrix ", MATRIXLBL, " cannot contain missing values.")); PRINT NOTE[LABEL=" " ROWNAME=" " COLNAME=" "] , MATRIX [LABEL = MATRIXLBL]; RETURN(CHECKER); END; END; RETURN(CHECKER); FINISH _TYPEMISSCHK; *************************************************************************; * *; * _SETOPT *; * *; * This module _SETOPT sets the options requested by the user. *; * *; * *; * INPUTS *; * *; * 1) Program supplied *; * *; * NEWOPTN, column of on/off switches (0 or 1) for options *; * NEWOPTNLBL, column of option names which label NEWOPTN *; * OPT_ON, column of user requested options to turn on *; * OPT_OFF, column of user requested options to turn off *; * *; * OUTPUTS *; * *; * NEWOPTN, column of switches with requested options turned on/off*; * *; *************************************************************************; START _SETOPT (NEWOPTN, NEWOPTNLBL, OPT_ON, OPT_OFF); DO CNT = 1 TO NROW(OPT_ON); NEWOPTN = NEWOPTN <> ( NEWOPTNLBL = OPT_ON[CNT,1] ); END; DO CNT = 1 TO NROW(OPT_OFF); NEWOPTN = NEWOPTN # ( NEWOPTNLBL ^= OPT_OFF[CNT,1] ); END; FINISH _SETOPT; *************************************************************************; * *; * _HLT *; * *; * This module calculates power for Hotelling-Lawley trace *; * based on the Pillai F approximation. HLT is the "population value" *; * Hotelling Lawley trace. F1 and DF2 are the hypothesis and *; * error degrees of freedom, OMEGA is the noncentrality parameter, and *; * FCRIT is the critical value from the F distribution. *; * *; * INPUTS *; * *; * 1) Program supplied *; * *; * A, rank of C matrix *; * B, rank of U matrix *; * S, minimum of A and B *; * N, total N *; * R, rank of X *; * EVAL, eigenvalues for H*INV(E) *; * ALPHASCALAR, size of test *; * OPTPOWERMAT2, options matrix specifying CL options *; * *; * 2) User supplied (or program default) - GLOBAL *; * *; * CLTYPE, N_EST, RANK_EST, ALPHA_CL, ALPHA_CU *; * *; * OUTPUTS *; * *; * POWER, power for Hotelling-Lawley trace & CL if requested *; * LBL, labels for entries in POWER *; * POWERWARN, vector of power calculation warning counts *; * *; *************************************************************************; START _HLT (POWER, LBL, POWERWARN, A, B, S, N, R, EVAL, ALPHASCALAR, MMETHODTEMP, OPTPOWERMAT2) GLOBAL (CLTYPE, N_EST, RANK_EST, ALPHA_CL, ALPHA_CU); IF (MMETHODTEMP[1,1] = 1) | (MMETHODTEMP[1,1] = 3) THEN DO; DF1 = A # B; DF2 = S # (N - R - B - 1) + 2; END; IF (MMETHODTEMP[1,1] = 2) | (MMETHODTEMP[1,1] = 4) THEN DO; DF1 = A # B; DF2 = (N - R) # (N - R) - (N - R) # (2 # B + 3) + B # (B + 3); DF2 = DF2 / ((N - R) # (A + B + 1) - (A + 2 # B + B # B - 1)); DF2 = 4 + (A # B + 2) # DF2; END; IF (DF2 <= 0) | (EVAL[1,1] = .) THEN DO; POWER = . ; POWERWARN[15,1] = POWERWARN[15,1] + 1; END; ELSE DO; IF (MMETHODTEMP[1,1] > 2) | (S = 1) THEN EVALT = EVAL # (N - R) / N; ELSE EVALT = EVAL; HLT = EVALT[+,]; IF (MMETHODTEMP[1,1] > 2) | (S = 1) THEN OMEGA = (N # S) # (HLT / S); ELSE OMEGA = (DF2) # (HLT / S); RUN _FINV (FCRIT, 1 - ALPHASCALAR, DF1, DF2); RUN _PROBF(PROB, FMETHOD, FCRIT, DF1, DF2, OMEGA); RUN _FWARN(POWERWARN, FMETHOD, 1); IF (FMETHOD = 4) & (PROB = 1) THEN POWER = ALPHASCALAR; ELSE POWER = 1 - PROB; END; LBL = {"POWER_HLT"}; IF CLTYPE >= 1 THEN DO; LBL = {"POWER_HLT_L" "POWER_HLT" "POWER_HLT_U"}; IF POWER = . THEN DO; POWER={. . .}; POWERWARN[16,1] = POWERWARN[16,1] + 1; END; ELSE DO; F_A = OMEGA / DF1; CALL _GLMMPCL(POWER_L, POWER_U, FMETHOD_L, FMETHOD_U, POWERWARN, F_A, ALPHASCALAR, DF1, N, DF2); POWER = POWER_L || POWER || POWER_U; END; END; IF OPTPOWERMAT2[1,1] = 1 THEN DO; IF CLTYPE >= 1 THEN DO; POWER = FMETHOD_L || FMETHOD || FMETHOD_U || POWER; LBL = "FMETHOD_HLT_L" || "FMETHOD_HLT" || "FMETHOD_HLT_U" || LBL; END; ELSE DO; POWER = FMETHOD || POWER; LBL = "FMETHOD_HLT" || LBL; END; END; FINISH _HLT; *************************************************************************; * *; * _PBT *; * *; * This module calculates power for Pillai-Bartlett trace based on the *; * F approx. method. V is the "population value" of PBT , *; * DF1 and DF2 are the hypothesis and error degrees of freedom, *; * OMEGA is the noncentrality parameter, and FCRIT is the *; * critical value from the F distribution. *; * *; * INPUTS *; * *; * 1) Program supplied *; * *; * A, rank of C matrix *; * B, rank of U matrix *; * S, minimum of A and B *; * N, total N *; * R, rank of X *; * EVAL, eigenvalues for H*INV(E) *; * ALPHASCALAR, size of test *; * OPTPOWERMAT2, options matrix specifying CL options *; * *; * 2) User supplied (or program default) - GLOBAL *; * *; * TOLERANCE, CLTYPE, N_EST, RANK_EST, ALPHA_CL, ALPHA_CU *; * *; * OUTPUTS *; * *; * POWER, power for Pillai-Bartlett trace & CL if requested *; * LBL, labels for entries in POWER *; * POWERWARN, vector of power calculation warning counts *; * *; *************************************************************************; START _PBT (POWER, LBL, POWERWARN, A, B, S, N, R, EVAL, ALPHASCALAR, MMETHODTEMP, OPTPOWERMAT2) GLOBAL(TOLERANCE, CLTYPE, N_EST, RANK_EST, ALPHA_CL, ALPHA_CU); FMETHOD = 0; IF (MMETHODTEMP[2,1] = 1) | (MMETHODTEMP[2,1] = 3) THEN DO; DF1 = A # B; DF2 = S # (N - R + S - B); END; IF (MMETHODTEMP[2,1] = 2) | (MMETHODTEMP[2,1] = 4) THEN DO; MU1= A # B / (N - R + A); FACTOR1 = (N - R + A - B) / (N - R + A - 1); FACTOR2 = (N - R) / (N - R + A + 2); VARIANCE = 2 # A # B # FACTOR1 # FACTOR2 / (N - R + A) ##2; MU2 = VARIANCE + MU1 ## 2; M1 = MU1 / S; M2 = MU2 / (S#S); DENOM = M2 - M1 # M1; DF1 = 2 # M1 # (M1 - M2) / DENOM; DF2 = 2 # (M1 - M2) # (1 - M1) / DENOM; END; IF (DF2 <= 0) | (EVAL[1,1] = .) THEN DO; POWER = .; POWERWARN[15,1] = POWERWARN[15,1] + 1; END; ELSE DO; IF (MMETHODTEMP[2,1] > 2) | (S = 1) THEN EVALT = EVAL # (N - R) / N; ELSE EVALT = EVAL; V = SUM(EVALT / (J(S,1,1) + EVALT)); IF (S - V) <= TOLERANCE THEN POWER=.; ELSE DO; IF (MMETHODTEMP[2,1] > 2) | (S = 1) THEN OMEGA = (N # S) # V / (S - V); ELSE OMEGA = (DF2) # V / (S - V); RUN _FINV(FCRIT, 1 - ALPHASCALAR, DF1, DF2); RUN _PROBF(PROB, FMETHOD, FCRIT, DF1, DF2, OMEGA); RUN _FWARN(POWERWARN, FMETHOD, 1); IF (FMETHOD = 4) & (PROB = 1) THEN POWER = ALPHASCALAR; ELSE POWER = 1 - PROB; END; END; LBL = {"POWER_PBT"}; IF CLTYPE >= 1 THEN DO; LBL = {"POWER_PBT_L" "POWER_PBT" "POWER_PBT_U"}; IF POWER = . THEN DO; POWER = {. . .}; POWERWARN[16,1] = POWERWARN[16,1] + 1; END; ELSE DO; F_A = OMEGA / DF1; CALL _GLMMPCL(POWER_L, POWER_U, FMETHOD_L, FMETHOD_U, POWERWARN, F_A, ALPHASCALAR, DF1, N, DF2); POWER = POWER_L || POWER || POWER_U; END; END; IF OPTPOWERMAT2[1,1] = 1 THEN DO; IF CLTYPE >= 1 THEN DO; POWER = FMETHOD_L || FMETHOD || FMETHOD_U || POWER; LBL = "FMETHOD_PBT_L" || "FMETHOD_PBT" || "FMETHOD_PBT_U" || LBL; END; ELSE DO; POWER = FMETHOD || POWER; LBL = "FMETHOD_PBT" || LBL; END; END; FINISH _PBT; *************************************************************************; * *; * _WLK *; * *; * This module calculates power for Wilk's Lambda based on *; * the F approx. method. W is the "population value" of Wilks` Lambda, *; * DF1 and DF2 are the hypothesis and error degrees of freedom, OMEGA *; * is the noncentrality parameter, and FCRIT is the critical value *; * from the F distribution. RM, RS, R1, and TEMP are intermediate *; * variables. *; * *; * INPUTS *; * *; * 1) Program supplied *; * *; * A, rank of C matrix *; * B, rank of U matrix *; * S, minimum of A and B *; * N, total N *; * R, rank of X *; * EVAL, eigenvalues for H*INV(E) *; * ALPHASCALAR, size of test *; * OPTPOWERMAT2, options matrix specifying CL options *; * *; * 2) User supplied (or program default) - GLOBAL *; * *; * CLTYPE, N_EST, RANK_EST, ALPHA_CL, ALPHA_CU *; * *; * OUTPUTS *; * *; * POWER, power for Wilk's Lambda & CL if requested *; * LBL, labels for entries in POWER *; * POWERWARN, vector of power calculation warning counts *; * *; *************************************************************************; START _WLK (POWER, LBL, POWERWARN, A, B, S, N, R, EVAL, ALPHASCALAR, MMETHODTEMP, OPTPOWERMAT2) GLOBAL (CLTYPE, N_EST, RANK_EST, ALPHA_CL, ALPHA_CU); FMETHOD = 0; DF1 = A # B; IF EVAL[1,1] = . THEN DO; W = . ; POWERWARN[15,1] = POWERWARN[15,1]; END; ELSE DO; IF (MMETHODTEMP[3,1] = 2) | (MMETHODTEMP[3,1] = 4) | (S = 1) THEN EVALT = EVAL # (N - R) / N; ELSE EVALT = EVAL; W = EXP(SUM(-LOG(J(S,1,1) + EVALT))); END; IF S = 1 THEN DO; DF2 = N - R - B + 1; RS = 1; TEMPW = W; END; ELSE DO; RM = N - R - (B - A + 1) / 2; RS = SQRT( (A # A # B # B - {4})/(A # A + B # B - {5}) ); R1 = (B # A - {2})/4; IF W = . THEN TEMPW = .; ELSE TEMPW = W ## (1 / RS); DF2 = (RM # RS) - 2 # R1; END; IF TEMPW = . THEN OMEGA = .; ELSE DO; IF (MMETHODTEMP[3,1] = 2) | (MMETHODTEMP[3,1] = 4) | (S = 1) THEN OMEGA = (N # RS) # (1 - TEMPW) / TEMPW; ELSE OMEGA = (DF2) # (1 - TEMPW) / TEMPW; END; IF (DF2 <= 0) | (W = .) | (OMEGA = .) THEN DO; POWER = .; POWERWARN[15,1] = POWERWARN[15,1] + 1; END; ELSE DO; RUN _FINV(FCRIT, 1 - ALPHASCALAR, DF1, DF2); RUN _PROBF(PROB, FMETHOD, FCRIT, DF1, DF2, OMEGA); RUN _FWARN(POWERWARN, FMETHOD, 1); IF (FMETHOD = 4) & (PROB = 1) THEN POWER = ALPHASCALAR; ELSE POWER = 1 - PROB; END; LBL = {"POWER_WLK"}; IF CLTYPE >= 1 THEN DO; LBL = {"POWER_WLK_L" "POWER_WLK" "POWER_WLK_U"}; IF POWER = . THEN DO; POWER = {. . .}; POWERWARN[16,1] = POWERWARN[16,1] + 1; END; ELSE DO; F_A = OMEGA / DF1; CALL _GLMMPCL(POWER_L, POWER_U, FMETHOD_L, FMETHOD_U, POWERWARN, F_A, ALPHASCALAR, DF1, N, DF2); POWER = POWER_L || POWER || POWER_U; END; END; IF OPTPOWERMAT2[1,1] = 1 THEN DO; IF CLTYPE >= 1 THEN DO; POWER = FMETHOD_L || FMETHOD || FMETHOD_U || POWER; LBL = "FMETHOD_WLK_L" || "FMETHOD_WLK" || "FMETHOD_WLK_U" || LBL; END; ELSE DO; POWER = FMETHOD || POWER; LBL = "FMETHOD_WLK" || LBL; END; END; FINISH _WLK; *************************************************************************; * *; * _SPECIAL *; * *; * This module performs two disparate tasks. For B=1 (UNIVARIATE *; * TEST), the powers are calculated more efficiently. For A=1 (SPECIAL *; * MULTIVARIATE CASE), exact multivariate powers are calculated. *; * Powers for the univariate tests require separate treatment. *; * DF1 & DF2 are the hypothesis and error degrees of freedom, *; * OMEGA is the noncentrality parameter, and FCRIT is the critical *; * value from the F distribution. *; * *; * INPUTS *; * *; * 1) Program supplied *; * *; * A, rank of C matrix *; * B, rank of U matrix *; * S, minimum of A and B *; * N, total N *; * R, rank of X *; * EVAL, eigenvalues for H*INV(E) *; * ALPHASCALAR, size of test *; * OPTPOWERMAT2, options matrix specifying CL options *; * *; * OUTPUTS *; * *; * POWER, power for special case when RANK(CBU)=1, & CL if req *; * LBL, labels for entries in POWER *; * POWERWARN, vector of power calculation warning counts *; * *; *************************************************************************; START _SPECIAL (POWER, LBL, POWERWARN, A, B, S, N, R, EVAL, ALPHASCALAR, OPTPOWERMAT2) GLOBAL (CLTYPE, N_EST, RANK_EST, ALPHA_CL, ALPHA_CU); FMETHOD = 0; DF1 = A # B; DF2 = N - R - B + 1; IF (DF2 <= 0) | (EVAL[1,1]=.) THEN DO; POWER = .; POWERWARN[15,1] = POWERWARN[15,1] + 1; END; ELSE DO; OMEGA = EVAL[1,1] # (N - R); RUN _FINV(FCRIT, 1 - ALPHASCALAR, DF1, DF2); RUN _PROBF(PROB, FMETHOD, FCRIT, DF1, DF2, OMEGA); RUN _FWARN(POWERWARN, FMETHOD, 1); IF (FMETHOD = 4) & (PROB = 1) THEN POWER = ALPHASCALAR; ELSE POWER = 1 - PROB; END; IF (A =1) & (B>1) THEN LBL = {"POWER_MULT"}; ELSE LBL = {"POWER"}; IF (A =1) & (B>1) THEN F_METHODLBL = {"F_METHOD_MULT"}; ELSE F_METHODLBL = {"F_METHOD"}; IF CLTYPE >= 1 THEN DO; LBL = CONCAT(LBL, "_L") || LBL || CONCAT(LBL, "_U") ; IF POWER = . THEN DO; POWER = {. . .}; POWERWARN[16,1] = POWERWARN[16,1] + 1; END; ELSE DO; F_A = OMEGA / DF1; CALL _GLMMPCL(POWER_L, POWER_U, FMETHOD_L, FMETHOD_U, POWERWARN, F_A, ALPHASCALAR, DF1, N, DF2); POWER = POWER_L || POWER || POWER_U; END; END; IF OPTPOWERMAT2[1,1] = 1 THEN DO; IF CLTYPE >= 1 THEN DO; POWER = FMETHOD_L || FMETHOD || FMETHOD_U || POWER; LBL = CONCAT(F_METHODLBL,"_L") || F_METHODLBL || CONCAT(F_METHODLBL,"_U") || LBL; END; ELSE DO; POWER = FMETHOD || POWER; LBL = F_METHODLBL || LBL; END; END; FINISH _SPECIAL; *************************************************************************; * *; * _FIRSTUNI *; * *; * Univariate STEP 1: *; * *; * This module produces matrices required for Geisser-Greenhouse, *; * Huynh-Feldt or uncorrected repeated measures power calculations. It *; * is the first step. Program uses approximations of expected values of *; * epsilon estimates due to Muller (1985), based on theorem of Fujikoshi *; * (1978). Program requires that U be orthonormal and orthogonal to a *; * columns of 1's. *; * *; * INPUTS *; * *; * 1) Program Supplied *; * *; * SIGMASTAR = U` * (SIGMA # SIGSCALTEMP) * U *; * B, rank of U *; * *; * 2) User supplied - GLOBAL *; * *; * TOLERANCE *; * *; * OUTPUTS *; * *; * D, number of distinct eigenvalues *; * MTP, multiplicities of eigenvalues *; * EPS, epsilon calculated from U`*SIGMA*U *; * DEIGVAL, first eigenvalue *; * SLAM1, sum of eigenvalues squared *; * SLAM2, sum of squared eigenvalues *; * SLAM3, sum of eigenvalues *; * *; *************************************************************************; START _FIRSTUNI (D, MTP, EPS, DEIGVAL, SLAM1, SLAM2, SLAM3, SIGMASTAR, B) GLOBAL(TOLERANCE); ***Get eigenvalues of covariance matrix associated with E. This is NOT; ***the USUAL sigma. This cov matrix is that of (Y-YHAT)*U, not of (Y-YHAT).; ***The covariance matrix is normalized to minimize numerical problems ; ESIG = SIGMASTAR / ( TRACE(SIGMASTAR) ); SEIGVAL = EIGVAL(ESIG); SLAM1 = SUM(SEIGVAL)##2; SLAM2 = SSQ(SEIGVAL); SLAM3 = SUM(SEIGVAL); EPS = SLAM1 / (B # SLAM2); ***Decide which eigenvalues are distinct; D = 1; ***Ends as number of distinct eigenvalues; MTP = 1; ***Ends as vector of multiplicities of eignvals; DEIGVAL = SEIGVAL[1,1]; DO CNT = 2 TO B; IF (DEIGVAL[D,1] - SEIGVAL[CNT,1]) > TOLERANCE THEN DO; D = D + 1; DEIGVAL = DEIGVAL // SEIGVAL[CNT,1]; MTP = MTP // {1}; END; ELSE MTP[D,1] = MTP[D,1] + 1; END; FINISH _FIRSTUNI; *************************************************************************; * *; * _HFEXEPS *; * *; * Univariate, HF STEP 2: *; * This function computes the approximate expected value of *; * the Huynh-Feldt estimate. *; * *; * FK = 1st deriv of FNCT of eigenvalues *; * FKK = 2nd deriv of FNCT of eigenvalues *; * For HF, FNCT is epsilon tilde *; * *; * INPUTS *; * *; * 1) Program supplied *; * *; * B, rank of U *; * N, total number of observations *; * R, rank of X *; * D, number of distinct eigenvalues *; * MTP, multiplicities of eigenvalues *; * EPS, epsilon calculated from U`*SIGMA*U *; * DEIGVAL, first eigenvalue *; * SLAM1, sum of eigenvalues squared *; * SLAM2, sum of squared eigenvalues *; * SLAM3, sum of eigenvalues *; * SIGMASTAREVAL, eigenvalues of SIGMASTAR=U`*SIGMA*U *; * SIGMASTAREVEC, eigenvectors of SIGMASTAR=U`*SIGMA*U *; * UMETHODTEMP, vector choosing approximation for epsilon in *; * GG and HF tests *; * *; * RETURNS *; * *; * EXEPS, expected value of epsilon estimator *; * *; *************************************************************************; START _HFEXEPS ( B, N, R, D, MTP, EPS, DEIGVAL, SLAM1, SLAM2, SLAM3, SIGMASTAREVAL, SIGMASTAREVEC, UMETHODTEMP); *** Compute approximate expected value of Huynh-Feldt estimate; H1 = N # SLAM1 - 2 # SLAM2; H2 = (N - R) # SLAM2 - SLAM1; DERH1 = J(D, 1, 2 # N # SLAM3) - 4 # DEIGVAL; DERH2 = 2 # (N - R) # DEIGVAL - J(D,1,2 # SQRT(SLAM1)); FK = DERH1 - H1 # DERH2 / H2; FK = FK / (B # H2); DER2H1 = J(D, 1, 2 # N - 4); DER2H2 = J(D, 1, 2 # (N - R) - 2); FKK = -DERH1 # DERH2 / H2 + DER2H1 - DERH1 # DERH2 / H2 + 2 # H1 # (DERH2 ## 2) / H2 ## 2 - H1 # DER2H2 / H2; FKK = FKK / (H2 # B); T1 = FKK # (DEIGVAL ## 2) # MTP; SUM1 = SUM(T1); IF D = 1 THEN SUM2 = 0; ELSE DO; T2 = FK # DEIGVAL # MTP; T3 = DEIGVAL # MTP; TM1 = T2 * T3`; T4 = DEIGVAL * J(1,D,1); TM2 = T4 - T4`; TM2INV = 1 / (TM2 + I(D)) - I(D); TM3 = TM1 # TM2INV; SUM2 = SUM(TM3); END; *** Define HF Approx E(.) for Method 0; E0EPSHF = H1 / (B # H2) + (SUM1 + SUM2) / (N - R); *** Computation of EXP(T1) and EXP(T2); DO CNT = 1 TO D; SEVAL = SEVAL // J(MTP[CNT],1,DEIGVAL[CNT]); END; NU = N - R; EXPT1 = (2 # NU # SLAM2) + ((NU # NU) # SLAM1); EXPT2 = ((NU # (NU + 1)) # SLAM2) + (NU # SUM(SEVAL * SEVAL`)); *** End computation of EXP(T1) and EXP(T2); ***For use with Method 1; NUM01 = (1 / B) # (((NU+1) # EXPT1) - (2 # EXPT2)); DEN01 = (NU # EXPT2) - EXPT1; *** Define HF Approx E(.) for Method 1; E1EPSHF = NUM01 / DEN01; *** Method 1 ; *** Select appropriate approximate value; IF UMETHODTEMP[1,1] = 1 THEN EXEPS = E0EPSHF; *** Select HF Method 1 - MB; IF UMETHODTEMP[1,1] = 2 THEN EXEPS = E1EPSHF; *** Select HF Method 2; RETURN (EXEPS); FINISH _HFEXEPS; *************************************************************************; * *; * _GGEXEPS *; * *; * Univariate, GG STEP 2: *; * This function computes the approximate expected value of the *; * Geisser-Greenhouse estimate. *; * *; * FK = 1st deriv of FNCT of eigenvalues *; * FKK = 2nd deriv of FNCT of eigenvalues *; * For GG FNCT is epsilon hat *; * *; * INPUTS *; * *; * 1) Program supplied *; * *; * B, rank of U *; * N, total number of observations *; * R, rank of X *; * D, number of distinct eigenvalues *; * MTP, multiplicities of eigenvalues *; * EPS, epsilon calculated from U`*SIGMA*U *; * DEIGVAL, first eigenvalue *; * SLAM1, sum of eigenvalues squared *; * SLAM2, sum of squared eigenvalues *; * SLAM3, sum of eigenvalues *; * SIGMASTAREVAL, Bx1 vector of eigenvalues of SIGMASTAR=U`*SIGMA*U*; * SIGMASTAREVEC, eigenvectors of SIGMASTAR=U`*SIGMA*U *; * UMETHODTEMP, vector choosing approximation for epsilon in *; * GG and HF tests *; * *; * RETURNS *; * *; * EXEPS, expected value of epsilon estimator *; * *; *************************************************************************; START _GGEXEPS ( B, N, R, D, MTP, EPS, DEIGVAL, SLAM1, SLAM2, SLAM3, SIGMASTAREVAL, SIGMASTAREVEC, UMETHODTEMP); FK = J(D,1,1) # 2 # SLAM3 / (SLAM2 # B) - 2 # DEIGVAL # SLAM1 / (B # SLAM2 ## 2); C0 = 1 - SLAM1 / SLAM2; C1 = -4 # SLAM3 / SLAM2; C2 = 4 # SLAM1 / SLAM2 ## 2; FKK = C0 # J(D,1,1) + C1 # DEIGVAL + C2 # DEIGVAL ## 2; FKK = 2 # FKK / (B # SLAM2); T1 = FKK # (DEIGVAL ## 2) # MTP; SUM1 = SUM(T1); IF D = 1 THEN SUM2 = 0; ELSE DO; T2 = FK # DEIGVAL # MTP; T3 = DEIGVAL # MTP; TM1 = T2 * T3`; T4 = DEIGVAL * J(1,D,1); TM2 = T4 - T4`; TM2INV = 1 / (TM2 + I(D)) - I(D); TM3 = TM1 # TM2INV; SUM2 = SUM(TM3); END; *** Define GG Approx E(.) for Method 0; E0EPSGG = EPS + (SUM1 + SUM2) / (N - R); ***Computation of EXP(T1) and EXP(T2); DO CNT = 1 TO D; SEVAL = SEVAL // J(MTP[CNT], 1, DEIGVAL[CNT]); END; NU = N - R; EXPT1 = (2 # NU # SLAM2) + ((NU # NU) # SLAM1); EXPT2 = ((NU # (NU + 1)) # SLAM2) + (NU # SUM(SEVAL * SEVAL`)); ***End computation of EXP(T1) and EXP(T2); ***Define GG Approx E(.) for Method 1; E1EPSGG = (1 / B) # (EXPT1 / EXPT2); *** Method 1 Approx E(.); *** Select appropriate approximate value; EXEPS = .; IF UMETHODTEMP[2,1] = 1 THEN EXEPS = E0EPSGG; *** Select GG Method 1 - MB; IF UMETHODTEMP[2,1] = 2 THEN EXEPS = E1EPSGG; *** Select GG Method 2; RETURN (EXEPS); FINISH _GGEXEPS; *************************************************************************; * *; * _LASTUNI *; * *; * Univariate STEP 3 *; * This module performs the final step for univariate repeated *; * measures power calculations. *; * *; * INPUTS *; * *; * 1) Program supplied *; * *; * EXEPS, expected value epsilon estimator *; * H, hypothesis sum of squares *; * E, error sum of squares *; * A, rank of C *; * B, rank of U *; * R, rank of X *; * N, total number of observations *; * EPS, epsilon calculated from U`*SIGMA*U *; * POWERCALC, indicates selected power calculation *; * OPTPOWERMAT2, options matrix specifying CL options *; * ALPHASCALAR, size of test *; * FIRSTUNI, indicates first requested univariate test *; * SIGMASTAREVAL, eigenvalues of SIGMASTAR=U`*SIGMA*U *; * SIGMASTAREVEC, eigenvectors of SIGMASTAR=U`*SIGMA*U *; * UCDFTEMP, vector of CDF form for Un, GG, HF, Box tests *; * UMETHODTEMP, vector choosing approximation for epsilon in *; * GG and HF tests *; * *; * 2) User supplied (or program default) - GLOBAL *; * *; * CLTYPE, N_EST, RANK_EST, ALPHA_CL, ALPHA_CU, ROUND *; * *; * OUTPUTS *; * *; * POWER, power for selected power calculation & CL if requested *; * LBL, labels for entries in POWER *; * POWERWARN, vector of power calculation warning counts *; * *; *************************************************************************; START _LASTUNI (POWER, LBL, POWERWARN, EXEPS, H, E, A, B, R, N, EPS, POWERCALC, OPTPOWERMAT2, ALPHASCALAR, FIRSTUNI, SIGMASTAREVAL, SIGMASTAREVEC, UCDFTEMP, UMETHODTEMP) GLOBAL (CLTYPE, N_EST, RANK_EST, ALPHA_CL, ALPHA_CU, ROUND, TOLERANCE, NONCENCL); *** POWERCALC= 5=>UN 6=>HF 7=>GG 8=>Box ***; OPTPOWERCALC1LBL = {"COLLAPSE" "HLT" "PBT" "WLK" "UN" "HF" "GG" "BOX"}`; LBL = CONCAT( {"POWER_"} , OPTPOWERCALC1LBL[POWERCALC]); CDFPOWERCALC = UCDFTEMP[ POWERCALC - {4} ]; POWER = .; FMETHOD=0; NUE = N - R; IF (EXEPS = .) | (NUE <= 0) THEN GO TO POWDONE; UNDF1 = B # A; UNDF2 = B # NUE ; NONCENCL=J(1,3,.); *** Create defaults - same for either SIGMA known or estimated ***; SIGSTAR=E/NUE; Q1 = TRACE(SIGSTAR); Q2 = TRACE(H); Q3 = Q1**2; Q4 = SIGSTAR[##]; Q5 = TRACE(SIGSTAR*H); LAMBAR = Q1/B; IF NCOL(N_EST) = 0 THEN DO; ***Enter loop to compute E1-E5 based on known SIGMA***; EPSN_NUM = Q3 + Q1 # Q2 # (2/A); EPSN_DEN = Q4 + Q5 # (2/A); EPSN = EPSN_NUM / (B # EPSN_DEN); E_1_2 = EXEPS; E_4 = EPS; E_3_5 = EPSN; *** MEST, (CDFPOWERCALC = 2,3,4 ***); IF (CDFPOWERCALC = 1) THEN E_3_5 = EPS; *** MB ***; END; ELSE IF NCOL(N_EST)>0 THEN DO; *** Enter loop to compute E1-E5 based on estimated SIGMA ***; NU_EST = N_EST-RANK_EST; IF NU_EST<=0 THEN DO; NOTE = "ERROR: Too few estimation df in LASTUNI. df = N_EST - RANK_EST < 0. "; PRINT NOTE[LABEL=" " ROWNAME=" " COLNAME=" "], NU_EST N_EST RANK_EST; STOP; END; EPSTILDE_R = ( (NU_EST+{1}) # Q3 - {2} # Q4 ) / ( B # (NU_EST # Q4 - Q3) ) ; *** Generalized approx unbiased ***; EPSTILDE_RM = MIN(EPSTILDE_R,{1}); MULT = NU_EST # NU_EST + NU_EST - {2}; EPSNHAT_NUM = NU_EST # (Q3 # NU_EST # (NU_EST+{1}) + Q1 # Q2 # ({2}#MULT/A) - Q4 # {2} # NU_EST) ; EPSNHAT_DEN = (Q4 # NU_EST # NU_EST + Q5 #(2#MULT/A) - Q3# NU_EST)# NU_EST ; EPSNHAT = EPSNHAT_NUM / (B # EPSNHAT_DEN); E_1_2 = EXEPS; IF (POWERCALC={6}) THEN E_1_2 = EPSTILDE_RM; *** for HF crit val ***; IF (POWERCALC={7}) THEN E_1_2 = EPS; *** for GG crit val ***; IF (CDFPOWERCALC = 1) THEN E_3_5 = EPS; *** MB ***; ELSE E_3_5 = EPSNHAT; E_4=EPS; CL1DF = B # NU_EST # E_4 / E_3_5; END; IF (E_1_2 < (1 / B)) & (E_1_2 ^= .) THEN DO; E_1_2 = {1} / B; POWERWARN[17,1] = POWERWARN[17,1] + 1; END; IF E_1_2 > 1 THEN DO; E_1_2 = 1; POWERWARN[18,1] = POWERWARN[18,1] + 1; END; *** Obtain critical value for power point estimate ***; OMEGA = E_3_5 # Q2 / LAMBAR; NONCENCL[1,2] = OMEGA; RUN _FINV(FCRIT, 1 - ALPHASCALAR, UNDF1#E_1_2, UNDF2#E_1_2); *** Compute power point estimate ***; *** (1) Muller, Edwards & Taylor 2002 CDF exact, Davies' algorithm; *** This order avoids duplicating code of 2 for 4 if POWER = . ; IF (CDFPOWERCALC = 3) | (CDFPOWERCALC = 4) THEN DO; DF1 = . ; DF2 = . ; FMETHOD=.; ACCURACY = 10 ## (-ROUND - 1); QWEIGHT = SIGMASTAREVAL // (-SIGMASTAREVAL # FCRIT # (UNDF1 / UNDF2)); QNUVEC = J(B,1,A) // J(B,1,N - R); DGOVER = DIAG(1 / SQRT(SIGMASTAREVAL)); FACTORI = SIGMASTAREVEC * DGOVER; OMEGSTAR = FACTORI` * H * FACTORI; QNONCEN = VECDIAG(OMEGSTAR) // J(B,1,0); CDFPOWR = _QPROB({0}, QWEIGHT, QNUVEC, QNONCEN, ACCURACY); IF CDFPOWR = . THEN POWERWARN[19,1] = POWERWARN[19,1] + 1; ELSE POWER = 1 - CDFPOWR; END; *** (2) Muller, Edwards & Taylor 2002 and Muller Barton 1989 CDF approx; *** UCDFTEMP[]=4 reverts to UCDFTEMP[]=2 if exact CDF fails; IF (CDFPOWERCALC = 1) | (CDFPOWERCALC = 2) | ( (CDFPOWERCALC = 4) & (POWER = .)) THEN DO; DF1 = UNDF1#E_3_5; DF2 = UNDF2#E_4; RUN _PROBF(PROB, FMETHOD, FCRIT, DF1, DF2, OMEGA); RUN _FWARN(POWERWARN, FMETHOD, 1); IF (FMETHOD = 4) & (PROB = 1) THEN POWER = ALPHASCALAR; ELSE POWER = 1 - PROB; END; *** Compute CL for power, if requested by user ***; IF CLTYPE = 2 THEN DO; *change from chi sq to F, and only change:); NOTE = "CLTYPE=2 for UNIREP awaiting implementation"; PRINT NOTE[LABEL=" " ROWNAME=" " COLNAME=" "]; STOP; END; IF CLTYPE =1 THEN DO; IF (CDFPOWERCALC >2) THEN DO; NOTE = "ERROR: Any use of Exact CDF incompatible with power CL so stop."; PRINT NOTE[LABEL=" " ROWNAME=" " COLNAME=" "]; STOP; END; *** Calculate lower bound for power ***; IF (ALPHA_CL <= TOLERANCE) THEN DO; PROB = 1 - ALPHASCALAR; FMETHOD_L = 5; NONCEN_L = .; NONCENCL[1,1] = .; END; ELSE DO; CHI_L = CINV(ALPHA_CL, CL1DF); NONCEN_L = OMEGA # (CHI_L / CL1DF ) ; NONCENCL[1,1] = NONCEN_L; RUN _PROBF(PROB, FMETHOD_L, FCRIT, DF1, DF2, NONCEN_L); RUN _FWARN(POWERWARN, FMETHOD_L, 2); END; IF (FMETHOD_L = 4) & (PROB = 1) THEN POWER_L = ALPHASCALAR; ELSE POWER_L = 1 - PROB; *** Calculate upper bound for power ***; IF (ALPHA_CU <= TOLERANCE) THEN DO; PROB = 0; FMETHOD_U = 5; NONCEN_U = .; NONCENCL[1,3] = .; END; ELSE DO; CHI_U = CINV(1-ALPHA_CU, CL1DF); NONCEN_U = OMEGA # (CHI_U / CL1DF ) ; NONCENCL[1,3] = NONCEN_U; RUN _PROBF(PROB, FMETHOD_U, FCRIT, DF1, DF2, NONCEN_U); RUN _FWARN(POWERWARN, FMETHOD_U, 2); END; IF (FMETHOD_U = 4) & (PROB = 1) THEN POWER_U = ALPHASCALAR; ELSE POWER_U = 1 - PROB; *** Assemble output vectors ***; LBL = COMPRESS(CONCAT(LBL,"_L") || LBL || CONCAT(LBL,"_U")); POWER = POWER_L || POWER || POWER_U; IF OPTPOWERMAT2[1,1] = 1 THEN DO; POWER = FMETHOD_L || FMETHOD_U || POWER; LBL = "FMETHOD_L" || "FMETHOD_U"|| LBL ; END; END; POWDONE: ; IF OPTPOWERMAT2[1,1] THEN DO; *** If FMETHOD requested in output ***; POWER = FMETHOD || POWER; LBL = CONCAT("FMETHOD_", OPTPOWERCALC1LBL[POWERCALC]) || LBL; END; IF OPTPOWERMAT2[2,1] THEN DO; ***If UCDF requested in output ***; POWER = CDFPOWERCALC || POWER; LBL = CONCAT("UCDF_", OPTPOWERCALC1LBL[POWERCALC]) || LBL; END; IF (POWERCALC = 6) & OPTPOWERMAT2[3,1] THEN DO; *** If UMETHOD requested in output & HF ***; POWER = UMETHODTEMP[1,1] || POWER; LBL = "UMETHOD_HF" || LBL; END; IF (POWERCALC = 7) & OPTPOWERMAT2[3,1] THEN DO; ***IF UMETHOD requested in output & GG ***; POWER = UMETHODTEMP[2,1] || POWER; LBL = "UMETHOD_GG" || LBL; END; POWER = EXEPS || POWER; IF POWERCALC = 5 THEN LBL = "EXEPS_UN" || LBL ; ELSE IF POWERCALC = 6 THEN LBL = "EXEPS_HF" || LBL ; ELSE IF POWERCALC = 7 THEN LBL = "EXEPS_GG" || LBL ; ELSE IF POWERCALC = 8 THEN LBL = "EXEPS_BOX" || LBL ; *** Include epsilon if lowest # <=> first UNIREP test encountered; IF FIRSTUNI = POWERCALC THEN DO; POWER = EPS || POWER; LBL = { "EPSILON"} || LBL; END; FINISH _LASTUNI; ***BELOW MODULE USED FOR UNIVARIATE AND MULTIREP TESTS ONLY***; *************************************************************************; * *; * _GLMMPCL *; * *; * This module computes confidence intervals for noncentrality and *; * power for a General Linear Hypothesis (GLH Ho:C*beta=theta0) test *; * in the General Linear Univariate Model (GLUM: y=X*beta+e, HILE *; * GAUSS), based on estimating the effect, error variance, or neither. *; * Methods from Taylor and Muller (1995). *; * *; * INPUTS *; * *; * 1) Program Defined *; * *; * F_A = MSH/MSE, the F value observed if BETAhat=BETA and *; * Sigmahat=Sigma, under the alternative hypothesis, with: *; * MSH=Mean Square Hypothesis (effect variance) *; * MSE=Mean Square Error (error variance) *; * NOTE: *; * F_A = (N2/N1)*F_EST and *; * MSH = (N2/N1)*MSH_EST, *; * with "_EST" indicating value which was observed *; * in sample 1 (source of estimates) *; * ALPHATEST, Significance level for target GLUM test *; * DFH, degrees of freedom for target GLH; *; * N2, *; * DFE2, Error df for target hypothesis *; * *; * 1) User supplied - required - GLOBAL *; * *; * CLTYPE =1 if Sigma estimated and Beta known *; * =2 if Sigma estimated and Beta estimated *; * N_EST, (scalar) # of observations in analysis which yielded *; * BETA and SIGMA estimates *; * RANK_EST, (scalar) design matrix rank in analysis which *; * yielded BETA and SIGMA estimates *; * *; * 2) User supplied - optional (program defaults supplied) - GLOBAL *; * *; * ALPHA_CL, Lower tail probability for confidence interval *; * ALPHA_CU, Upper tail probability for confidence interval *; * *; * OUTPUTS *; * *; * POWER_L, power confidence interval lower bound *; * POWER_U, power confidence interval upper bound *; * FMETHOD_L, Method used to calculate probability from F CDF *; * used in lower confidence limits power calculation *; * FMETHOD_U, Method used to calculate probability from F CDF *; * used in lower confidence limits power calculation *; * POWERWARN, vector of power calculation warning counts *; * *; *************************************************************************; START _GLMMPCL (POWER_L, POWER_U, FMETHOD_L, FMETHOD_U, POWERWARN, F_A, ALPHATEST, DFH, N2, DFE2) GLOBAL (CLTYPE, N_EST, RANK_EST, ALPHA_CL, ALPHA_CU, TOLERANCE); ***1) Set outputs to missing or zero; ***NONCEN_E = .; NONCEN_L = .; NONCEN_U = . ; ***POWER_E = .; POWER_L = .; POWER_U = .; ***2) input validity checked in main program; ***3) Calculate noncentrality; DFE1 = N_EST - RANK_EST; NONCEN_E = DFH # F_A ; RUN _FINV(FCRIT, 1 - ALPHATEST, DFH, DFE2); ***4) Calculate lower bound for noncentrality; ***CLTYPE (formerly LIMCASE) = 0 is an old option for BETA and SIGMA known; ***IF (CLTYPE = 0) THEN DO; *** NONCEN_L = NONCEN_E; *** GO TO FIVE; *** END; IF (ALPHA_CL <= TOLERANCE) THEN DO; NONCEN_L = 0; GO TO FIVE; END; IF (CLTYPE = 1) THEN DO; *Sigma estimated and Beta fixed; CHI_L = CINV(ALPHA_CL, DFE1); NONCEN_L = (CHI_L / DFE1) # NONCEN_E; GO TO FIVE; END; IF (CLTYPE = 2) THEN DO; *Sigma and Beta estimated; RUN _FINV(BOUND_L, 1 - ALPHA_CL, DFH, DFE1); IF F_A <= BOUND_L THEN NONCEN_L = 0; ELSE NONCEN_L = FNONCT(F_A, DFH, DFE1, 1 - ALPHA_CL); END; FIVE: ; ***5) Calculate lower bound for power; IF (ALPHA_CL <= TOLERANCE) THEN DO; PROB = 1 - ALPHATEST; FMETHOD_L = 5; END; ELSE DO; RUN _PROBF(PROB, FMETHOD_L, FCRIT, DFH, DFE2, NONCEN_L); RUN _FWARN(POWERWARN, FMETHOD_L, 2); END; IF (FMETHOD_L = 4) & (PROB = 1) THEN POWER_L = ALPHATEST; ELSE POWER_L = 1 - PROB; ***6) Calculate upper bound for noncentrality; ***CLTYPE = 0 is an old option for BETA and SIGMA known; ***IF (CLTYPE = 0) THEN DO; *** NONCEN_U = NONCEN_E; *** GO TO SEVEN; *** END; IF (ALPHA_CU <= TOLERANCE) THEN DO; NONCEN_U = .I ; GO TO SEVEN; END; IF (CLTYPE = 1) THEN DO; *Sigma estimated and Beta fixed; CHI_U = CINV(1 - ALPHA_CU, DFE1); NONCEN_U = (CHI_U / DFE1) # NONCEN_E; GO TO SEVEN; END; IF (CLTYPE = 2) THEN DO; *Sigma and Beta estimated; RUN _FINV(BOUND_U, ALPHA_CU, DFH, DFE1); IF F_A <= BOUND_U THEN NONCEN_U = 0; ELSE NONCEN_U = FNONCT(F_A, DFH, DFE1, ALPHA_CU); END; SEVEN: ; ***7) Calculate upper bound for power; IF (ALPHA_CU <= TOLERANCE) THEN DO; PROB = 0; FMETHOD_U = 5; END; ELSE DO; RUN _PROBF(PROB, FMETHOD_U, FCRIT, DFH, DFE2, NONCEN_U); RUN _FWARN(POWERWARN, FMETHOD_U, 3); END; IF (FMETHOD_U = 4) & (PROB = 1) THEN POWER_U = ALPHATEST; ELSE POWER_U = 1 - PROB; ***8) Warning for conservative confidence interval; IF (CLTYPE > 1) & (N2 ^= N_EST) THEN DO; IF (ALPHA_CL > 0) & (NONCEN_L = 0) THEN POWERWARN[5,1] = POWERWARN[5,1] + 1; IF (ALPHA_CL = 0) & (NONCEN_U = 0) THEN POWERWARN[10,1] = POWERWARN[10,1] + 1; END; FINISH _GLMMPCL; *************************************************************************; * *; * _PROBF *; * *; * This module is an updated version of the _PROBFW module. *; * _PROBF calculates Pr(FCRIT < F(df1,df2,noncen)) using one of four *; * methods. The first, most common method uses the cumulative *; * distribution function of the non-central F. If the CDF method will *; * fail, the second method uses the Tiku approximation to the non- *; * central F (Johnson and Kotz). In situations where the TIKU will *; * fail or be inaccurate, we use a Normal approximation (from _PROBF). *; * *; * INPUTS *; * *; * FCRIT, Critical value of F distribution under null hypothesis *; * DF1, Numerator (hypothesis) degrees of freedom *; * DF2, Denominator (error) degrees of freedom; *; * NONCEN, Noncentrality parameter *; * *; * OUTPUTS *; * *; * PROB, Probability that variable distributed F(DF1, DF2, NONCEN) *; * will take a value <= FCRIT *; * FMETHOD, *; * =1, CDF function (no approximation) *; * =2, Tiku approximation (best approximation) *; * =3, Normal approximation, |Z-score| < 6 (worst *; * approximation) *; * =4, Normal approximation, |Z-score| > 6 (approximation *; * but power is almost certainly zero or one) *; * =5, Power missing *; * *; *************************************************************************; START _PROBF (PROB, FMETHOD, FCRIT, DF1, DF2, NONCEN); *** Check for missing FCRIT; IF FCRIT = . THEN DO; PROB = .; FMETHOD = 5; END; ELSE DO; *** Method 1: Use the CDF function; IF (DF1 < 10##4.4) & (DF2 < 10##5.4) & (NONCEN < 10##6.4) THEN DO; PROB = CDF("F", FCRIT, DF1, DF2, NONCEN); FMETHOD = 1; END; ELSE IF (DF2 < 10) & (DF1 < 10##6) & (NONCEN <= 10##6) THEN DO; PROB = CDF("F", FCRIT, DF1, DF2, NONCEN); FMETHOD = 1; END; *** Method 2: Use Tiku approximation; ELSE IF (DF1 < 10##9.2) & (DF2 < 10##9.2) & (NONCEN < 10##6.4) & (DF1 >= 1) & (DF2 > 10##0.6) THEN DO; FMETHOD = 2; H_TIKU = 2 * (DF1 + NONCEN) ## 3 + 3 * (DF1 + NONCEN) * (DF1 + 2 * NONCEN) *(DF2 - 2) + (DF1 + 3 * NONCEN) * (DF2 - 2) ## 2; K_TIKU = (DF1 + NONCEN) ## 2 + (DF2 - 2) * (DF1 + 2 * NONCEN); DF1_TIKU = FLOOR( .5 * (DF2 - 2)*( SQRT (H_TIKU ## 2 / (H_TIKU ## 2 - 4 * K_TIKU ## 3)) - 1 ) ); C_TIKU = (DF1_TIKU / DF1) / (2 * DF1_TIKU + DF2 - 2) * (H_TIKU / K_TIKU); B_TIKU = - DF2 / (DF2 - 2) * (C_TIKU - 1 - NONCEN / DF1); FCRIT_TIKU = (FCRIT - B_TIKU) / C_TIKU; PROB = CDF("F", FCRIT_TIKU, DF1_TIKU, DF2, 0); END; *** Methods 3 and 4: Calculate the Z-score; ELSE DO; P1 = 1 / 3; P2 = -2; P3 = 1 / 2; P4 = 2 / 3; ARG1 = ((DF1 * FCRIT) / (DF1 + NONCEN)); ARG2 = (2 / 9) * (DF1 + (2 * NONCEN)) * ((DF1 + NONCEN) ## P2); ARG3 = (2 / 9) * (1 / DF2); NUMZ = (ARG1 ## P1) - (ARG3 * (ARG1 ## P1)) - (1 - ARG2); DENZ = ((ARG2 + ARG3 * (ARG1 ## P4))) ## P3; ZSCORE = NUMZ / DENZ; *** Method 3: |Z| < 6, so use the normal distribution; IF ABS(ZSCORE) < 6 THEN DO; FMETHOD = 3; PROB = CDF('NORMAL',ZSCORE); END; *** Method 4: |Z| > 6, so power is 0 or 1; ELSE DO; FMETHOD = 4; IF (ZSCORE < -6) THEN PROB = 0; IF (ZSCORE > 6) THEN PROB = 1; END; END; END; FINISH _PROBF; *************************************************************************; * *; * _FINV *; * *; * This module returns the critical value from a central F(DF1, DF2) *; * distribution for a given signficance level alpha. It screens the *; * FINV function for numerator DF greater than 1*10^7.6 or denominator *; * DF greater than 1*10^9.4. For large degrees of freedom, it returns *; * a missing value. *; * *; * INPUTS *; * *; * ALPHA, Signficance level *; * DF1, Numerator degrees of freedom *; * DF2, Denominator degrees of freedom *; * *; * OUTPUTS *; * *; * FCRIT, Critical value from the probability that a variable *; * distributed F(DF1,DF2) <= FCRIT is equal to (1-ALPHA) *; * *; *************************************************************************; START _FINV (FCRIT, ALPHA, DF1, DF2); IF (DF1 > 10**7.6) | (DF2 > 10**9.4) | (DF1 < 0) | (DF2 < 0) THEN FCRIT = .; ELSE FCRIT = FINV(ALPHA, DF1, DF2); FINISH _FINV; *************************************************************************; * *; * _FWARN *; * *; * This module updates the appropriate warning counts in POWERWARN based *; * on the value of FMETHOD produced by the PROBF module (used to *; * calculate powers). *; * *; * INPUTS *; * *; * 1) Program supplied *; * *; * FMETHOD, Method used to calculate probability from F CDF *; * CL, =1 if calculation of power of a test *; * =2 if calculation of lower CL for power of a test *; * =3 if calculation of upper CL for power of a test *; * *; * OUTPUTS *; * *; * POWERWARN, vector of power calculation warning counts *; * *; *************************************************************************; START _FWARN (POWERWARN, FMETHOD, CL); IF CL = 1 THEN DO; IF FMETHOD = 2 THEN POWERWARN[1,1] = POWERWARN[1,1] + 1; ELSE IF FMETHOD = 3 THEN POWERWARN[2,1] = POWERWARN[2,1] + 1; ELSE IF FMETHOD = 4 THEN POWERWARN[3,1] = POWERWARN[3,1] + 1; ELSE IF FMETHOD = 5 THEN POWERWARN[4,1] = POWERWARN[4,1] + 1; END; ELSE IF CL = 2 THEN DO; IF FMETHOD = 2 THEN POWERWARN[6,1] = POWERWARN[6,1] + 1; ELSE IF FMETHOD = 3 THEN POWERWARN[7,1] = POWERWARN[7,1] + 1; ELSE IF FMETHOD = 4 THEN POWERWARN[8,1] = POWERWARN[8,1] + 1; ELSE IF FMETHOD = 5 THEN POWERWARN[9,1] = POWERWARN[9,1] + 1; END; ELSE IF CL = 3 THEN DO; IF FMETHOD = 2 THEN POWERWARN[11,1] = POWERWARN[11,1] + 1; ELSE IF FMETHOD = 3 THEN POWERWARN[12,1] = POWERWARN[12,1] + 1; ELSE IF FMETHOD = 4 THEN POWERWARN[13,1] = POWERWARN[13,1] + 1; ELSE IF FMETHOD = 5 THEN POWERWARN[14,1] = POWERWARN[14,1] + 1; END; FINISH _FWARN; *************************************************************************; * *; * _SASDS *; * *; * Creates SAS dataset if requested. *; * *; * INPUTS *; * *; * 1) User supplied (or program default) *; * *; * DSNAME, { "libref" "dataset name" "default library" } *; * *; * The "default library" is optional. If omitted WORK *; * is used. Default {WORK DODFAULT}. *; * *; * 2) Program supplied *; * *; * _HOLDPOW, matrix of power values and optional output. *; * _HOLDNM, matrix of variable names for the SAS dataset. *; * OPTWARNPRINT=1 if warnings to be printed and 0 otherwise. *; * *; *************************************************************************; START _SASDS(_HOLDPOW, _HOLDNM, DSNAME,OPTWARNPRINT); IF NCOL(DSNAME) = 2 THEN DSNAME = DSNAME || {WORK}; IF NCOL(DSNAME) = 3 THEN DSNAME = DSNAME; LIB = DSNAME[1,1]; DATASET = DSNAME[1,2]; DEFAULT = DSNAME[1,3]; *** Reset default library to libref; RESET NONAME DEFLIB = LIB; LISTDS = DATASETS(LIB); ENDIT=0; NUMDS = NROW(LISTDS); ***Check to see if _PWRDTMP or user-specified DATASET already exists; IF NUMDS > 0 THEN DO CNT = 1 TO NUMDS; IF LISTDS[CNT,1] = "_PWRDTMP" THEN DO; ENDIT=1; NOTE = CONCAT("WARNING 16: The program uses an ", "intermediate dataset called _PWRDTMP. This ", "dataset already exists in the requested ", "library. New datasets cannot be created."); IF OPTWARNPRINT THEN PRINT NOTE[LABEL=" " ROWNAME=" " COLNAME=" "]; END; IF LISTDS[CNT,1] = DATASET THEN DO; NOTE1 = COMPBL(CONCAT("Dataset ", LIB, ".", DATASET, " already exists. The default ", "PWRDT### will be created instead. ")); NOTE2 = COMPBL(CONCAT("To overwrite dataset ", LIB, ".", DATASET, " specify the statement ", "CALL DELETE(", LIB, ", ", DATASET, "); prior to the RUN ", "POWER; statement.")); NOTE = CONCAT(NOTE1, NOTE2); IF OPTWARNPRINT THEN PRINT ,, NOTE; DATASET = "DODFAULT"; END; END; *** Set default dataset name if required; NEWNUM = 0; IF DATASET = "DODFAULT" THEN NEWNUM = 1; IF DATASET = "DODFAULT" & NUMDS > 0 THEN DO; DO CNT2 = 1 TO NUMDS; ***Are any PWRDT### datasets in the library ?; IF SUBSTR(LISTDS[CNT2,1],1,5) = "PWRDT" THEN LISTPDS = LISTPDS // LISTDS[CNT2,1]; END; ***What numbers do PWRDT### datasets have? Set number for new DS.**; IF NROW(LISTPDS) > 0 THEN DO; PDSNUMS = NUM(SUBSTR(LISTPDS,6,3)); MAXNUM = MAX(PDSNUMS); NEWNUM = MAXNUM + 1; END; ***Maximum number of PWRDT### datasets is 999**; IF NEWNUM > 999 THEN DO; ENDIT = 1; NOTE = CONCAT("There are already 999 PWRDT### datasets. ", "No more can be created."); IF OPTWARNPRINT THEN PRINT NOTE[LABEL=" " ROWNAME=" " COLNAME=" "]; END; END; *** New default name; IF (DATASET = "DODFAULT" & (1 <= NEWNUM <= 999) & ENDIT ^= 1) THEN DATASET = COMPRESS(CONCAT("PWRDT",CHAR(NEWNUM))); *** Create intermediate dataset called _PWRDTMP; IF ENDIT ^= 1 THEN DO; CREATE _PWRDTMP VAR _HOLDNM; APPEND FROM _HOLDPOW; CLOSE _PWRDTMP; ***Change name of intermediate dataset to user specified or default; CALL RENAME(LIB,"_PWRDTMP",DATASET); NAMED = COMPRESS(CONCAT(LIB,".",DATASET)); NOTE = COMPBL(CONCAT("The results have been stored in SAS file ", NAMED)); IF OPTWARNPRINT THEN PRINT NOTE[LABEL=" " ROWNAME=" " COLNAME=" "]; END; *** Reset to original default library; RESET NAME DEFLIB = DEFAULT; FINISH _SASDS; *************************************************************************; * *; * _QPROB *; * *; * This function returns the CDF of a weighted sum of independent *; * chi squares via Davies' algorithm. *; * *; * INPUTS *; * *; * 1) Program supplied *; * *; * Q, point at which CDF is evaluated: Prob{Q <= q} *; * Note: Q is often zero if start from ratio of positively *; * weighted forms *; * LAMBDA, Nx1 vector of coefficients *; * NU, Nx1 vector of degrees of freedom *; * OMEGA, Nx1 vector of noncentralities *; * ACCURACY, maximum error in probability allowed *; * *; * RETURNS *; * *; * PROB, Prob{Q <= q}, the CDF *; * This is set to zero if any problems occur *; * *; *************************************************************************; START _QPROB ( Q, LAMBDA, NU, OMEGA, ACCURACY); N = NROW(LAMBDA); IF MAX(NCOL(LAMBDA), NCOL(NU), NCOL(OMEGA)) > 1 THEN RETURN(.); IF (NROW(NU) ^= N) | (NROW(OMEGA) ^= N) THEN RETURN(.); LIMIT1 = 50000; ***Steps allowed; SIGMAC = 0; ***Coefficient on Gaussian; PRNTPROB = "OFF"; ERRORCHK = "OFF"; RUN _AS(PROB, TRACE, ICOUNT, IFAULT, N, LIMIT1, LAMBDA, SIGMAC, Q, ACCURACY, OMEGA, NU, PRNTPROB, ERRORCHK); IF IFAULT > 0 THEN DO; * PRINT PROB LAMBDA NU OMEGA IFAULT TRACE ICOUNT; RETURN(.); END; RETURN(PROB); FINISH _QPROB; *******Following modules, _AS through _CFE, needed for _QPROB module*****; *************************************************************************; * *; * _AS *; * *; * Computes distribution of a linear combination of non-central chi- *; * squared random variables. Taken from Algorithm AS 155, Applied *; * Statistics (1980), Vol 29, No 3. *; * *; * INPUTS *; * *; * 1) Program supplied *; * *; * IRR, Number of chi-squared terms in the sum *; * LIM1, Maximum number of integration terms *; * ALB, IRRx1 vector of constant multipliers *; * SIGMA, Coefficient of normal term *; * CC, Point at which the distribution function should be *; * evaluated *; * ACC, Error bound *; * ANC, Vector of noncentrality parameters *; * N, Vector of degrees of freedom *; * PRINTPROB, ="ON" or ="OFF", "ON" prints QF *; * ERRORCHK, ="ON" or ="OFF", "ON" prints ICOUNT, TRACE, IFAULT *; * *; * OUTPUTS *; * *; * QF, Probability that the quadratic form is less than CC *; * TRACE, 7x1 vector of variables that indicate the performance *; * of the procedure *; * TRACE[1] = Absolute value sum *; * TRACE[2] = Total number of integration terms *; * TRACE[3] = Number of integrations *; * TRACE[4] = Integration interval in main integration *; * TRACE[5] = Truncation point in initial integration *; * TRACE[6] = Standard deviation of convergence factor term *; * TRACE[7] = Number of cycles to locate integration parameters *; * IFAULT, Output fault indicator *; * =0 No error *; * =1 Requested accuracy could not be obtained *; * =2 Round-off error possibly significant *; * =3 Invalid parameters *; * =4 Unable to location integration parameters *; * ICOUNT, Number of times the function was called *; * *; *************************************************************************; START _AS (QF, TRACE, ICOUNT, IFAULT, IRR, LIM1, ALB, SIGMA, CC, ACC, ANC, N, PRNTPROB, ERRORCHK); ***Defines constants; ALN28 = LOG(2)/8; PI = 2#ARCOS(0); NDTSRT = "TRUE"; FAIL = "FALSE"; ***Initialize variables; IFAULT = 0; ICOUNT = 0; AINTL = 0; ERSM=0; QF = -1; TRACE = J(7,1,0); ***Produce local copies of some variables; C = CC; IR = IRR; LIM = LIM1; ACC1 = ACC; *** AMEAN, Scalar representing the expected value of Q; *** SD, Scalar representing the squared deviation of Q- the second moment; *** ALMAX, Maximum of the constants; *** ALMIN, Minimum of the constants; XLIM = LIM; SIGSQ = SIGMA**2; SD = SIGSQ; ALMAX = 0; ALMIN = 0; AMEAN = 0; J = 1; L20: IF (^(J<=IR)) THEN GOTO L60; NJ = N[J]; ALJ = ALB[J]; ANCJ = ANC[J]; IF (^((NJ<0) | (ANCJ<0))) THEN GOTO L30; IFAULT = 3; GOTO L260; L30: SD = SD + (ALJ**2)*(2*NJ + 4*ANCJ); AMEAN = AMEAN + ALJ*(NJ + ANCJ); IF (^(ALMAX < ALJ)) THEN GOTO L40; ALMAX = ALJ; GOTO L50; L40: IF (^(ALMIN > ALJ)) THEN GOTO L50; ALMIN = ALJ; L50: J = J + 1; GOTO L20; L60: IF (^(SD = 0)) THEN GOTO L80; IF (^(C = 0)) THEN GOTO L70; QF = 1; GOTO L260; L70: QF = 0; GOTO L260; L80: IF (^((ALMIN=0) & (ALMAX= 0) & (SIGMA= 0))) THEN GOTO L90; IFAULT = 3; GOTO L260; L90: SD = SQRT(SD); IF (^(ALMAX < -ALMIN)) THEN GOTO L100; ALMX = -ALMIN; GOTO L110; L100: ALMX = ALMAX; *** Define starting values for modules _FINDU and _CTFF; L110: UTX = 16# INV(SD); UP = 4.5 # INV(SD); UN = -1#UP; *** Calculate the Truncation point without any convergence factor; CALL _FINDU(UTX, N, ALB, ANC, (.5#ACC1), LIM, ICOUNT, SIGSQ, IR); *** Does convergence factor help?; IF ^((C ^= 0) & (ALMX > .07 * SD)) THEN GOTO L130; CALL _CFE(FAIL, CFE, N, ALB, ANC, ITH, C, LIM, ICOUNT, NDTSRT, IR); TAUSQ = .25 * ACC1 / CFE; IF (FAIL = "FALSE") THEN GOTO L120; FAIL = "FALSE"; GOTO L130; L120: IF (^(_TRUNCN(N, ALB, ANC, UTX, TAUSQ, LIM, ICOUNT, SIGSQ, IR) < 0.2 * ACC1)) THEN GOTO L130; SIGSQ = SIGSQ + TAUSQ; CALL _FINDU(UTX, N, ALB, ANC, 0.25 * ACC1, LIM, ICOUNT, SIGSQ, IR); TRACE[6] = SQRT(TAUSQ); L130: TRACE[5] = UTX; ACC1 = .5 * ACC1; *** Find 'range' of distribution; *** Quit if outside of this; L140: CALL _CTFF(UP, FCTFF, N, ALB, ANC, ACC1, AMEAN, ALMIN, ALMAX, LIM, ICOUNT, SIGSQ, IR); D1 = FCTFF - C; IF (^(D1 < 0)) THEN GOTO L150; QF = 1; GOTO L260; L150: CALL _CTFF(UN, FCTFF, N, ALB, ANC, ACC1, AMEAN, ALMIN, ALMAX, LIM, ICOUNT, SIGSQ, IR); D2 = C - FCTFF; IF (^(D2 < 0)) THEN GOTO L160; QF = 0; GOTO L260; *** Find integration interval; L160: IF (^(D1 > D2)) THEN GOTO L170; AINTV = D1; GOTO L180; L170: AINTV = D2; L180: AINTV = 2 * (2 # ARCOS(0)) / AINTV; *** Calculate number of terms required for main and auxillary integrations; XNT = UTX / AINTV; XNTM = 3.0 / SQRT(ACC1); IF (^(XNT > XNTM*1.5)) THEN GOTO L220; IF (^(XNTM > XLIM)) THEN GOTO L190; IFAULT = 1; GOTO L260; *** Parameters for auxiliary integration; L190: NTM = _IROUND(XNTM); AINTV1 = UTX / XNTM; X = 2*(2 # ARCOS(0)) / AINTV1; IF (^(X <= ABS(C))) THEN GOTO L200; GOTO L220; *** Calculate convergence factor; L200: CALL _CFE(FAIL, CX1, N, ALB, ANC, ITH, C - X, LIM, ICOUNT, NDTSRT, IR); CALL _CFE(FAIL, CX2, N, ALB, ANC, ITH, C + X, LIM, ICOUNT, NDTSRT, IR); TAUSQ = CX1 + CX2; TAUSQ = .33 * ACC1 / (1.1*TAUSQ); IF (FAIL = "FALSE") THEN GOTO L210; GOTO L220; L210: ACC1 = 0.67 * ACC1; *** Auxiliary integration; CALL _INTEGR (AINTL, ERSM, N, ALB, ANC, NTM, AINTV1, TAUSQ, "FALSE", C, SIGSQ, IR); XLIM = XLIM - XNTM; SIGSQ = SIGSQ + TAUSQ; TRACE[3] = TRACE[3] + 1; TRACE[2] = TRACE[2] + NTM + 1; *** Find truncation point with new convergence factor; CALL _FINDU(UTX, N, ALB, ANC, 0.25 * ACC1, LIM, ICOUNT, SIGSQ, IR); ACC1 = 0.75 * ACC1; GOTO L140; *** Main integration; L220: TRACE[4] = AINTV; IF (^(XNT > XLIM)) THEN GOTO L230; IFAULT = 1; GOTO L260; L230: NT = _IROUND(XNT); CALL _INTEGR(AINTL, ERSM, N, ALB, ANC, NT, AINTV, 0, "TRUE", C, SIGSQ, IR); TRACE[3] = TRACE[3] + 1; TRACE[2] = TRACE[2] + NT + 1; QF = .5 - AINTL; TRACE[1] = ERSM; UP = ERSM; *** Test whether round-off error could be significant; *** Allow for radix 8 or 16 machines; X = UP + (ACC / 10.0); J = 1; L240: IF (^(J <= 8)) THEN GOTO L260; IF (^(J * X = J*UP)) THEN GOTO L250; IFAULT = 2; L250: J = J*2; GOTO L240; L260: TRACE[7] = ICOUNT; *** Print output; IF PRNTPROB= "ON" THEN DO; PRINT "The probability that the quadratic form is less than C is"; PRINT QF; END; *** Print analysis of run; IF ERRORCHK = "ON" THEN DO; PRINT "TRACE(7) Output array of 7 variables that ", " indicate the performance of the ", " procedure. ", " 1: absolute value sum ", " 2: total number of integration terms ", " 3: number of integrations ", " 4: integration interval in main ", " integration. ", " 5: truncation point in initial ", " integration ", " 6: standard deviation of convergence ", " factor term ", " 7: number of cycles to locate ", " integration parameters ", ,,TRACE; PRINT / " IFAULT Output fault indicator: ", " 0: No error ", " 1: Requested accuracy could not ", " be obtained. ", " 2: Round-off error possibly ", " significant ", " 3: Invalid parameters ", " 4: Unable to locate integration ", " parameters ", ,,IFAULT; PRINT " ICOUNT Number of times the evaluation was called ", ,,ICOUNT; END; FINISH _AS; *************************************************************************; * *; * _IROUND *; * *; * This function module rounds numbers to the nearest integer. *; * *; * INPUT *; * *; * X, Number to round *; * *; * RETURNS *; * *; * ANSWER, Rounded value *; * *; *************************************************************************; START _IROUND (X); IF (X >= 0) THEN DO; ANSWER = INT(X + .5); END; IF (X < 0) THEN DO; ANSWER = INT(X - .5); END; RETURN (ANSWER); FINISH _IROUND; *************************************************************************; * *; * _COUNTR *; * *; * This module counts number of calls to _ERRBD, _TRUNCN, and _CFE. *; * *; * INPUT *; * *; * LIM, Maximum number of integration terms *; * *; * OUTPUT *; * *; * ICOUNT, count of number of calls *; * *; *************************************************************************; START _COUNTR (ICOUNT, LIM); IF ICOUNT > LIM THEN DO; PRINT "QF: CANNOT LOCATE INTEGRATION PARAMETERS"; ABORT; END; ELSE ICOUNT = ICOUNT + 1; FINISH _COUNTR; *************************************************************************; * *; * _ALOG1 *; * *; * For a number X, this function computes either LN(1 + X) *; * or LN(1 + X) - X . *; * *; * INPUTS *; * *; * X, Number on which to perform computation *; * FIRST, ="TRUE", when this module is called for a matrix for *; * the first time *; * ="FALSE", otherwise *; * *; * RETURNS *; * *; * LN(1 + X), if FIRST = "TRUE" *; * LN(1 + X) - X, if FIRST = "FALSE" *; * *; *************************************************************************; START _ALOG1 ( X, FIRST); IF (^(ABS(X) > .1)) THEN GOTO L20; IF (FIRST = "FALSE") THEN GOTO L10; ALOG1 = LOG(1 + X); GOTO L70; L10: ALOG1 = LOG(1 + X) - X; GOTO L70; L20: Y = X / (2 + X); TERM = 2 * Y ** 3; AK = 3.0; IF (FIRST = "FALSE") THEN GOTO L30; S = 2; GOTO L40; L30: S = -X; L40: S = S * Y; Y = Y ** 2; S1 = S + TERM / AK; L50: IF (^(S1 ^= S)) THEN GOTO L60; AK = AK + 2; TERM = TERM * Y; S = S1; S1 = S + TERM / AK; GOTO L50; L60: ALOG1 = S; L70: RETURN (ALOG1); FINISH _ALOG1; *************************************************************************; * *; * _EXP1 *; * *; * This function computes e^X for X > -706. *; * *; * INPUT *; * X, scalar *; * *; * RETURNS *; * *; * 0, if X <= -706 *; * e^X, if X > -706 *; * *; *************************************************************************; START _EXP1 (X); IF X <= -706 THEN RETURN(0); ELSE RETURN(EXP(X)); FINISH _EXP1; *************************************************************************; * *; * _ORDER *; * *; * This module finds the ranks of absolute values of elements of ALB. *; * Ties are ranked arbitrarily, e.g., the matrix {2,2} is ranked {1,2}. *; * *; * INPUT *; * *; * ALB, IRx1 vector of constant multipliers *; * *; * OUTPUTS *; * *; * ITH, Vector of ranks of absolute values of ALB *; * NDTSRT, ="FALSE" if this module has been run *; * *; *************************************************************************; START _ORDER (ITH, NDTSRT, ALB); ITH = RANK(ABS(ALB)); NDTSRT = "FALSE"; FINISH _ORDER; *************************************************************************; * *; * _ERRBD *; * *; * This module finds bound on tail probability using moment-generating *; * function *; * *; * INPUTS *; * *; * N, Vector of degrees of freedom *; * ALB, IRx1 vector of constant multipliers *; * ANC, Vector of noncentrality parameters *; * UU, *; * LIM, Maximum number of integration terms *; * ICOUNT, Count of number of times this module is called *; * SIGSQ, square of SIGMA, the coefficient of normal term *; *; * IR, Number of chi-squared terms in the sum *; * *; * *; * OUTPUTS *; * *; * ERRBD, Bound on tail probability *; * CX, *; * *; *************************************************************************; START _ERRBD (ERRBD, CX, N, ALB, ANC, UU, LIM, ICOUNT, SIGSQ, IR); CALL _COUNTR(ICOUNT, LIM); U = UU; CONST = U * SIGSQ; SUM1 = U * CONST; U = 2 * U; J = IR; L10: IF J <= 0 THEN GOTO L20; NJ = N[J]; ALJ = ALB[J]; ANCJ = ANC[J]; X = U # ALJ; Y = 1 - X; CONST = CONST + ALJ # (ANCJ / Y + NJ) / Y; SUM1 = SUM1 + ANCJ # ((X / Y) ## 2); SUM1 = SUM1 + NJ # ((X ## 2) / Y + _ALOG1(-X, "FALSE")); J = J-1; GO TO L10; L20: ERRBD = _EXP1(-.5 * SUM1); CX = CONST; FINISH _ERRBD; *************************************************************************; * *; * _CTFF *; * *; * This module finds CTF so that: *; * If UPN > 0: P(QF > CTFF) < ACCX *; * Otherwise: P(QF < CTFF) < ACCX *; * *; * INPUTS *; * *; * N, Vector of degrees of freedom *; * ALB, IRx1 vector of constant multipliers *; * ANC, Vector of noncentrality parameters *; * ACCX, Error bound *; * AMEAN, Scalar representing the expected value of the QF *; * ALMIN, Minimum of the constant multipliers *; * ALMAX, Maximum of the constant multipliers *; * LIM, Maximum number of integration terms *; * ICOUNT, Count of number of times this module is called *; * SIGSQ, square of SIGMA, the coefficient of normal term *; * IR, Number of chi-squared terms in the sum *; * *; * OUTPUTS *; * *; * UPN, *; * FCTFF, *; * *; *************************************************************************; START _CTFF (UPN, FCTFF, N, ALB, ANC, ACCX, AMEAN, ALMIN, ALMAX, LIM, ICOUNT, SIGSQ, IR); U2 = UPN; U1 = 0; C1 = AMEAN; IF (^(U2 > 0)) THEN GOTO L10; RB = ALMAX; GOTO L20; L10: RB = ALMIN; L20: RB = 2 * RB; U = U2 / (1 + U2 * RB); L30: CALL _ERRBD(ERRBD, C2, N, ALB, ANC, U, LIM, ICOUNT, SIGSQ, IR); IF ^(ERRBD > ACCX) THEN GOTO L40; U1 = U2; C1 = C2; U2 = 2 * U2; U = U2 / (1 + U2 * RB); GOTO L30; L40: U = (C1 - AMEAN) / (C2 - AMEAN); L50: IF (^(U < 0.9)) THEN GOTO L80; U = (U1 + U2) / 2; CALL _ERRBD(ERRBD, CONST, N, ALB, ANC, U/(1 + U * RB), LIM, ICOUNT, SIGSQ, IR); IF (^(ERRBD > ACCX)) THEN GOTO L60; U1 = U; C1 = CONST; GOTO L70; L60: U2 = U; C2 = CONST; L70: U = (C1 - AMEAN) / (C2 - AMEAN); GOTO L50; L80: FCTFF = C2; UPN = U2; FINISH _CTFF; *************************************************************************; * *; * _TRUNCN *; * *; * This function bounds integration error due to truncation at U. *; * *; * INPUTS *; * *; * N, Vector of degrees of freedom *; * ALB, IRx1 vector of constant multipliers *; * ANC, Vector of noncentrality parameters *; * UU, *; * TAUSQ, *; * LIM, Maximum number of integration terms *; * ICOUNT, Count of number of times this module is called *; * SIGSQ, square of SIGMA, the coefficient of normal term *; * IR, Number of chi-squared terms in the sum *; * *; * RETURNS *; * *; * TRUNCN, integration error *; * *; *************************************************************************; START _TRUNCN( N, ALB, ANC, UU, TAUSQ, LIM, ICOUNT, SIGSQ, IR); CALL _COUNTR(ICOUNT, LIM); U = UU; SUM1 = 0; PROD2 = 0; PROD3 = 0; NS = 0; SUM2 = (SIGSQ + TAUSQ)*U**2; PROD1 = 2*SUM2; U = 2*U; J = 1; L10: IF (^(J <= IR)) THEN GOTO L40; ALJ = ALB[J]; ANCJ = ANC[J]; NJ = N[J]; X = (U * ALJ) ** 2; SUM1 = SUM1 + ANCJ * X / (1 + X); IF (^(X > 1)) THEN GOTO L20; PROD2 = PROD2 + NJ * LOG(X); PROD3 = PROD3 + NJ * _ALOG1(X,"TRUE"); NS = NS + NJ; GOTO L30; L20: PROD1 = PROD1 + NJ * _ALOG1(X,"TRUE"); L30: J = J + 1; GOTO L10; L40: SUM1 = .5 * SUM1; PROD2 = PROD1 + PROD2; PROD3 = PROD1 + PROD3; X = (_EXP1(-SUM1 - 0.25 * PROD2)) / (2 # ARCOS(0)); Y = (_EXP1(-SUM1 - 0.25 * PROD3)) / (2 # ARCOS(0)); IF (^(NS = 0)) THEN GOTO L50; ERR1 = 1; GOTO L60; L50: ERR1 = X * 2 / NS; L60: IF (^(PROD3 > 1)) THEN GOTO L70; ERR2 = 2.5 * Y; GOTO L80; L70: ERR2 = 1; L80: IF (^(ERR2 ACCX. *; * *; * INPUTS *; * *; * N, Vector of degrees of freedom *; * ALB, IRx1 vector of constant multipliers *; * ANC, Vector of noncentrality parameters *; * ACCX, *; * LIM, Maximum number of integration terms *; * ICOUNT, Count of number of times this module is called *; * SIGSQ, square of SIGMA, the coefficient of normal term *; * IR, Number of chi-squared terms in the sum *; * *; * OUTPUTS *; * *; * UTX, *; * *; *************************************************************************; START _FINDU (UTX, N, ALB, ANC, ACCX, LIM, ICOUNT, SIGSQ, IR); DIVIS = {2.0, 1.4, 1.2, 1.1}; UT = UTX; U = UT / 4; IF (^(_TRUNCN(N, ALB, ANC, U, 0, LIM, ICOUNT, SIGSQ, IR) > ACCX)) THEN GOTO L20; U = UT; L10: IF (^(_TRUNCN(N, ALB, ANC, U, 0, LIM, ICOUNT, SIGSQ, IR) > ACCX)) THEN GOTO L40; UT = UT * 4; U = UT; GOTO L10; L20: UT = U; U = U / 4; L30: IF (^(_TRUNCN(N, ALB, ANC, U, 0, LIM, ICOUNT, SIGSQ, IR) <= ACCX)) THEN GOTO L40; UT = U; U = U / 4; GOTO L30; L40: DO I = 1 TO 4; U = UT / DIVIS[I]; IF (^(_TRUNCN(N, ALB, ANC, U, 0, LIM, ICOUNT, SIGSQ, IR) <= ACCX)) THEN GOTO L50; UT = U; END; L50: UTX = UT; FINISH _FINDU; *************************************************************************; * *; * _INTEGR *; * *; * INPUTS *; * *; * N, Vector of degrees of freedom *; * ALB, IRx1 vector of constant multipliers *; * ANC, Vector of noncentrality parameters *; * NTERM, Number of terms in integration *; * AINTRV, *; * TAUSQ, *; * MAIN, *; * C, Point at which the distribution function should be *; * evaluated *; * SIGSQ, square of SIGMA, the coefficient of normal term *; * IR, Number of chi-squared terms in the sum *; * *; * OUTPUTS *; * *; * AINTL, *; * ERSM, *; * *; *************************************************************************; START _INTEGR (AINTL, ERSM, N, ALB, ANC, NTERM, AINTRV, TAUSQ, MAIN, C, SIGSQ, IR); PI = 2 # ARCOS(0); AINPI = AINTRV/PI; K = NTERM; L10: IF (^(K >= 0)) THEN GOTO L50; U = (K + .5) * AINTRV; SUM1 = -2 * U * C; SUM2 = ABS(SUM1); SUM3 = -.5 * SIGSQ * (U ** 2); J = IR; L20: IF (^(J>0)) THEN GOTO L30; NJ = N[J,1]; X = 2 * ALB[J,1] * U; Y = X ** 2; SUM3 = SUM3 - 0.25 * NJ * _ALOG1(Y, "TRUE"); Y = ANC[J,1] * X / (1 + Y); Z = NJ * ATAN(X) + Y; SUM1 = SUM1 + Z; SUM2 = SUM2 + ABS(Z); SUM3 = SUM3 - .5*X*Y; J = J - 1; GOTO L20; L30: X = AINPI * _EXP1(SUM3) / U; IF (MAIN = "TRUE") THEN GOTO L40; X = X * (1 - _EXP1(-.5 * TAUSQ * U**2)); L40: SUM1 = SIN(.5 * SUM1) * X; SUM2 = .5 * SUM2 * X; AINTL = AINTL + SUM1; ERSM = ERSM + SUM2; K = K - 1; GOTO L10; L50: ; FINISH _INTEGR; *************************************************************************; * *; * _CFE *; * *; * This module computes the coefficient of TAUSQ in error when *; * the convergence factor of Exp(-0.5 * TAUSQ * U**2) is used when DF *; * is evaluated at X. *; * *; * INPUTS *; * *; * N, Vector of degrees of freedom *; * ALB, IRx1 vector of constant multipliers *; * ANC, Vector of noncentrality parameters *; * ITH, Vector of ranks of absolute values of ALB *; * X, *; * LIM, Maximum number of integration terms *; * ICOUNT, Count of number of times this module is called *; * NDRSRT, ="TRUE" if _ORDER module has not been run *; * ="FALSE" if _ORDER module has been run *; * IR, Number of chi-squared terms in the sum *; * *; * OUTPUTS *; * *; * FAIL, ="TRUE" if module produces unreasonable values *; * FCFE, Coefficient of TAUSQ *; * *; *************************************************************************; START _CFE (FAIL, FCFE, N, ALB, ANC, ITH, X, LIM, ICOUNT, NDTSRT, IR); CALL _COUNTR(ICOUNT, LIM); IF (NDTSRT = "FALSE") THEN GOTO L10; CALL _ORDER(ITH, NDTSRT, ALB); L10: AXL = ABS(X); SXL = X / ABS(X); SUM1 = 0; J = IR; L20: IF (^(J > 0)) THEN GOTO L70; IT = ITH[J]; IF (^(ALB[IT] * SXL > 0)) THEN GOTO L60; ALJ = ABS(ALB[IT]); AXL1 = AXL - ALJ*(N[IT] + ANC[IT]); ALN28 = LOG(2) / 8; AXL2 = ALJ / ALN28; IF (^(AXL1 > AXL2)) THEN GOTO L30; AXL = AXL1; GOTO L60; L30: IF (^(AXL > AXL2)) THEN GOTO L40; AXL = AXL2; L40: SUM1 = (AXL - AXL1) / ALJ; K = J - 1; L50: IF (^(K > 0)) THEN GOTO L70; ITK = ITH[K]; SUM1 = SUM1 + (N[ITK] + ANC[ITK]); K = K - 1; GOTO L50; L60: J = J - 1; GOTO L20; L70: IF (^(SUM1 > 100)) THEN GOTO L80; FCFE = 1; FAIL = "TRUE"; GOTO L90; L80: FCFE = 2 ** (SUM1 / 4) / ((2 # ARCOS(0)) * AXL ** 2); L90: ; FINISH _CFE; *** Modules NAMELIST, UTREND, UMEAN, UPOLY1, UPOLY2 taken from ZUTIL in LINMOD; *************************************************************************; * *; * NAMELIST *; * *; * This function generates row of names, STEMlow to STEMhigh, by BY. *; *; * *; * INPUTS *; * *; * STEM, 1x1 character string to from prefix in names *; * LOW, 1x1 scalar at which to start naming suffixes *; * HIGH, 1x1 scalar at which to end naming suffixes *; *; * BY, 1x1 scalar by which to increment naming suffixes *; * Require 0 <= LOW <= HIGH and 1 <= BY *; * *; * RETURNS *; * *; * LIST, list of names *; * *; *************************************************************************; START NAMELIST ( STEM, LOW, HIGH, BY); IF (NROW(LOW )^=1) | (NCOL(LOW )^=1) | (NROW(HIGH)^=1) | (NCOL(HIGH)^=1) | (NROW(BY )^=1) | (NCOL(BY )^=1) THEN RETURN("?"); ARGS = LOW // HIGH // BY; IF MAX(ABS(ROUND(ARGS, 1) - ARGS)) > 0 THEN RETURN("?"); IF (HIGH < LOW) | (LOW < 0) | (BY < 1) < 0 THEN RETURN("?"); ***IML does not create DO(1,1,1) ; ***Code below avoids limitation; IF HIGH = LOW THEN VALUES = LOW; ELSE VALUES = DO(LOW, HIGH, BY); LIST = COMPRESS(CONCAT(STEM, CHAR(VALUES ))); RETURN(LIST); FINISH NAMELIST; *************************************************************************; * *; * UMEAN *; * *; * This module creates a Px1 U averager matrix. *; * *; * INPUTS *; * *; * P, Number of outcome variables *; * *; * RETURNS *; * *; * UAVE, Px1 matrix *; * *; *************************************************************************; START UMEAN (P); UAVE = J(P, 1, 1/P); RETURN(UAVE); FINISH UMEAN; *************************************************************************; * *; * UPOLY1 *; * *; * This module creates a U contrast matrix with orthogonal polynomial *; * coding for one within subject factor. *; * *; * INPUTS *; * *; * VALUES, Set of numeric treatment levels/values *; * NAME, Character string providing stem of labels for trends *; * *; * OUTPUTS *; * *; * U, Matrix of orthonormal polynomial coefficient columns *; * ULBL, Matrix of labels for columns of U *; * *; *************************************************************************; START UPOLY1 (VALUES, NAME, U, ULBL ); IF (NROW(VALUES) = 0) | (NROW(NAME) = 0) THEN DO; NOTE = CONCAT("ERROR 66: One or more inputs to UPOLY2 has ", "not been assigned a value."); PRINT NOTE[LABEL=" " ROWNAME=" " COLNAME=" "]; ABORT; END; IF ((NCOL(VALUES) > 1) & (NROW(VALUES) > 1)) THEN DO; NOTE = CONCAT("ERROR 67: VALUES has more than one column ", "and more than one row. It must be a row or ", "column vector."); PRINT NOTE[LABEL=" " ROWNAME=" " COLNAME=" "]; ABORT; END; IF (NCOL(NAME) > 1) | (NROW(NAME) > 1) THEN DO; NOTE = CONCAT("ERROR 68: NAME has more than one column or more ", "than one row. It must be a scalar (1x1)."); PRINT NOTE[LABEL=" " ROWNAME=" " COLNAME=" "]; ABORT; END; IF NCOL(VALUES) > 1 THEN V = VALUES`; ELSE V = VALUES ; IF (TYPE(VALUES)^="N") | (TYPE(NAME)^="C") THEN DO; NOTE = CONCAT("ERROR 69: NAMES not character or VALUES not ", "numeric."); PRINT NOTE[LABEL=" " ROWNAME=" " COLNAME=" "]; ABORT; END; P = NROW(V); P_1 = P - 1; IF (P = 1) THEN DO; NOTE = "ERROR 70: VALUES must have at least 2 elements."; PRINT NOTE[LABEL=" " ROWNAME=" " COLNAME=" "]; ABORT; END; ***All polynomials, 0 through p1-1; MEAN = SUM(V) / P; CENTER = V - MEAN; CENTSCL = CENTER / SQRT( CENTER`*CENTER ); POLY = ORPOL(CENTSCL); ***All trends; ZEROTREND = POLY[,1]; U = POLY[,2:P]; ULBL = NAMELIST(NAME,1,P_1,1); FINISH UPOLY1; *************************************************************************; * *; * UPOLY2 *; * *; * This module creates a U contrast matrix with orthogonal polynomial *; * coding for two within subject factors. It assumes Factor 1, with *; * levels VALUES1, varies slowly and that Factor 2, with levels VALUES2, *; * varies rapidly. *; * *; * INPUTS *; * *; * VALUES1, set of numeric treatment levels/ values for Factor 1 *; * NAME1, character string providing stem of labels for Factor 1 *; * VALUES2, set of numeric treatment levels/values for Factor 2 *; * NAME2, character string providing stem of labels for Factor 2 *; * *; * OUTPUTS *; * *; * U1, Matrix of orthonormal polynomial coefficient columns *; * corresponding to the main effect of Factor 1 *; * U1LBL, Matrix of labels for columns of U1 *; * U2, Matrix of orthonormal polynomial coefficient columns *; * corresponding to the main effect of Factor 2 *; * U2LBL, Matrix of labels for columns of U2 *; * U12, Matrix of orthonormal polynomial coefficient columns *; * corresponding the Factor 1 by Factor 2 interactions *; * U12LBL, Matrix of labels for columns of U12 *; * *; *************************************************************************; START UPOLY2(VALUES1, NAME1, VALUES2, NAME2, U1, U1LBL, U2, U2LBL , U12, U12LBL); IF (NROW(VALUES1) = 0) | (NROW(VALUES2) = 0) | (NROW(NAME1) = 0) | (NROW(NAME2) = 0) THEN DO; NOTE = CONCAT("ERROR 71: One or more inputs to UPOLY2 has ", "not been assigned a value."); PRINT NOTE[LABEL=" " ROWNAME=" " COLNAME=" "]; ABORT; END; IF ((NCOL(VALUES1) > 1) & (NROW(VALUES1) > 1)) | ((NCOL(VALUES2) > 1) & (NROW(VALUES2) > 1)) THEN DO; NOTE = CONCAT("ERROR 72: VALUES1 or VALUES2 has more than one ", "column and more than one row. Both must be a row ", "or column vector"); PRINT NOTE[LABEL=" " ROWNAME=" " COLNAME=" "]; ABORT; END; IF (NCOL(NAME1) > 1) | (NROW(NAME1) > 1) | (NCOL(NAME2) > 1) | (NROW(NAME2) > 1) THEN DO; NOTE = CONCAT("ERROR 73: One or both of NAME1 and NAME2 has ", "more than one column or more than one row. Both ", "must be a scalar (1x1)."); PRINT NOTE[LABEL=" " ROWNAME=" " COLNAME=" "]; ABORT; END; IF NCOL(VALUES1) > 1 THEN V1 = VALUES1`; ELSE V1 = VALUES1 ; IF NCOL(VALUES2) > 1 THEN V2 = VALUES2`; ELSE V2 = VALUES2; IF (TYPE(VALUES1)^="N") | (TYPE(VALUES2)^="N")| (TYPE(NAME1)^="C") | (TYPE(NAME2)^="C") THEN DO; NOTE = CONCAT("ERROR 74: NAMES1 or NAMES2 not character or ", "VALUES1 or VALUES2 not numeric."); PRINT NOTE[LABEL=" " ROWNAME=" " COLNAME=" "]; ABORT; END; P1 = NROW(V1); P1_1 = P1 - 1; P2 = NROW(V2); P2_1 = P2 - 1; IF (P1 = 1) | (P2 = 1) THEN DO; NOTE = CONCAT("ERROR 75: VALUES1 and VALUES2 must both have ", "at least 2 elements."); PRINT NOTE[LABEL=" " ROWNAME=" " COLNAME=" "]; ABORT; END; ***All Factor 1 polynomials, 0 through p1-1; MEAN1 = SUM(V1) / P1; CENTER1 = V1 - MEAN1; CENTSCL1 = CENTER1 / SQRT( CENTER1`*CENTER1 ); POLY1 = ORPOL(CENTSCL1); ***All Factor 1 trends; ZEROTREND1 = POLY1[,1]; TREND1 = POLY1[,2:P1]; ***All Factor 2 polynomials, 0 through p2-1; MEAN2 = SUM(V2) / P2; CENTER2 = V2 - MEAN2; CENTSCL2 = CENTER2 / SQRT( CENTER2`*CENTER2 ); POLY2 = ORPOL(CENTSCL2); ***All Factor 2 trends; ZEROTREND2 = POLY2[,1]; TREND2 = POLY2[,2:P2]; ***Grand mean/zero order trend; U0 = ZEROTREND1 @ ZEROTREND2; NMOUT0 = "GRAND MEAN"; ***Main effect Factor 1; U1 = TREND1 @ ZEROTREND2; U1LBL = NAMELIST(NAME1,1,P1_1,1); ***Main effect Factor 2; U2 = ZEROTREND1 @ TREND2; U2LBL = NAMELIST(NAME2,1,P2_1,1); ***Interactions; U12 = TREND1 @ TREND2; DO I = 1 TO P1_1; DO J = 1 TO P2_1; U12LBL = U12LBL || COMPRESS(CONCAT(U1LBL[I], U2LBL[J])); END; END; FINISH UPOLY2; *************************************************************************; * *; * UPOLY3 *; * *; * This module creates a U contrast matrix with orthogonal polynomial *; * coding for three within subject factors. It assumes Factor 1, with *; * levels VALUES1, varies most slowly and that Factor 3, with levels *; * VALUES3, varies most rapidly. *; * *; * INPUTS *; * *; * VALUES1, set of numeric treatment levels/ values for Factor 1 *; * NAME1, character string providing stem of labels for Factor 1 *; * VALUES2, set of numeric treatment levels/values for Factor 2 *; * NAME2, character string providing stem of labels for Factor 2 *; * VALUES3, set of numeric treatment levels/ values for Factor 3 *; * NAME3, character string providing stem of labels for Factor 3 *; * *; * OUTPUTS *; * *; * U1, Matrix of orthonormal polynomial coefficient columns *; * corresponding to the main effect of Factor 1 *; * U1LBL, Matrix of labels for columns of U1 *; * U2, Matrix of orthonormal polynomial coefficient columns *; * corresponding to the main effect of Factor 2 *; * U2LBL, Matrix of labels for columns of U2 *; * U3, Matrix of orthonormal polynomial coefficient columns *; * corresponding to the main effect of Factor 3 *; * U3LBL, Matrix of labels for columns of U3 *; * U12, Matrix of orthonormal polynomial coefficient columns *; * corresponding the Factor 1 by Factor 2 interactions *; * U12LBL, Matrix of labels for columns of U12 *; * U13, Matrix of orthonormal polynomial coefficient columns *; * corresponding the Factor 1 by Factor 3 interactions *; * U13LBL, Matrix of labels for columns of U13 *; * U23, Matrix of orthonormal polynomial coefficient columns *; * corresponding the Factor 2 by Factor 3 interactions *; * U23LBL, Matrix of labels for columns of U23 *; * U123, Matrix of orthonormal polynomial coefficient columsn *; * corresponding to the Factor 1 by Factor 2 by Factor 3 *; * interactions *; * U123LBL, Matrix of labels for the columns of U123 *; * *; *************************************************************************; START UPOLY3(VALUES1, NAME1, VALUES2, NAME2, VALUES3, NAME3, U1, U1LBL, U2, U2LBL, U3, U3LBL, U12, U12LBL, U13, U13LBL, U23, U23LBL, U123, U123LBL); IF (NROW(VALUES1) = 0)| (NROW(VALUES2) = 0)| (NROW(VALUES3) = 0)| (NROW(NAME1) = 0) | (NROW(NAME2) = 0) | (NROW(NAME3) = 0) THEN DO; NOTE = CONCAT("ERROR 76: One or more inputs to UPOLY3 has not been ", "assigned a value."); PRINT NOTE[LABEL=" " ROWNAME=" " COLNAME=" "]; ABORT; END; IF ((NCOL(VALUES1) > 1) & (NROW(VALUES1) > 1)) | ((NCOL(VALUES2) > 1) & (NROW(VALUES2) > 1)) | ((NCOL(VALUES3) > 1) & (NROW(VALUES3) > 1)) THEN DO; NOTE = CONCAT("ERROR 77: One or more of VALUES1, VALUES2, and ", "VALUES3 has more than one column and more than one ", "row. Each must be a row or column vector"); PRINT NOTE[LABEL=" " ROWNAME=" " COLNAME=" "]; ABORT; END; IF (NCOL(NAME1) > 1) | (NROW(NAME1) > 1) | (NCOL(NAME2) > 1) | (NROW(NAME2) > 1) | (NCOL(NAME3) > 1) | (NROW(NAME3) > 1) THEN DO; NOTE = CONCAT("ERROR 78: One or more of NAME1, NAME2, and NAME3 ", "has more than one column or more than one row. Each ", "must be a scalar (1x1)."); PRINT NOTE[LABEL=" " ROWNAME=" " COLNAME=" "]; ABORT; END; IF NCOL(VALUES1) > 1 THEN V1 = VALUES1`; ELSE V1 = VALUES1 ; IF NCOL(VALUES2) > 1 THEN V2 = VALUES2`; ELSE V2 = VALUES2; IF NCOL(VALUES3) > 1 THEN V3 = VALUES3`; ELSE V3 = VALUES3; IF (TYPE(VALUES1)^="N") | (TYPE(VALUES2)^="N")| (TYPE(VALUES3)^="N") | (TYPE(NAME1)^="C") | (TYPE(NAME2)^="C") | (TYPE(NAME3)^="C") THEN DO; NOTE = CONCAT("ERROR 79: VALUES1, VALUES2, or VALUES3 not ", "numeric or NAME1, NAME2, or NAME3 not character."); PRINT NOTE[LABEL=" " ROWNAME=" " COLNAME=" "]; ABORT; END; P1 = NROW(V1); P1_1 = P1 - 1; P2 = NROW(V2); P2_1 = P2 - 1; P3 = NROW(V3); P3_1 = P3 - 1; IF (P1 = 1) | (P2 = 1) | (P3 = 1) THEN DO; NOTE = CONCAT("ERROR 80: VALUES1, VALUES2, and VALUES3 must ", "each have at least 2 elements."); PRINT NOTE[LABEL=" " ROWNAME=" " COLNAME=" "]; ABORT; END; ***All Factor 1 polynomials, 0 through p1-1; MEAN1 = SUM(V1) / P1; CENTER1 = V1 - MEAN1; CENTSCL1 = CENTER1 / SQRT( CENTER1`*CENTER1 ); POLY1 = ORPOL(CENTSCL1); ***All Factor 1 trends; ZEROTREND1 = POLY1[,1]; TREND1 = POLY1[,2:P1]; ***All Factor 2 polynomials, 0 through p2-1; MEAN2 = SUM(V2) / P2; CENTER2 = V2 - MEAN2; CENTSCL2 = CENTER2 / SQRT( CENTER2`*CENTER2 ); POLY2 = ORPOL(CENTSCL2); ***All Factor 2 trends; ZEROTREND2 = POLY2[,1]; TREND2 = POLY2[,2:P2]; ***All Factor 3 polynomials, 0 through p3-1; MEAN3 = SUM(V3) / P3; CENTER3 = V3 - MEAN3; CENTSCL3 = CENTER3 / SQRT( CENTER3`*CENTER3 ); POLY3 = ORPOL(CENTSCL3); ***All Factor 3 trends; ZEROTREND3 = POLY3[,1]; TREND3 = POLY3[,2:P3]; ***Grand mean/zero order trend; U0 = ZEROTREND1 @ ZEROTREND2 @ ZEROTREND3; NMOUT0 = "GRAND MEAN"; ***Main effect Factor 1; U1 = TREND1 @ ZEROTREND2 @ ZEROTREND3; U1LBL = NAMELIST(NAME1,1,P1_1,1); ***Main effect Factor 2; U2 = ZEROTREND1 @ TREND2 @ ZEROTREND3; U2LBL = NAMELIST(NAME2,1,P2_1,1); ***Main effect Factor 3; U3 = ZEROTREND1 @ ZEROTREND2 @ TREND3; U3LBL = NAMELIST(NAME3,1,P3_1,1); ***2 and 3 Interactions; U12 = TREND1 @ TREND2 @ ZEROTREND3; U13 = TREND1 @ ZEROTREND2 @ TREND3; U23 = ZEROTREND1 @ TREND2 @ TREND3; U123 = TREND1 @ TREND2 @ TREND3; DO I = 1 TO P1_1; DO J = 1 TO P2_1; U12LBL = U12LBL || COMPRESS(CONCAT(U1LBL[I], U2LBL[J])); DO K = 1 TO P3_1; IF J = 1 THEN U13LBL = U13LBL || COMPRESS(CONCAT(U1LBL[I], U3LBL[K])); IF I = 1 THEN U23LBL = U23LBL || COMPRESS(CONCAT(U2LBL[J], U3LBL[K])); U123LBL = U123LBL || COMPRESS(CONCAT(U12LBL[I*J], U3LBL[K])); END; END; END; FINISH UPOLY3;