OPTIONS NOCENTER FORMCHAR = "|----|+|---+=";; TITLE1 'Bock (1983) data - Fit of Models'; DATA ONE; INPUT ID Depress Intcpt Linear LinChang Order OrdLin OrdLCha; DATALINES; 1 1 1 -2.5 -0.5 0 0 0 1 1 1 -1.5 0.0 0 0 0 1 1 1 -0.5 0.5 0 0 0 1 1 1 0.5 0.5 0 0 0 1 1 1 1.5 0.0 0 0 0 1 1 1 2.5 -0.5 0 0 0 .......................................... more data .......................................... 75 5 1 -2.5 -0.5 1 -2.5 -0.5 75 5 1 -1.5 0 1 -1.5 0 75 5 1 -0.5 0.5 1 -0.5 0.5 75 5 1 0.5 0.5 1 0.5 0.5 75 5 1 1.5 0 1 1.5 0 75 4 1 2.5 -0.5 1 2.5 -0.5 ; /* get means across time for the two groups */ proc means; class order linear; var depress; run; /* get correlations & sds of variables across time */ DATA t1;SET one; IF linear=-2.5; dep1 = depress; DATA t2;SET one; IF linear=-1.5; dep2 = depress; DATA t3;SET one; IF linear=-0.5; dep3 = depress; DATA t4;SET one; IF linear= 0.5; dep4 = depress; DATA t5;SET one; IF linear= 1.5; dep5 = depress; DATA t6;SET one; IF linear= 2.5; dep6 = depress; DATA comp (KEEP=id dep1-dep6); MERGE t1 t2 t3 t4 t5 t6; BY id; PROC CORR ; VAR DEP1-DEP6; RUN; ---------- output ---------- Bock (1983) data - Fit of Models The MEANS Procedure Analysis Variable : Depress N Order Linear Obs N Mean Std Dev Minimum Maximum -------------------------------------------------------------------------------------------- 0 -2.5 46 46 3.7608696 1.3529356 1.0000000 6.0000000 -1.5 46 46 3.4565217 1.3777038 1.0000000 6.0000000 -0.5 46 46 3.1086957 1.4640622 1.0000000 6.0000000 0.5 46 46 2.8913043 1.5524641 1.0000000 6.0000000 1.5 46 46 2.8043478 1.6140985 1.0000000 6.0000000 2.5 46 46 2.7391304 1.6657969 1.0000000 6.0000000 1 -2.5 29 29 4.7241379 0.9597824 3.0000000 6.0000000 -1.5 29 29 4.6206897 1.1152768 2.0000000 6.0000000 -0.5 29 29 4.5517241 1.1827996 2.0000000 6.0000000 0.5 29 29 4.4482759 1.2126184 2.0000000 6.0000000 1.5 29 29 4.2068966 1.3464055 1.0000000 6.0000000 2.5 29 29 3.8965517 1.3717775 1.0000000 6.0000000 -------------------------------------------------------------------------------------------- Bock (1983) data - Fit of Models The CORR Procedure 6 Variables: dep1 dep2 dep3 dep4 dep5 dep6 Simple Statistics Variable N Mean Std Dev Sum Minimum Maximum dep1 75 4.13333 1.29795 310.00000 1.00000 6.00000 dep2 75 3.90667 1.39665 293.00000 1.00000 6.00000 dep3 75 3.66667 1.52753 275.00000 1.00000 6.00000 dep4 75 3.49333 1.61390 262.00000 1.00000 6.00000 dep5 75 3.34667 1.65622 251.00000 1.00000 6.00000 dep6 75 3.18667 1.64968 239.00000 1.00000 6.00000 Pearson Correlation Coefficients, N = 75 Prob > |r| under H0: Rho=0 dep1 dep2 dep3 dep4 dep5 dep6 dep1 1.00000 0.90896 0.74520 0.68424 0.59426 0.60040 <.0001 <.0001 <.0001 <.0001 <.0001 dep2 0.90896 1.00000 0.86568 0.81806 0.70353 0.67629 <.0001 <.0001 <.0001 <.0001 <.0001 dep3 0.74520 0.86568 1.00000 0.91176 0.77808 0.72217 <.0001 <.0001 <.0001 <.0001 <.0001 dep4 0.68424 0.81806 0.91176 1.00000 0.87550 0.84303 <.0001 <.0001 <.0001 <.0001 <.0001 dep5 0.59426 0.70353 0.77808 0.87550 1.00000 0.95530 <.0001 <.0001 <.0001 <.0001 <.0001 dep6 0.60040 0.67629 0.72217 0.84303 0.95530 1.00000 <.0001 <.0001 <.0001 <.0001 <.0001 -------------- IML code -------------- PROC IML; /* 3 random effects and CI errors */ xt0 = { 1 1 1 1 1 1 , -2.5 -1.5 -0.5 0.5 1.5 2.5 , -0.5 0.0 0.5 0.5 0.0 -0.5 , 0 0 0 0 0 0 , 0 0 0 0 0 0 , 0 0 0 0 0 0 }; x0 = T(xt0); xt1 = { 1 1 1 1 1 1 , -2.5 -1.5 -0.5 0.5 1.5 2.5 , -0.5 0.0 0.5 0.5 0.0 -0.5 , 1 1 1 1 1 1 , -2.5 -1.5 -0.5 0.5 1.5 2.5 , -0.5 0.0 0.5 0.5 0.0 -0.5}; x1 = T(xt1); zt = { 1 1 1 1 1 1 , -2.5 -1.5 -0.5 0.5 1.5 2.5 , -0.5 0.0 0.5 0.5 0.0 -0.5}; z = T(zt); beta = { 3.127 , -0.208 , -0.250 , 1.281 , 0.051 , 0.440 }; var_ups = { 1.463 0.099 0.143 , 0.099 0.079 -0.001 , 0.143 -0.001 0.413}; var_err = { 0.149 }; means0 = T(x0*beta); means1 = T(x1*beta); varcov = z*var_ups*t(z) + var_err*I(6); sd = sqrt(vecdiag(varcov)); sdt = T(sd); s = 1/sd; corr = diag(s)*varcov*diag(s); PRINT '3 random effects model - CI errors'; PRINT 'Estimated Means across time - group=0', means0 [FORMAT=10.4]; PRINT 'Estimated Means across time - group=1', means1 [FORMAT=10.4]; PRINT 'Estimated SDs across time', sdt [FORMAT=10.4]; PRINT 'Estimated Correlations', corr [FORMAT=10.4]; /* 1 random effect and Toeplitz errors with 3 lags */ xt0 = { 1 1 1 1 1 1 , -2.5 -1.5 -0.5 0.5 1.5 2.5 , -0.5 0.0 0.5 0.5 0.0 -0.5 , 0 0 0 0 0 0 , 0 0 0 0 0 0 , 0 0 0 0 0 0 }; x0 = T(xt0); xt1 = { 1 1 1 1 1 1 , -2.5 -1.5 -0.5 0.5 1.5 2.5 , -0.5 0.0 0.5 0.5 0.0 -0.5 , 1 1 1 1 1 1 , -2.5 -1.5 -0.5 0.5 1.5 2.5 , -0.5 0.0 0.5 0.5 0.0 -0.5}; x1 = T(xt1); zt = { 1 1 1 1 1 1 }; z = T(zt); beta = { 3.126 , -0.200 , -0.252 , 1.281 , 0.024 , 0.438 }; var_ups = { 1.062 }; var_err = { 0.860 }; rho = { 1.000 , 0.738 , 0.489 , 0.214 }; omega = j(6,6,0); do i = 1 to 6; do j = 1 to 6; diagind = 1 + ABS(i-j); if diagind>4 then omega[i,j] = 0; if diagind<=4 then omega[i,j] = rho[diagind,1]; end; end; means0 = T(x0*beta); means1 = T(x1*beta); varcov = z*var_ups*t(z) + var_err*omega; sd = sqrt(vecdiag(varcov)); sdt = T(sd); s = 1/sd; corr = diag(s)*varcov*diag(s); PRINT '1 random effect model - Toeplitz(4) errors'; PRINT 'Estimated Means across time - group=0', means0 [FORMAT=10.4]; PRINT 'Estimated Means across time - group=1', means1 [FORMAT=10.4]; PRINT 'Estimated SDs across time', sdt [FORMAT=10.4]; PRINT 'Estimated Correlations', corr [FORMAT=10.4]; PRINT 'Autocorrelated Error Structure', omega [FORMAT=10.4]; ---------------- IML output ---------------- Bock (1983) data - Fit of Models 3 random effects model - CI errors Estimated Means across time - group=0 MEANS0 3.7720 3.4390 3.1060 2.8980 2.8150 2.7320 Estimated Means across time - group=1 MEANS1 4.7055 4.6435 4.5815 4.4245 4.1725 3.9205 Estimated SDs across time SDT 1.2524 1.2218 1.3340 1.4059 1.4446 1.6011 Estimated Correlations CORR 1.0000 0.8437 0.6958 0.6046 0.5511 0.4635 0.8437 1.0000 0.8568 0.8017 0.7282 0.6101 0.6958 0.8568 1.0000 0.9009 0.8165 0.6824 0.6046 0.8017 0.9009 1.0000 0.8819 0.7795 0.5511 0.7282 0.8165 0.8819 1.0000 0.9013 0.4635 0.6101 0.6824 0.7795 0.9013 1.0000 1 random effect model - Toeplitz(4) errors Estimated Means across time - group=0 MEANS0 3.7520 3.4260 3.1000 2.9000 2.8260 2.7520 Estimated Means across time - group=1 MEANS1 4.7540 4.6710 4.5880 4.4120 4.1430 3.8740 Estimated SDs across time SDT 1.3864 1.3864 1.3864 1.3864 1.3864 1.3864 Estimated Correlations CORR 1.0000 0.8828 0.7714 0.6483 0.5525 0.5525 0.8828 1.0000 0.8828 0.7714 0.6483 0.5525 0.7714 0.8828 1.0000 0.8828 0.7714 0.6483 0.6483 0.7714 0.8828 1.0000 0.8828 0.7714 0.5525 0.6483 0.7714 0.8828 1.0000 0.8828 0.5525 0.5525 0.6483 0.7714 0.8828 1.0000 Autocorrelated Error Structure OMEGA 1.0000 0.7380 0.4890 0.2140 0.0000 0.0000 0.7380 1.0000 0.7380 0.4890 0.2140 0.0000 0.4890 0.7380 1.0000 0.7380 0.4890 0.2140 0.2140 0.4890 0.7380 1.0000 0.7380 0.4890 0.0000 0.2140 0.4890 0.7380 1.0000 0.7380 0.0000 0.0000 0.2140 0.4890 0.7380 1.0000