Thanks for the help.  Just to confirm a modelling issue, what I'm trying to
recreate is the ML version of a hierarchical Bayesian analysis, in which
each class exhibits a main effect, and each observed covariate (all members
of some underlying class of unobserved effects) are random effects.

Here's what the hierarchical Bayesian model looks like (in BUGS model
language):

model {
  for (i in 1:N) {
    Y[i] ~ dnorm(mu[i], tau)
    mu[i] ~ intercept + inprod(covBeta[], covariates[i,])
  }

  intercept ~ dnorm(0, 1.0E-03)

  for (k in 1:numCov) {
    covBeta[k] ~ dnorm(classBeta[cov2class[k]],
                        classTau[cov2class[k]])
  }

  for (j in 1:numClass) {
    classBeta[j] ~ dnorm(0, 1.0E-03)
    classTau[j] <- pow(classSigma[j], -2)
    classSigma[j] ~ dunif(0, 1000)
  }

  tau <- pow(sigma, -2)
  sigma <- dunif(0, 1000)
}

Is this lmer formula equivalent in nature? (disregarding the impact of the
hyperpriors):

Y ~ 1 + cov1 + (0+cov1 | class1) + cov2 + (0+cov2 | class1) + ... + cov100 +
(0+cov100 | class8)

Or should this be written as:

Y ~ 1 + class1 + (0+cov1 | class1) + (0+cov2|class1) + ... class8 + ... +
(0+cov100 | class8)

Thanks again,

-Aaron


On Tue, Aug 19, 2008 at 8:04 AM, Douglas Bates <bates@stat.wisc.edu> wrote:

> On Tue, Aug 19, 2008 at 1:59 AM, Martin Maechler
> <maechler@stat.math.ethz.ch> wrote:
> >>>>>> "AM" == Aaron Mackey <ajmackey@gmail.com>
> >>>>>>     on Mon, 18 Aug 2008 15:41:20 -0400 writes:
> >
> >    AM> Dear Dr. Bates, Thanks for the info about the SIG on
> >    AM> mixed models.
> >
> >    AM> A followup question: if I did not want a generalized
> >    AM> variance-covariance matrix, and instead wanted to assume
> >    AM> independence between grouped covariates, then instead of
> >    AM> (3n-n^2)/2 variance-covariance parameters per grouping
> >    AM> (of size n), I'd have only n variance parameters, yes?
> >    AM> Would that be specified (x1|fac) + (x2|fac) instead of
> >    AM> (x1 + x2 | fac)?
> >
> > yes.
>
> For lmer I think you would want to write that as (0+x1|fac) + (0+x2|fac) +
> ...
>
> By default a linear model formula includes the intercept term so that
> (x1|fac) is equivalent to (1+x1|fac).  If you do not suppress the
> implicit intercept you end up with multiple terms including it.
>
> >    AM> And either way I formulate it, my syntax question
> >    AM> remains: is there a way to write the formula that is
> >    AM> generic for a given covariate matrix of unspecified
> >    AM> dimension?  Said another way, how would you write a lmer
> >    AM> "wrapper" function that was given a matrix of
> >    AM> covariates?  How would you "unpack" the columns of that
> >    AM> matrix to generate the formula required for lmer()?
> >    AM> Would you build up a formula string programmatically,
> >    AM> and then "eval" it?
> >
> > yes / "kind of" :
> >
> > I would build the formula programmatically
> > but then use  as.formula(.)  or possibly in some cases  parse().
> >
> > Package 'sfsmisc' has a utility  wrapFormula() that I had need
> > for at some time:
> >
> > ----------------------------------------
> >
> > install.packages("sfsmisc")
> > library(sfsmisc)
> >
> > myF <- wrapFormula(Fertility ~ . , data = swiss)
> > myF
> > ## Fertility ~ s(Agriculture) + s(Examination) + s(Education) +
> > ##             s(Catholic) + s(Infant.Mortality)
> >
> > ------------------------------------------
> >
> > and you can look at the short definition of  wrapFormula
> > to see how I use gsub(.) on character strings to construct the
> > new formula.
> >
> > Martin Maechler, ETH Zurich
> >
> >
> >    AM> Thanks again,
> >
> >    AM> -Aaron
> >
> >    AM> On Mon, Aug 18, 2008 at 3:21 PM, Douglas Bates
> >    AM> <bates@stat.wisc.edu> wrote:
> >
> >    >> On Mon, Aug 18, 2008 at 11:20 AM, Aaron Mackey
> >    >> <ajmackey@gmail.com> wrote: > I have a fairly large
> >    >> model:
> >    >>
> >    >> >> length(Y) > [1] 3051 >> dim(covariates) > [1] 3051 211
> >    >>
> >    >> > All of these 211 covariates need to be nested
> >    >> hierarchically within a > grouping "class", of which
> >    >> there are 8.  I have an accessory vector, " > cov2class"
> >    >> that specifies the mapping between covariates and the 8
> >    >> classes.
> >    >>
> >    >> > Now, I understand I can break all this information up
> >    >> into individual > vectors (cov1, cov2, ..., cov211,
> >    >> class1, class2, ..., class8), and do > something like
> >    >> this:
> >    >>
> >    >> > model <- lmer(Y ~ 1 + cov1 + cov2 + ... + cov211 + >
> >    >> (cov1 + cov2 + ... | class1) + > (...) + > (... + cov210
> >    >> + cov211 | class8)
> >    >>
> >    >> > But I'd like keep things syntactically simpler, and use
> >    >> the covariates > and cov2class > variables directly.  I
> >    >> haven't been able to find the right syntactic sugar > to
> >    >> get this done.
> >    >>
> >    >> Before considering the syntax you should consider the
> >    >> complexity of the model.  Even with 3000+ observations
> >    >> estimating fixed effects for over 200 covariates is a
> >    >> risky business.  To follow that up by incorporating
> >    >> dozens of covariates in a vector-valued random effects
> >    >> term with a general variance covariance matrix is not
> >    >> likely to be productive.  When estimating variances and
> >    >> covariances of random effects the complexity of the
> >    >> estimation problem increases according to the square of
> >    >> the dimension of each random-effects vector.  For a term
> >    >> of the form (1|fac) there is one variance-covariance
> >    >> parameter to estimate.  For a term of the form (1+x|fac)
> >    >> there are 3 such parameters to estimate, (1+x1+x2|fac)
> >    >> requires 6, etc.  Assuming that your 211 covariates are
> >    >> more-or-less equally distributed among the 8 classes you
> >    >> would have about 27 per class which means that you have
> >    >> 378 variance-covariance parameters to estimate for each
> >    >> of your 8 classes.  Even with 3000+ observations it would
> >    >> be optimistic to expect to estimate that many
> >    >> variance-covariance parameters.
> >    >>
> >    >> You may want to move your question to the
> >    >> R-SIG-Mixed-Models mailing list and follow up on the
> >    >> model specification there.
> >    >>
> >
> >        [[alternative HTML version deleted]]
> >
> >    AM> _______________________________________________
> >    AM> R-sig-mixed-models@r-project.org mailing list
> >    AM> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
> >
> > _______________________________________________
> > R-sig-mixed-models@r-project.org mailing list
> > https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
> >
>

	[[alternative HTML version deleted]]


