Skip to content
Prev 9794 / 20628 Next

2 correlated random effects with quadrature?

On Wed, 13 Mar 2013, Ben Bolker wrote:

            
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