[R-sig-ME] 2 correlated random effects with quadrature?

[Joshua Wiley]

Hi Ross,

I do not know any functionality for GHQ, but you already acknowledged
you were okay with slow, at which point, you can probably do better
with MCMC.  If you're not feeling like rolling your own, the MCMCglmm
package makes it quite easy.

Here is a (not particularly sensible) example:

require(MCMCglmm)
set.seed(1234)
dat <- mtcars[sample(1:32, 1000, replace = TRUE), ]
dat <- within(dat, {
  qsec <- scale(qsec)
  hp <- scale(hp)
  mpg <- scale(mpg)
  disp <- scale(disp)
})
dat$ID <- factor(rep(letters, length.out = 1000))
dat$cyl <- factor(dat$cyl)
# set seed and estimate model
set.seed(10)
m <- MCMCglmm(vs ~ qsec + mpg + drat, random = ~ us(1 + mpg):ID,
family = "categorical",
  data = dat, prior = list(
  B = list(mu = c(0, 0, 0, 0), V = diag(4) * 1e2),
  R = list(V = 1, fix = 1),
  G = list(G1 = list(V = diag(2), nu = 2))), pr=TRUE,
  nitt = 55000, thin = 20, burnin = 5000, verbose=FALSE)

# print a summary
summary(m)


 Iterations = 5001:54981
 Thinning interval  = 20
 Sample size  = 2500
 DIC: 22.87558
 G-structure:  ~us(1 + mpg):ID
                           post.mean l-95% CI u-95% CI eff.samp
(Intercept):(Intercept).ID   0.75374   0.1228    1.754   217.21
mpg:(Intercept).ID          -0.03048  -1.1675    0.922   108.96
(Intercept):mpg.ID          -0.03048  -1.1675    0.922   108.96
mpg:mpg.ID                   0.97473   0.1383    2.644    41.16
 R-structure:  ~units
      post.mean l-95% CI u-95% CI eff.samp
units         1        1        1        0
 Location effects: vs ~ qsec + mpg + drat
            post.mean l-95% CI u-95% CI eff.samp  pMCMC
(Intercept)    2.4890  -8.4259  11.9981   40.374  0.652
qsec          18.9466  13.3597  23.8444    2.672 <4e-04 ***
mpg            8.2778   5.6246  11.2228    9.550 <4e-04 ***
drat          -0.3849  -2.9258   2.3827   40.419  0.800
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1


Now of course I did not run enough iterations, and probably have a low
thinning interval, you may want better priors, and if you are looking
to compare to standard logistic models, don't forget to rescale the
estimates (note the non zero variance used to improve mixing).

Still, it is relatively straightforward, highly stable given enough
time, flexible to have many effects, and gives you posteriors which
are great for inference.

Cheers,

Josh



On Mon, Mar 11, 2013 at 1:25 PM, Ross Boylan <ross at biostat.ucsf.edu> wrote:
> Is there a way to fit generalized linear mixed model with 2 correlated
> random effects in R, using quadrature?  At the moment, I'm only
> concerned with binary outcomes.
>
> When I try glmer from lme4 with the quadrature argument I get
> Error: AGQ only defined for a single scalar random-effects term
>
> Yes, I know 2 dimensional quadrature is slow.
>
> Ross Boylaln
>
> _______________________________________________
> R-sig-mixed-models at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models



-- 
Joshua Wiley
Ph.D. Student, Health Psychology
University of California, Los Angeles
http://joshuawiley.com/
Senior Analyst - Elkhart Group Ltd.
http://elkhartgroup.com