Whoops, that should have been       Yhat = b0 + b1*Var1 + b2*Var2 + theta + delta*Var2 or       Yhat = (b0 + theta) + b1*Var1 + (b2 + delta)*Var2 I left off Var2 at the end of the first equation (the product delta*Var2). Does this result in the same predicted values? ... View more

I don't think your model is what you believe it to be. It is fine to have Var1 as fixed and Var2 as random. But you appear to be operating under the assumption that variables specified on the RANDOM statement are assumed to have a distribution with nonzero mean. For a model written as       Y{ij} = X{ij}*beta + gamma{i} + epsilon{ij} where i indexes store and j indexes an observation within the i-th store, it appears that you believe that the random effect gamma{i} has distribution       gamma{i} ~ N(mu, V) That is, it appears that you believe that there is a nonzero expectation for gamma{i} that is implicit in the model. Because you believe that there is a nonzero expectation for gamma{i}, you have specified NOINT on the MODEL statement (assuming that the random intercept would have nonzero expectation) and have not specified Var2 on the MODEL statement. That is NOT how the MIXED procedure operates. Effects specified on the RANDOM statement are assumed to have distribution       gamma{i} ~ N(0, V) To be even more explicit, the code which you have specified estimates the following model:       Sales{ij} = b1*Var1{ij} + theta{i} + delta{i}*Var2{ij} + epsilon{ij}       theta{i} ~ N(0, V{theta})       delta{i} ~ N(0, V{delta})       cov(theta{i},delta{i}) = rho*sqrt(V{theta}*V{delta}) If you remove the NOINT option from the MODEL statement and include both Var1 and Var2 on the MODEL statement, then you would fit what is almost certainly a better model (and definitely NOT A WRONG model) where you have       theta{i} ~ N(mu{theta}, V{theta})       delta{i} ~ N(mu{delta}, V{delta})       cov(theta{i},delta{i}) = rho*sqrt(V{theta}*V{delta}) Now, as to your original question about the MIXED procedure returning values which are different from the values which you compute, I don't see anything in the code you supplied which indicates how you constructed predicted values either from the MIXED procedure or in a data step (using parameter estimates from the MIXED procedure). You don't show either an OUTP or OUTPM option on the MODEL statement. You also don't show any ODS OUPUT statements for capturing the fixed effect estimates (ODS OUTPUT solutionF=...) or for capturing the random effect estimates (ODS OUTPUT solutionR=...). Without seeing exactly what you requested and exactly how you computed predicted values, one can only speculate on the reason for the differences. My guess is that you requested the marginal expectation (using OUTPM=) and then constructed conditional predicted values (computing b1*Var1 + theta + delta*Var2) OR you requested the conditional predicted values (using OUTP=) and then constructed marginal predicted values (computing b1*Var1). You have to be consistent about whether you do or do not include the random effects when constructing predicted values. Finally, I would note that your PROC MIXED code did not name your subject effect (Store_number) on a CLASS statement. That is OK if your data are sorted by Store_number. But if your data are not sorted by Store_number, then you either need to first sort by Store_number or name Store_number on a CLASS statement. To address all of the issues which I have indicated above, I would revise your model as PROC MIXED DATA = mylib.data;   class store_number;   model sales = Var1 Var2 / solution cl outp=Preds_RE outpm=Preds_FE;   random intercept Var2 / type=un subject=Store_Number solution cl gcorr; run; Note that the data set Preds_RE represent predictions which include the random effects and is obtained as       Yhat = b0 + b1*Var1 + b2*Var2 + theta + delta while Preds_FE represent predictions which include only the fixed effects:       Yhat = b0 + b1*Var1 + b2*Var2 ... View more

Susan, Thanks for posting the example. I had tried using the NOPROFILE option previously without success. What I did not have is the RANDOM statement specifying the residual with a unique subject specification. It is nice to know that it is possible to obtain a test of the residual variance using the COVTEST statement. However, it is not the most convenient syntax. ... View more

In addition to the COVTEST statement of the GLIMMIX procedure (which does not seem to function when the only variance parameter is a residual variance), the UNIVARIATE procedure will construct 1-tailed and 2-tailed confidence limits for the variance of a response which is assumed to be normally distributed. That does not directly test a null hypothesis about the variance being a specified value. However, if you construct an alpha=0.05 CI and the CI includes the value under the null hypothesis, then you would not reject the null hypothesis at alpha=0.05. If your null hypothesis has a two-tailed alternative (and you can assume that the response is normally distributed), then you can download from SAS a macro which will compute the p-value for the specified null. See http://support.sas.com/kb/25/024.html. The macro will not test a null hypothesis with a 1-sided alternative. ... View more

Your SAS code is fine. However, the (default) method of estimating the random effects which is employed by the GLIMMIX procedure (essentially a Penalized Quasi-Likelihood estimator) is known to produce random effect variance estimates which are biased toward zero. The GLLAMM program employs Gaussian quadrature to estimate the variances of the random effects. Gaussian quadrature does not suffer the same variance bias as does PQL. If you have only a few students per school (say, 10), then it is possible to obtain quadrature-based estimates of the variances of the random effects from SAS because such models can be written with schools as the subject and then specifying as the random effects an intercept and student within school where student within school is numbered from 1 to n{school} (max(n{school})=10). I would be happy to clarify how to obtain quadrature-based variance estimates if you have such a limited number of students per school. But my guess is that you have more than 10 students/school. ... View more

Yes, that would be appropriate. ... View more

Data step code could certainly construct the frequencies which you need. proc sort data=mydata;   by facilityID SSN; run; data freqs;   set mydata;   by facilityID ssn;   if first.facilityID then do;     NRestr1=0;        /* Number of times restraint 1 used */     NPatRestr1=0;   /*Number of patients for whom restraint 1 was used */     NRestr2=0;        /* Number of times restraint 2 used */     NPatRestr2=0;   /*Number of patients for whom restraint 2 was used */     ...   end;   if first.SSN then do;     PatRestr1=0;   /* Restraint 1 ever used on patient */     PatRestr2=0;   /* Restraint 2 ever used on patient */     ...   end;   NRestr1 + Restr1;   /* Increment total use of restraint 1 */   NRestr2 + Restr2;   /* Increment total use of restraint 2 */   ...   PatRestr1 = max(PatRestr1, Restr1);   /* Update indicator of patient ever using restraint 1 */   PatRestr2 = max(PatRestr2, Restr2);   /* Update indicator of patient ever using restraint 2 */   ...   if last.SSN then do;     NPatRestr1 + PatRestr1;   /* Increment count of patients using restraint 1 */     NPatRestr2 + PatRestr2;   /* Increment count of patients using restraint 2 */     ...   end;   if last.FacilityID then output;   keep facilityID NRestr1 NpatRestr1 NRestr2 NpatRestr2; run; proc print data=freqs; run; ... View more

I don't believe that there is any direct way to convert parameter estimates obtained from the MIANALYZE procedure into standardized parameter estimates. However, you could use a process as described below to obtain standardized parameter estimates: 1) impute data using PROC MI 2) standardize the response and predictors to have mean 0 and variance 1 using PROC STANDARD 3) fit the regression model to each of the standardized imputation data sets outputting parameter estimates and the covariance matrix of the parameter estimates. Any of a number of procedures could be used for this step. 4) Combine standardized regression coefficients from each imputation set into a single estimate of the standardized regression coefficients using PROC MI ... View more

You have to specify the HOLD= option on the PARMS statement. I assume that you have specified an unstructured covariance matrix for the intercept, slope, and slope^2 and that the terms in your random model are named in that order. I also assume that conditional on the random effects, the residuals are IID (that you do not have a REPEATED statement modifying the residual variance structure). If the above is correct, then the covariance of slope and slope^2 random effects will be the 6-th and last random effects covariance parameter in the table of covariance parameters, and that there is one parameter (residual variance) that appears in the table after the random effect covariance estimates. For covariance estimates , , ..., , (where explicit values should be assumed for each of these terms), you could write your MIXED code as: proc mixed data=mydata;   class subjectID ...;   model response = ...;   random intercept slope slope2 / subject=subjectID type=un;   parms ( ) ( ) ( ) ( ) ( ) (0) ( ) /             hold=6; run; ... View more

Well, I just Googled "ake carpenter phreg" (without the quotes - the quotes here serve only to indicate the terms that were in the Google search). Low and behold, right on top was a a paper with the title "Extending the use of PROC PHREG in survival analysis" which was presented at a SUG. ... View more