Skip to content

Standard errors of variances; negative variance

16 messages · Ben Bolker, Thierry Onkelinx, Sorkin, John +6 more

Original
Will Hopkins
# 
I've just joined this list to get answers to a few questions. I would have
searched the archive before posting, but there seems to be no way of
searching except via quarterly summaries. I searched the last four quarters
without success.

I'm a SAS user, but a few years ago I tried the mixed model in R, with the
help of an R user (Alice Sweeting). At that time, the lme package did not
provide standard errors for the variances, nor did it allow negative
variance. Have these limitations been addressed?  As I recall, Alice found
some code that provided standard errors, but it gave values different from
those of the mixed model in SAS (Proc Mixed). I therefore opted to stay with
SAS, because allowing for negative variance for random effects other than
residuals and estimating uncertainties in variances are both fundamental to
mixed modeling, in my view. I also became fluent with SAS coding over the
years and did not want to make the effort with R coding. Interestingly, SAS
introduced a free cloud version called SAS Studio, evidently modeled on R
Studio, to try to win back customers!

Will
Ben Bolker
# 
Hi, Will,

Googling

  site:https://stat.ethz.ch/pipermail/r-sig-mixed-models/ "negative 
variance"

claims to get you 67 results (although only 11 that Google considers 
worth displaying by default),

   You can filter by date - there are only two hits since 2020:

https://www.google.com/search?q=site%3Ahttps%3A%2F%2Fstat.ethz.ch%2Fpipermail%2Fr-sig-mixed-models%2F+%22negative+variance%22&client=firefox-b-d&tbs=cdr%3A1%2Ccd_min%3A01-01-2020%2Ccd_max%3A&tbm= 


   and these both look like false positives (they include "negative" and 
"variance" but not "negative variance" ...)

   The short answer is that I am not aware of any mixed effect package 
in R that will allow you to return negative values. You can see my 
answer here: 
https://stat.ethz.ch/pipermail/r-sig-mixed-models/2018q1/026437.html ...


   As for uncertainties in variances - the merDeriv package (dating back 
to 2017) will give you the standard errors of variances and covariances 
(although again, note that there's a theoretical issue here - the 
standard errors are often extremely poor summaries of the uncertainty or 
RE variances, the original authors of lme4 would certainly prefer that 
you use profile confidence intervals to summarize the uncertainty ...)

   That's what I know -- perhaps someone else knows about a package that 
allows for negative variances ... ??
On 2023-06-19 6:13 p.m., Will Hopkins wrote:
Thierry Onkelinx
# 
Dear all,

A negative variance seems odd to me. Isn't a variance positive by
definition? How would you interpret a random intercept with a negative
variance?

Best regards,

ir. Thierry Onkelinx
Statisticus / Statistician

Vlaamse Overheid / Government of Flanders
INSTITUUT VOOR NATUUR- EN BOSONDERZOEK / RESEARCH INSTITUTE FOR NATURE AND
FOREST
Team Biometrie & Kwaliteitszorg / Team Biometrics & Quality Assurance
thierry.onkelinx at inbo.be
Havenlaan 88 bus 73, 1000 Brussel
www.inbo.be

///////////////////////////////////////////////////////////////////////////////////////////
To call in the statistician after the experiment is done may be no more
than asking him to perform a post-mortem examination: he may be able to say
what the experiment died of. ~ Sir Ronald Aylmer Fisher
The plural of anecdote is not data. ~ Roger Brinner
The combination of some data and an aching desire for an answer does not
ensure that a reasonable answer can be extracted from a given body of data.
~ John Tukey
///////////////////////////////////////////////////////////////////////////////////////////

<https://www.inbo.be>


Op di 20 jun 2023 om 01:36 schreef Ben Bolker <bbolker at gmail.com>:

  
  


  


      

    

    

      
      

        


  

        

      Sorkin, John
    

    

      
      #
    

  


  

      
Thierry,
I have the same question that you posed, how can a variance, a measure that is squared be negative? Will, please enlighten us!
John
John David Sorkin M.D., Ph.D.
Professor of Medicine
Chief, Biostatistics and Informatics
University of Maryland School of Medicine Division of Gerontology and Geriatric Medicinet
Baltimore VA Medical Center
10 North Greene Street<x-apple-data-detectors://12>
GRECC<x-apple-data-detectors://12> (BT/18/GR)
Baltimore, MD 21201-1524<x-apple-data-detectors://13/0>
(Phone) 410-605-711<tel:410-605-7119>9
(Fax) 410-605-7913<tel:410-605-7913> (Please call phone number above prior to faxing)

        
On Jun 20, 2023, at 3:45 AM, Thierry Onkelinx via R-sig-mixed-models <r-sig-mixed-models at r-project.org> wrote:

        
?Dear all,

A negative variance seems odd to me. Isn't a variance positive by
definition? How would you interpret a random intercept with a negative
variance?

Best regards,

ir. Thierry Onkelinx
Statisticus / Statistician

Vlaamse Overheid / Government of Flanders
INSTITUUT VOOR NATUUR- EN BOSONDERZOEK / RESEARCH INSTITUTE FOR NATURE AND
FOREST
Team Biometrie & Kwaliteitszorg / Team Biometrics & Quality Assurance
thierry.onkelinx at inbo.be
Havenlaan 88 bus 73, 1000 Brussel
https://nam11.safelinks.protection.outlook.com/?url=http%3A%2F%2Fwww.inbo.be%2F&data=05%7C01%7Cjsorkin%40som.umaryland.edu%7C2c0d1047636b43c2bfcc08db71625f82%7C717009a620de461a88940312a395cac9%7C0%7C0%7C638228439560019013%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000%7C%7C%7C&sdata=sbP2Xemwg0iHjvxcc6zR8vtPbZf5ZoFCv1BIFIeDdBE%3D&reserved=0

///////////////////////////////////////////////////////////////////////////////////////////
To call in the statistician after the experiment is done may be no more
than asking him to perform a post-mortem examination: he may be able to say
what the experiment died of. ~ Sir Ronald Aylmer Fisher
The plural of anecdote is not data. ~ Roger Brinner
The combination of some data and an aching desire for an answer does not
ensure that a reasonable answer can be extracted from a given body of data.
~ John Tukey
///////////////////////////////////////////////////////////////////////////////////////////

<https://nam11.safelinks.protection.outlook.com/?url=https%3A%2F%2Fwww.inbo.be%2F&data=05%7C01%7Cjsorkin%40som.umaryland.edu%7C2c0d1047636b43c2bfcc08db71625f82%7C717009a620de461a88940312a395cac9%7C0%7C0%7C638228439560019013%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000%7C%7C%7C&sdata=I5HJ05o0Bd4ly0GL0Smc1iYm3cuZMc2IGCcxy8OWn%2B8%3D&reserved=0>


Op di 20 jun 2023 om 01:36 schreef Ben Bolker <bbolker at gmail.com>:

  Hi, Will,

Googling

 site:https://nam11.safelinks.protection.outlook.com/?url=https%3A%2F%2Fstat.ethz.ch%2Fpipermail%2Fr-sig-mixed-models%2F&data=05%7C01%7Cjsorkin%40som.umaryland.edu%7C2c0d1047636b43c2bfcc08db71625f82%7C717009a620de461a88940312a395cac9%7C0%7C0%7C638228439560019013%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000%7C%7C%7C&sdata=vw%2FnMxnvA7L9A0qYuJfHa9suA2YBTnZ%2BqxPltyR3G%2FU%3D&reserved=0 "negative
variance"

claims to get you 67 results (although only 11 that Google considers
worth displaying by default),

  You can filter by date - there are only two hits since 2020:


https://nam11.safelinks.protection.outlook.com/?url=https%3A%2F%2Fwww.google.com%2Fsearch%3Fq%3Dsite%253Ahttps%253A%252F%252Fstat.ethz.ch%252Fpipermail%252Fr-sig-mixed-models%252F%2B%2522negative%2Bvariance%2522%26client%3Dfirefox-b-d%26tbs%3Dcdr%253A1%252Ccd_min%253A01-01-2020%252Ccd_max%253A%26tbm%3D&data=05%7C01%7Cjsorkin%40som.umaryland.edu%7C2c0d1047636b43c2bfcc08db71625f82%7C717009a620de461a88940312a395cac9%7C0%7C0%7C638228439560019013%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000%7C%7C%7C&sdata=xbWrjRwClEMYFKVkncDFgZga6SPNcmsPvbNs4kpUQjY%3D&reserved=0


  and these both look like false positives (they include "negative" and
"variance" but not "negative variance" ...)

  The short answer is that I am not aware of any mixed effect package
in R that will allow you to return negative values. You can see my
answer here:
https://nam11.safelinks.protection.outlook.com/?url=https%3A%2F%2Fstat.ethz.ch%2Fpipermail%2Fr-sig-mixed-models%2F2018q1%2F026437.html&data=05%7C01%7Cjsorkin%40som.umaryland.edu%7C2c0d1047636b43c2bfcc08db71625f82%7C717009a620de461a88940312a395cac9%7C0%7C0%7C638228439560019013%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000%7C%7C%7C&sdata=2E5S51uIXiPEsuH11f693MEOLYFb2dh5e39t2GSLor8%3D&reserved=0 ...


  As for uncertainties in variances - the merDeriv package (dating back
to 2017) will give you the standard errors of variances and covariances
(although again, note that there's a theoretical issue here - the
standard errors are often extremely poor summaries of the uncertainty or
RE variances, the original authors of lme4 would certainly prefer that
you use profile confidence intervals to summarize the uncertainty ...)

  That's what I know -- perhaps someone else knows about a package that
allows for negative variances ... ??

        
On 2023-06-19 6:13 p.m., Will Hopkins wrote:

        
I've just joined this list to get answers to a few questions. I would
have
searched the archive before posting, but there seems to be no way of
searching except via quarterly summaries. I searched the last four
quarters
without success.

I'm a SAS user, but a few years ago I tried the mixed model in R, with
the
help of an R user (Alice Sweeting). At that time, the lme package did not
provide standard errors for the variances, nor did it allow negative
variance. Have these limitations been addressed?  As I recall, Alice
found
some code that provided standard errors, but it gave values different
from
those of the mixed model in SAS (Proc Mixed). I therefore opted to stay
with
SAS, because allowing for negative variance for random effects other than
residuals and estimating uncertainties in variances are both fundamental
to
mixed modeling, in my view. I also became fluent with SAS coding over the
years and did not want to make the effort with R coding. Interestingly,
SAS
introduced a free cloud version called SAS Studio, evidently modeled on R
Studio, to try to win back customers!

Will

_______________________________________________
R-sig-mixed-models at r-project.org mailing list
https://nam11.safelinks.protection.outlook.com/?url=https%3A%2F%2Fstat.ethz.ch%2Fmailman%2Flistinfo%2Fr-sig-mixed-models&data=05%7C01%7Cjsorkin%40som.umaryland.edu%7C2c0d1047636b43c2bfcc08db71625f82%7C717009a620de461a88940312a395cac9%7C0%7C0%7C638228439560019013%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000%7C%7C%7C&sdata=cPlRiZIOr3%2B4lUfOPHm26RrP4peEw9ZAkjtm5enO4Uw%3D&reserved=0

_______________________________________________
R-sig-mixed-models at r-project.org mailing list
https://nam11.safelinks.protection.outlook.com/?url=https%3A%2F%2Fstat.ethz.ch%2Fmailman%2Flistinfo%2Fr-sig-mixed-models&data=05%7C01%7Cjsorkin%40som.umaryland.edu%7C2c0d1047636b43c2bfcc08db71625f82%7C717009a620de461a88940312a395cac9%7C0%7C0%7C638228439560019013%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000%7C%7C%7C&sdata=cPlRiZIOr3%2B4lUfOPHm26RrP4peEw9ZAkjtm5enO4Uw%3D&reserved=0



_______________________________________________
R-sig-mixed-models at r-project.org mailing list
https://nam11.safelinks.protection.outlook.com/?url=https%3A%2F%2Fstat.ethz.ch%2Fmailman%2Flistinfo%2Fr-sig-mixed-models&data=05%7C01%7Cjsorkin%40som.umaryland.edu%7C2c0d1047636b43c2bfcc08db71625f82%7C717009a620de461a88940312a395cac9%7C0%7C0%7C638228439560019013%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000%7C%7C%7C&sdata=cPlRiZIOr3%2B4lUfOPHm26RrP4peEw9ZAkjtm5enO4Uw%3D&reserved=0

  
  


  


      

    

    

      
      

        


  

        

      Ben Bolker
    

    

      
      #
    

  


  

      
From Molenberghs and Verbeke 2011:

In the first view, where the focus is entirely on the resulting marginal 
model (1.2), negative values for ?2 are perfectly acceptable (Nelder, 
1954; Verbeke and Molenberghs, 2000, Section 5.6.2), because this merely 
corresponds to the occurrence of negative within-cluster correlation ? = 
?2/(?2 + ?2). In such a case, the only requirement is that ?2 + ?2 > 
0andVi = ?2 Jni + ?2 Ini is a positive definite, marginal covariance 
matrix. Further discussions on negative variance components, and their 
implications for inferences, can be found in Thompson (1962), Searle et 
al. (1996), Verbeke and Molenberghs (2003) and Molenberghs and Verbeke 
(2007). In the second view, when the link between marginal model (1.2) 
and its generating hierarchical model (1.1) is preserved, thereby 
including the concept of random effects bi and perhaps even requiring 
inferences about them, it has been considered imperative to restrict ?2 
to non-negative values.

   Molenberghs, Geert, and Geert Verbeke. 2011. ?A Note on a 
Hierarchical Interpretation for Negative Variance Components.? 
Statistical Modelling 11 (5): 389?408. 
https://doi.org/10.1177/1471082X1001100501.

    Since lme4 deals with the hierarchical model, I don't think that 
negative variances are going to work.

        
On 2023-06-20 9:21 a.m., Sorkin, John wrote:

            

      
      
      
    

        

      
      
    

            

      
      
      
    

        

      
      
    

            

      
      
      
    

        

      
      
    

        

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      James Pustejovsky
    

    

      
      #
    

  


  

      
Regarding standard errors for variance components: Will, I'm not sure if
you were asking about lme4::lmer() or nlme::lme(). If you're not familiar
with the latter, you might find it interesting. It is in some respects more
flexible than lmer(), such as providing model components for different
level-1 error structures and correlation structures, and so you might find
it more comparable to SAS. If you're using nlme::lme(), the lmeInfo package
provides a variance-covariance matrix for the variance component parameters
(based on the inverse expected information or average information):
https://jepusto.github.io/lmeInfo/
That said, I agree with the other folks who've suggested that likelihood
ratio tests and profile likelihood CIs might be a better choice for
inference on variance components.

Regarding negative variances, Ben's post is instructive. In addition, in
some specific problems related to meta-analysis, we've found that allowing
for negative variance components can improve the performance of score and
likelihood ratio tests involving other model parameters. My working theory
here is that bounding variances at zero means that the log likelihood is no
longer smooth, which seems to  muck up the behavior or tests that involve
the first and second derivatives of the log likelihood.

James

        
On Tue, Jun 20, 2023 at 8:37?AM Ben Bolker <bbolker at gmail.com> wrote:

        

            

      
      
      
    

        

      
      
    

            

      
      
      
    

        

      
      
    

            

      
      
      
    

        

      
      
    

            

      
      
      
    

        

      
      
    

        

  
  


  


      

    

    

      
      

        


  

        

      Will Hopkins
    

    

      
      #
    

  


  

      
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

  
  


  


      

    

    

      
      

        


  

        

      Ben Bolker
    

    

      
      #
    

  


  

      
A few *short* responses below.

        
On 2023-06-20 8:08 p.m., Will Hopkins wrote:

            

      
      
      
    

        
not sure what the evidence is that merDeriv's results are 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.

            

      
      
      
    

        
Interestingly, the default in SAS is positive-only variance, and if you 
ask f

            

      
      
      
    

        
I would argue that an equally good (possibly even better?) way to 
test the maximum plausible size of the effect is computing an upper 
*likelihood profile* confidence limit.  This should not require allowing 
a negative estimate of the variance to get good coverage ... ??  (I can 
certainly see why a Wald-based confidence interval would fail if you 
didn't allow negative variances.  How does SAS compute these intervals? 
Are they assumed to be symmetric on the variance scale or on the 
log-variance scale?)

            

      
      
      
    

        
This argument could presumably be turned around to say that we 
shouldn't rely on a method that only considers the marginal variance 
structure (i.e., the context in which negative variances make the most 
sense), but that we want to consider the 'true' structure of the 
variance (in which case we should probably be considering 
compound-symmetric variance structures ...)

            

      
      
      
    

        
An alternative approach to this would be to model the 
heteroscedasticity directly, which can be done with (e.g.) lme and 
glmmTMB, although not at present in lme4 - no need to make up latent 
variables with negative variances ...

   cheers
     Ben Bolker

            

      
      
      
    

        

      
      
    

  
  


  


      

    

    

      
      

        


  

        

      stevedrd m@iii@g oii y@hoo@com
    

    

      
      #
    

  


  

      
Just sayin'
Part of the "negative variance" issue in this particular case could be due to using change scores as the dependent variable. Using the "pre" value as a covariate adjuster of the post values should go a long way toward alleviating this issue. See Frank Harrell's website/blog for better reasoning than I am able to do.
Steve Denham

        
On Tuesday, June 20, 2023 at 08:59:43 PM EDT, Ben Bolker <bbolker at gmail.com> wrote:

        
? A few *short* responses below.

        
On 2023-06-20 8:08 p.m., Will Hopkins wrote:

            

      
      
      
    

        
? not sure what the evidence is that merDeriv's results are 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.

            

      
      
      
    

        
Interestingly, the default in SAS is positive-only variance, and if you 
ask f

            

      
      
      
    

        
? I would argue that an equally good (possibly even better?) way to 
test the maximum plausible size of the effect is computing an upper 
*likelihood profile* confidence limit.? This should not require allowing 
a negative estimate of the variance to get good coverage ... ??? (I can 
certainly see why a Wald-based confidence interval would fail if you 
didn't allow negative variances.? How does SAS compute these intervals? 
Are they assumed to be symmetric on the variance scale or on the 
log-variance scale?)

            

      
      
      
    

        
? This argument could presumably be turned around to say that we 
shouldn't rely on a method that only considers the marginal variance 
structure (i.e., the context in which negative variances make the most 
sense), but that we want to consider the 'true' structure of the 
variance (in which case we should probably be considering 
compound-symmetric variance structures ...)

            

      
      
      
    

        
wi

            

      
      
      
    

        
? An alternative approach to this would be to model the 
heteroscedasticity directly, which can be done with (e.g.) lme and 
glmmTMB, although not at present in lme4 - no need to make up latent 
variables with negative variances ...

? cheers
? ? Ben Bolker

            

      
      
      
    

        

      
      
    

        
_______________________________________________
R-sig-mixed-models at r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models

  
  


  


      

    

    

      
      

        


  

        

      John Maindonald
    

    

      
      #
    

  


  

      
Dear  All,
As with the square root of -1, negative variances are a convenient 
computational fiction.  The square root of -1 is an invitation to move
to a 2-dimensional word, where they make perfect sense.

A negative variance estimate, in those systems that allow it, is a useful 
warning that one needs to think and move beyond the one-dimensional 
representation in which the variances have been set.
John Maindonald
Statistics Research Associates, Wellington NZ.

            

      
      
      
    

        

      
      
    

  
  


  


      

    

    

      
      

        


  

        

      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 at gmail.com> 
Sent: Wednesday, June 21, 2023 12:08 PM
To: r-sig-mixed-models at 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

  
  


  


      

    

    

      
      

        


  

        

      John Maindonald
    

    

      
      #
    

  


  

      
Another way to look at the issue is that variance components are just
that ? components that cannot necessarily be interpreted as variances.
They cannot even necessarily be interpreted in a particularly meaningful 
way as variances, even if the values come out as positive.
John Maindonald
Statistics Research Associates, Wellington NZ.

Dear  All,
As with the square root of -1, negative variances are a convenient 
computational fiction.  The square root of -1 is an invitation to move
to a 2-dimensional word, where they make perfect sense.

A negative variance estimate, in those systems that allow it, is a useful 
warning that one needs to think and move beyond the one-dimensional 
representation in which the variances have been set.
John Maindonald
Statistics Research Associates, Wellington NZ.

            

      
      
      
    

        

      
      
    

  
  


  


      

    

    

      
      

        


  

        

      J.D. Haltigan
    

    

      
      #
    

  


  

      
This negative variance issue crops up in factor analysis oftentimes with
Heywood cases etc. One approach that has often been invoked in this
literature is to set them to zero, though it has been some time since I
have delved deep into this literature.

On Wed, Jun 21, 2023 at 9:16?PM John H Maindonald <jhmaindonald at gmail.com>
wrote:

            

      
      
      
    

        

      
      
    

            

      
      
      
    

        

      
      
    

            

      
      
      
    

        

      
      
    

            

      
      
      
    

        

      
      
    

            

      
      
      
    

        

      
      
    

        

  
  


  


      

    

    

      
      

        


  

        

      Ben Pelzer
    

    

      
      #
    

  


  

      
Hi Will,

I was wondering if your "strange variance" result is basically comparable
to the following simple R script. Data are simulated, with no higher level
at all. For each "subject" there is only one single observation. Trying to
estimate both the residual variance SE and also a random intercept variance
with lmer (after switching off all mechanisms which prevent a user to
estimate such a "strange" model) and lme, leads to quite arbitrary results
for the two variance estimates. Obviously, estimating two quantities when
there exists only one is problematic.... so why would I like to estimate
such a model? But maybe your situation is different, a reproducible example
(data + sas code) might be enlightening. Best regards, Ben P.

y <- rnorm(100, 0, 1)
subject <- 1:100

library(lme4)
m1 <- lmer(y ~ (1|subject),
           control=lmerControl(check.nobs.vs.nlev = "ignore",
                               check.nobs.vs.rankZ = "ignore",
                               check.nobs.vs.nRE="ignore"))
summary(m1)

library(nlme)
m2 <- lme(y ~ 1, random=~1|subject)
summary(m2)

        
On Thu, 22 Jun 2023 at 00:25, Will Hopkins <willthekiwi at gmail.com> wrote:

        

            

      
      
      
    

        

      
      
    

        

  
  


  


      

    

    

      
      

        


  

        

      Will Hopkins
    

    

      
      #
    

  


  

      
In response to John Maindonald's latest post... John, you suggested that the variance components "cannot necessarily be interpreted as variances". I don?t think so. For any real data, variance components can and should be interpreted as the square of standard deviations, which have a direct interpretation in terms of differences between and/or changes within subjects and/or settings. You only have to simulate data to realize that: all the random effects are generated randomly with a mean of zero and a given SD. (Admittedly the data have to be generated with positive SDs. I'm not sure how to simulate data directly with a negative SD, although I can do it indirectly by adding random variability in the "wrong" or unexpected place.) To allow for sampling variation in the estimation of the SD, you have to allow negative variance. When you do that, you estimate the simulated SD without bias, and the confidence intervals have correct coverage. If you don't allow negative variance, you will get bias and unrealistic upper confidence limits, with small sample sizes anyway. Negative point estimates and negative lower confidence limits are realistic in the sense that they properly allow for sampling uncertainty, and in any case they may be real, in the sense that the extra variance could be in a different place from what you expected and modeled.

in response to Ben Pelzer's latest post... Ben, you can't get the extra variance with only one observation per subject in one group. You can get it with two observations on the same subject in one group; so for a group of subjects measured in two trials, without necessarily anything else happening, you can estimate the variance representing true differences between subjects with random int/subject=SubjectID and two residual variances with repeated/group=Trial. And with the same data and exactly the same results, you can estimate one residual variance and the difference between the residuals as extra variance, either on the first trial or the second trial, and it will positive or negative variance, depending on which residual was greater with the repeated/group=Trial analysis. (The variances add up exactly, but the predicted values differ very slightly between the models where the extra variance is positive and negative.) Here is the SAS code for doing all that, written with macro variables to allow you to change values of means, SDs and sample size. Remember that sample size has to large to get good precision for the SDs. I've set the sample size to 100. If you make it 10, there's huge uncertainty, and sometimes it gives infinite likelihood (although you can probably overcome that with suitable starting values for the covparms). Note that I have not added extra variance representing individual responses in Trial 2 (as shown with %let ir=0;), so half the simulations will produce negative extra variance in Trial 2, and you can play with estimating the extra variance in either of the trials to see that you get the same extra variance differing only in sign, and the same standard error.

*Simple test-retest reliability, allowing for different residuals in the two trials
 and individual responses in the second trial;
*Modify the program as shown in the proc mixed step to estimate one residual
 and extra variance one or other trial;

%let mean=100; *true grand mean in the first test;
%let sd=10; *true between-subject SD;
%let tem=5; *typical (standard) error of measurement;
%let deltamean=2; *change in mean test2-test1;
%let ir=0; *SD for individual responses (in Trial 2);
%let ss=100; *sample size;

*generate the data;
data dat1;
do SubjectID=1 to &ss;
  TruePerform=&mean+&sd*rannor(0);
  Trial=1;
  IndivResp=0;
  ObsvdPerform=TruePerform+&tem*rannor(0)+IndivResp;
  output;
  Trial=2;
  IndivResp=&ir*rannor(0);
  ObsvdPerform=TruePerform+&tem*rannor(0)+&deltamean+IndivResp;
  output;
  end;

title "Raw data, first 10 observations";
proc print data=dat1(obs=10);
format TruePerform ObsvdPerform IndivResp 5.1;
run;

title "Simple stats";
proc means maxdec=1;
var TruePerform ObsvdPerform IndivResp;
class Trial;
run;

*Set up dummy variables for estimating extra variance in Trial 1 or Trial 2;
data dat2;
set dat1;
xVarTrial1=0;
if Trial=1 then xVarTrial1=1;
xVarTrial2=0;
if Trial=2 then xVarTrial2=1;

ods select none;
title "Rely mixed model, ssize=&ss";
title2 "True mean=&mean, true between SD=&sd, TEM=&tem, delta mean=&deltamean";
proc mixed data=dat2 covtest cl alpha=0.1 nobound;
class SubjectID Trial;
model ObsvdPerform=Trial/noint ddfm=sat solution cl alpha=0.1 outp=pred alphap=0.1;
random int/subject=SubjectID s cl alpha=0.1; *This estimates true differences between subjects;
repeated/group=Trial; *This estimates different residuals in Trial 1 and 2;
*Star off the above two lines and unstar one or other (but not both) of the following two lines...;
*...to estimate one residual and extra variance on either Trial 1 or Trial 2;
*random int xVarTrial1/subject=SubjectID s cl alpha=0.1; *extra error on Trial 1;
*random int xVarTrial2/subject=SubjectID s cl alpha=0.1; *or extra variance on Trial 1;
estimate "Trial 2-1" Trial -1 1/cl alpha=0.1;
ods output classlevels=clev;
ods output covparms=cov;
ods output solutionf=solf;
ods output solutionr=solr;
ods output estimates=est;
run;
ods select all;

title3 "Class levels";
proc print data=clev;
run;

title3 "Covparms";
proc print data=cov;
run;

title3 "Predicteds and residuals, first 10 observations";
proc print data=pred(obs=10);
run;

title3 "Fixed-effects solution";
proc print data=solf;
run;

title3 "Random-effects solution, first 10 observations";
proc print data=solr(obs=10);
run;

title3 "Estimates";
proc print data=est;
run;

  
  


  


      

    

    

      
      

        


  

        

      John Maindonald
    

    

      
      #
    

  


  

      
What I am saying is that if in a randomized  block design blocks are 
chosen to be overly heterogeneous, the within blocks component 
of variance will be exaggerated, and the between blocks component
reduced accordingly. A consequence is that the variance components
apply only to the particular layout of blocks and of plots within blocks
used for that particular experiment.  

There may in such cases be ways to use the data, together with what 
one knows of the field layout, to get a wider perspective on what the 
data may have to say. One is then moving outside of the confines of 
the analysis that gave the problematic ?variances?, whether negative 
or just distorted.  These sorts of issues carry over in different ways to
different contexts.

John Maindonald
Statistics Research Associates, Wellington NZ.