Skip to content

Another case of -1.0 correlation of random effects

28 messages · Ken Knoblauch, Ben Bolker, Jarrod Hadfield +4 more

Original

Messages 1–25 of 28

Kevin E. Thorpe
# 
Hello.

I know this has come up a couple times recently, but I'm still not sure 
what to do about it in my data.  Note that my sessionInfo() will be at 
the bottom.

My data come from a crossover trial and are balanced.

 > str(gluc)
'data.frame':	96 obs. of  4 variables:
  $ Subject  : int  1 2 3 5 6 7 10 11 12 13 ...
  $ Treatment: Factor w/ 2 levels "Barley","Oat": 1 1 1 1 1 1 1 1 1 1 ...
  $ Dose     : int  8 8 8 8 8 8 8 8 8 8 ...
  $ iAUC     : num  110 256 129 207 244 ...

 > xtabs(~Treatment+Dose,data=gluc)
          Dose
Treatment  0  2  4  8
    Barley 12 12 12 12
    Oat    12 12 12 12

I plot the data (attached as gluc.pdf, if it comes through).

 From the plot, I think I want to fit the model as:

lmer(iAUC~Treatment+Dose+(Treatment|Subject)+(Dose|Subject),data=gluc)

It could possibly be argued that the (Treatment|Subject) part is not 
needed.  When I fit this, I got -1.0 correlation within the Dose random 
effects.  To simplify, I will fit a simpler model, since the issue persists.

 > lmer(iAUC~Dose+(Dose|Subject),data=gluc,subset=Treatment=="Oat")
Linear mixed model fit by REML
Formula: iAUC ~ Dose + (Dose | Subject)
    Data: gluc
  Subset: Treatment == "Oat"
    AIC   BIC logLik deviance REMLdev
  562.6 573.9 -275.3    563.1   550.6
Random effects:
  Groups   Name        Variance Std.Dev. Corr
  Subject  (Intercept) 8274.324 90.9633
           Dose          16.214  4.0266  -1.000
  Residual             4862.319 69.7303
Number of obs: 48, groups: Subject, 12

Fixed effects:
             Estimate Std. Error t value
(Intercept)  309.352     30.539  10.130
Dose         -14.424      3.596  -4.012

Correlation of Fixed Effects:
      (Intr)
Dose -0.647

Now, a plot created by (and attached as lmlist.pdf):

plot(confint(lmList(iAUC~Dose|Subject,data=gluc,subset=Treatment=="Oat"),pooled=TRUE),order=1)

shows (I think) a strong negative correlation between the intercept and 
slope random effects for Dose.

So, I would appreciate some advice on how I might specify these random 
effects correctly.

One last thing I tried.  If I treat Dose as a factor (which might be 
reasonable) rather than numeric, I don't get any -1.0 correlations.

 > lmer(iAUC~dose+(dose|Subject),data=gluc,subset=Treatment=="Oat")
Linear mixed model fit by REML
Formula: iAUC ~ dose + (dose | Subject)
    Data: gluc
  Subset: Treatment == "Oat"
    AIC   BIC logLik deviance REMLdev
  545.2 573.3 -257.6      547   515.2
Random effects:
  Groups   Name        Variance Std.Dev. Corr
  Subject  (Intercept)  7509.9   86.660
           dose2       11993.0  109.513  -0.321
           dose4        6399.5   79.997   0.043  0.873
           dose8        6051.7   77.793  -0.743  0.433  0.306
  Residual              1206.4   34.733
Number of obs: 48, groups: Subject, 12

Fixed effects:
             Estimate Std. Error t value
(Intercept)  293.567     26.951  10.893
dose2          6.692     34.648   0.193
dose4        -39.975     27.099  -1.475
dose8       -105.517     26.558  -3.973

Correlation of Fixed Effects:
       (Intr) dose2  dose4
dose2 -0.380
dose4 -0.103  0.786
dose8 -0.724  0.443  0.360

Thanks in advance and here is my sessionInfo().

 > sessionInfo()
R version 2.10.1 Patched (2009-12-29 r50852)
i686-pc-linux-gnu

locale:
  [1] LC_CTYPE=en_US       LC_NUMERIC=C         LC_TIME=en_US
  [4] LC_COLLATE=C         LC_MONETARY=C        LC_MESSAGES=en_US
  [7] LC_PAPER=en_US       LC_NAME=C            LC_ADDRESS=C
[10] LC_TELEPHONE=C       LC_MEASUREMENT=en_US LC_IDENTIFICATION=C

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base

other attached packages:
[1] lme4_0.999375-32   Matrix_0.999375-33 lattice_0.17-26

loaded via a namespace (and not attached):
[1] grid_2.10.1
Ben Bolker
#
Kevin E. Thorpe wrote:
It's not entirely clear whether Subjects are unique within Treatments
(OK, maybe we should say "nested") or whether the same Subject can have
both Treatments  -- what would xtabs(~Subject+Dose+Treatment) look like?

  *If* each Subject appears in only one Treatment (i.e., 5 doses per
Subject, each Subject is either Oat or Barley, 12 Subjects per
Treatment), which is what I initially took for granted, then it makes
perfect sense to me that you can't specify Treatment|Subject, because
Treatment does not vary within Subject -- there doesn't seem to be any
way you can recover information about how the effect of Treatment varies
among Subject.

  Since your plots do look reasonably linear within Subject,

iAUC~Treatment*Dose+(Dose|Subject)

 seems best to me.
For what it's worth, this example closely parallels the Orthodont
example in the nlme package ...
Kevin E. Thorpe
#
Ben Bolker wrote:
It is a true crossover.  There are 12 patients and each patient gets 
both treatments and all four doses.

I will look at the Orthodont example too.
Ken Knoblauch
# 
Kevin E. Thorpe <kevin.thorpe at ...> writes:
Shouldn't you make Subject into a factor?

Ken
Ben Bolker
#
Ken Knoblauch wrote:
It would make the plot a little bit prettier but I don't think it
matters in this case because variable that appears as a grouping
variable (i.e. on the right of the | ) is automatically treated as a
factor?  I think?

  Since it is really a crossover trial, it would seem reasonable in
principle to have the (Treatment|Subject) random effect in there as
well. I'm not sure what to do about the -1 correlation: it seems the
choices (not necessarily in order) are (1) throw up your hands and say
there's not enough data to estimate independently; (2) try WinBUGS,
possibly with slightly informative priors; (3) try using lme4a to create
profiles of the parameters and see if you can figure out what's happening.
Kevin E. Thorpe
#
Ken Knoblauch wrote:
I just tried that and it had no effect on the results.
Viechtbauer Wolfgang (STAT)
# 
Maybe I am totally off here, but wouldn't it help if you make what is currently Dose = 0 equal to Dose = -8 and then have what is currently Dose = 8 be equal to Dose = 0? This should help to decrease the correlation between the intercepts and the slopes.

Best,

--
Wolfgang Viechtbauer                        http://www.wvbauer.com/
Department of Methodology and Statistics    Tel: +31 (43) 388-2277
School for Public Health and Primary Care   Office Location:
Maastricht University, P.O. Box 616         Room B2.01 (second floor)
6200 MD Maastricht, The Netherlands         Debyeplein 1 (Randwyck)


----Original Message----
From: r-sig-mixed-models-bounces at r-project.org
[mailto:r-sig-mixed-models-bounces at r-project.org] On Behalf Of Kevin E.
Thorpe Sent: Friday, April 09, 2010 13:04 To:
r-sig-mixed-models at r-project.org Subject: [R-sig-ME] Another case of
-1.0 correlation of random effects
Jarrod Hadfield
# 
Hi,

As a solution to suggestion 2: an inverse-Wishart prior with null  
scale matrix and degree of belief parameter equal to 3 (for a 2X2  
case) is non-informative for the correlation.  As a quick stab you  
could try:

prior<-list(R=list(V=1, nu=0), G=list(G1=list(V=diag(2)*1e-6, nu=3)))

m1<-MCMCglmm(iAUC~Dose, random=~us(Dose):Subject, prior=prior,  
data=gluc)

plot(posterior.cor(m1$VCV[,1:4]))


NOTE: I do not advocate this prior as it does not imply non- 
informativeness for the (co)variances (amongst other things). It is  
just useful to know the uncertainty regarding the correlation.

Jarrod
On 9 Apr 2010, at 14:03, Ben Bolker wrote:

            

      
      
      
    

        

      
      
    

        

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Ken Knoblauch
    

    

      
      #
    

  


  

      
Yes. It was a knee-jerk reaction, recalling a trap that
I fell into before with glm.  But, lmer seems to handle
this smoothly for the grouping factor.  Thanks.

Ken


Quoting Ben Bolker <bolker at ufl.edu>:

            

      
      
      
    

        

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Kevin E. Thorpe
    

    

      
      #
    

  


  

      
Ben Bolker wrote:

            

      
      
      
    

        
Let's see.  I wish (1) was an option.  (2) would be promising if my 
knowledge of BUGS and Bayesian methods filled more than a thimble. 
Thanks to Jarrod for his suggestion in response to this.  I'll take a 
look at that too.  Option (3) is probably worth a go too.

Aside from the fact that the Dose variable are the actual doses and not 
categories, and we all know not to categorize continuous variables, what 
are your thoughts on treating Dose as a factor (since it seems to behave)?

Thanks all for taking the time to provide your suggestions.

Kevin

  
    

      
      
    

  


  


      

    

    
2 days later

    

      
      

        


  

        

      Kevin E. Thorpe
    

    

      
      #
    

  


  

      
Kevin E. Thorpe wrote:

            

      
      
      
    

        
Regarding lme4a: how do I obtain it?  I guess that comes down to, what 
is the repository to give to install.packages()?  Does it require a 
different Matrix package than the one I have, which is, 0.999375-33 and 
if so, how do I not break my current lme4/Matrix combination?

By the way, the problems with these data get stranger.  In a different 
outcome from the same trial the following results from a model fitting 
attempt.

lmer(iAUC~Treatment+Dose+(Treatment|Subject)+(Dose|Subject),data=insulin)
Linear mixed model fit by REML
Formula: iAUC ~ Treatment + Dose + (Treatment | Subject) + (Dose | Subject)
    Data: insulin
   AIC  BIC logLik deviance REMLdev
  1956 1982   -968     1983    1936
Random effects:
  Groups   Name         Variance   Std.Dev.   Corr
  Subject  (Intercept)  4.2678e-02 2.0659e-01
           TreatmentOat 1.6115e+07 4.0144e+03 0.000
  Subject  (Intercept)  3.0430e+08 1.7444e+04
           Dose         1.5173e+06 1.2318e+03 -1.000
  Residual              3.1907e+07 5.6486e+03
Number of obs: 96, groups: Subject, 12

Fixed effects:
              Estimate Std. Error t value
(Intercept)   40142.4     5146.7   7.800
TreatmentOat   1340.3     1634.7   0.820
Dose          -2675.1      405.5  -6.597

Correlation of Fixed Effects:
             (Intr) TrtmnO
TreatmentOt -0.079
Dose        -0.922  0.000

As you can see, I get a 0 correlation within one set of random effects 
and -1.0 in the other.  Also, the fact the fixed effects estimates are 
huge makes me suspicious.

Note that if I drop the treatment portions and fit the Dose model to 
only one treatment, the correlation is again -1.0 and the fixed effects 
are similar.

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Douglas Bates
    

    

      
      #
    

  


  

      
On Mon, Apr 12, 2010 at 8:00 AM, Kevin E. Thorpe

        
<kevin.thorpe at utoronto.ca> wrote:

            

      
      
      
    

        
The repository is http://R-forge.R-project.org but be aware that lme4a
is under active development and not guaranteed against breakage.  It
would be inadvisable to rely on functions and classes in that package
to persist.

            

      
      
      
    

        
And did you notice that you have an Intercept term by Subject in there
twice?  It is not surprising that the parameter estimates are
unstable.  You will need to rethink the model.

            

      
      
      
    

        

      
      
    

  
  


  


      

    

    

      
      

        


  

        

      Kevin E. Thorpe
    

    

      
      #
    

  


  

      
Douglas Bates wrote:

            

      
      
      
    

        
Thanks Doug.  Any suggestions?  Note that the instability is present 
without the random effects for treatment in the model too.

            

      
      
      
    

        

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Viechtbauer Wolfgang (STAT)
    

    

      
      #
    

  


  

      
Dear Kevin,

I am curious to see what happens if you fit the same model with lme and if you play around with the optimizer used. Also, try changing the definition of the intercept. In particular:

library(nlme)
res <- lme(iAUC~Treatment+Dose, random=~Dose|Subject, data=insulin)
summary(res)

res <- lme(iAUC~Treatment+Dose, random=~Dose|Subject, data=insulin,
           control=list(msVerbose=T, opt="nlm"))
summary(res)

insulin$Dose <- insulin$Dose - 8

res <- lme(iAUC~Treatment+Dose, random=~Dose|Subject, data=insulin)
summary(res)

res <- lme(iAUC~Treatment+Dose, random=~Dose|Subject, data=insulin,
           control=list(msVerbose=T, opt="nlm"))
summary(res)

Best,

--
Wolfgang Viechtbauer                        http://www.wvbauer.com/
Department of Methodology and Statistics    Tel: +31 (43) 388-2277
School for Public Health and Primary Care   Office Location:
Maastricht University, P.O. Box 616         Room B2.01 (second floor)
6200 MD Maastricht, The Netherlands         Debyeplein 1 (Randwyck)


----Original Message----
From: r-sig-mixed-models-bounces at r-project.org
[mailto:r-sig-mixed-models-bounces at r-project.org] On Behalf Of Kevin E.
Thorpe Sent: Monday, April 12, 2010 15:23 To: Douglas Bates
Cc: r-sig-mixed-models at r-project.org
Subject: Re: [R-sig-ME] Another case of -1.0 correlation of random
effects

            

      
      
      
    

        

      
      
    

  
  


  


      

    

    

      
      

        


  

        

      Douglas Bates
    

    

      
      #
    

  


  

      
lme won't return a correlation of exactly -1 but that is a deficiency
of lme, not lmer.

Mixed-effects software should be able to handle the case of a singular
variance-covariance matrix for the random effects.  Because most
numerical methods, including those in lme, iterate with respect to the
precision matrix (inverse of the variance-covariance matrix) they
usually declare convergence at a non-singular matrix.  That is not an
advantage.  If the MLEs or REML estimates are at a singular covariance
the user should learn this.

On Mon, Apr 12, 2010 at 8:51 AM, Viechtbauer Wolfgang (STAT)

        
<Wolfgang.Viechtbauer at stat.unimaas.nl> wrote:

            

      
      
      
    

        

      
      
    

            

      
      
      
    

        

      
      
    

  
  


  


      

    

    

      
      

        


  

        

      Kevin E. Thorpe
    

    

      
      #
    

  


  

      
Douglas Bates wrote:

            

      
      
      
    

        
I tried to get lme4a but am having problems.

 > install.packages("lme4a", repos="http://R-Forge.R-project.org")
Warning message:
In getDependencies(pkgs, dependencies, available, lib) :
   package 'lme4a' is not available

Also, if I try and download the package source directly from:

http://r-forge.r-project.org/src/contrib/lme4a_0.999375-48.tar.gz

which comes from the "R Packages" tab of the lme4 project pages I get 
"PAGE NOT FOUND."

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Ben Bolker
    

    

      
      #
    

  


  

      
Something (don't know what) is currently broken in the lme4a build
process on r-forge, which doesn't necessarily mean that the package
itself is broken (it passes tests and installs fine on my system).  You
have two choices:

 (1) follow the instructions on r-forge to download the full SVN
repository; build the package from source. (You will need SVN tools as
well as the full suite of R-building tools -- compiler etc etc.)

 (2) I can send/post a binary build of the current SVN version easily
(which will install if you have the necessary tools to build from source
-- this gets you out of needing the SVN tools above, but you still need
compiler etc). I have shipped a version off to the R-project windows
cross-building service, so I could also send/post a Windows binary if
you like. Can you send sessionInfo() output?

        
Kevin E. Thorpe wrote:

            

      
      
      
    

        

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Douglas Bates
    

    

      
      #
    

  


  

      
On Tue, Apr 13, 2010 at 7:57 AM, Kevin E. Thorpe

        
<kevin.thorpe at utoronto.ca> wrote:

            

      
      
      
    

        
lme4a is under active development and I will often commit things from
work then update at home, add some more code, and commit again.   My
standard for commits is that I can run

R CMD INSTALL lme4a

on an Ubuntu 64-bit system, which I run both at home and at work.  I'm
not surprised that it doesn't compile and pass checks everywhere.
Things will need to stabilize with a week because I have project
reports to read and grade starting next Tuesday and I won't have time
for further development.

  
  


  


      

    

    

      
      

        


  

        

      Kevin E. Thorpe
    

    

      
      #
    

  


  

      
Ben Bolker wrote:

            

      
      
      
    

        
I posted my sessionInfo() in my OP, but I'll give it below again.  In 
any case, I run on Linux and always build R from source, so I have the 
compiler tools.  I also have subversion.

 > sessionInfo()
R version 2.10.1 Patched (2009-12-29 r50852)
i686-pc-linux-gnu

locale:
  [1] LC_CTYPE=en_US       LC_NUMERIC=C         LC_TIME=en_US
  [4] LC_COLLATE=C         LC_MONETARY=C        LC_MESSAGES=en_US
  [7] LC_PAPER=en_US       LC_NAME=C            LC_ADDRESS=C
[10] LC_TELEPHONE=C       LC_MEASUREMENT=en_US LC_IDENTIFICATION=C

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base

other attached packages:
[1] lme4_0.999375-32   Matrix_0.999375-33 lattice_0.17-26

loaded via a namespace (and not attached):
[1] grid_2.10.1  tools_2.10.1

            

      
      
      
    

        

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Ben Bolker
    

    

      
      #
    

  


  

      
In that case, just go to http://r-forge.r-project.org/scm/?group_id=60
and follow the instructions:

svn checkout svn://svn.r-forge.r-project.org/svnroot/lme4

then navigate to the lme4a directory and build/install it ...

        
Kevin E. Thorpe wrote:

            

      
      
      
    

        

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Kevin E. Thorpe
    

    

      
      #
    

  


  

      
Ben Bolker wrote:

            

      
      
      
    

        
It's not my day.  I go to the directory where the lme4a tree is located 
and type "R CMD build lme4a" and get an error when the process gets to:

g++ -I/usr/local/lib/R/include  -I/usr/local/include 
-I"/usr/local/lib/R/library/Matrix/include" 
-I"/usr/local/lib/R/library/Rcpp/include"   -fpic  -g -O2 -c 
glmFamily.cpp -o glmFamily.o

In file included from glmFamily.cpp:1:
glmFamily.h:5:18: error: Rcpp.h: No such file or directory
In file included from glmFamily.cpp:1:
[errors continue]

I looked at lme4a/src/glmFamily.h and find

#include <Rcpp.h>
#include <Rmath.h>

Now, I have both of these files in my installation, but not in the 
locations referenced with the -I switches.  My Rcpp.h is found in 
/usr/local/lib/R/library/Rcpp/lib (not include) and my Rmath.h is
found in /usr/local/lib/R/include which doesn't even show up in an -I 
switch.

So, I don't know how these -I switches get generated, but there is a 
problem somewhere, either with my installation or something else.

My installation was a standard build (configure; make; make install).

Yes, I fully understand what "unstable" means.  Hopefully these errors 
are helpful.

Kevin

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Douglas Bates
    

    

      
      #
    

  


  

      
On Tue, Apr 13, 2010 at 10:19 AM, Kevin E. Thorpe

        
<kevin.thorpe at utoronto.ca> wrote:

            

      
      
      
    

        
The DESCRIPTION file says that the lme4a package is linking to Matrix
and Rcpp so you need to have both of those packages installed before
you can compile lme4a.

  
  


  


      

    

    

      
      

        


  

        

      Kevin E. Thorpe
    

    

      
      #
    

  


  

      
Douglas Bates wrote:

            

      
      
      
    

        
They are both installed.  From the DESCRIPTION files on my system:

Package: Rcpp
Title: Rcpp R/C++ interface package
Version: 0.7.11
Date: $Date: 2010-03-26 15:15:51 -0500 (Fri, 26 Mar 2010) $

and

Package: Matrix
Version: 0.999375-38
Date: 2010-03-31

  
    

      
      
    

  


  


      

    

    
2 days later

    

      
      

        


  

        

      Kevin E. Thorpe
    

    

      
      #
    

  


  

      
Kevin E. Thorpe wrote:

            

      
      
      
    

        
Okay, I now have lme4a installed and I get an error message when I do 
(note: this is the same model from my OP):

 > gluc.lmer1a <- 
lmer(iAUC~Dose+(Dose|Subject),data=gluc,subset=Treatment=="Oat",REML=FALSE)

 > gluc.lmer1a
Linear mixed model fit by maximum likelihood  ['lmer']
Formula: iAUC ~ Dose + (Dose | Subject)
    Data: gluc
  Subset: Treatment == "Oat"
      AIC      BIC   logLik deviance
    575.1    586.3   -281.6    563.1

Random effects:
  Groups   Name        Variance Std.Dev. Corr
  Subject  (Intercept) 7492.19  86.557
           Dose          14.68   3.831   -1.000
  Residual             4727.27  68.755
Number of obs: 48, groups: Subject, 12

Fixed effects:
             Estimate Std. Error t value
(Intercept)  309.352     29.338  10.544
Dose         -14.424      3.533  -4.083

Correlation of Fixed Effects:
      (Intr)
Dose -0.647

 > pr1 <- profile(gluc.lmer1a at env)  ## using @env base on other threads
Error in x[ndat + (1L:deg) - deg] :
   only 0's may be mixed with negative subscripts

Is this because I'm trying to profile a model that profile() cannot 
handle yet, or does it indicate there really are serious problems with 
my model?

I'm at a loss as to how determine what is really going on with these data.

Kevin