Wow. This is a very information dense post, and I don't think I can address all of it, so let's try some pieces. First, nesting in SAS. You need to look at the documentation for MIXED or GLM that covers this. The parameterization of A nested in B and A crossed with B is the same. That is why it makes no difference in the rcorr. Here is what the online manual says: Nested Effects Nested effects are generated in the same manner as crossed effects. Hence, the design columns generated by the following two statements are the same (but the ordering of the columns is different): model Y=A B(A); model Y=A A*B; The nesting operator in PROC MIXED is more a notational convenience than an operation distinct from crossing. Nested effects are typically characterized by the property that the nested variables never appear as main effects. The order of the variables within nesting parentheses is made to correspond to the order of these variables in the CLASS statement. The order of the columns is such that variables outside the parentheses index faster than those inside the parentheses, and the rightmost nested variables index faster than the leftmost variables (Table 56.18). Table 56.18 Example of Nested Effects Data (I deleted the table as it doesn't come across very well) Note that nested effects are often distinguished from interaction effects by the implied randomization structure of the design. That is, they usually indicate random effects within a fixed-effects framework. The fact that random effects can be modeled directly in the RANDOM statement might make the specification of nested effects in the MODEL statement unnecessary. The other part that I want to address is the doubly repeated measures. This assumes that the correlation of residuals amongst days is the same, no matter which period we have. I was thinking that observations on the close time points would be more highly related, and probably similar, so that an ar(1) or cs structure might fit. The separated time points (i.e., periods) would have a different relationship. If you want to deal with all at once without those assumptions, you have unequal spacing and the spatial power method goes after that with the assumption that the errors are "constantly correlated" as in ar(1), but that the points used to estimate that correlation are not equally spaced. A read of the TYPE= section of the PROC MIXED documentation might point you at the right kind of structure. Good luck on all of this. Steve Denham ... View more

Aha. Now I understand part of what is going on. Think about the repeated statement: repeated day(period)/subject=animal(farm) type=cs; and how SAS will parameterize this statement. The parameterization is the same as: repeated day*period/subject=animal*farm type=cs; This makes sense to me that you would get the same results if you respecify as you mention. Try adding the rcorr option after the slash, and take a look at the correlations under your various approaches. They should be identical, as they are restatements of the same approach. Some alternate approaches might be: repeated day*period/subject=animal*farm type=csh repeated day*period/subject=animal*farm type=sp(pow)(time) where time is actual elapsed days post study start or even: repeated period day/subject=animal*farm type=UN@AR(1) which would model an unstructured relationship between periods, but a common autoregressive relationship within periods. Note also that for the latter two structures, you should include animal*farm in the random statement, so that between animal variability is correctly modeled. Good luck. Steve Denham ... View more

I have some thoughts, but I don't view any of them as a real solution. They might give some insight into what's happening if you give them a try. First, what happens to the estimates if you lift the constraint on R to the fixed values? Does the model converge without problems, and if so, how much change do you see in those parameters. Second, have you tried type=FA0(2) as a covariance structure? It is constrained to be non-negative definite, so maybe something good would happen. Third, along the same lines, change over to PROC GLIMMIX, and try type=chol, for a Cholesky root parameterization, this would give at least a positive semi-definite matrix. Fourth, now that you are in GLIMMIX, try at least a couple of different optimization methods. I don't know that any of these will succeed as I'm not familiar with the data, but since I've long since gone through the effort of pulling out my hair on problems like this, I thought I'd give it a shot as to what I might try. Good luck. Steve Denham Associate Director, Biostatistics MPI Research, Inc. ... View more

I can think of at least two ways. First, do not specify it as a parameter, but instead write the model with the "forced" value in place. Second, let it be a parameter, but use the BOUNDS statement to restrict the value. Good luck. Steve Denham ... View more

This would be my preference as well. The algorithm could also prune branches that are obvious dead ends as well, as in the exact algorithms in FREQ, MULTTEST and LOGISTIC. But this bangs up against theoretical computer science, in an interesting way. We have traded a very tedious polynomial time process for a non-polynomial time one (a classical stopping problem). That doesn't mean it will take longer--it's just interesting. ... View more

Even the smaller problem is not tractable, if we are to examine all possible combinations. Let's assume that you have access to a petaflop speed machine, and it takes 10 floating point operations to characterize a single combination. That means we could generate 10^14 combinations in one second. Using Ramanujan's approximation for log n!, I get something like 10^2106 years to generate all possible 900 record subsets. So let's not look at all possibles, it is just too daunting. Instead, let's see how many subsets we would have to look at to have a 50% chance of finding one that has at least 5 equal values. First, using a hypergeometric distribution, the probability that any single subset of 900 has at least 5 equal values is 4.13*10^-9. Now we consider a geometric distribution, apply some approximations, and see that we need to generate about 121 million successive subsets to have a roughly 50% chance of finding one that has at least five equal values. That might be doable in SAS in a reasonable amount of time The original problem is still intractable. The probability of at least 5 equal values in this situation is on the order of 10^-40, so we are looking at about 2.6*10^40 subsets to have a roughly 50% chance of finding one. On the petaflop machine, this is still 10^25 seconds or 3.2*10^17 years. ... View more

Anne, What criteria are you going to use to decide a ranking of the subsets? Is it critical to determine the absolute best on that scale, or would any of, say, 100 subsets that were approximately equal be sufficient? And you are absolutely right that Binomial(158 792, 5 000 ) is a very big number--you are in the cosmological realm with that. And that is what makes me wonder why you would need the "best" subset, rather than a properly stratified random sample. Good luck. Steve Denham ... View more

Hey Peter, Look at this (code added to your original): proc format; value racfmt 1 = 'Black' 2 = 'White' 3 = 'Latino'; run; data today; input catv1 catv2 catv3 @@; dv = catv1 * 3 + catv2 * 5 + catv3 + rannor(123); format catv3 racfmt.; datalines; 0 1 2 0 1 1 0 1 3 1 0 1 0 0 3 0 1 3 0 1 2 1 0 3 1 0 1 0 1 2 0 0 3 1 1 2 1 1 1 0 1 3 1 1 1 1 0 3 0 1 3 1 1 2 1 0 3 1 0 1 0 1 2 1 1 3 0 0 2 1 0 1 1 1 3 0 0 1 0 0 3 0 0 3 0 0 2 1 0 3 1 0 1 0 1 2 1 1 3 ; run; title 'Version with all on CLASS statement'; proc glm data = today; class catv1 catv2 catv3; model dv = catv1 catv2 catv3; lsmeans catv1 catv2 catv3; run; title 'Version with only race on CLASS statement'; proc glm data = today; class catv3; model dv = catv1 catv2 catv3; lsmeans catv3; run; proc means data=today; var dv catv1 catv2 catv3; run; title 'Version with only race on CLASS statement, with at=0.5'; title2 'Results same as all on CLASS statement'; proc glm data = today; class catv3; model dv = catv1 catv2 catv3; lsmeans catv3/at (catv1 catv2)=(0.5 0.5); run; title 'Version with only race on CLASS statement, with at= '; title2 'Results same as only race on CLASS statement'; proc glm data = today; class catv3; model dv = catv1 catv2 catv3; lsmeans catv3/at (catv1 catv2)=(0.48484848485 0.5151515151515); run; So the difference is in the solution for the OLS equations. The first calculates lsmeans with equal weighting by class membership (which was the whole point of Searle, Speed and Milliken (1980), I think), while the second calculates lsmeans at the mean value. I would bet that the large differences in your real data arise from substantial differences in class size. Good luck. ... View more

David, That approach looks pretty good, but there might be an even better approach. Is there any reason to believe that observations within a station and site are correlated from year to year? If so, you may want to expand the model to something like: class = site reach year station model environmental data = site|reach|year random site/subject=station repeated year/subject=station*site type=un Good luck. Steve Denham ... View more