That's the drawback to the PL approach--it is hard to get metrics to compare models, as the pseudo-data is not the same under the two models, thus even looking at the -2*log pseudo-likelihood values is problematic.  Now if it were fit using integration methods, you could look at information criteria to come up with a decision--one that incorporates both model complexity and information completeness. Have you tried METHOD=LAPLACE in the PROC GLIMMIX statement?  It may run into memory problems, or execution time problems, or convergence problems... To compare the two models, you will need to have the RANDOM statements in subject= format. So, 1: random intercept hospital/subject=physician; and 2: random intercept/subject=hospital; random intercept/subject=physician; These two models compute variance components due to physician and physician-hospital in the first case, and due to hospital and physician in the second case (no correlation).  I guess you could throw in a third that only considers physician-hospital dyads as: 3: random intercept/subject=physician*hospital; See if any of this works with the Laplacian integration method. Steve Denham ... View more