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

[David Duffy]

On Wed, 13 Mar 2013, Ben Bolker wrote:

> Ross Boylan <ross at ...> writes:
>
>>
>> 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
>
>  I don't know offhand of an R package that will do this.

If you're happy with a probit link, it is pretty straightforward to hook 
mvtnorm up to a maximizer, especially if you have just the one problem to 
fit.  Your groups can't be too large unfortunately.

If you look in 
http://genepi.qimr.edu.au/staff/davidD/Sib-pair/Src/sib-pair.R

under varcomp.mft(), there is an example of two*N crossed correlated 
random effects in a standard genetics variance components model:

lik.mft.fam <- function(vc.data, m, va, vq, nthresh) {
   ln <- function(x) ifelse(x>0, log(x), 0)
   ndata <- length(vc.data$y)
   thresh <- matrix(NA, nr=ndata, nc=(nthresh+2))
   thresh[,1] <- -Inf
   thresh[,ncol(thresh)] <- Inf
   thresh[, -c(1, ncol(thresh))] <- vc.data$covariates %*% matrix(m, nr=1)
   S <- va*vc.data$rel + vq*vc.data$ibd + diag(1-va-vq,nrow(vc.data$rel))
   res <- ln(pmvnorm(lower=thresh[cbind(1:ndata,vc.data$y+1)],
                     upper=thresh[cbind(1:ndata,vc.data$y+2)],
                     corr=S))

where va and vq are the variance components for two classes of random 
effects (additive genetic and QTL, one each for each individual 
observation), and rel and ibd are the correlation matrices for the random 
effects (which we can specify using the pedigree structure and 
observed sharing of genetic markers).  This likelihood is maximized 
using optim().  For high dimensional integration, you need to tweak abseps 
and releps of the Quasi-Monte-Carlo algorithm.  It can do up to 1000 
dimensions, but your likelihood can bounce around a fair bit (because it 
is MC), which upsets your maximizer, but you'll get the right answer 
eventually.

Cheers, David Duffy.




| David Duffy (MBBS PhD)                                         ,-_|\
| email: davidD at qimr.edu.au  ph: INT+61+7+3362-0217 fax: -0101  /     *
| Epidemiology Unit, Queensland Institute of Medical Research   \_,-._/
| 300 Herston Rd, Brisbane, Queensland 4029, Australia  GPG 4D0B994A v