/* Edps/Psych/Stat 587 C.J. Anderson */ * Program to simmulate data with difference sample sizes * to demonstrate Effect on estimates of random effects ; options ls=80; data sample2; gamma00=12.0; gamma11=2.0; do id=1 to 160; z1=rannor(-1); Uo= 3*z1; do i=1 to 2; z2=rannor(-1); x =rannor(-1); R1=6*z2; y2 = gamma00 + gamma11*x + Uo + R1; output; end; end; run; proc mixed data=sample2 noclprint covtest method=ML ic; title 'N=160, n_j=2'; class id; model y2 = x /solution outpred=pred2 outpredm=pm2; random intercept / subject=id type=un solution cl; ods output SolutionR=RanUs2; run; proc sort data=RanUs2; by Estimate; run; data stdranU2; set RanUs2; stdU = estimate/StdErrPred; run; data tmp2; set pred2; predij=pred; residij=resid; lowij=lower; upij= upper; stdij=StdErrPred; keep predij residij lowij upij stdij id; run; proc sort data=tmp2; by id; proc sort data=pm2; by id; data predict; merge tmp2 pm2; by id; run; goptions reset=(axis, legend, pattern, symbol, title, footnote) norotate hpos=0 vpos=0 htext=2 ftext=Swiss ctext= target= gaccess= gsfmode= ; goptions device=WIN ctext=blue graphrc interpol=join; symbol1 c=DEFAULT ci=blue v=NONE line=3; symbol2 c=DEFAULT ci=BLACK v=NONE ; axis1 color=black width=2.0 label=(height=2 'Group Centered SES'); axis2 color=black width=2.0 label=(angle=90 height=2 'Predicted Mean'); axis3 color=black width=2.0 label=(angle=90 height=2 'Predicted School Mean'); legend1 frame; proc gplot data=WORK.PREDICT; plot Pred *x / haxis=axis1 vaxis=axis2 frame ; title 'Overall Mean Regression Line: Y= gamma00 + gamma10*X'; run; proc gplot data=WORK.PREDICT; plot predij*x=id / haxis=axis1 vaxis=axis3 frame nolegend; title1 height=2 'Group regression lines'; title2 height=1.5 ' Y = gamma00 + gamma10*X + U_j '; run; data newid; do N=1 to 160; newid = N; output; end; run; proc sort data=RanUs2; by estimate; data RanUxx2; merge RanUs2 newid; run; goptions reset=(axis, legend, pattern, symbol, title, footnote) norotate hpos=0 vpos=0 htext=2 ftext=swiss ctext= target= gaccess= gsfmode= ; goptions device=WIN ctext=black graphrc interpol=join ; symbol1 value=none interpol=join color=black; symbol2 value=plus interpol=none color=black; symbol3 value=plus interpol=none color=black ; symbol4 value=none interpol=join color=black line=3 ci=green; symbol5 value=none interpol=join color=black line=1 ci=red; symbol6 value=none interpol=join color=black line=2 ci=blue; axis1 color=black width=2.0 length=60 pct label=(angle=90 'BLUPS & 95% CI') order=-15 to 15 by 5; axis2 color=black width=2.0 length=60 pct label=('Marco Units ordered by BLUPS'); axis3 color=black width=2.0 ; axis4 color=black width=2.0 length=60 pct label=(angle=90 'OLS & Beta0j'); proc gplot data=WORK.RANUxx2 gout=fig1nj; plot (Estimate Lower Upper ) * newid / overlay haxis=axis2 vaxis=axis1 frame; title1 'N = 160, n_j = 2'; title2 '95% confidence intervals'; run; legend1 position=(inside left top) frame cshadow=grey label=none; goptions reset=(axis, legend, pattern, symbol, title, footnote) norotate hpos=0 vpos=0 htext= ftext=swiss ctext=black target= gaccess= gsfmode= ; goptions device=WIN ctext=black graphrc interpol=join; pattern1 value=SOLID color=black;pattern2 value=SOLID color=black; pattern3 value=SOLID color=black;pattern4 value=SOLID color=black; pattern5 value=SOLID color=black;pattern6 value=SOLID color=black; pattern7 value=SOLID color=black;pattern8 value=SOLID color=black; pattern9 value=SOLID color=black;pattern10 value=SOLID color=black; pattern11 value=SOLID color=black;pattern12 value=SOLID color=black; pattern13 value=SOLID color=black;pattern14 value=SOLID color=black; pattern15 value=SOLID color=black;pattern16 value=SOLID color=black; pattern17 value=SOLID color=black;pattern18 value=SOLID color=black; pattern19 value=SOLID color=black;pattern20 value=SOLID color=black; axis1 length=80 pct color=black width=2.0 label=("Empirical Bayes Estimate of Uo") ; axis21 length=80 pct color=black width=2.0 label=("Sigma^2=1") value=none ; axis3 length=80 pct color=black width=2.0 label=("Sigma^2=9") value=none ; axis5 length=80 pct color=black width=2.0 label=("Sigma^2=25") value=none ; axis10 length=80 pct color=black width=2.0 label=("Sigma^2=100") value=none ; axis2 length=70 pct color=black width=2.0 label=(" ") major=none minor=none value=none; *axis3 color=black width=2.0 ; proc gchart data=WORK.RANUxx2 gout=fig2njaa; VBAR Estimate / maxis=axis1 raxis=axis2 frame type=FREQ patternid=midpoint; title 'N = 160, n_j = 2'; run; ****************************************** n_j = 5 *********************************; data sample5; gamma00=12.0; gamma11=2.0; do id=1 to 160; z1=rannor(-1); Uo= 3*z1; do i=1 to 5; z2=rannor(-1); x =rannor(-1); R1=6*z2; y5 = gamma00 + gamma11*x + Uo + R1; output; end; end; run; proc mixed data=sample5 noclprint covtest method=ML ic; title 'N=160, n_j=5'; class id; model y5 = x /solution outpred=pred5 outpredm=pm5; random intercept / subject=id type=un solution cl; ods output SolutionR=RanUs5; run; proc sort data=RanUs5; by Estimate; run; data stdranU5; set RanUs5; stdU = estimate/StdErrPred; run; data tmp5; set pred5; predij=pred; residij=resid; lowij=lower; upij= upper; stdij=StdErrPred; keep predij residij lowij upij stdij id; run; proc sort data=tmp5; by id; proc sort data=pm5; by id; data predict; merge tmp5 pm5; by id; run; goptions reset=(axis, legend, pattern, symbol, title, footnote) norotate hpos=0 vpos=0 htext=2 ftext=Swiss ctext= target= gaccess= gsfmode= ; goptions device=WIN ctext=blue graphrc interpol=join; symbol1 c=DEFAULT ci=blue v=NONE line=3; symbol2 c=DEFAULT ci=BLACK v=NONE ; axis1 color=black width=2.0 label=(height=2 'Group Centered SES'); axis2 color=black width=2.0 label=(angle=90 height=2 'Predicted Mean'); axis3 color=black width=2.0 label=(angle=90 height=2 'Predicted School Mean'); legend1 frame; proc gplot data=WORK.PREDICT; plot Pred *x / haxis=axis1 vaxis=axis2 frame ; title 'Overall Mean Regression Line: Y= gamma00 + gamma10*X'; run; proc gplot data=WORK.PREDICT; plot predij*x=id / haxis=axis1 vaxis=axis3 frame nolegend; title1 height=2 'Group regression lines'; title2 height=1.5 ' Y = gamma00 + gamma10*X + U_j '; run; data newid; do N=1 to 160; newid = N; output; end; run; proc sort data=RanUs5; by estimate; data RanUxx5; merge RanUs5 newid; run; goptions reset=(axis, legend, pattern, symbol, title, footnote) norotate hpos=0 vpos=0 htext=2 ftext=swiss ctext= target= gaccess= gsfmode= ; goptions device=WIN ctext=black graphrc interpol=join ; symbol1 value=none interpol=join color=black; symbol2 value=plus interpol=none color=black; symbol3 value=plus interpol=none color=black ; symbol4 value=none interpol=join color=black line=3 ci=green; symbol5 value=none interpol=join color=black line=1 ci=red; symbol6 value=none interpol=join color=black line=2 ci=blue; axis1 color=black width=2.0 length=60 pct label=(angle=90 'BLUPS & 95% CI') order=-15 to 15 by 5; axis2 color=black width=2.0 length=60 pct label=('Marco Units ordered by BLUPS'); axis3 color=black width=2.0 ; axis4 color=black width=2.0 length=60 pct label=(angle=90 'OLS & Beta0j'); proc gplot data=WORK.RANUxx5 gout=fig1nj; plot (Estimate Lower Upper ) * newid / overlay haxis=axis2 vaxis=axis1 frame; title1 'N = 160, n_j = 5'; title2 '95% confidence intervals'; run; legend1 position=(inside left top) frame cshadow=grey label=none; goptions reset=(axis, legend, pattern, symbol, title, footnote) norotate hpos=0 vpos=0 htext= ftext=swiss ctext=black target= gaccess= gsfmode= ; goptions device=WIN ctext=black graphrc interpol=join; pattern1 value=SOLID color=black;pattern2 value=SOLID color=black; pattern3 value=SOLID color=black;pattern4 value=SOLID color=black; pattern5 value=SOLID color=black;pattern6 value=SOLID color=black; pattern7 value=SOLID color=black;pattern8 value=SOLID color=black; pattern9 value=SOLID color=black;pattern10 value=SOLID color=black; pattern11 value=SOLID color=black;pattern12 value=SOLID color=black; pattern13 value=SOLID color=black;pattern14 value=SOLID color=black; pattern15 value=SOLID color=black;pattern16 value=SOLID color=black; pattern17 value=SOLID color=black;pattern18 value=SOLID color=black; pattern19 value=SOLID color=black;pattern20 value=SOLID color=black; axis1 length=80 pct color=black width=2.0 label=("Empirical Bayes Estimate of Uo") ; axis21 length=80 pct color=black width=2.0 label=("Sigma^2=1") value=none ; axis3 length=80 pct color=black width=2.0 label=("Sigma^2=9") value=none ; axis5 length=80 pct color=black width=2.0 label=("Sigma^2=25") value=none ; axis10 length=80 pct color=black width=2.0 label=("Sigma^2=100") value=none ; axis2 length=70 pct color=black width=2.0 label=(" ") major=none minor=none value=none; *axis3 color=black width=2.0 ; proc gchart data=WORK.RANUxx5 gout=fig2njaa; VBAR Estimate / maxis=axis1 raxis=axis2 frame type=FREQ patternid=midpoint; title 'N = 160, n_j = 5'; run; ************************************ n_j = 10 *********************************; data sample10; gamma00=12.0; gamma11=2.0; do id=1 to 160; z1=rannor(-1); Uo= 3*z1; do i=1 to 10; z2=rannor(-1); x =rannor(-1); R1=6*z2; y10 = gamma00 + gamma11*x + Uo + R1; output; end; end; run; proc mixed data=sample10 noclprint covtest method=ML ic; title 'N=160, n_j=10'; class id; model y10 = x /solution outpred=pred10 outpredm=pm10; random intercept / subject=id type=un solution cl; ods output SolutionR=RanUs10; run; proc sort data=RanUs10; by Estimate; run; data stdranU10; set RanUs10; stdU = estimate/StdErrPred; run; data tmp10; set pred10; predij=pred; residij=resid; lowij=lower; upij= upper; stdij=StdErrPred; keep predij residij lowij upij stdij id; run; proc sort data=tmp10; by id; proc sort data=pm10; by id; data predict; merge tmp10 pm10; by id; run; goptions reset=(axis, legend, pattern, symbol, title, footnote) norotate hpos=0 vpos=0 htext=2 ftext=Swiss ctext= target= gaccess= gsfmode= ; goptions device=WIN ctext=blue graphrc interpol=join; symbol1 c=DEFAULT ci=blue v=NONE line=3; symbol2 c=DEFAULT ci=BLACK v=NONE ; axis1 color=black width=2.0 label=(height=2 'Group Centered SES'); axis2 color=black width=2.0 label=(angle=90 height=2 'Predicted Mean'); axis3 color=black width=2.0 label=(angle=90 height=2 'Predicted School Mean'); legend1 frame; proc gplot data=WORK.PREDICT; plot Pred *x / haxis=axis1 vaxis=axis2 frame ; title 'Overall Mean Regression Line: Y= gamma00 + gamma10*X'; run; proc gplot data=WORK.PREDICT; plot predij*x=id / haxis=axis1 vaxis=axis3 frame nolegend; title1 height=2 'Group regression lines'; title2 height=1.5 ' Y = gamma00 + gamma10*X + U_j '; run; data newid; do N=1 to 160; newid = N; output; end; run; proc sort data=RanUs10; by estimate; data RanUxx10; merge RanUs10 newid; run; goptions reset=(axis, legend, pattern, symbol, title, footnote) norotate hpos=0 vpos=0 htext=2 ftext=swiss ctext= target= gaccess= gsfmode= ; goptions device=WIN ctext=black graphrc interpol=join ; symbol1 value=none interpol=join color=black; symbol2 value=plus interpol=none color=black; symbol3 value=plus interpol=none color=black ; symbol4 value=none interpol=join color=black line=3 ci=green; symbol5 value=none interpol=join color=black line=1 ci=red; symbol6 value=none interpol=join color=black line=2 ci=blue; axis1 color=black width=2.0 length=60 pct label=(angle=90 'BLUPS & 95% CI') order=-15 to 15 by 5; axis2 color=black width=2.0 length=60 pct label=('Marco Units ordered by BLUPS'); axis3 color=black width=2.0 ; axis4 color=black width=2.0 length=60 pct label=(angle=90 'OLS & Beta0j'); proc gplot data=WORK.RANUxx10 gout=fig1nj; plot (Estimate Lower Upper ) * newid / overlay haxis=axis2 vaxis=axis1 frame; title1 'N = 160, n_j = 10'; title2 '95% confidence intervals'; run; legend1 position=(inside left top) frame cshadow=grey label=none; goptions reset=(axis, legend, pattern, symbol, title, footnote) norotate hpos=0 vpos=0 htext= ftext=swiss ctext=black target= gaccess= gsfmode= ; goptions device=WIN ctext=black graphrc interpol=join; pattern1 value=SOLID color=black;pattern2 value=SOLID color=black; pattern3 value=SOLID color=black;pattern4 value=SOLID color=black; pattern5 value=SOLID color=black;pattern6 value=SOLID color=black; pattern7 value=SOLID color=black;pattern8 value=SOLID color=black; pattern9 value=SOLID color=black;pattern10 value=SOLID color=black; pattern11 value=SOLID color=black;pattern12 value=SOLID color=black; pattern13 value=SOLID color=black;pattern14 value=SOLID color=black; pattern15 value=SOLID color=black;pattern16 value=SOLID color=black; pattern17 value=SOLID color=black;pattern18 value=SOLID color=black; pattern19 value=SOLID color=black;pattern20 value=SOLID color=black; axis1 length=80 pct color=black width=2.0 label=("Empirical Bayes Estimate of Uo") ; axis21 length=80 pct color=black width=2.0 label=("Sigma^2=1") value=none ; axis3 length=80 pct color=black width=2.0 label=("Sigma^2=9") value=none ; axis5 length=80 pct color=black width=2.0 label=("Sigma^2=25") value=none ; axis10 length=80 pct color=black width=2.0 label=("Sigma^2=100") value=none ; axis2 length=70 pct color=black width=2.0 label=(" ") major=none minor=none value=none; *axis3 color=black width=2.0 ; proc gchart data=WORK.RANUxx10 gout=fig2njaa; VBAR Estimate / maxis=axis1 raxis=axis2 frame type=FREQ patternid=midpoint; title 'N = 160, n_j = 10'; run; ************************************ n_j = 100 *********************************; data sample100; gamma00=12.0; gamma11=2.0; do id=1 to 160; z1=rannor(-1); Uo= 3*z1; do i=1 to 100; z2=rannor(-1); x =rannor(-1); R1=6*z2; y100 = gamma00 + gamma11*x + Uo + R1; output; end; end; run; proc mixed data=sample100 noclprint covtest method=ML ic; title 'N=160, n_j=100'; class id; model y100 = x /solution outpred=pred100 outpredm=pm100; random intercept / subject=id type=un solution cl; ods output SolutionR=RanUs100; run; proc sort data=RanUs100; by Estimate; run; data stdranU100; set RanUs100; stdU = estimate/StdErrPred; run; data tmp100; set pred100; predij=pred; residij=resid; lowij=lower; upij= upper; stdij=StdErrPred; keep predij residij lowij upij stdij id; run; proc sort data=tmp100; by id; proc sort data=pm100; by id; data predict; merge tmp100 pm100; by id; run; goptions reset=(axis, legend, pattern, symbol, title, footnote) norotate hpos=0 vpos=0 htext=2 ftext=Swiss ctext= target= gaccess= gsfmode= ; goptions device=WIN ctext=blue graphrc interpol=join; symbol1 c=DEFAULT ci=blue v=NONE line=3; symbol2 c=DEFAULT ci=BLACK v=NONE ; axis1 color=black width=2.0 label=(height=2 'Group Centered SES'); axis2 color=black width=2.0 label=(angle=90 height=2 'Predicted Mean'); axis3 color=black width=2.0 label=(angle=90 height=2 'Predicted School Mean'); legend1 frame; proc gplot data=WORK.PREDICT; plot Pred *x / haxis=axis1 vaxis=axis2 frame ; title 'Overall Mean Regression Line: Y= gamma00 + gamma10*X'; run; proc gplot data=WORK.PREDICT; plot predij*x=id / haxis=axis1 vaxis=axis3 frame nolegend; title1 height=2 'Group regression lines'; title2 height=1.5 ' Y = gamma00 + gamma10*X + U_j '; run; data newid; do N=1 to 160; newid = N; output; end; run; proc sort data=RanUs100; by estimate; data RanUxx100; merge RanUs100 newid; run; goptions reset=(axis, legend, pattern, symbol, title, footnote) norotate hpos=0 vpos=0 htext=2 ftext=swiss ctext= target= gaccess= gsfmode= ; goptions device=WIN ctext=black graphrc interpol=join ; symbol1 value=none interpol=join color=black; symbol2 value=plus interpol=none color=black; symbol3 value=plus interpol=none color=black ; symbol4 value=none interpol=join color=black line=3 ci=green; symbol5 value=none interpol=join color=black line=1 ci=red; symbol6 value=none interpol=join color=black line=2 ci=blue; axis1 color=black width=2.0 length=60 pct label=(angle=90 'BLUPS & 95% CI') order=-15 to 15 by 5 ; axis2 color=black width=2.0 length=60 pct label=('Marco Units ordered by BLUPS'); axis3 color=black width=2.0 ; axis4 color=black width=2.0 length=60 pct label=(angle=90 'OLS & Beta0j'); proc gplot data=WORK.RANUxx100 gout=fig1nj; plot (Estimate Lower Upper ) * newid / overlay haxis=axis2 vaxis=axis1 frame; title1 'N = 160, n_j = 100'; title2 '95% confidence intervals'; run; legend1 position=(inside left top) frame cshadow=grey label=none; goptions reset=(axis, legend, pattern, symbol, title, footnote) norotate hpos=0 vpos=0 htext= ftext=swiss ctext=black target= gaccess= gsfmode= ; goptions device=WIN ctext=black graphrc interpol=join; pattern1 value=SOLID color=black;pattern2 value=SOLID color=black; pattern3 value=SOLID color=black;pattern4 value=SOLID color=black; pattern5 value=SOLID color=black;pattern6 value=SOLID color=black; pattern7 value=SOLID color=black;pattern8 value=SOLID color=black; pattern9 value=SOLID color=black;pattern10 value=SOLID color=black; pattern11 value=SOLID color=black;pattern12 value=SOLID color=black; pattern13 value=SOLID color=black;pattern14 value=SOLID color=black; pattern15 value=SOLID color=black;pattern16 value=SOLID color=black; pattern17 value=SOLID color=black;pattern18 value=SOLID color=black; pattern19 value=SOLID color=black;pattern20 value=SOLID color=black; axis1 length=80 pct color=black width=2.0 label=("Empirical Bayes Estimate of Uo") ; axis21 length=80 pct color=black width=2.0 label=("Sigma^2=1") value=none ; axis3 length=80 pct color=black width=2.0 label=("Sigma^2=9") value=none ; axis5 length=80 pct color=black width=2.0 label=("Sigma^2=25") value=none ; axis10 length=80 pct color=black width=2.0 label=("Sigma^2=100") value=none ; axis2 length=70 pct color=black width=2.0 label=(" ") major=none minor=none value=none; proc gchart data=WORK.RANUxx100 gout=fig2njaa; VBAR Estimate / maxis=axis1 raxis=axis2 frame type=FREQ patternid=midpoint; title 'N = 160, n_j = 100'; run; proc greplay nofs igout=fig1nj gout=fit1 tc=my1cat; tdef t1x2 /* Upper Left Panel */ 1 / llx= 0 ulx= 0 urx=50 lrx=50 lly=50 uly=100 ury=100 lry=50 /* Upper Right Panel */ 2 / llx =50 ulx=50 urx=100 lrx=100 lly=50 uly=100 ury=100 lry=50 /* Lower Left Panel */ 3 / llx=0 ulx=0 urx=50 lrx=50 lly=0 uly=50 ury=50 lry=0 /* Lower Right Panel */ 4 / llx= 50 ulx= 50 urx=100 lrx=100 lly=0 uly=50 ury=50 lry=0 ; template t1x2; tplay 1:1 2:2 3:3 4:4; run; proc greplay nofs igout=fig2njaa gout=fit2 tc=my2cat; tdef t1x2 /* Upper Left Panel */ 1 / llx= 0 ulx= 0 urx=50 lrx=50 lly=50 uly=100 ury=100 lry=50 /* Upper Right Panel */ 2 / llx =50 ulx=50 urx=100 lrx=100 lly=50 uly=100 ury=100 lry=50 /* Lower Left Panel */ 3 / llx=0 ulx=0 urx=50 lrx=50 lly=0 uly=50 ury=50 lry=0 /* Lower Right Panel */ 4 / llx= 50 ulx= 50 urx=100 lrx=100 lly=0 uly=50 ury=50 lry=0 ; template t1x2; tplay 1:1 2:2 3:3 4:4; run; proc mixed data=sample2 noclprint covtest method=ML ic; title 'N=160, n_j=2'; class id; model y2 = x /solution outpred=pred2 outpredm=pm2 ddfm=sather; random intercept / subject=id type=un solution cl; ods output SolutionR=RanUs2; run; proc mixed data=sample5 noclprint covtest method=ML ic; title 'N=160, n_j=5'; class id; model y5 = x /solution outpred=pred2 outpredm=pm2 ddfm=sather; random intercept / subject=id type=un solution cl; ods output SolutionR=RanUs2; run; proc mixed data=sample10 noclprint covtest method=ML ic; title 'N=160, n_j=10'; class id; model y10 = x /solution outpred=pred2 outpredm=pm2 ddfm=sather; random intercept / subject=id type=un solution cl; ods output SolutionR=RanUs2; run; proc mixed data=sample100 noclprint covtest method=ML ic; title 'N=160, n_j=100'; class id; model y100 = x /solution outpred=pred2 outpredm=pm2 ddfm=sather; random intercept / subject=id type=un solution cl; ods output SolutionR=RanUs2; run;