R/MCMCgrm.R

In gap: Genetic Analysis Package

#' Mixed modeling with genetic relationship matrices
#'
#' @param model statistical model.
#' @param prior a list of priors for parameters in the model above.
#' @param data a data.frame containing outcome and covariates.
#' @param GRM a relationship matrix.
#' @param eps a small number added to the diagonal of the a nonpositive definite GRM.
#' @param n.thin thinning parameter in the MCMC.
#' @param n.burnin the number of burn-in's.
#' @param n.iter the number of iterations.
#' @param ... other options as appropriate for MCMCglmm.
#'
#' @details
#' Mixed modeling with genomic relationship matrix. This is appropriate with relationship
#' matrix derived from family structures or unrelated individuals based on whole genome data.
#'
#' The function was created to address a number of issues involving mixed modelling with
#' family data or population sample with whole genome data. First, the implementaiton
#' will shed light on the uncertainty involved with polygenic effect in that posterior
#' distributions can be obtained. Second, while the model can be used with the MCMCglmm
#' package there is often issues with the specification of pedigree structures but this
#' is less of a problem with genetic relationship matrices. We can use established
#' algorithms to generate kinship or genomic relationship matrix as input to the MCMCglmm
#' function. Third, it is more intuitive to specify function arguments in line with other
#' packages such as R2OpenBUGS, R2jags or glmmBUGS. In addition, our experiences of tuning
#' the model would help to reset the input and default values.
#'
#' @export
#' @return
#' The returned value is an object as generated by MCMCglmm.
#'
#' @references
#' \insertRef{hadfield10}{gap}
#'
#' @examples
#' \dontrun{
#' ### with kinship
#' 
#' # library(kinship) 
#' # fam <- with(l51,makefamid(id,fid,mid))
#' # s <-with(l51, makekinship(fam, id, fid, mid))
#' # K <- as.matrix(s)*2   
#'
#' ### with gap
#'
#' s <- kin.morgan(l51)
#' K <- with(s,kin.matrix*2)
#' prior <- list(R=list(V=1, nu=0.002), G=list(G1=list(V=1, nu=0.002)))
#' m <- MCMCgrm(qt~1,prior,l51,K)
#' save(m,file="l51.m")
#' pdf("l51.pdf")
#' plot(m)
#' dev.off()
#'
#' # A real analysis on bats
#' ## data
#' bianfu.GRM <- read.table("bianfu.GRM.txt", header = TRUE)
#' bianfu.GRM[1:5,1:6]
#' Data <- read.table(file = "PHONE.txt", header = TRUE, 
#'                    colClasses=c(rep("factor",3),rep("numeric",7)))
#' ## MCMCgrm
#' library("MCMCglmm")
#' GRM <- as.matrix(bianfu.GRM[,-1])
#' colnames(GRM) <- rownames(GRM) <- bianfu.GRM[,1]
#' library(gap)
#' names(Data)[1] <- "id"
#' prior <- list(G = list(G1 = list(V = 1, nu = 0.002)), R = list(V = 1, nu = 0.002))
#' model1.1 <- MCMCgrm(WEIGTHT ~ 1, prior, Data, GRM, n.burnin=100, n.iter=1000, verbose=FALSE)
#' ## an alternative
#' names(Data)[1] <- "animal"
#' N <- nrow(Data)
#' i <- rep(1:N, rep(N, N))
#' j <- rep(1:N, N)
#' s <- Matrix::spMatrix(N, N, i, j, as.vector(GRM))
#' Ginv <- Matrix::solve(s)
#' class(Ginv) <- "dgCMatrix"
#' rownames(Ginv) <- Ginv@Dimnames[[1]] <- with(Data, animal)
#' model1.2 <- MCMCglmm(WEIGTHT ~ 1, random= ~ animal, data = Data,
#'   ginverse=list(animal=Ginv), prior = prior, burnin=100, nitt=1000, verbose=FALSE)
#' ## without missing data
#' model1.3 <- MCMCglmm(Peak_Freq ~ WEIGTHT, random = ~ animal, 
#'   data=subset(Data,!is.na(Peak_Freq)&!is.na(WEIGTHT)), 
#'   ginverse=list(animal=Ginv), prior = prior, burnin=100, nitt=1000, verbose=FALSE)
#' }
#'
#' @author Jing Hua Zhao
#' @keywords htest

MCMCgrm <- function(model,prior,data,GRM,eps=0,n.thin=10,n.burnin=3000,n.iter=13000,...)
{
  for(p in c("Matrix", "MCMCglmm")) {
     if (length(grep(paste("^package:", p, "$", sep=""), search())) == 0) {
        if (!requireNamespace(p, quietly = TRUE))
        warning(paste("MCMCgrm needs package `", p, "' to be fully functional; please install", sep=""))
     }
  }
  N <- dim(data)[1]
  GRM <- GRM + diag(eps,N)
  i <- rep(1:N,rep(N,N))
  j <- rep(1:N,N)
  s <- Matrix::spMatrix(N,N,i,j,as.vector(GRM))
  Ginv<-Matrix::solve(s)
  class(Ginv) <- "dgCMatrix"
  rownames(Ginv) <- Ginv@Dimnames[[1]] <- with(data,id)
  m <- MCMCglmm::MCMCglmm(as.formula(model), random=~id, ginverse=list(id=Ginv), data=data,
                prior=prior, thin=n.thin, burnin=n.burnin, nitt=n.iter, ...)
  summary(m)
  return(m)
}

# 8-1-2016 MRC-Epid JHZ