Another case of -1.0 correlation of random effects

On Mon, Apr 12, 2010 at 8:00 AM, Kevin E. Thorpe
<kevin.thorpe at utoronto.ca> wrote:

Kevin E. Thorpe wrote:

Ben Bolker wrote:

Ken Knoblauch wrote:

Kevin E. Thorpe <kevin.thorpe at ...> writes:

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 ...

clip>

Shouldn't you make Subject into a factor?

Ken

 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.

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

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?

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.

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

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.

 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.