It is very unusual to define a G matrix based on continuous variables and no subject. This would not be your usual random-coefficient model. Your model is likely overparameterized with that many random effects. Without a subject= specification for the a's, the fixed and random terms are competing to explain the same variability. These comments apply to either version of your code. One of the warning signs is that you are getting a missing (.) for the estimated SE for your estimated variance-covariance parameters. What warnings or notes are you getting in your Log? This does not preclude using these parameters, but I would worry about this. It is also true that your second model is identical to the first, and you are simply asking for additional output. This should not affect your model fit, obviously since the G option only kicks in when deciding what to display. I am sure the problems go back to my initial comments. ... View more

You need to use the solution option: model Running= Treatment*Day / solution ; This will show you the parameters being estimated, and the order. ... View more

These are totally different hypotheses being tested. With      Treat Day Treat*Day, test of the last term is strictly for the interaction (not for main effects). But if you don't have main effects, the test of      Treat*Day is for everything (i.e., a test of 'group' where group corresponds to a combination of treat, day, and the interaction). The significant result you now get makes perfect sense. It does not necessarily mean you have an interaction. I think you should get the book "SAS for Mixed Models, 2nd edition" by Littell et al. (2006). ... View more

I am bothered whey you say the F test for fixed effects are giving you garbage. What do you mean? If you are using the means model as we described, then the F test is for the effect of 'group', where group is the combination of treatment and time. That is, are there differences among the groups? You should be getting realistic numbers for F and the df. If not, you are doing something wrong. Contrasts can get very tricky with interaction means. You have to be careful about how the levels are sorted, based on the order in the class statement and the values of the codes. You can go very wrong, and I can't tell. But if there are not 14 days for each treatment, there would be less than the 5*14 values. All of this assumes that the data you have are unbiased measures of the actual random variable. ... View more

With missing cells, there are different approaches one could take, and no simple guidance for you (even assuming that you have missing at random). There will be detractors and proponents for each possibility. If there is good evidence that you do not have an interaction, then you could just model the main effect(s). model  = Treatment Day; lsmeans Treatment Day / diff; Results would be biased if you really have an interaction (but fail to detect it because of poor power). If you don't have many missing cells at the end, you could just end the data set at the last time where you have at least one observation in each treatment-time cell. But if this means giving up too many times for some treatments, then this would not be a great approach. One could use multiple imputation to fill in the cells m times, and use the methodology for imputations. Or you could use the full data set and the means model (just treatment*days, no main effect), as I described above. You could use some contrasts for treatment, but this will be tricky to code if you are not used to contrast syntax in mixed. I can't get into this. Given that the missing data are not random, this raises all kinds of other issues. There is a big literature on this for mixed models, but you don't have any easy solutions. You probably will need to work with a statistician to work this out.  I am sure I forgot some other approaches. ... View more

model Treatment*Day; only gives Treatment*Day, no main effect of Treatment and no main effect of Day. This is often called a means model in comparison to an effect model (what you coded). If you have combinations treatment and day with no observations, you cannot have the main effect lsmeans (they are not estimable, by definition, if you also have an interaction in the model). ... View more

You can fit models with missing data, but you cannot obtain least squares means for main effects when interaction cells are completely empty (when you have interactions in the model). There are different approaches to handle this. One common way is to fit the 'means model'. Just use the interaction term: model  = Treatment*Day; lsmeans Treatment*Day / diff; You can write contrast or estimate statements for selected contrasts of interest. ... View more

See my response to "Output Standard Deviation Rather ..." which I posted a few minutes ago. It is on the same topic. ... View more

There is very little information in the GLIMMIX documentation on this. My previous posting from ages ago basically tells you what is given in the documentation. In essence, for the fixed effects, GLIMMIX apparently scales (standardizes) all the X variables (including interactions)  internally; this would make the optimization methods more stable (avoiding extremely large or small values). This scaling involves dividing each X by the square root of the sum of squares of the relevant X variables.This is not the usual standardization in linear models where the standard deviation of the Y is also utilized. After model fitting, 'unscaling' is done to get the parameters back on the scale of the original X (this is what is displayed). If you put in the STDCOEF option, you will get the scaled parameter estimates. You can read about this in the help documentation for GLIMMIX. Go to STDCOEF (under MODEL), and then go to NOCENTER option to see the equations. Centering means subtracting mean X (it can be done with or without this). I don't think you should overly interpret the standardized parameter estimates. These are not the same as those you get in linear models (they are not unitless because the Y units are still embedded in the parameters, based on the documentation). However, since  the response variable is the same for all the X predictors, then I suppose that one can roughly compare the magnitude of the standardized parameter estimates. But I would be careful even doing that. ... View more

Regarding your last comment, although you can use the data as you coded, you should not do that. With GLMMs (what you have), there can be large biases in fixed effect parameter estimates (not relevant for GLMs, as in PROC LOGISTIC). See SAS for Mixed Models, 2nd edtion (2006), or the new book by Walter Stroup. The larger the n (number of surgeries per doctor), the less the bias when there are random effects. The only time there is a really high bias with a GLMM is with the situation you have (equivalent to n=1). So be careful. Plus, you really should not do a goodness of fit test with the data in your form. The ICC only makes sense when your data are in the y/n format. The ICC roughly gives the average correlation of binary status of the individuals within the lowest sampling unit. But with your approach, there is one individual within this sampling unit (the success of the surgery). A correlation of one value makes no sense. I will have to look up a citation on other versions of ICC. I did not say that department or doctor does not have an effect. I gave you code to properly test for a random department effect (you did not test for this, based on your post). But I did say it is hard to test for department when there are only two departments. The results you show indicate that doctor within department has an effect. This is really the same as an interaction of doctor and department. Most of the effect of doctor(department) could be from the department or the doctors. You can consider departments directly, but you can only consider doctors within departments. One last point for now. You mention random error. For you, doctor(department) is the random error, there is no other random error. You talked about variability, with an ANOVA mindset. I understand what you are trying to get at, but with GLMM, one must move away form ideas like explained variability. THere is a different mindset to fully appreciate the jump from LM to GLM to GLMM. Stroup's new book on generalized linear mixed models explains this in detail. Good luck. ... View more

One slight mistake in my previous message. If you use the y/n format, then use dist=binomial, not dist=binary. If you use the y/n format for the data, then it makes more sense to look at the conditional Pearson chi-squared / df for goodness of fit. ... View more

Since you have multiple observations per doctor, I think you should reformat your data file in the y/n format. That is, each record should correspond to a doctor (within a department). Give the number of successful surgeries (y), the number of surgeries for that doctor (n), treatment, department, and anything else. Then the model would be: proc glimmix data = ... method = quad;      class department doctor_id treatment;      model y/n= treatment / dist = binary oddsratio solution;      random int / solution subject = doctor_id(department); run; or random int doctor_id / subject=department; if you want to try the hierarchy. Things are usually better behaved, and the conditional model (what you have by using random effects) has much better properties when n (number of surgeries) is greater than 1. Your coding makes n=1 since an observatin can be either 0 or 1; this can lead to large biases for GLMMs. With these changes, your wish to calculate the intra-cluster correlation may be more appropriate. You are taking the approach that there is an underlying standardized logistic distribution with your correlation formula (where 3.29 is pi^2/3). There are other approaches to this based on different assumptions. ... View more

I don't have time to fully assess your design, but will give you some pointers on your code, and the interpretation of the results. The message about the G matrix is really a note or warning, not an error (usually). In your case, it simply means that your between-doctor variance is 0 (4.4E-13 means 4.4*10^-13, which is basically 0). I would take out the doctor term, as you did, because your results indicate no substantial doctor effect.In fact, it does NOT make sense to include this term. As you describe the design, doctors are nested within departments. There is no real meaning to doctor as a stand-alone term in the model. If you labeled the doctors uniquely across the two departments (for instance, that there cannot be a doctor "1" in each department; perhaps doctors 1-10 in department A and doctors 11-20 in department B), then you could use either random int / solution subject = doctor_id(department); or random int / solution subject = doctor_id; It would be the same thing in terms of the model for the random effect. This would not be true if there are doctors 1-10 in A and 1-10 in B (for instance). Basically, with nesting, there is an interaction, a consequence of your particular design. With your chosen model, you can't address this (they are tied together). However, you could fit your model with random effects for department and doctor within department: random int doctor_id / solution subject = department; which gives you two random effects (department and doctor nested within department). This can tell you about the effect of department. But with only two departments, the variance estimate for department will be very imprecise. Not sure how much you can interpret it (it may have a large SE). You could get a 0 for one of the variances. ... View more

This has come up before. Do some searching of this website. You should read chapter 12 in SAS for Mixed Models (2006) by Littell et al. and chapter 16 in Generalized Linear Mixed Models (2012) by Walter Stroup. You should not use GLMPOWER or POWER procedures when there are random effects. ... View more