* THIS PROGRAM PERFORMS A RANDOM-INTERCEPTS LINEAR REGRESSION ANALYSIS ; * AS DISCUSSED IN HEDEKER, GIBBONS, & FLAY "RANDOM-EFFECTS REGRESSION ; * MODELS FOR CLUSTERED DATA: WITH AN EXAMPLE FROM SMOKING PREVENTION RESEARCH"; * JOURNAL OF CONSULTING AND CLINICAL PSYCHOLOGY (1994); * THIS PROGRAM IS WRITTEN UTILZING THE IML MATRIX ROUTINES OF SAS; * USING MAXIMUM MARGINAL LIKELIHOOD AND EMPIRICAL BAYES ESTIMATION; * ANY COMMENTS OR QUESTIONS CAN BE DIRECTED TO DON HEDEKER U41098@UICVM ; * PLEASE FEEL FREE TO REDISTRIBUTE THIS MACRO ; * NOTE THAT IN THE PROGRAM BELOW - THE DATA AND SAS COMMANDS ARE ; * INCLUDED AS 1 FILE - THIS IS ONLY FOR CONVENIENCE IN DISTRIBUTION ; * THE DATA CAN BE EXTRACTED INTO A SEPARATE FILE AND THEN READ IN ; * BY THE PROGRAM FILE - SEE THE SAS MANUAL FOR INSTRUCTIONS ON THIS ; * ; * this macro was first written in SPSS by Don Hedeker (10/91) ; * and then ported over to SAS by Chitra Shaligram (2/94); * ; OPTIONS LS=80 NODATE; TITLE 'ANALYSIS OF RIESBY DATA'; DATA ONE; INPUT ID HamD Intcpt Week Endog EndWeek ; * ; * ; LABEL HamD = 'Hamilton Depression Score' Endog = 'Endogenous Diagnosis'; * ; CARDS; 101 26 1 0 0 0 101 22 1 1 0 0 101 18 1 2 0 0 101 7 1 3 0 0 101 4 1 4 0 0 101 3 1 5 0 0 ..... MORE DATA ...... 361 30 1 0 1 0 361 22 1 1 1 1 361 11 1 2 1 2 361 8 1 3 1 3 361 7 1 4 1 4 361 19 1 5 1 5 ; * ; PROC FORMAT; VALUE ENDOG 0='NonEndog' 1='Endog' ; * ; * ; * RUN PROC MIXED FIRST TO VERIFY THAT THE RESULTS; * THAT ARE OBTAINED BY THE IML MACRO ARE CORRECT ; * ; PROC MIXED METHOD=ML; CLASS ID; MODEL HAMD = WEEK ENDOG ENDWEEK / SOLUTION ; RANDOM INTERCEPT / SUB=ID TYPE=UN G SOLUTION; * ; * ; * for the Random-effects program the data ; * needs to be sorted by the cluster ID ; * and observations with missing data on either ; * the dependent variable, the cluster ID ; * or any of the independent variables ; * need to be excluded ; * ; * ; DATA TWO; SET ONE; IF ID =. OR HamD =. OR Endog =. OR Week =. OR EndWeek =. THEN DELETE; PROC SORT DATA=TWO; BY ID; * ; * ; * ; * 2-LEVEL RANDOM REGRESSION MODEL (with one random effect) ; * ; * N = number of clusters ; * NTOT = number of total observations ; * NI = vector of size N with number of observations per cluster ; * ; * NOTE THAT THE ONLY LINES BELOW THAT NEED TO BE CHANGED ARE ; * USE TWO VAR{ } ; * ID = ; * Y = ; * X = ; * YNAME = ; * XNAMES = ; * ; * ; PROC IML; USE TWO VAR{ID HamD Intcpt Week Endog EndWeek}; READ ALL INTO REGDAT; ID = REGDAT[,1]; Y = REGDAT[,2]; X = REGDAT[,3:6]; YNAME = 'HamD'; XNAMES = {'INTERCEPT' 'Week' 'Endog' 'EndWeek'}; NTOT = NROW(ID); PRINT 'RANDOM-EFFECTS REGRESSION ANALYSIS '; PRINT ' '; PRINT 'TOTAL NUMBER OF OBSERVATIONS '; PRINT NTOT; N = 1; IDOLD = 0; COUNT = 0; NITEMP = J(NTOT,1,0); DO I = 1 TO NTOT; IDTEMP = ID[I]; IF IDOLD = 0 THEN IDOLD = IDTEMP; IF IDTEMP = IDOLD THEN COUNT = COUNT + 1; ELSE DO; NITEMP[N] = COUNT; N = N + 1; COUNT = 1; IDOLD = IDTEMP; END; IF I = NTOT THEN NITEMP[N] = COUNT; END; PRINT 'NUMBER OF CLUSTERS '; PRINT N; NI = J(N,1,0); DO K = 1 TO N; NI[K] = NITEMP[K]; END; TRANSNI = T(NI); PRINT 'CLUSTER SAMPLE SIZES '; PRINT TRANSNI; /* COMPUTE STARTING VALUES */ BETA = INV(X` * X) * (X`* Y); YHAT = X * BETA; RESID = Y - YHAT; SSE = RESID` * RESID; ERVAR = (1.0/NTOT) # SSE; CLUSTVAR = 0.10 # ERVAR; /* SET CONVERGENCE CRITERIA AND START ITERATIONS */ MAXCORR=1; CONVERGE = 0.0001; IT = 0; LN2PI = LOG( 2 * 3.141592654); P = NCOL(X); VARDER = J( 2,1,0); VARINF = J( 2,2,0); ALPHA = J(N,1,0); ALPHAVAR = J(N,1,0); DO WHILE ( MAXCORR > CONVERGE); IT = IT + 1; /* PRINT ' ITERATION = '; PRINT IT; */ /* COMPUTE THE INTRACLASS CORRELATION */ INTCLASS = CLUSTVAR / ( CLUSTVAR + ERVAR ); /* COMPUTE THE MODEL DEVIATES (WITHOUT THE CLUSTER EFFECTS) */ U = Y - X * BETA; /* COMPUTE EAP ESTIMATES */ IS = 0; IE = 0; DO I = 1 TO N; NII = NI [I] ; IS = IE + 1; IE = IS - 1 + NII; SB = (INTCLASS # NII) / ( 1 + (NII - 1) # INTCLASS ) ; V = SUM (U [IS : IE,]) ; ALPHA [I] = (SB / NII) # V ; ALPHAVAR [I] = (1.0 - SB) # CLUSTVAR ; END; /* COMPUTE THE MODEL DEVIATES ( WITH THE CLUSTER EFFECTS ) */ I=1; I2=0; DO J = 1 TO NTOT; I2 = I2 + 1; U[J] = U[J] - ALPHA[I]; IF (J ^= NTOT & I2 = NI[I]) THEN DO; I2 = 0; I = I + 1; END; END; /* COMPUTE THE FIRST DERIVATIVES */ XX = T(X); BETADER = (1.0 / ERVAR) # (XX * U) ; SSA= ALPHA` * ALPHA ; SUMALPHA = SUM(ALPHAVAR) ; VARDER[1] = .5 # (1.0 / CLUSTVAR) # (SSA + SUMALPHA - N # CLUSTVAR); SSU = U` * U; VARDER[2] = .5 # (1.0/ ERVAR) # (SSU + (TRANSNI * ALPHAVAR) - (NTOT # ERVAR)); /* COMPUTE THE INFORMATION MATRIX */ IS=0; IE=0; BETAINF = J(P,P,0); DO K = 1 TO N; NII = NI [K] ; IS = IE + 1 ; IE = IS - 1 + NII; XTEMP = X [ IS : IE ,] ; ONEVEC = J(NII,1,1); SSXT = XTEMP` * XTEMP; TONE = T(ONEVEC) * XTEMP; SSO = TONE` * TONE; BETAINF = BETAINF + (ERVAR # SSXT) - (ALPHAVAR[K] # SSO); END; BETAINF = (1.0 / (ERVAR * ERVAR)) # BETAINF; DIFVAR = CLUSTVAR - ALPHAVAR; SSVAR = DIFVAR` * DIFVAR; VARINF[1,1] = 0.5 # (1.0/ ( CLUSTVAR ** 2)) # SSVAR; VARINF[2,2] = 0.5 # (1.0/ ( ERVAR ** 2)) # ( TRANSNI * (ERVAR - ALPHAVAR) ## 2); VARINF[2,1] = 0.5 # (1.0 / ERVAR) # (1.0 / CLUSTVAR ) # ( TRANSNI * (ALPHAVAR ## 2 )); VARINF[1,2] = VARINF[2,1]; /* COMPUTE THE CORRECTIONS */ BETACORR = INV(BETAINF) * BETADER; VARCORR = INV(VARINF) * VARDER; /* PRINT ' BETA CORRECTIONS '; PRINT BETACORR; PRINT ' VARIANCE CORRECTIONS '; PRINT VARCORR; */ BETA = BETA + BETACORR; PI = LOG(CLUSTVAR) + VARCORR[1] ; CLUSTVAR = EXP(PI); IF CLUSTVAR < 0.0000001 THEN CLUSTVAR = 0.0000001; TAU = LOG(ERVAR) + VARCORR[2]; ERVAR = EXP (TAU); /* PRINT ' REGRESSION COEFFICIENTS '; PRINT BETA; PRINT 'ERROR VARIANCE '; PRINT ERVAR; PRINT ' CLUSTER VARIANCE'; PRINT CLUSTVAR; */ LOGLIK = - 0.5 # ( SUM(NI) # (LN2PI + LOG(ERVAR)) + N # LOG(CLUSTVAR) - SUM(LOG(ALPHAVAR)) + ((ALPHA` * ALPHA) / CLUSTVAR) + ( U` * U) / ERVAR ) ; /* PRINT ' LOG-LIKELIHOOD '; PRINT LOGLIK; */ MAXCORRB = MAX(ABS(BETACORR)); MAXCORRV = ABS(VARCORR[2]); IF ( CLUSTVAR > 0.0000001 & ABS(VARCORR[1]) > ABS(VARCORR[2])) THEN MAXCORRV = ABS(VARCORR[1]) ; MAXCORR = MAXCORRV; IF(MAXCORRB > MAXCORRV) THEN MAXCORR = MAXCORRB; END; /* COMPUTE STANDARD ERRORS, Z STATISTICS, AND P-VALUES FOR ALL PARAMETERS AND PUT THESE VECTORS AS COLUMN VECTORS IN BETARES & VARRES MATRICES */ BETARES = J(P,4,0); BETARES[1:P,1] = BETA; BETARES[1:P,2] = SQRT(VECDIAG(INV(BETAINF))); BETARES[1:P,3] = BETARES[1:P,1] / BETARES[1:P,2]; BETARES[1:P,4] = 2 * PROBNORM( - ABS(BETARES[1:P,3])); /* THE INFORMATION MATRIX NEEDS TO BE PUT IN TERMS OF THE RAW PARAMETERS */ VARRES = J(2,4,0) ; VARRES[1,1] = CLUSTVAR; VARRES[2,1] = ERVAR; VARINF[1,1] = VARINF[1,1] / (CLUSTVAR * CLUSTVAR); VARINF[2,2] = VARINF[2,2] / (ERVAR * ERVAR) ; VARINF[1,2] = VARINF[1,2] / (ERVAR * CLUSTVAR); VARINF[2,1] = VARINF[2,1] / (ERVAR * CLUSTVAR) ; VARRES[1:2,2] = SQRT(VECDIAG(INV(VARINF))); VARRES[1:2,3] = VARRES[1:2,1] / VARRES[1:2,2] ; VARRES[1:2,4] = PROBNORM ( - ABS(VARRES[1:2,3])); INTCLASS = CLUSTVAR / (CLUSTVAR + ERVAR) ; PRINT ' ' /; PRINT '*******************'; PRINT ' FINAL RESULTS '; PRINT '*******************'; PRINT ' '; PRINT 'DEPENDENT VARIABLE:' YNAME; PRINT ' '; PRINT 'LOG-LIKELIHOOD :' LOGLIK; PRINT ' '; PRINT 'REGRESSION PARAMETERS'; MATTRIB BETARES ROWNAME = (XNAMES) COLNAME = ({ESTIMATE SE Z 'P-VAL(2)'}) LABEL = ' VARIABLE ' FORMAT = 12.6; PRINT BETARES; PRINT ' '; PRINT ' '; PRINT 'VARIANCE PARAMETERS'; MATTRIB VARRES ROWNAME = ({CLUSTER ERROR}) COLNAME = ({ESTIMATE SE Z 'P-VAL(1)'}) LABEL = ' ' FORMAT = 12.6; PRINT VARRES; PRINT 'INTRACLASS CORRELATION'; PRINT INTCLASS; PRINT ' ' /; PRINT 'CLUSTER EMPIRICAL BAYES (EAP) EFFECTS'; PRINT ALPHA; ALPHASD = SQRT(ALPHAVAR); PRINT 'POSTERIOR STANDARD DEVIATIONS'; PRINT ALPHASD; ANALYSIS OF RIESBY DATA The MIXED Procedure Class Level Information Class Levels Values ID 66 101 103 104 105 106 107 108 113 114 115 117 118 120 121 123 302 303 304 305 308 309 310 311 312 313 315 316 318 319 322 327 328 331 333 334 335 337 338 339 344 345 346 347 348 349 350 351 352 353 354 355 357 360 361 501 502 504 505 507 603 604 606 607 608 609 610 ML Estimation Iteration History Iteration Evaluations Objective Criterion 0 1 1701.2349747 1 2 1592.9333232 0.00000000 Convergence criteria met. G Matrix Effect ID Row COL1 INTERCEPT 101 1 15.28474645 Covariance Parameter Estimates (MLE) Cov Parm Subject Estimate UN(1,1) ID 15.28474645 Residual 19.03509762 Model Fitting Information for HAMD Description Value Observations 375.0000 Log Likelihood -1141.07 Akaike's Information Criterion -1143.07 Schwarz's Bayesian Criterion -1147.00 -2 Log Likelihood 2282.137 Null Model LRT Chi-Square 108.3017 Null Model LRT DF 1.0000 Null Model LRT P-Value 0.0000 Solution for Fixed Effects Effect Estimate Std Error DF t Pr > |t| INTERCEPT 22.44162996 0.93906232 64 23.90 0.0001 WEEK -2.35184269 0.19801780 307 -11.88 0.0001 ENDOG 1.99287252 1.26341082 307 1.58 0.1157 ENDWEEK -0.04417597 0.27058959 307 -0.16 0.8704 Solution for Random Effects Effect ID Estimate SE Pred DF t Pr > |t| INTERCEPT 101 -2.67372835 1.75091167 307 -1.53 0.1278 INTERCEPT 103 3.39911909 1.75091167 307 1.94 0.0531 INTERCEPT 104 -1.33419562 1.72474670 307 -0.77 0.4398 INTERCEPT 105 -1.56957427 1.75091167 307 -0.90 0.3707 INTERCEPT 106 1.40707121 1.83780108 307 0.77 0.4445 INTERCEPT 107 -2.05213081 1.83692129 307 -1.12 0.2648 INTERCEPT 108 -4.37061934 1.72474670 307 -2.53 0.0118 INTERCEPT 113 3.25205574 1.86108640 307 1.75 0.0816 INTERCEPT 114 -3.51148529 1.86126730 307 -1.89 0.0602 INTERCEPT 115 -1.79849080 1.83814735 307 -0.98 0.3286 INTERCEPT 117 -5.19873490 1.72474670 307 -3.01 0.0028 INTERCEPT 118 0.12293394 1.83814735 307 0.07 0.9467 INTERCEPT 120 1.05279167 1.75091167 307 0.60 0.5481 INTERCEPT 121 0.91477241 1.75091167 307 0.52 0.6017 INTERCEPT 123 -4.32995947 1.75091167 307 -2.47 0.0139 INTERCEPT 302 -3.54250378 1.72474670 307 -2.05 0.0408 INTERCEPT 303 -2.12165131 1.75091167 307 -1.21 0.2265 INTERCEPT 304 3.13507850 1.83651589 307 1.71 0.0888 INTERCEPT 305 -4.32995947 1.75091167 307 -2.47 0.0139 INTERCEPT 308 0.22467611 1.75091167 307 0.13 0.8980 INTERCEPT 309 1.46684945 1.75091167 307 0.84 0.4028 INTERCEPT 310 -4.99767793 1.83780108 307 -2.72 0.0069 INTERCEPT 311 -1.09141844 1.83692129 307 -0.59 0.5528 INTERCEPT 312 -2.18213826 1.83712916 307 -1.19 0.2358 INTERCEPT 313 -2.67372835 1.75091167 307 -1.53 0.1278 INTERCEPT 315 -2.43577827 1.83780108 307 -1.33 0.1860 INTERCEPT 316 5.98082515 1.72474670 307 3.47 0.0006 INTERCEPT 318 -1.05815710 1.72474670 307 -0.61 0.5400 INTERCEPT 319 -2.71438822 1.72474670 307 -1.57 0.1166 INTERCEPT 322 6.69098924 1.83780108 307 3.64 0.0003 INTERCEPT 327 -0.74145871 1.75091167 307 -0.42 0.6722 INTERCEPT 328 7.95375467 1.75091167 307 4.54 0.0001 INTERCEPT 331 -0.18938167 1.75091167 307 -0.11 0.9139 INTERCEPT 333 3.12308057 1.75091167 307 1.78 0.0755 INTERCEPT 334 0.05857225 1.85957384 307 0.03 0.9749 INTERCEPT 335 -0.32740093 1.75091167 307 -0.19 0.8518 INTERCEPT 337 2.57100353 1.75091167 307 1.47 0.1430 INTERCEPT 338 3.39911909 1.75091167 307 1.94 0.0531 INTERCEPT 339 -2.02201953 1.83712916 307 -1.10 0.2719 INTERCEPT 344 -1.47825335 1.83814735 307 -0.80 0.4219 INTERCEPT 345 6.98761985 1.75091167 307 3.99 0.0001 INTERCEPT 346 2.11628588 1.72474670 307 1.23 0.2208 INTERCEPT 347 -5.31791538 1.83780108 307 -2.89 0.0041 INTERCEPT 348 -1.70759353 1.75091167 307 -0.98 0.3302 INTERCEPT 349 -3.54250378 1.72474670 307 -2.05 0.0408 INTERCEPT 350 0.22467611 1.75091167 307 0.13 0.8980 INTERCEPT 351 1.70222810 1.72474670 307 0.99 0.3244 INTERCEPT 352 -0.09202228 1.72474670 307 -0.05 0.9575 INTERCEPT 353 1.42618958 1.72474670 307 0.83 0.4089 INTERCEPT 354 5.24992069 1.83780108 307 2.86 0.0046 INTERCEPT 355 -1.74825340 1.72474670 307 -1.01 0.3116 INTERCEPT 357 3.49647848 1.72474670 307 2.03 0.0435 INTERCEPT 360 7.63705627 1.72474670 307 4.43 0.0001 INTERCEPT 361 -1.88627266 1.72474670 307 -1.09 0.2750 INTERCEPT 501 5.42874811 1.72474670 307 3.15 0.0018 INTERCEPT 502 -5.61279268 1.72474670 307 -3.25 0.0013 INTERCEPT 504 -1.15551649 1.75091167 307 -0.66 0.5098 INTERCEPT 505 -6.26222911 1.75091167 307 -3.58 0.0004 INTERCEPT 507 5.01469033 1.72474670 307 2.91 0.0039 INTERCEPT 603 1.74288797 1.75091167 307 1.00 0.3203 INTERCEPT 604 -1.55150611 1.86108640 307 -0.83 0.4051 INTERCEPT 606 -4.23260008 1.72474670 307 -2.45 0.0147 INTERCEPT 607 5.42874811 1.72474670 307 3.15 0.0018 INTERCEPT 608 -3.22580539 1.75091167 307 -1.84 0.0664 INTERCEPT 609 -1.31813462 1.83814735 307 -0.72 0.4739 INTERCEPT 610 5.18975771 1.98081923 307 2.62 0.0092 Tests of Fixed Effects Source NDF DDF Type III F Pr > F WEEK 1 307 141.06 0.0001 ENDOG 1 307 2.49 0.1157 ENDWEEK 1 307 0.03 0.8704 ANALYSIS OF RIESBY DATA RANDOM-EFFECTS REGRESSION ANALYSIS TOTAL NUMBER OF OBSERVATIONS NTOT 375 NUMBER OF CLUSTERS N 66 CLUSTER SAMPLE SIZES TRANSNI 6 6 6 6 5 5 6 : 5 5 5 6 5 6 6 : 6 6 6 5 6 6 6 : 5 5 5 6 5 6 6 : 6 5 6 6 6 6 5 : 6 6 6 5 5 6 6 : 5 6 6 6 6 6 6 : 5 6 6 6 6 6 6 : 6 6 6 6 5 6 6 : 6 5 4 ******************* FINAL RESULTS ******************* YNAME DEPENDENT VARIABLE: HamD LOGLIK LOG-LIKELIHOOD : -1141.069 REGRESSION PARAMETERS VARIABLE ESTIMATE SE Z P-VAL(2) INTERCEPT 22.441627 0.939095 23.897078 0.000000 Week -2.351842 0.198015 -11.877101 0.000000 Endog 1.992873 1.263454 1.577322 0.114722 EndWeek -0.044176 0.270586 -0.163260 0.870313 VARIANCE PARAMETERS ESTIMATE SE Z P-VAL(1) CLUSTER 15.286630 3.260563 4.688341 0.000001 ERROR 19.034726 1.630563 11.673714 0.000000 INTRACLASS CORRELATION INTCLASS 0.445397 CLUSTER EMPIRICAL BAYES (EAP) EFFECTS ALPHA -2.673799 3.3992229 -1.334247 -1.569613 1.4071144 -2.052206 -4.370758 3.2521782 -3.511605 -1.798572 -5.198897 0.1229173 1.052828 0.9148048 -4.330078 -3.542618 -2.121706 3.1351726 -4.330078 0.2246886 1.4668977 -4.997848 -1.091462 -2.182228 -2.673799 -2.435863 5.9809846 -1.0582 -2.714479 6.6912088 -0.741474 7.9539896 -0.189381 3.1231765 0.0585785 -0.327404 2.5710836 3.3992229 -2.022104 -1.478323 6.9878269 2.1163341 -5.318097 -1.707637 -3.542618 0.2246886 1.7022644 -0.092038 1.4262179 5.2500922 -1.748316 3.4965664 7.6372633 -1.88634 5.4288916 -5.612967 -1.155544 -6.262403 5.014822 1.7429442 -1.551544 -4.232735 5.4288916 -3.225892 -1.318199 5.1899523 POSTERIOR STANDARD DEVIATIONS ALPHASD 1.6208657 1.6208657 1.6208657 1.6208657 1.7458192 1.7458192 1.6208657 1.7458192 1.7458192 1.7458192 1.6208657 1.7458192 1.6208657 1.6208657 1.6208657 1.6208657 1.6208657 1.7458192 1.6208657 1.6208657 1.6208657 1.7458192 1.7458192 1.7458192 1.6208657 1.7458192 1.6208657 1.6208657 1.6208657 1.7458192 1.6208657 1.6208657 1.6208657 1.6208657 1.7458192 1.6208657 1.6208657 1.6208657 1.7458192 1.7458192 1.6208657 1.6208657 1.7458192 1.6208657 1.6208657 1.6208657 1.6208657 1.6208657 1.6208657 1.7458192 1.6208657 1.6208657 1.6208657 1.6208657 1.6208657 1.6208657 1.6208657 1.6208657 1.6208657 1.6208657 1.7458192 1.6208657 1.6208657 1.6208657 1.7458192 1.9049856