[R-sig-ME] Strange mcmcsamp issue

[Douglas Bates]

On Thu, Feb 21, 2008 at 6:43 PM, Jarrett Byrnes <jebyrnes at ucdavis.edu> wrote:
> I'm attempting to pull out the simple effects from a mixed model with
>  two crossed treatments.  The model structure is such that
>
>  a.lmer<-lmer(response ~ trta*trtb+(1|pot))
>
>  In the experiment, I have an array of pots.  Each pot has a type of
>  treatment A applied to it.  Within the pot, there are two types of
>  treatment B that are applied, one on either side.  I am using a mixed
>  model as I wanted to account for non-independence within a pot.
>
>  There is an interaction between a and b, but I want to look at the 95%
>  credible intervals of the simple effects to see which treatment
>  combinations overlap 0, are greater than 0, or are less than 0.  While
>  mcmcsamp works great on this object I am unclear on how to then
>  combine parameter values and error to get this interval.
>
>  So, I attempted a model such as the following
>
>  a.lmer<-lmer(response ~ trta*trtb+(1|pot))

I would suggest trying to fit that model in the development version of
the lme4 package, available via

install.packages("lme4", repos = "http://r-forge.r-project.org")

requesting verbose output.  That is, fit it as

a.lmer <- lmer(response ~ trta * trtb + (1|pot), verbose = TRUE)

For each iteration of the optimizer the iteration number, the value of
the deviance and the current estimate of the relative standard
deviation of the random effects (i.e. the ratio of the standard
deviation of the random effects and the standard deviation of the
residuals) are printed.  Check if the relative standard deviation is
being driven to zero.

>  Which yielded the following error:
>   Leading minor of order 15 in downdated X'X is not positive definite


>  Thinking that this might be an intercept issue, I fit the following
>  model:
>
>  a.lmer<-lmer(response ~ trta*trtb+0+(1|pot))
>
>  This fit just fine.  summary() showed me a table of parameter values
>  that seemed about what I would expect, although the correlation of
>  fixed effects matrix was populated largely by 0's However,
>
>  a.mcmcsamp<-(a.lmer, 1000)
>  yielded the following error
>
>  Error: Omega[[1]] is not positive definite
>  Error in t(.Call(mer_MCMCsamp, object, saveb, n, trans, verbose,
>  deviance)) :
>    error in evaluating the argument 'x' in selecting a method for
>  function 't'
>
>  However, if I try for roughly 30 or fewer replicates, everything works
>  just fine.
>
>  Even more strange, when I next looked at a.lmer using summary() all of
>  the error values for parameters had become 0, and the matrix for
>  correlation of fixed effects was filled with NaNs.  This strikes me as
>  rather odd.
>
>  1) Perhaps this has been fixed in later releases - I'm working off of
>  lme4 version 0.99875-9 on R 2.6.2.  Should I try these instead?
>
>  2) Or, am I going about this attempt to get simple effect estimates
>  all wrong?  Is it possible to use the output from the first model with
>  mcmcsamp to get estimates of the simple effects?
>
>  Thanks for any advice you might have!
>
>  _______________________________________________
>  R-sig-mixed-models at r-project.org mailing list
>  https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>