[R-sig-ME] Standard errors of variances; negative variance

[Will Hopkins]

Someone emailed me personally to point out a mistake in my SAS random-effects code for individual responses. I had random xVarExpt*SubjectID, which is correct, but it is equivalent to random xVarExpt/subject=SubjectID, not intercept/subject=SubjectID. Sorry about that. It's distressingly easy to make such simple errors, but hopefully they become obvious when you run the program and get silly results, infinite likelihood, or failed to converge. 

If you model original scores rather than change scores, the fixed effect is Group*Trial and there are two random effects: SubjectID xVarExpt*SubjectID (or intercept xVarExpt/subject=SubjectID) along with the residual. By the way, with a big-enough sample size, you can estimate individual responses with this model for the original scores and data for only one group. The fixed effect is then simply Trial. Of course, there is a lot of uncertainty, depending on the sample size and the relative magnitudes of the between-subject SD and error-of-measurement (residual) SD. It still amazes me that Proc Mixed can partition out three variances from two SDs on the same subjects. If you don't estimate individual responses directly with a dummy variable, you can estimate a different residual variance in the pre- and post-test. The individual-response variance is the difference between the residual variances. Does lme have the option for specifying different random effects and/or residuals in different groups? In SAS its simply done with the /group= syntax (here repeated/group=Trial).

Stevedrd has raised the issue of adjusting for the pretest value. Yes, that improves precision by accounting for regression to the mean with change scores, when error of measurement is substantial relative to the between-subject SD in the pretest. I routinely include the pretest when modeling change scores; the fixed effects are then Group Group*PreTest. You thereby estimate the individual responses with improved precision, but it has no effect on their expected value, unless the pretest really does modify the treatment effect (which sometimes happens, especially with human performance). The resulting individual responses are those not explained by the pretest and by any other effect modifiers you put in the fixed-effects model. You can't include the pretest when you model original scores, unless you make the post-test value the dependent and adjust for the pretest; I don't normally do it that way, because modeling change scores is easier for people to understand and to visualize in scatterplots, IMHO.

In response to Ben Bolker... I said merDeriv's results may be untrustworthy, because you said previously that "the standard errors are often extremely poor summaries of the uncertainty or RE variances". I thought your comment was about the standard errors produced by merDeriv, not any standard errors for variances produced by mixed models. I can't remember what package Alice found for producing standard errors with lme, but it was around 2017 and so it may have been merDeriv. As I said in my first message, the SEs weren't the same as SAS's, so I assumed "extremely poor" referred to merDeriv. As for the standard errors being extremely poor summaries in general, I don't know where that assessment comes from. In my experience with SAS they are excellent summaries. For simple models you get exactly what you would expect for the degrees of freedom, and the confidence intervals have acceptable coverage. SAS assumes a normal sampling distribution for the variance, centered on the observed variance, not the log of the variance.

In response to Ben Peizer... The change score is the dependent variable in a controlled trial, with one observation per subject in the control group and one per subject in the experimental group. If the SD of the change scores in the experimental group is less than that in the control group, and you specify a model that estimates extra variance in the experimental group with a dummy variable, you will get negative variance. If instead you specify extra variance in the control group, you will get equal and opposite variance. The SE is the same. Why bother to specify extra variance in the experimental group, if you can see that the variance in that group is less than that in the control group?  Because there probably are real individual responses in the experimental group, and you want to see how big they could be, by estimating their upper confidence limit.  The observed SD of change scores in the experimental group may be less than that in the control group simply because of sampling variation, especially for small sample sizes. 

I like John Maindonald's comparison with square root of -1. I wonder if it is possible to work with SDs as complex numbers. I can't see how at the moment.

Will

-----Original Message-----
From: Will Hopkins <willthekiwi using gmail.com> 
Sent: Wednesday, June 21, 2023 12:08 PM
To: r-sig-mixed-models using r-project.org
Subject: RE: [R-sig-ME] Standard errors of variances; negative variance

Thanks heaps for the seven replies from five people on this list. In summary, it looks like standard errors could be derived for variances in R: with asreml-R (which you have to pay for?); possibly with nlme; with merDeriv, but they may be untrustworthy; and with Bayesian packages (but I'd be specifying weakly informative priors that in the limit have no effect on the posterior, which might not work in some packages). Negative variance seems to be off the table, however. Ben, James and Timothee have indicated scenarios where negative variance made some sense. 

I had a particular reason for seeking advice on this list; see the bottom of my message. First, here's my justification for negative variance and standard errors of variances.

My favorite example is the variance representing individual responses to a treatment, derived as the variance of the change scores in an experimental group minus the variance of the change scores in a control group. In a mixed model, the dependent is the change score, there is a fixed effect for group (to estimate the difference in the mean changes), other fixed-effect modifiers if you want them, and a random effect for the interaction of a dummy variable (1 in the experimental group, 0 in the control group) with subject identity. I usually call the dummy variable something like xVarExpt, to indicate extra variance in the experimental change scores, and the code in SAS is random xVarExpt*SubjectID; or random intercept/subject=SubjectID;. And of course there is the residual. 

Usually the variance of the change scores in the experimental group is greater than that in the control group, owing to individual responses. But with small sample sizes, the variance on the experimental group could be less than that in the control, simply because of sampling variation. You will therefore get a negative variance point estimate. Importantly, its upper confidence limit will tell you how big and positive the individual responses could be. The confidence interval for variances estimated by Proc Mixed in SAS have good coverage in my experience with many simulations over the years. (The coverage isn't quite so good in meta-analytic mixed models when the study-estimate sample sizes are small, but the coverage is generally conservative, i.e., what you specify as a 90%CI has greater than 90% coverage, as I recall.) The coverage works well ONLY when you allow negative variance for the random effects. Interestingly, the default in SAS is positive-only variance, and if you ask for confidence limits, they are calculated from a standard error and the chi-squared distribution. The resulting confidence limits are unrealistic: you can get enormous upper limits. When you allow negative variance, they are realistic; quoting from the on-line documentation, "Wald Z-scores and normal quantiles are used to construct the limits". (See  https://documentation.sas.com/doc/en/pgmsascdc/9.4_3.3/statug/statug_mixed_syntax01.htm There is no further information about how the standard errors are calculated and no mention of profile confidence limits.) The confidence limits for the residual(s) are always calculated from the chi-squared distribution, unless you specify otherwise, and these are always realistic and have correct coverage, in my experience.

The variance in the experimental group could truly be less than that in the control group, either because the treatment compresses the responses towards some upper or lower limit, or less likely because there are individual responses to the control treatment that exceed those to the experimental treatment. Either way negative variance is meaningful. You can turn the dummy variable around if you like, but you don't have to.

When I get a negative variance, or a negative lower confidence limit, I have no hesitation in changing the sign, taking the square root, and calling it a negative SD. It just simply means that the individual responses, or more generally whatever differences or variability are represented by the random effect, could be reduced differences or variability. I know a negative SD for the intercept random effect is a hard one to swallow, but it's better to see that and think again about your model and/or what might really be happening than to simply let it be set to zero. 

I would like to emphasize just how important the random effects are in the analysis of repeated measurement or other clustered data. That's where the individual differences and individual responses are found, which are crucial in this age of personalized medicine. Hence you need to make inferences about the random effects; they are not simply nuisance parameters that you have to include to make sure the uncertainty in the mean effects are trustworthy, then forget about. I make inferences about magnitudes of the random effects (Bayesian or hypothesis testing against smallest importants) assuming a normal distribution (chi-squared for residuals) and smallest importants for SDs that are half those for mean effects.

OK, enough of the hobby horse. I wanted to find out whether R allows for negative variances and standard errors of variances, because I am working with a colleague on a manuscript focusing on the correct way to analyze, report and interpret errors of measurement in various kinds of reliability study. Example: for some measures, the first measurement in a series of repeated measurements can have more error, which is easily specified and estimated with a dummy variable, as explained above for individual responses. Another example: the possibility of additional error arising between time-separated clusters of repeated measurements needs to be specified as a random effect that can have negative variance. We've simulated data and analyzed them all with SAS Studio, and we hope to publish the programs as supplementary material. We're looking for a co-author/collaborator who had experience of measurement studies and skill with mixed models in R to provide the same simulations and analyses with R. It's contingent upon estimating negative variance and standard errors of variances, both of which require more than lme or nlme, by the look of it. Maybe it's a good opportunity for someone to get involved and figure out the easiest way to get the same estimates and uncertainties as in SAS, because that would or should be useful for mixed modeling more generally with R.

Will