/* Modeling of n=23 schools from Kreft & de Leeuw */ title 'Baseline Model (Null/Empty HLM)'; proc mixed data=school23c method=ml covtest noclprint ic; class school ; model math = / solution ddfm=satterth; random intercept / subject=school type=un; run; title 'Preliminary Linear Mixed Model'; proc mixed data=school23c method=ml covtest noclprint ic; class school gender whitec public; model math = homew cses gender whitec cses*whitec meanses public / solution ddfm=satterth; random intercept homew / subject=school type=un; run; /******* Other random *********/ title 'cses'; proc mixed data=school23c method=ml covtest noclprint ic; class school gender whitec public; model math = homew cses whitec cses*whitec meanses public / solution ddfm=satterth; random intercept cses / subject=school type=un; run; title 'gender'; proc mixed data=school23c method=ml covtest noclprint ic; class school gender whitec public; model math = homew cses whitec cses*whitec meanses public / solution ddfm=satterth; random gender / subject=school type=un; run; title 'whitec'; proc mixed data=school23c method=ml covtest noclprint ic; class school gender whitec public; model math = homew cses whitec cses*whitec meanses public / solution ddfm=satterth; random whitec / subject=school type=un; run; /************** Reduced model to test U_1j: HO: tau_11 = tau_10 =0 *********/ title 'Reduced Preliminary Linear Mixed Model'; proc mixed data=school23c method=ml covtest noclprint ic; class school gender whitec public; model math = homew cses gender whitec cses*whitec meanses public / solution ddfm=satterth; random intercept / subject=school type=un; run; title 'Test for random cses'; data mixpval; n2lnlike0=3673.6011; n2lnlike1=3598.2708; LR = n2lnlike0-n2lnlike1; df2=2; df1=1; p1 = 1 - cdf('chisquare',LR,df1); p2 = 1 - cdf('chisquare',LR,df2); mix = .5*(p1+p2); run; proc print; run; /************* Remove meanSES from preliminary model & observed PUBLIC *******/ title 'No mean SES'; proc mixed data=school23c method=ml covtest noclprint ic; class school gender whitec public; model math = homew cses gender whitec cses*whitec public / solution ddfm=satterth; random intercept homew / subject=school type=un; run; /******** Using constrast to test fixed effects */ title 'Preliminary Linear Mixed Model'; proc mixed data=school23c method=ml covtest noclprint ic; class school gender whitec public; model math = homew cses gender whitec cses*whitec meanses public / solution ddfm=satterth; random intercept homew / subject=school type=un; contrast 'gender & sector' public 1 -1, gender 1 -1; run; /********* Drop public & gender ****/ title 'Reduced Model1'; proc mixed data=school23c method=ml covtest noclprint ic; class school whitec; model math = homew cses whitec cses*whitec meanses / solution ddfm=satterth; random intercept homew / subject=school type=un; run; data pval; LR = 0.5015; df=2; pvalue = 1 - cdf('chisquare',LR,df); proc print; run; title 'Reduced Model = Final Model '; ods html; ods graphics on; proc mixed data=school23c method=ml covtest noclprint ic namelen=200; class school white; model math = homew cses white cses*white meanses / solution ddfm=satterth; random intercept homew / subject=school type=un solution g; ods output SolutionR=Us CovParms=cov G=gmat ModelInfo=mod SolutionF=solf; run; %hlmrsq(CovParms=cov,GMatrix=gmat,ModelInfo=mod,SolutionF=solf); run; ods graphics off; ods html close; title 'Reduced Model = Final Model '; ods html; ods graphics on; proc mixed data=school23c method=ml covtest noclprint ic namelen=200; class school white; model math = homew cses white meanses / solution ddfm=satterth; random intercept homew / subject=school type=un solution g; ods output SolutionR=Us CovParms=cov G=gmat ModelInfo=mod SolutionF=solf; run; %hlmrsq(CovParms=cov,GMatrix=gmat,ModelInfo=mod,SolutionF=solf); run; ods graphics off; ods html close; /******** Take a look at Uhats ********/ data U0 U1; set Us; lower=estimate- 1.96*StdErrPred; upper=estimate+ 1.96*StdErrPred; if Effect='Intercept' then output U0; else output U1; run; data integers; do i=1 to 23; output; end; run; proc sort data=U0; by estimate; proc sort data=U1; by estimate; data U0; merge integers U0; U0=estimate; run; data U1; merge integers U1; u1=estimate; run; title1 'Estimated U_o and their 95% Confidence intervals'; title2 'Ordered by their values of U_o'; proc sgplot data=U0;; scatter x=i y=Estimate / YERRORLOWER=lower YERRORUPPER=upper; run; title1 'Estimated U_1 and their 95% Confidence intervals'; title2 'Ordered by their values of U_1'; proc sgplot data=U1;; scatter x=i y=Estimate / YERRORLOWER=lower YERRORUPPER=upper; run; /******** Histrograms won't look good with only N=23 *******/ proc sgplot data=u0; histogram Estimate; run; proc sgplot data=u1; histogram Estimate; run; /******** Compare U0 and U1 ***********/ proc sort data=u0; by school; proc sort data=u1; by school; data U0andU1; merge u0 u1; by school; run; title 'U_1 versus U_0'; proc sgplot data=U0andU1; scatter x=u0 y=u1; run;