Steve, thanks heaps for the detailed reply. I haven't followed those links yet, because they may not be relevant. What I have done instead is write a simulation (shown below) for a controlled trial in which there are individual responses in an experimental group and in which the change scores in the control and experimental groups are analyzed with two simple mixed models: the first estimates extra variance in the experimental group, representing the individual responses in that group; the second estimates extra variance in the control group, in which case the estimate of the variance is expected to be negative, equal and opposite to the estimate in the first model. It all hangs together nicely. Both models give exactly the same values for the mean effect of the treatment, the residual variance, and the random-effect variance (with a change of sign). So I conclude there is no problem with allowing negative variance when the data need it, and we should allow negative variance anyway, to get unbiased estimates and compatibility limits. (Note that, if the estimate of the variance is positive, allowing only positive variance makes no difference to the estimate and its standard error, but the compatibility limits provided by SAS are unrealistic, because they are based on an unrealistic assumption about the sampling distribution of the variance. The assumption of normality when variance is allowed to be negative produces correct coverage of the interval. Bias would arise when only positive variance is allowed, because negative values of variance arising from sampling variation--or indeed if the true value is negative because of truly less variance in the experimental group--are set to zero.) However, there is a bit of a problem with the random-effect solution when the corresponding variance is negative: SAS doesn't want to give it a standard error (and therefore compatibility limits). This problem occurred occasionally with some but not all values of the random-effect solution in my generalized mixed model with much more data and more complex fixed- and random-effects models, but it occurs every time with every value of the solution in this simulation. I figured that, because of the way the random-effect solution is added to the other effects to get predicted values, there would have to be a negative correlation between the random-effect solution and the residuals, when the variance is negative (the second model), because the residual variance is the variance in the experimental group, so the random-effect solution in the control group has to add to and reduce the values of the residuals in the control group to give the lower SD in the control group. Yes, the correlation was negative alright, but unexpectedly it was perfect, -1.00. In the first model, I didn't expect a correlation, but it was also perfect, this time positive, 1.00. These correlations must arise from the way the variance is partitioned between the residuals and the random effect. I conclude that, at least in this simple individual-responses data and model, the random-effect solution for a negative variance has to be interpreted in the light of the residuals and any other random effects that it is competing with. When the true extra variance is in the experimental group, but the extra variance is estimated in the control group (the second model), the random-effect solution for the negative variance represents individual responses that reduce the residual variance (which represents the variance of the experimental group) rather than increase it, as in the first model. So the random-effect solution for negative variance does not have any immediate useful interpretation, and the lack of compatibility limits is even more problematic for the interpretation of the values. Not shown is the estimation of the variance representing individual responses derived by squaring the SD of the random-effect solution and multiplying by a degrees-of-freedom correction factor, but I did it in a spreadsheet, and the correspondence is good, so the mismatch between the two estimates is just a degrees-of-freedom issue. I am still undecided about whether I should use the random-effect solution as is, or whether I should inflate the values so that their simple SD squared is the same as the covparm variance. Will /*
Simulation of a controlled trial, with individual responses in an exptal group.
Analysis of the change scores is done with a general linear mixed model allowing for extra variance, first in the exptal group, then in the control group.
*/
%let SDdelta=1; *SD of change scores in control group (error of measurement is 1/sqrt(2) times this value);
%let MeanResp=5; *mean response to a treatment in exptal group;
%let IndResp=2; *SD of individual responses in exptal group;
%let GrpSampSize=10; *sample size in control and exptal groups;
%let alpha=0.1; *alpha for compatibility limits;
%let seed=2; *set to 0 for a new random sample every time;
data dat1;
do SubjectID=1 to &GrpSampSize;
Group="Control";
IndResp=0;
Y=&SDdelta*rannor(&seed);
xVarCont=1; *dummy variable to estimate extra variance representing individual responses in control group;
xVarExpt=0; *dummy variable to estimate extra variance representing individual responses in exptal group;
output;
end;
do SubjectID=&GrpSampSize+1 to 2*&GrpSampSize;
Group="Exptal";
IndResp=&IndResp*rannor(&seed);
Y=&SDdelta*rannor(&seed)+&MeanResp+IndResp;
xVarCont=0;
xVarExpt=1;
output;
end;
format _numeric_ 5.2;
title "The data"; proc print;
run;
proc means maxdec=2;
class Group;
run;
title1 "Individual responses estimated in Exptal group";
ods select none;
proc mixed data=dat1 covtest cl nobound alpha=0.1;
class SubjectID Group;
model Y=Group/noint outp=pred ddfm=sat alphap=&alpha ddfm=sat;
random xVarExpt/subject=SubjectID s cl alpha=α
estimate "Treatment effect" Group -1 1/cl alpha=α
ods output covparms=cov;
ods output solutionr=solr;
ods output estimates=est;
run;
ods select all;
title2 "Treatment mean effect (Expt-Control); expected value = &MeanResp";
proc print data=est;
run;
title2 "Covparms, expected values: xVarExpt = &IndResp**2; Residual = &SDdelta**2";
data cov1;
set cov;
DegFree=2*Zvalue**2;
proc print data=cov1;
run;
title2 "Random-effect solution for xVarExpt";
proc print data=solr;
run;
title2 "Predicted and residual values";
proc print data=pred;
run;
title2 "Merge residuals and random-effect solution and get correlation";
data residrand;
merge pred solr(rename=(Alpha=AlphaEst Lower=LowerEst Upper=UpperEst));
by SubjectID;
proc print;
var SubjectID Group IndResp Resid Estimate LowerEst UpperEst AlphaEst;
run;
proc corr;
var Resid;
with Estimate;
by Group;
run;
title1 "Individual responses estimated in Control group";
ods select none;
proc mixed data=dat1 covtest cl nobound alpha=0.1;
class SubjectID Group;
model Y=Group/noint outp=pred ddfm=sat alphap=&alpha ddfm=sat;
random xVarCont/subject=SubjectID s cl alpha=α
estimate "Treatment effect" Group -1 1/cl alpha=α
ods output covparms=cov;
ods output solutionr=solr;
ods output estimates=est;
run;
ods select all;
title2 "Treatment mean effect (Expt-Control); expected value = &MeanResp";
proc print data=est;
run;
title2 "Covparms, expected values: xVarCont = -&IndResp**2; Residual = &SDdelta**2+&IndResp**2";
data cov1;
set cov;
DegFree=2*Zvalue**2;
proc print data=cov1;
run;
title2 "Random-effect solution for xVarCont";
proc print data=solr;
run;
title2 "Predicted and residual values";
proc print data=pred;
run;
title2 "Merge residuals and random-effect solution and get correlation";
data residrand;
merge pred solr(rename=(Alpha=AlphaEst Lower=LowerEst Upper=UpperEst));
by SubjectID;
proc print;
var SubjectID Group IndResp Resid Estimate LowerEst UpperEst AlphaEst;
run;
proc corr;
var Resid;
with Estimate;
by Group;
run;
... View more
