R/EM.MIMv.R

In QTLEMM: QTL EM Algorithm Mapping and Hotspots Detection

#' EM Algorithm for QTL MIM with Asymptotic Variance-Covariance Matrix
#'
#' Expectation-maximization algorithm for QTL multiple interval mapping.
#' It can obtain the asymptotic variance-covariance matrix of the result
#' from the EM algorithm and the approximate solution of variances of
#' parameters.
#'
#' @param QTL matrix. A q*2 matrix contains the QTL information, where the
#' row dimension 'q' represents the number of QTLs in the chromosomes. The
#' first column labels the chromosomes where the QTLs are located, and the
#' second column labels the positions of QTLs (in morgan (M) or centimorgan
#' (cM)).
#' @param marker matrix. A k*2 matrix contains the marker information,
#' where the row dimension 'k' represents the number of markers in the
#' chromosomes. The first column labels the chromosomes where the markers
#' are located, and the second column labels the positions of markers (in
#' morgan (M) or centimorgan (cM)). It's important to note that chromosomes
#' and positions must be sorted in order.
#' @param geno matrix. A n*k matrix contains the genotypes of k markers
#' for n individuals. The marker genotypes of P1 homozygote (MM),
#' heterozygote (Mm), and P2 homozygote (mm) are coded as 2, 1, and 0,
#' respectively, with NA indicating missing values.
#' @param D.matrix matrix. The design matrix of QTL effects is a g*p
#' matrix, where g is the number of possible QTL genotypes, and p is the
#' number of effects considered in the MIM model. The design matrix can
#' be easily generated by the function D.make().
#' @param cp.matrix matrix. The conditional probability matrix is an
#' n*g matrix, where n is the number of genotyped individuals, and g is
#' the number of possible genotypes of QTLs. If cp.matrix=NULL, the
#' function will calculate the conditional probability matrix, and markers
#' whose positions are the same as QTLs will be skipped.
#' @param y vector. A vector with n elements that contain the phenotype
#' values of individuals.
#' @param type character. The population type of the dataset. Includes
#' backcross (type="BC"), advanced intercross population (type="AI"), and
#' recombinant inbred population (type="RI"). The default value is "RI".
#' @param ng integer. The generation number of the population type. For
#' instance, in a BC1 population where type="BC", ng=1; in an AI F3
#' population where type="AI", ng=3.
#' @param cM logical. Specify the unit of marker position. If cM=TRUE, it
#' denotes centimorgan; if cM=FALSE, it denotes morgan.
#' @param E.vector0 vector. The initial value for QTL effects. The
#' number of elements corresponds to the column dimension of the design
#' matrix. If E.vector0=NULL, the initial value for all effects will be
#' set to 0.
#' @param X matrix. The design matrix of the fixed factors except for
#' QTL effects. It is an n*k matrix, where n is the number of
#' individuals, and k is the number of fixed factors. If X=NULL,
#' the matrix will be an n*1 matrix where all elements are 1.
#' @param beta0 vector. The initial value for effects of the fixed
#' factors. The number of elements corresponds to the column dimension
#' of the fixed factor design matrix.  If beta0=NULL, the initial value
#' will be set to the average of y.
#' @param variance0 numeric. The initial value for variance. If
#' variance0=NULL, the initial value will be the variance of
#' phenotype values.
#' @param crit numeric. The convergence criterion of EM algorithm.
#' The E and M steps will iterate until a convergence criterion is met.
#' It must be a value between 0 and 1.
#' @param stop numeric. The stopping criterion of EM algorithm. The E and
#' M steps will halt when the iteration number reaches the stopping
#' criterion, treating the algorithm as having failed to converge.
#' @param conv logical. If set to False, it will disregard the failure to
#' converge and output the last result obtained during the EM algorithm
#' before reaching the stopping criterion.
#' @param var.pos logical. Determines whether the variance of the position
#' of QTLs will be considered in the asymptotic variance-covariance matrix
#' or not.
#' @param console logical. Determines whether the process of the algorithm
#' will be displayed in the R console or not.
#' @param IMresult list. The data list of the output from IM.search(). The
#' required parameters for this function will be extracted from the data list.
#' @param MIMresult list. The data list of the output from MIM.search() or
#' MIM.points(). The required parameters for this function will be extracted
#' from the data list.
#'
#' @return
#' \item{E.vector}{The QTL effects are calculated by the EM algorithm.}
#' \item{beta}{The effects of the fixed factors are calculated by the EM
#' algorithm.}
#' \item{variance}{The error variance is calculated by the EM algorithm.}
#' \item{PI.matrix}{The posterior probabilities matrix after the
#' process of the EM algorithm.}
#' \item{log.likelihood}{The log-likelihood value of this model.}
#' \item{LRT}{The LRT statistic of this model.}
#' \item{R2}{The coefficient of determination of this model. This
#' can be used as an estimate of heritability.}
#' \item{y.hat}{The fitted values of trait values are calculated by
#' the estimated values from the EM algorithm.}
#' \item{iteration.number}{The iteration number of the EM algorithm.}
#' \item{avc.matrix}{The asymptotic variance-covariance matrix contains
#' position of QTLs, QTL effects, variance, and fix effects.}
#' \item{EMvar}{The asymptotic approximate variances include the position
#' of QTLs, QTL effects, variance, and fixed effects. Parameters for which
#' the approximate variance cannot be calculated will be shown as 0. The
#' approximate variance of the position of QTLs is calculated in morgan.}
#'
#' @note
#'
#' When IMresult and MIMresult are entered simultaneously, only IMresult will
#' be processed.
#'
#' If an error occurs, please check whether the original output data of
#' IM.search(), MIM.result(), or MIM.points() is used.
#'
#' @export
#'
#' @references
#'
#' KAO, C.-H. and Z.-B. ZENG 1997 General formulas for obtaining the maximum
#' likelihood estimates and the asymptotic variance-covariance matrix in QTL
#' mapping when using the EM algorithm. Biometrics 53, 653-665. <doi: 10.2307/2533965.>
#'
#' KAO, C.-H., Z.-B. ZENG and R. D. TEASDALE 1999 Multiple interval mapping
#' for Quantitative Trait Loci. Genetics 152: 1203-1216. <doi: 10.1093/genetics/152.3.1203>
#'
#' @seealso
#' \code{\link[QTLEMM]{D.make}}
#' \code{\link[QTLEMM]{Q.make}}
#' \code{\link[QTLEMM]{EM.MIM}}
#' \code{\link[QTLEMM]{IM.search}}
#' \code{\link[QTLEMM]{MIM.search}}
#' \code{\link[QTLEMM]{MIM.points}}
#'
#' @examples
#' # load the example data
#' load(system.file("extdata", "exampledata.RDATA", package = "QTLEMM"))
#'
#' # run and result
#' D.matrix <- D.make(3, type = "RI", aa = c(1, 3, 2, 3), dd = c(1, 2, 1, 3), ad = c(1, 2, 2, 3))
#' result <- EM.MIMv(QTL, marker, geno, D.matrix, cp.matrix = NULL, y)
#' result$EMvar
EM.MIMv <- function(QTL = NULL, marker = NULL, geno = NULL, D.matrix = NULL, cp.matrix = NULL, y = NULL,
                    type = "RI", ng = 2, cM = TRUE, E.vector0 = NULL, X = NULL, beta0 = NULL,
                    variance0 = NULL, crit = 10^-5, stop = 1000, conv = TRUE, var.pos = TRUE,
                    console = interactive(), IMresult = NULL, MIMresult = NULL){

  if(!is.null(IMresult)){
    check1 <- names(IMresult) == c("effect", "LRT.threshold", "detect.QTL", "model", "inputdata")
    if(!is.list(IMresult) | sum(check1) != 5){
      stop("IMresult data error, please input all of the original output data of IM.search().", call. = FALSE)
    } else {
      check2 <- names(IMresult$inputdata) == c("marker", "geno", "y", "yu", "sele.g", "type", "ng", "cM", "d.eff" )
      if(!is.list(IMresult$inputdata) | sum(check2) != 9){
        stop("IMresult data error, please input all of the original output data of IM.search().", call. = FALSE)
      } else {
        if(length(IMresult$inputdata$sele.g) != 0){
          if(IMresult$inputdata$sele.g != "n"){
            stop("IMresult data error, this function does not support selective genotyping data.", call. = FALSE)
          }
        }
        if(is.data.frame(IMresult$detect.QTL) | is.matrix(IMresult$detect.QTL)){
          QTL <- as.matrix(IMresult$detect.QTL[,1:2])
        } else {stop("IMresult data error, please input all of the original output data of IM.search().", call. = FALSE)}
        marker <- IMresult$inputdata$marker
        geno <- IMresult$inputdata$geno
        y <- IMresult$inputdata$y
        type <- IMresult$inputdata$type
        ng <- IMresult$inputdata$ng
        cM <- IMresult$inputdata$cM
        if(is.null(D.matrix)){
          D.matrix <- D.make(nrow(QTL), type = type)
          if(length(IMresult$inputdata$d.eff) != 0){
            if(IMresult$inputdata$d.eff[1] == 0){
              D.matrix <- D.make(nrow(QTL), type = type, d = 0)
            }
          }
        }
      }
    }
  } else if(!is.null(MIMresult)){
    check1 <- names(MIMresult) == c("effect", "QTL.best", "effect.best", "model", "inputdata")
    if(!is.list(MIMresult) | sum(check1) != 5){
      stop("MIMresult data error, please input all of the original output data of MIM.search()/MIM.points().", call. = FALSE)
    } else {
      check2 <- names(MIMresult$inputdata) == c("marker", "geno", "y", "yu", "sele.g", "type", "ng", "cM", "D.matrix" )
      if(!is.list(MIMresult$inputdata) | sum(check2) != 9){
        stop("MIMresult data error, please input all of the original output data of MIM.search()/MIM.points().", call. = FALSE)
      } else {
        if(length(MIMresult$inputdata$sele.g) != 0){
          if(length(MIMresult$inputdata$sele.g) != "n"){
            stop("MIMresult data error, this function does not support selective genotyping data.", call. = FALSE)
          }
        }
        if(is.data.frame(MIMresult$QTL.best) | is.matrix(MIMresult$QTL.best)){
          QTL <- as.matrix(MIMresult$QTL.best[,1:2])
        } else {stop("MIMresult data error, please input all of the original output data of MIM.search()/MIM.points().", call. = FALSE)}
        marker <- MIMresult$inputdata$marker
        geno <- MIMresult$inputdata$geno
        y <- MIMresult$inputdata$y
        type <- MIMresult$inputdata$type
        ng <- MIMresult$inputdata$ng
        cM <- MIMresult$inputdata$cM
        if(is.null(D.matrix)){
          D.matrix <- MIMresult$inputdata$D.matrix
        }
      }
    }
  }

  if(is.null(QTL) | is.null(marker) | is.null(geno) | is.null(D.matrix) | is.null(y)){
    stop("Input data is missing. The argument QTL, marker, geno, D.matrix, and y must be input.", call. = FALSE)
  }

  genotest <- table(c(geno))
  datatry <- try(geno*geno, silent=TRUE)
  if(class(datatry)[1] == "try-error" | FALSE %in% (names(genotest) %in% c(0, 1, 2))  | length(dim(geno)) != 2){
    stop("Genotype data error, please cheak your genotype data.", call. = FALSE)
  }

  marker <- as.matrix(marker)
  markertest <- c(ncol(marker) != 2, NA %in% marker, marker[,1] != sort(marker[,1]), nrow(marker) != ncol(geno))
  datatry <- try(marker*marker, silent=TRUE)
  if(class(datatry)[1] == "try-error" | TRUE %in% markertest){
    stop("Marker data error, or the number of marker does not match the genotype data.", call. = FALSE)
  }

  QTL <- as.matrix(QTL)
  datatry <- try(QTL*QTL, silent=TRUE)
  if(class(datatry)[1] == "try-error" | ncol(QTL) != 2 | NA %in% QTL | max(QTL[, 1]) > max(marker[, 1])){
    stop("QTL data error, please cheak your QTL data.", call. = FALSE)
  }

  for(i in 1:nrow(QTL)){
    ch0 <- QTL[i, 1]
    q0 <- QTL[i, 2]
    if(!ch0 %in% marker[, 1]){
      stop("The specified QTL is not in the position that can estimate by the marker data.", call. = FALSE)
    }
    if(q0 > max(marker[marker[, 1] == ch0, 2]) | q0 < min(marker[marker[, 1] == ch0, 2])){
      stop("The specified QTL is not in the position that can estimate by the marker data.", call. = FALSE)
    }
  }

  if(!type[1] %in% c("AI","RI","BC") | length(type) > 1){
    stop("Parameter type error, please input AI, RI, or BC.", call. = FALSE)
  }

  if(!is.numeric(ng) | length(ng) > 1 | min(ng) < 1){
    stop("Parameter ng error, please input a positive integer.", call. = FALSE)
  }
  ng <- round(ng)

  if(is.null(D.matrix) | is.null(y)){
    stop("Input data is missing, please cheak and fix", call. = FALSE)
  }

  datatry <- try(y%*%t(y), silent=TRUE)
  if(class(datatry)[1] == "try-error" | length(y) < 2){
    stop("Input data error, please check your input data.", call. = FALSE)
  }

  if(!is.numeric(crit) | length(crit) > 1 | min(crit) <= 0 | max(crit) >= 1){
    stop("Parameter crit error, please input a positive number between 0 and 1.", call. = FALSE)
  }

  if(!is.numeric(stop) | length(stop) > 1 | min(crit) <= 0){
    stop = 1000
    warning("Parameter stop error, adjust to 1000.")
  }

  if(!is.null(cp.matrix)){
    datatry <- try(y%*%cp.matrix%*%D.matrix, silent=TRUE)
    if(class(datatry)[1] == "try-error" | NA %in% D.matrix | NA %in% cp.matrix){
      stop("Input data error, please check your input data.", call. = FALSE)
    }
  } else {
    datatry <- try(y%*%t(y), silent=TRUE)
    datatry1 <- try(D.matrix*D.matrix, silent=TRUE)
    if(class(datatry)[1] == "try-error" | class(datatry1)[1] == "try-error" | NA %in% D.matrix | length(y) < 2){
      stop("Input data error, please check your input data.", call. = FALSE)
    }
    cp.matrix <- Q.make(QTL, marker, geno, type = type, ng = ng, interval = TRUE)$cp.matrix
  }

  Y <- c(y)
  ind <- length(Y)
  g <- nrow(D.matrix)
  eff <- ncol(D.matrix)

  Y[is.na(Y)] <- mean(Y,na.rm = TRUE)

  if(is.null(E.vector0)){E.vector0 <- rep(0, eff)}
  if(is.null(X)){
    X <- matrix(1, ind, 1)
    colnames(X) <- "mu"
  } else if (is.vector(X)){
    X <- matrix(X, length(X), 1)
  }
  if(is.null(beta0)){
    beta0 <- matrix(rep(mean(Y), ncol(X)), ncol(X), 1)
  } else if (is.numeric(beta)){
    beta0 <- matrix(rep(beta, ncol(X)), ncol(X), 1)
  }
  if(is.null(variance0)){variance0 <- stats::var(c(Y))}
  if(!console[1] %in% c(0,1) | length(console) > 1){console <- TRUE}
  if(!conv[1] %in% c(0,1) | length(conv) > 1){conv <- TRUE}

  datatry <- try(D.matrix%*%E.vector0, silent=TRUE)
  if(class(datatry)[1] == "try-error" | NA %in% E.vector0){
    stop("Parameter E.vector0 error, please check and fix.", call. = FALSE)
  }

  datatry <- try(Y%*%X%*%beta0, silent=TRUE)
  if(class(datatry)[1] == "try-error" | NA %in% X | NA %in% beta0){
    stop("Parameter X or bata0 error, please check and fix.", call. = FALSE)
  }

  if(!is.numeric(variance0) | length(variance0) > 1 | min(variance0) < 0){
    stop("Parameter variance0 error, please input a positive number.", call. = FALSE)
  }

  if(!is.numeric(crit) | length(crit) > 1 | min(crit) <= 0 | max(crit) >= 1){
    stop("Parameter crit error, please input a positive number between 0 and 1.", call. = FALSE)
  }

  if(!is.numeric(stop) | length(stop) > 1 | min(crit) <= 0){
    stop = 1000
    warning("Parameter stop error, adjust to 1000.")
  }

  EM <- EM.MIM0(D.matrix, cp.matrix, Y, E.vector0, X, beta0, variance0, crit, stop, conv, console)
  PI <- EM$PI.matrix
  sigma2 <- EM$variance
  E <- EM$E.vector
  beta <- EM$beta
  ind <- nrow(geno)

  ncd <- ncol(D.matrix)
  nx <- ncol(X)

  if(is.null(colnames(X))){colnames(X) <- paste("X", 1:nx, sep = "")}

  cr <- QTL[, 1]
  QTL <- QTL[, 2]
  if(cM){
    QTL <- QTL/100
    marker[, 2] <- marker[, 2]/100
  }

  marker <- cbind(marker, 1:nrow(marker))

  cr0 <- marker[marker[, 1]==cr[1], ]

  if(QTL[1] == cr0[1, 2]){
    marker1 <- cr0[1, ]
    marker2 <- cr0[2, ]
    marker0 <- as.numeric(c(marker1[3], marker2[3]))

    dMQ <- as.numeric(QTL[1]-marker1[2])
    dMN <- as.numeric(marker2[2]-marker1[2])
  } else {
    marker1 <- cr0[max(which(cr0[, 2] < QTL[1])), ]
    marker2 <- cr0[min(which(cr0[, 2] > QTL[1])), ]
    marker0 <- as.numeric(c(marker1[3], marker2[3]))

    dMQ <- as.numeric(QTL[1]-marker1[2])
    dMN <- as.numeric(marker2[2]-marker1[2])
  }

  nd <- 3
  if(type == "BC"){
    nd <- 2
    funQ <- function(dMN, dMQ, ng){
      R <- (1-exp(-2*dMN))/2
      rMQ <- (1-exp(-2*dMQ))/2
      P <- rMQ/R

      M11 <- (1-R)/2
      M10 <- R/2
      M01 <- R/2
      M00 <- (1-R)/2

      Q111 <- (1-R)/2
      Q110 <- R/2-P*R/2
      Q101 <- 0
      Q011 <- P*R/2
      Q001 <- R/2-P*R/2
      Q010 <- 0
      Q100 <- P*R/2
      Q000 <- (1-R)/2

      Q111_1 <- Q101_1 <- Q010_1 <- Q000_1 <- 0
      Q110_1 <- Q001_1 <- -R/2
      Q011_1 <- Q100_1 <- R/2

      lnQ111_1 <- lnQ101_1 <- lnQ010_1 <- lnQ000_1 <- 0
      lnQ110_1 <- lnQ001_1 <- -1/(1-P)
      lnQ011_1 <- lnQ100_1 <- 1/P

      Q111_2 <- Q110_2 <- Q101_2 <- Q011_2 <- Q001_2 <- Q010_2 <- Q100_2 <- Q000_2 <- 0

      lnQ111_2 <- lnQ101_2 <- lnQ010_2 <- lnQ000_2 <- 0
      lnQ110_2 <- lnQ001_2 <- -1/(1-P)^2
      lnQ011_2 <- lnQ100_2 <- -1/P^2


      if(ng > 1){
        for(i in 2:ng){
          M11[i] <- M11[i-1]+1/2*M10[i-1]+1/2*M01[i-1]+(1-R)/2*M00[i-1]
          M10[i] <- 1/2*M10[i-1]+R/2*M00[i-1]
          M01[i] <- 1/2*M01[i-1]+R/2*M00[i-1]
          M00[i] <- (1-R)/2*M00[i-1]

          Q111[i] <- Q111[i-1]+1/2*Q110[i-1]+1/2*Q101[i-1]+1/2*Q011[i-1]+
            (1-P*R)/2*Q001[i-1]+(1-R+P*R)/2*Q100[i-1]+(1-R)/2*Q000[i-1]
          Q110[i] <- 1/2*Q110[i-1]+(R-P*R)/2*Q100[i-1]+(R-P*R)/2*Q000[i-1]
          Q101[i] <- 1/2*Q101[i-1]+P*R/2*Q001[i-1]+(R-P*R)/2*Q100[i-1]
          Q011[i] <- 1/2*Q011[i-1]+P*R/2*Q001[i-1]+P*R/2*Q000[i-1]
          Q001[i] <- (1-P*R)/2*Q001[i-1]+(R-P*R)/2*Q000[i-1]
          Q010[i] <- 0
          Q100[i] <- (1-R+P*R)/2*Q100[i-1]+P*R/2*Q000[i-1]
          Q000[i] <- (1-R)/2*Q000[i-1]

          Q111_1[i] <- Q111_1[i-1]+1/2*Q110_1[i-1]+1/2*Q101_1[i-1]+1/2*Q011_1[i-1]+
            (1-P*R)/2*Q001_1[i-1]-R/2*Q001[i-1]+(1-R+P*R)/2*Q100_1[i-1]+R/2*Q100[i-1]+(1-R)/2*Q000_1[i-1]
          Q110_1[i] <- 1/2*Q110_1[i-1]+(R-P*R)/2*Q100_1[i-1]-R/2*Q100[i-1]+(R-P*R)/2*Q000_1[i-1]-R/2*Q000[i-1]
          Q101_1[i] <- 1/2*Q101_1[i-1]+P*R/2*Q001_1[i-1]+R/2*Q001[i-1]+(R-P*R)/2*Q100_1[i-1]-R/2*Q100[i-1]
          Q011_1[i] <- 1/2*Q011_1[i-1]+P*R/2*Q001_1[i-1]+R/2*Q001[i-1]+P*R/2*Q000_1[i-1]+R/2*Q000[i-1]
          Q001_1[i] <- (1-P*R)/2*Q001_1[i-1]-R/2*Q001[i-1]+(R-P*R)/2*Q000_1[i-1]-R/2*Q000[i-1]
          Q100_1[i] <- (1-R+P*R)/2*Q100_1[i-1]+R/2*Q100[i-1]+P*R/2*Q000_1[i-1]+R/2*Q000[i-1]
          Q000_1[i] <- (1-R)/2*Q000_1[i-1]

          lnQ111_1[i] <- Q111_1[i]/Q111[i]
          lnQ110_1[i] <- Q110_1[i]/Q110[i]
          lnQ101_1[i] <- Q101_1[i]/Q101[i]
          lnQ011_1[i] <- Q011_1[i]/Q011[i]
          lnQ001_1[i] <- Q001_1[i]/Q001[i]
          lnQ010_1[i] <- 0
          lnQ100_1[i] <- Q100_1[i]/Q100[i]
          lnQ000_1[i] <- Q000_1[i]/Q000[i]

          Q111_2[i] <- Q111_2[i-1]+1/2*Q110_2[i-1]+1/2*Q101_2[i-1]+1/2*Q011_2[i-1]+
            (1-P*R)/2*Q001_2[i-1]-R*Q001_1[i-1]+(1-R+P*R)/2*Q100_2[i-1]+R*Q100_1[i-1]+(1-R)/2*Q000_2[i-1]
          Q110_2[i] <- 1/2*Q110_2[i-1]+(R-P*R)/2*Q100_2[i-1]-R*Q100_1[i-1]+(R-P*R)/2*Q000_2[i-1]-R*Q000_1[i-1]
          Q101_2[i] <- 1/2*Q101_2[i-1]+P*R/2*Q001_2[i-1]+R*Q001_1[i-1]+(R-P*R)/2*Q100_2[i-1]-R*Q100_1[i-1]
          Q011_2[i] <- 1/2*Q011_2[i-1]+P*R/2*Q001_2[i-1]+R*Q001_1[i-1]+P*R/2*Q000_2[i-1]+R*Q000_1[i-1]
          Q001_2[i] <- (1-P*R)/2*Q001_2[i-1]-R*Q001_1[i-1]+(R-P*R)/2*Q100_2[i-1]-R*Q100_1[i-1]
          Q100_2[i] <- (1-R+P*R)/2*Q100_2[i-1]+R*Q100_1[i-1]+P*R/2*Q000_2[i-1]+R*Q000_1[i-1]
          Q000_2[i] <- (1-R)/2*Q000_2[i-1]

          lnQ111_2[i] <- (Q111_2[i]*Q111[i]-Q111_1[i]*Q111_1[i])/Q111[i]^2
          lnQ110_2[i] <- (Q110_2[i]*Q110[i]-Q110_1[i]*Q110_1[i])/Q110[i]^2
          lnQ101_2[i] <- (Q101_2[i]*Q101[i]-Q101_1[i]*Q101_1[i])/Q101[i]^2
          lnQ011_2[i] <- (Q011_2[i]*Q011[i]-Q011_1[i]*Q011_1[i])/Q011[i]^2
          lnQ001_2[i] <- (Q001_2[i]*Q001[i]-Q001_1[i]*Q001_1[i])/Q001[i]^2
          lnQ010_2[i] <- 0
          lnQ100_2[i] <- (Q100_2[i]*Q100[i]-Q100_1[i]*Q100_1[i])/Q100[i]^2
          lnQ000_2[i] <- (Q000_2[i]*Q000[i]-Q000_1[i]*Q000_1[i])/Q000[i]^2
        }
      }

      MN <- c(M11[ng], M10[ng], M01[ng], M00[ng])
      MQN <- c(Q111[ng], Q110[ng], Q101[ng], Q011[ng], Q001[ng], Q010[ng], Q100[ng], Q000[ng])

      lnMQN_1 <- c(lnQ111_1[ng], lnQ110_1[ng], lnQ101_1[ng], lnQ011_1[ng],
                   lnQ001_1[ng], lnQ010_1[ng], lnQ100_1[ng], lnQ000_1[ng])

      lnMQN_2 <- c(lnQ111_2[ng], lnQ110_2[ng], lnQ101_2[ng], lnQ011_2[ng],
                   lnQ001_2[ng], lnQ010_2[ng], lnQ100_2[ng], lnQ000_2[ng])

      Q11_1 <- c(Q111_1[ng], Q110_1[ng], Q101_1[ng], Q011_1[ng], Q001_1[ng], Q010_1[ng], Q100_1[ng], Q000_1[ng])

      Q11_2 <- c(Q111_2[ng], Q110_2[ng], Q101_2[ng], Q011_2[ng], Q001_2[ng], Q010_2[ng], Q100_2[ng], Q000_2[ng])

      Q0 <- c(MQN[1]/MN[1], MQN[3]/MN[1], MQN[2]/MN[2], MQN[7]/MN[2],
              MQN[4]/MN[3], MQN[5]/MN[3], MQN[6]/MN[4], MQN[8]/MN[4])
      Q1 <- matrix(Q0, 4, 2, byrow = TRUE)

      lnQ0_1 <- c(lnMQN_1[1], lnMQN_1[3], lnMQN_1[2], lnMQN_1[7], lnMQN_1[4], lnMQN_1[5], lnMQN_1[6], lnMQN_1[8])
      lnQ1_1 <- matrix(lnQ0_1, 4, 2, byrow = TRUE)
      lnQ1_1[is.na(lnQ1_1)] <- 0
      lnQ1_1[abs(lnQ1_1) == Inf] <- 0

      lnQ0_2 <- c(lnMQN_2[1], lnMQN_2[3], lnMQN_2[2], lnMQN_2[7], lnMQN_2[4], lnMQN_2[5], lnMQN_2[6], lnMQN_2[8])
      lnQ1_2 <- matrix(lnQ0_2, 4, 2, byrow = TRUE)
      lnQ1_2[is.na(lnQ1_2)] <- 0
      lnQ1_2[abs(lnQ1_2) == Inf] <- 0

      Q0_1 <- c(Q11_1[1]/MN[1], Q11_1[3]/MN[1], Q11_1[2]/MN[2], Q11_1[7]/MN[2],
                Q11_1[4]/MN[3], Q11_1[5]/MN[3], Q11_1[6]/MN[4], Q11_1[8]/MN[4])
      Q1_1 <- matrix(Q0_1, 4, 2, byrow = TRUE)
      Q1_1[is.na(Q1_1)] <- 0
      Q1_1[abs(Q1_1) == Inf] <- 0

      Q0_2 <- c(Q11_2[1]/MN[1], Q11_2[3]/MN[1], Q11_2[2]/MN[2], Q11_2[7]/MN[2],
                Q11_2[4]/MN[3], Q11_2[5]/MN[3], Q11_2[6]/MN[4], Q11_2[8]/MN[4])
      Q1_2 <- matrix(Q0_2, 4, 2, byrow = TRUE)
      Q1_2[is.na(Q1_2)] <- 0
      Q1_2[abs(Q1_2) == Inf] <- 0

      colnames(Q1) <- colnames(lnQ1_1) <- colnames(lnQ1_2) <- colnames(Q1_1) <- colnames(Q1_2) <- c("QQ", "Qq")
      rownames(Q1) <- rownames(lnQ1_1) <- rownames(lnQ1_2) <- rownames(Q1_1) <- rownames(Q1_2) <- c(22, 21, 12, 11)

      re <- list(Q1, lnQ1_1, lnQ1_2, Q1_1, Q1_2)
      return(re)
    }
  } else if (type =="RI"){
    funQ <- function(dMN, dMQ, ng){
      r <- (1-exp(-2*dMN))/2
      rMQ <- (1-exp(-2*dMQ))/2
      P <- rMQ/r
      c0 <- (1-r)^2+r^2

      TMN <- matrix(1:25, nrow = 5)
      TMN[1, ] <- c(4, 0, 2, (1-r)^2, r^2)/4
      TMN[2, ] <- c(0, 4, 2, r^2, (1-r)^2)/4
      TMN[3, ] <- c(0, 0, 1, r*(1-r), r*(1-r))/2
      TMN[4, ] <- c(0, 0, 0, (1-r)^2, r^2)/2
      TMN[5, ] <- c(0, 0, 0, r^2, (1-r)^2)/2
      Freq <- c(0, 0, 0, 1, 0)
      n.Iteration <- ng-1
      ZygoteFreq <- matrix(rep(0, (9*ng-9)), nrow = 9)
      for (i in 1:n.Iteration) {
        Freq <- TMN %*% Freq
        ZygoteFreq[, i] <- c(Freq[c(1, 3, 2, 3)], sum(Freq[4:5]),
                             Freq[c(3, 2, 3, 1)])
      }
      MN <- ZygoteFreq[, n.Iteration]

      TMQN <- TMQN_1 <- TMQN_2 <- matrix(rep(0, 20^2), nrow = 20)
      TMQN[1, ] <- c(4, 1, 0, 1, (1-r+P*r)^2, (r-P*r)^2, 0, 0, 0, 0, 1, (1-P*r)^2, (P*r)^2, 0, (1-r)^2,
                     r^2, (1-r)^2,  0, (r-P*r)^2, (P*r)^2)/4
      TMQN[2, ] <- c(0, 1, 0, 0, (r-P*r)*(1-r+P*r), (r-P*r)*(1-r+P*r), 0, 0, 0, 0, 0, (P*r)*(1-P*r), (P*r)*(1-P*r),
                     0, 0, 0, 0, 0, (r-P*r)*(P*r), (r-P*r)*(P*r))/2
      TMQN[3, ] <- c(0, 1, 4, 0, (r-P*r)^2, (1-r+P*r)^2, 1, 0, 0, 0, 0, (P*r)^2, (1-P*r)^2, 1, (1-r)^2,
                     r^2,  0, (1-r)^2, (P*r)^2, (r-P*r)^2)/4
      TMQN[4, ] <- c(0, 0, 0, 1, (r-P*r)*(1-r+P*r), (r-P*r)*(1-r+P*r), rep(0, 8), r*(1-r), r*(1-r),
                     (r-P*r)*(1-r), 0, (r-P*r)*(1-r), 0)/2
      TMQN[5, ] <- c(0, 0, 0, 0, (1-r+P*r)^2, (r-P*r)^2, rep(0, 10), (P*r)*(1-r),
                     0, 0, (P*r)*(1-r))/2
      TMQN[6, ] <- c(0, 0, 0, 0, (r-P*r)^2, (1-r+P*r)^2, rep(0, 10), 0,
                     (P*r)*(1-r), (P*r)*(1-r), 0)/2
      TMQN[7, ] <- c(0, 0, 0, 0, (r-P*r)*(1-r+P*r), (r-P*r)*(1-r+P*r), 1, rep(0, 7), r*(1-r), r*(1-r),
                     0, (r-P*r)*(1-r), 0, (r-P*r)*(1-r))/2
      TMQN[8, ] <- c(0, 0, 0, 1, (r-P*r)^2, (1-r+P*r)^2, 0, 4, 1, 0, 0, (1-P*r)^2, (P*r)^2, 1, r^2,
                     (1-r)^2, (r-P*r)^2, (P*r)^2, (1-r)^2,  0)/4
      TMQN[9, ] <- c(0, 0, 0, 0, (r-P*r)*(1-r+P*r), (r-P*r)*(1-r+P*r), 0, 0, 1, 0, 0, (P*r)*(1-P*r),
                     (P*r)*(1-P*r), 0, 0, 0, (r-P*r)*(P*r), (r-P*r)*(P*r), 0, 0)/2
      TMQN[10, ] <- c(0, 0, 0, 0, (1-r+P*r)^2, (r-P*r)^2, 1, 0, 1, 4, 1, (P*r)^2, (1-P*r)^2, 0, r^2,
                      (1-r)^2, (P*r)^2, (r-P*r)^2,  0, (1-r)^2)/4
      TMQN[11, ] <- c(rep(0, 10), 1, (P*r)*(1-P*r), (P*r)*(1-P*r), 0, r*(1-r), r*(1-r),
                      (P*r)*(1-r), 0, 0, (P*r)*(1-r))/2
      TMQN[12, ] <- c(rep(0, 11), (1-P*r)^2, (P*r)^2, 0, 0, 0, (r-P*r)*(1-r),
                      0, (r-P*r)*(1-r), 0)/2
      TMQN[13, ] <- c(rep(0, 11), (P*r)^2, (1-P*r)^2, 0, 0, 0, 0, (r-P*r)*(1-r),
                      0, (r-P*r)*(1-r))/2
      TMQN[14, ] <- c(rep(0, 11), (P*r)*(1-P*r), (P*r)*(1-P*r), 1, r*(1-r), r*(1-r), 0,
                      (P*r)*(1-r), (P*r)*(1-r), 0)/2
      TMQN[15, ] <- c(rep(0, 14), (1-r)^2, r^2, 0, 0, (r-P*r)*(P*r), (r-P*r)*(P*r))/2
      TMQN[16, ] <- c(rep(0, 14), r^2, (1-r)^2, (r-P*r)*(P*r), (r-P*r)*(P*r), 0, 0)/2
      TMQN[17, ] <- c(rep(0, 16), (1-r)^2,  0, (r-P*r)^2, (P*r)^2)/2
      TMQN[18, ] <- c(rep(0, 16),  0, (1-r)^2, (P*r)^2, (r-P*r)^2)/2
      TMQN[19, ] <- c(rep(0, 16), (r-P*r)^2, (P*r)^2, (1-r)^2,  0)/2
      TMQN[20, ] <- c(rep(0, 16), (P*r)^2, (r-P*r)^2,  0, (1-r)^2)/2


      TMQN_1[1, ] <- c(0, 0, 0, 0, 2*r*(1-r+P*r), -2*r^2*(1-P), 0, 0, 0, 0, 0, -2*r*(1-P*r), 2*P*r^2, 0, 0,
                       0, 0,  0, -2*r^2*(1-P), 2*P*r^2)/4
      TMQN_1[2, ] <- c(0, 0, 0, 0, -r*(1-2*r+2*P*r), -r*(1-2*r+2*P*r), 0, 0, 0, 0, 0, r*(1-2*P*r), r*(1-2*P*r),
                       0, 0, 0, 0, 0, r^2*(1-2*P), r^2*(1-2*P))/2
      TMQN_1[3, ] <- c(0, 0, 0, 0, -2*r^2*(1-P), 2*r*(1-r+P*r), 0, 0, 0, 0, 0, 2*P*r^2, -2*r*(1-P*r), 0, 0,
                       0,  0, 0, 2*P*r^2, -2*r^2*(1-P))/4
      TMQN_1[4, ] <- c(0, 0, 0, 0, -r*(1-2*r+2*P*r), -r*(1-2*r+2*P*r), rep(0, 8),  0,  0,
                       -r*(1-r), 0, -r*(1-r), 0)/2
      TMQN_1[5, ] <- c(0, 0, 0, 0, 2*r*(1-r+P*r), -2*r^2*(1-P), rep(0, 10), r*(1-r),
                       0, 0, r*(1-r))/2
      TMQN_1[6, ] <- c(0, 0, 0, 0, -2*r^2*(1-P), 2*r*(1-r+P*r), rep(0, 10), 0,
                       r*(1-r), r*(1-r), 0)/2
      TMQN_1[7, ] <- c(0, 0, 0, 0, -r*(1-2*r+2*P*r), -r*(1-2*r+2*P*r), 0, rep(0, 7),  0,  0,
                       0, -r*(1-r), 0, -r*(1-r))/2
      TMQN_1[8, ] <- c(0, 0, 0, 0, -2*r^2*(1-P), 2*r*(1-r+P*r), 0, 0, 0, 0, 0, -2*r*(1-P*r), 2*P*r^2, 0, 0,
                       0, -2*r^2*(1-P), 2*P*r^2, 0,  0)/4
      TMQN_1[9, ] <- c(0, 0, 0, 0, -r*(1-2*r+2*P*r), -r*(1-2*r+2*P*r), 0, 0, 0, 0, 0, r*(1-2*P*r),
                       r*(1-2*P*r), 0, 0, 0, r^2*(1-2*P), r^2*(1-2*P), 0, 0)/2
      TMQN_1[10, ] <- c(0, 0, 0, 0, 2*r*(1-r+P*r), -2*r^2*(1-P), 0, 0, 0, 0, 0, 2*P*r^2, -2*r*(1-P*r), 0, 0,
                        0, 2*P*r^2, -2*r^2*(1-P),  0, 0)/4
      TMQN_1[11, ] <- c(rep(0, 10), 0, r*(1-2*P*r), r*(1-2*P*r), 0,  0,  0, r*(1-r), 0, 0, r*(1-r))/2
      TMQN_1[12, ] <- c(rep(0, 11), -2*r*(1-P*r), 2*P*r^2, 0, 0, 0, -r*(1-r), 0, -r*(1-r), 0)/2
      TMQN_1[13, ] <- c(rep(0, 11), 2*P*r^2, -2*r*(1-P*r), 0, 0, 0, 0, -r*(1-r), 0, -r*(1-r))/2
      TMQN_1[14, ] <- c(rep(0, 11), r*(1-2*P*r), r*(1-2*P*r), 0,  0,  0, 0, r*(1-r), r*(1-r), 0)/2
      TMQN_1[15, ] <- c(rep(0, 14), 0, 0, 0, 0, r^2*(1-2*P), r^2*(1-2*P))/2
      TMQN_1[16, ] <- c(rep(0, 14), 0, 0, r^2*(1-2*P), r^2*(1-2*P), 0, 0)/2
      TMQN_1[17, ] <- c(rep(0, 16), 0,  0, -2*r^2*(1-P), 2*P*r^2)/2
      TMQN_1[18, ] <- c(rep(0, 16),  0, 0, 2*P*r^2, -2*r^2*(1-P))/2
      TMQN_1[19, ] <- c(rep(0, 16), -2*r^2*(1-P), 2*P*r^2, 0,  0)/2
      TMQN_1[20, ] <- c(rep(0, 16), 2*P*r^2, -2*r^2*(1-P),  0, 0)/2


      TMQN_2[1, ] <- c(0, 0, 0, 0, 2*r^2, 2*r^2, 0, 0, 0, 0, 0, 2*r^2, 2*r^2, 0, 0, 0, 0,  0, 2*r^2, 2*r^2)/4
      TMQN_2[2, ] <- c(0, 0, 0, 0, -2*r^2, -2*r^2, 0, 0, 0, 0, 0, -2*r^2, -2*r^2, 0, 0, 0, 0, 0, -2*r^2, -2*r^2)/2
      TMQN_2[3, ] <- c(0, 0, 0, 0, 2*r^2, 2*r^2, 0, 0, 0, 0, 0, 2*r^2, 2*r^2, 0, 0, 0,  0, 0, 2*r^2, 2*r^2)/4
      TMQN_2[4, ] <- c(0, 0, 0, 0, -2*r^2, -2*r^2, rep(0, 8),  0,  0, 0, 0, 0, 0)/2
      TMQN_2[5, ] <- c(0, 0, 0, 0, 2*r^2, 2*r^2, rep(0, 10), 0, 0, 0, 0)/2
      TMQN_2[6, ] <- c(0, 0, 0, 0, 2*r^2, 2*r^2, rep(0, 10), 0, 0, 0, 0)/2
      TMQN_2[7, ] <- c(0, 0, 0, 0, -2*r^2, -2*r^2, 0, rep(0, 7),  0,  0, 0, 0, 0, 0)/2
      TMQN_2[8, ] <- c(0, 0, 0, 0, 2*r^2, 2*r^2, 0, 0, 0, 0, 0, 2*r^2, 2*r^2, 0, 0, 0, 2*r^2, 2*r^2, 0,  0)/4
      TMQN_2[9, ] <- c(0, 0, 0, 0, -2*r^2, -2*r^2, 0, 0, 0, 0, 0, -2*r^2, -2*r^2, 0, 0, 0, -2*r^2, -2*r^2, 0, 0)/2
      TMQN_2[10, ] <- c(0, 0, 0, 0, 2*r^2, 2*r^2, 0, 0, 0, 0, 0, 2*r^2, 2*r^2, 0, 0, 0, 2*r^2, 2*r^2,  0, 0)/4
      TMQN_2[11, ] <- c(rep(0, 10), 0, -2*r^2, -2*r^2, 0,  0,  0, 0, 0, 0, 0)/2
      TMQN_2[12, ] <- c(rep(0, 11), 2*r^2, 2*r^2, 0, 0, 0, 0, 0, 0, 0)/2
      TMQN_2[13, ] <- c(rep(0, 11), 2*r^2, 2*r^2, 0, 0, 0, 0, 0, 0, 0)/2
      TMQN_2[14, ] <- c(rep(0, 11), -2*r^2, -2*r^2, 0,  0,  0, 0, 0, 0, 0)/2
      TMQN_2[15, ] <- c(rep(0, 14), 0, 0, 0, 0, -2*r^2, -2*r^2)/2
      TMQN_2[16, ] <- c(rep(0, 14), 0, 0, -2*r^2, -2*r^2, 0, 0)/2
      TMQN_2[17, ] <- c(rep(0, 16), 0,  0, 2*r^2, 2*r^2)/2
      TMQN_2[18, ] <- c(rep(0, 16),  0, 0, 2*r^2, 2*r^2)/2
      TMQN_2[19, ] <- c(rep(0, 16), 2*r^2, 2*r^2, 0,  0)/2
      TMQN_2[20, ] <- c(rep(0, 16), 2*r^2, 2*r^2,  0, 0)/2

      FMQN0 <- c(rep(0, 16), 1, 0, 0, 0)
      n.Iteration <- ng - 1
      ZMQN <- ZMQN_1 <- ZMQN_2 <- matrix(rep(0, (20*ng - 20)), nrow = 20)

      FMQN <- TMQN%*%FMQN0
      FMQN_1 <- TMQN_1%*%FMQN0
      FMQN_2 <- TMQN_2%*%FMQN0
      ZMQN[, 1] <- FMQN0 <- FMQN
      ZMQN_1[, 1] <- FMQN_1
      ZMQN_2[, 1] <- FMQN_2

      if(ng > 2){
        for (i in 2:n.Iteration) {
          FMQN <- TMQN%*%FMQN
          FMQN_11 <- TMQN_1%*%FMQN0+TMQN%*%FMQN_1
          FMQN_22 <- TMQN_2%*%FMQN0+2*TMQN_1%*%FMQN_1+TMQN%*%FMQN_2
          ZMQN[, i] <- FMQN0 <- FMQN
          ZMQN_1[, i] <- FMQN_1 <- FMQN_11
          ZMQN_2[, i] <- FMQN_2 <- FMQN_22
        }
      }
      A <- ZMQN[, n.Iteration]
      B <- ZMQN_1[, n.Iteration]
      D <- ZMQN_2[, n.Iteration]

      q0 <- c(A[1:3], A[4], A[5]+A[6], A[7], A[8:10],
              A[11], A[12]+A[13], A[14],
              A[15]+A[16], sum(A[17:20]), A[15]+A[16],
              A[14], A[12]+A[13], A[11],
              A[10:8], A[7], A[5]+A[6], A[4], A[3:1])
      M0 <- MN[rep(1:9, each = 3)]
      Q0 <- q0/M0
      Q1 <- matrix(Q0, 9, 3, byrow = TRUE)

      Q0_1 <- c(B[1:3], B[4], B[5]+B[6], B[7], B[8:10],
                B[11], B[12]+B[13], B[14],
                B[15]+B[16], sum(B[17:20]), B[15]+B[16],
                B[14], B[12]+B[13], B[11],
                B[10:8], B[7], B[5]+B[6], B[4], B[3:1])
      Q0_2 <- c(D[1:3], D[4], D[5]+D[6], D[7], D[8:10],
                D[11], D[12]+D[13], D[14],
                D[15]+D[16], sum(D[17:20]), D[15]+D[16],
                D[14], D[12]+D[13], D[11],
                D[10:8], D[7], D[5]+D[6], D[4], D[3:1])

      lnQ0_1 <- Q0_1/q0
      lnQ1_1 <- matrix(lnQ0_1, 9, 3, byrow = TRUE)
      lnQ1_1[is.na(lnQ1_1)] <- 0
      lnQ1_1[abs(lnQ1_1) == Inf] <- 0

      lnQ0_2 <- (Q0_2*q0-Q0_1*Q0_1)/(q0)^2
      lnQ1_2 <- matrix(lnQ0_2, 9, 3, byrow = TRUE)
      lnQ1_2[is.na(lnQ1_2)] <- 0
      lnQ1_2[abs(lnQ1_2) == Inf] <- 0

      Q1_1 <- matrix(Q0_1/M0, 9, 3, byrow = TRUE)
      Q1_1[is.na(Q1_1)] <- 0
      Q1_1[abs(Q1_1) == Inf] <- 0

      Q1_2 <- matrix(Q0_2/M0, 9, 3, byrow = TRUE)
      Q1_2[is.na(Q1_2)] <- 0
      Q1_2[abs(Q1_2) == Inf] <- 0

      colnames(Q1) <- colnames(lnQ1_1) <- colnames(lnQ1_2) <- colnames(Q1_1) <- colnames(Q1_2) <- c("QQ", "Qq", "qq")
      rownames(Q1) <- rownames(lnQ1_1) <- rownames(lnQ1_2) <- rownames(Q1_1) <- rownames(Q1_2) <- c(22, 21, 20, 12, 11, 10, "02", "01", "00")

      re <- list(Q1, lnQ1_1, lnQ1_2, Q1_1, Q1_2)
      return(re)
    }
  } else {
    funQ <- function(dMN, dMQ, ng){
      r <- (1-exp(-2*dMN))/2
      rMQ <- (1-exp(-2*dMQ))/2
      P <- rMQ/r

      ZygoteFreq <- matrix(rep(0, (9*ng-9)), nrow = 9)
      GameteFreq <- c((1-r)/2, r/2)
      C <- GameteFreq[1]^2
      D <- GameteFreq[2]^2
      E <- 2*GameteFreq[1]*GameteFreq[2]
      F0 <- 2*GameteFreq[1]^2
      G <- 2*GameteFreq[2]^2
      ZygoteFreq[, 1] <- c(C, E, D, E, F0+G, E, D, E, C)
      if (ng > 2) {
        for (i in 3:(ng)) {
          GameteFreq <- (1-r)*GameteFreq+r/4
          C <- GameteFreq[1]^2
          D <- GameteFreq[2]^2
          E <- 2*GameteFreq[1]*GameteFreq[2]
          F0 <- 2*GameteFreq[1]^2
          G <- 2*GameteFreq[2]^2
          ZygoteFreq[, (i-1)] <- c(C, E, D, E, F0+G, E, D, E, C)
        }
      }
      MN <- ZygoteFreq[, (ng-1)]

      ZMQN <- ZMQN_1 <- ZMQN_2 <- matrix(rep(0, (20*ng-20)), nrow = 20)
      N <- (1-r)/2
      R <- (r-P*r)/2
      B <- 0
      L <- P*r/2
      COEF <- rep(2, 20)
      COEF[c(1, 3, 8, 10)] <- rep(1, 4)
      GS0 <- c(N, N, B, N, N, R, B, R, R, L, N, N, B, R, N, R, N, B, R, L)
      GE0 <- c(N, B, B, R, L, B, L, R, L, L, L, R, L, B, B, L, N, B, R, L)
      ZMQN[, 1] <- COEF*GS0*GE0
      N1 <- 0
      R1 <- -r/2
      B1 <- 0
      L1 <- r/2
      GS1 <- c(N1, N1, B1, N1, N1, R1, B1, R1, R1, L1, N1, N1, B1, R1, N1, R1, N1, B1, R1, L1)
      GE1 <- c(N1, B1, B1, R1, L1, B1, L1, R1, L1, L1, L1, R1, L1, B1, B1, L1, N1, B1, R1, L1)
      N2 <- R2 <- B2 <- L2 <- 0
      ZMQN_1[, 1] <- COEF*(GS1*GE0+GS0*GE1)
      ZMQN_2[, 1] <- COEF*(2*GS1*GE1)

      if (ng > 2) {
        for (i in 3:ng) {
          F23 <- c(N+L, B+R, B+R, N+L)
          F13 <- rep(c(N+B, L+R), 2)
          F12 <- rep(c(N+R, L+B), each = 2)
          F23_1 <- c(N1+L1, B1+R1, B1+R1, N1+L1)
          F12_1 <- rep(c(N1+R1, L1+B1), each = 2)
          F23_2 <- c(N2+L2, B2+R2, B2+R2, N2+L2)
          F12_2 <- rep(c(N2+R2, L2+B2), each = 2)
          NGF <- (1-r)*c(N, R, B, L)+0.5*(P*r)*F23+0*F13+0.5*(r-P*r)*F12
          NGF1 <- (1-r)*c(N1, R1, B1, L1)+0.5*r*F23+0.5*(P*r)*F23_1+0*F13-0.5*r*F12+0.5*(r-P*r)*F12_1
          NGF2 <- (1-r)*c(N2, R2, B2, L2)+r*F23_1+0.5*(P*r)*F23_2+0-r*F12_1+0.5*(r-P*r)*F12_2

          N <- NGF[1]
          R <- NGF[2]
          B <- NGF[3]
          L <- NGF[4]
          GS <- c(N, N, B, N, N, R, B, R, R, L, N, N, B, R, N, R, N, B, R, L)
          GE <- c(N, B, B, R, L, B, L, R, L, L, L, R, L, B, B, L, N, B, R, L)

          N1 <- NGF1[1]
          R1 <- NGF1[2]
          B1 <- NGF1[3]
          L1 <- NGF1[4]
          GS1 <- c(N1, N1, B1, N1, N1, R1, B1, R1, R1, L1, N1, N1, B1, R1, N1, R1, N1, B1, R1, L1)
          GE1 <- c(N1, B1, B1, R1, L1, B1, L1, R1, L1, L1, L1, R1, L1, B1, B1, L1, N1, B1, R1, L1)

          N2 <- NGF2[1]
          R2 <- NGF2[2]
          B2 <- NGF2[3]
          L2 <- NGF2[4]
          GS2 <- c(N2, N2, B2, N2, N2, R2, B2, R2, R2, L2, N2, N2, B2, R2, N2, R2, N2, B2, R2, L2)
          GE2 <- c(N2, B2, B2, R2, L2, B2, L2, R2, L2, L2, L2, R2, L2, B2, B2, L2, N2, B2, R2, L2)

          ZMQN[, (i-1)] <- COEF*GS*GE
          ZMQN_1[, (i-1)] <- COEF*(GS1*GE+GS*GE1)
          ZMQN_2[, 1] <- COEF*(GS2*GE+2*GS1*GE1+GS*GE2)
        }
      }
      A <- ZMQN[, (ng-1)]
      B <- ZMQN_1[, (ng-1)]
      D <- ZMQN_2[, (ng-1)]

      q0 <- c(A[1:3], A[4], A[5]+A[6], A[7], A[8:10],
              A[11], A[12]+A[13], A[14],
              A[15]+A[16], sum(A[17:20]), A[15]+A[16],
              A[14], A[12]+A[13], A[11],
              A[10:8], A[7], A[5]+A[6], A[4], A[3:1])
      M0 <- MN[rep(1:9, each = 3)]
      Q0 <- q0/M0
      Q1 <- matrix(Q0, 9, 3, byrow = TRUE)

      Q0_1 <- c(B[1:3], B[4], B[5]+B[6], B[7], B[8:10],
                B[11], B[12]+B[13], B[14],
                B[15]+B[16], sum(B[17:20]), B[15]+B[16],
                B[14], B[12]+B[13], B[11],
                B[10:8], B[7], B[5]+B[6], B[4], B[3:1])
      Q0_2 <- c(D[1:3], D[4], D[5]+D[6], D[7], D[8:10],
                D[11], D[12]+D[13], D[14],
                D[15]+D[16], sum(D[17:20]), D[15]+D[16],
                D[14], D[12]+D[13], D[11],
                D[10:8], D[7], D[5]+D[6], D[4], D[3:1])

      lnQ0_1 <- Q0_1/q0
      lnQ1_1 <- matrix(lnQ0_1, 9, 3, byrow = TRUE)
      lnQ1_1[is.na(lnQ1_1)] <- 0
      lnQ1_1[abs(lnQ1_1) == Inf] <- 0

      lnQ0_2 <- (Q0_2*q0-Q0_1*Q0_1)/(q0)^2
      lnQ1_2 <- matrix(lnQ0_2, 9, 3, byrow = TRUE)
      lnQ1_2[is.na(lnQ1_2)] <- 0
      lnQ1_2[abs(lnQ1_2) == Inf] <- 0

      Q1_1 <- matrix(Q0_1/M0, 9, 3, byrow = TRUE)
      Q1_1[is.na(Q1_1)] <- 0
      Q1_1[abs(Q1_1) == Inf] <- 0

      Q1_2 <- matrix(Q0_2/M0, 9, 3, byrow = TRUE)
      Q1_2[is.na(Q1_2)] <- 0
      Q1_2[abs(Q1_2) == Inf] <- 0

      colnames(Q1) <- colnames(lnQ1_1) <- colnames(lnQ1_2) <- colnames(Q1_1) <- colnames(Q1_2) <- c("QQ", "Qq", "qq")
      rownames(Q1) <- rownames(lnQ1_1) <- rownames(lnQ1_2) <- rownames(Q1_1) <- rownames(Q1_2) <- c(22, 21, 20, 12, 11, 10, "02", "01", "00")

      re <- list(Q1, lnQ1_1, lnQ1_2, Q1_1, Q1_2)
      return(re)
    }
  }

  Q0 <- funQ(dMN, dMQ, ng)
  Q1 <- list(Q0[[1]])
  lnQ1_1 <- list(Q0[[2]])
  lnQ1_2 <- list(Q0[[3]])
  Q1_1 <- list(Q0[[4]])
  Q1_2 <- list(Q0[[5]])

  m <- 1
  nQTL <- length(QTL)
  if(nQTL > 1){
    for(m in 2:nQTL){
      cr0 <- marker[marker[, 1] == cr[m], ]

      if(QTL[m] == cr0[1, 2]){
        marker1 <- cr0[1, ]
        marker2 <- cr0[2, ]

        dMQ <- as.numeric(QTL[m]-marker1[2])
        dMN <- as.numeric(marker2[2]-marker1[2])
      } else {
        marker1 <- cr0[max(which(cr0[, 2] < QTL[m])), ]
        marker2 <- cr0[min(which(cr0[, 2] > QTL[m])), ]
        marker0 <- as.numeric(c(marker0, marker1[3], marker2[3]))

        dMQ <- as.numeric(QTL[m]-marker1[2])
        dMN <- as.numeric(marker2[2]-marker1[2])
      }

      Q0 <- funQ(dMN, dMQ, ng)
      Q1[[m]] <- Q0[[1]]
      lnQ1_1[[m]] <- Q0[[2]]
      lnQ1_2[[m]] <- Q0[[3]]
      Q1_1[[m]] <- Q0[[4]]
      Q1_2[[m]] <- Q0[[5]]
    }
  }

  ngt <- nd^nQTL
  red.genotype <- geno[, marker0]

  invector <- rep(1,ind)
  iqvector <- rep(1,ngt)
  T.matrix <- (y-X%*%beta)%x%t(iqvector)-invector%x%t(D.matrix%*%E)
  S.matrix <- T.matrix^2/(2*sigma2)-0.5
  TPI <- T.matrix*PI
  SPII <- (S.matrix*PI)%*%iqvector

  mnd0 <- matrix(0, ind, nd)
  colnames(mnd0) <- 2:(3-nd)

  namePI <- matrix(unlist(strsplit(colnames(PI), split = "")), ngt, nQTL, byrow = TRUE)

  D2 <- gtools::permutations(nd, nQTL, 2:(3-nd), FALSE, TRUE)
  qname <- apply(D2, 1, paste, collapse = "")
  D3 <- -D2+3

  P1I0 <- P2I0 <- list()
  for(i in 1:nQTL){
    P1I0[[i]] <- P2I0[[i]] <- mnd0
  }

  for(j in 1:ind){
    geno.j <- red.genotype[j, ]
    for(k in 1:nQTL){
      g0 <- c(geno.j)[(k*2-1):(k*2)]
      if(!is.na(g0[1]) & !is.na(g0[2])){
        g0 <- paste(g0[1], g0[2], sep = "")
        p11 <- lnQ1_1[[k]]
        p11 <- p11[row.names(p11) == g0]
        p22 <- lnQ1_2[[k]]
        p22 <- p22[row.names(p22) == g0]
      } else if (!is.na(g0[1]) & is.na(g0[2])){
        a <- row.names(Q1[[k]])
        a <- as.numeric(unlist(strsplit(a,split = "", fixed = TRUE)) )
        a <- t(matrix(a, 2))
        a0 <- Q1[[k]][a[, 1] == as.numeric(g0[1]), ]
        a0 <- matrix(unlist(a0), nrow(a0), ncol(a0))
        b0 <- Q1_1[[k]][a[, 1] == as.numeric(g0[1]), ]
        b0 <- matrix(unlist(b0), nrow(b0), ncol(b0))
        p11 <- apply(b0, 2, sum)/apply(a0, 2, sum)
        d0 <- Q1_2[[k]][a[, 1] == as.numeric(g0[1]), ]
        d0 <- matrix(unlist(d0), nrow(d0), ncol(d0))
        p22 <- apply(d0, 2, sum)/apply(a0, 2, sum)
      } else if (is.na(g0[1]) & !is.na(g0[2])){
        a <- row.names(Q1[[k]])
        a <- as.numeric(unlist(strsplit(a,split = "", fixed = TRUE)) )
        a <- t(matrix(a, 2))
        a0 <- Q1[[k]][a[, 2] == as.numeric(g0[2]), ]
        a0 <- matrix(unlist(a0), nrow(a0), ncol(a0))
        b0 <- Q1_1[[k]][a[, 2] == as.numeric(g0[2]), ]
        b0 <- matrix(unlist(b0), nrow(b0), ncol(b0))
        p11 <- apply(b0, 2, sum)/apply(a0, 2, sum)
        d0 <- Q1_2[[k]][a[, 2] == as.numeric(g0[2]), ]
        d0 <- matrix(unlist(d0), nrow(d0), ncol(d0))
        p22 <- apply(d0, 2, sum)/apply(a0, 2, sum)
      } else {
        a0 <- Q1[[k]]
        a0 <- matrix(unlist(a0), nrow(a0), ncol(a0))
        b0 <- Q1_1[[k]]
        b0 <- matrix(unlist(b0), nrow(b0), ncol(b0))
        p11 <- apply(b0, 2, sum)/apply(a0, 2, sum)
        d0 <- Q1_2[[k]]
        d0 <- matrix(unlist(d0), nrow(d0), ncol(d0))
        p22 <- apply(d0, 2, sum)/apply(a0, 2, sum)
      }
      p11[is.na(p11)] <- 0
      p22[is.na(p22)] <- 0
      P1I0[[k]][j, ] <- p11
      P2I0[[k]][j, ] <- p22
    }
  }

  P1I <- P2I <- list()
  for(i in 1:nQTL){
    p10 <- p20 <- matrix(0, ind, ngt)
    np10 <- colnames(P1I0[[i]])
    for(j in 1:nd){
      p10[, namePI[, i] == np10[j]] <- P1I0[[i]][, j]
      p20[, namePI[, i] == np10[j]] <- P2I0[[i]][, j]
    }
    P1I[[i]] <- p10
    P2I[[i]] <- p20
  }

  if(var.pos){
    IOC <- IOM <- matrix(0, ncd+nQTL+nx+1, ncd+nQTL+nx+1)
    for(i in 1:nQTL){
      IOC[i, i] <- -sigma2*t(invector)%*%(P2I[[i]]*PI)%*%iqvector
    }

    for(i in 1:nQTL){
      P1IPI <- P1I[[i]]*PI
      P1IPII <- (P1IPI)%*%iqvector
      for(j in i:nQTL){
        IOM[i, j] <- IOM[j, i] <- t(invector)%*%(P1I[[j]]*P1IPI)%*%iqvector+
          sum(((P1I[[j]]*PI)%*%iqvector)%*%t(P1IPII))-sum(((P1I[[j]]*PI)%*%iqvector)*P1IPII)
      }

      P1TPIi <- P1I[[i]]*TPI
      for(j in 1:ncd){
        Di <- D.matrix[, j]
        IOM[i, nQTL+j] <- IOM[nQTL+j, i] <- (t(invector)%*%P1TPIi%*%Di+sum(P1IPII%*%t(TPI%*%Di))-sum(P1IPII*(TPI%*%Di)))/sigma2
      }

      IOM[i, nQTL+ncd+1] <- IOM[nQTL+ncd+1, i] <- (t(invector)%*%(S.matrix*P1IPI)%*%iqvector+sum(P1IPII%*%t(SPII))-sum(P1IPII*SPII))/sigma2
      for(j in 1:nx){
        Xi <- X[, j]
        qTPIX <- t(iqvector)%*%t(TPI)*Xi
        IOM[i, nQTL+ncd+1+j] <- IOM[nQTL+ncd+1+j, i] <- (t(iqvector)%*%t(P1TPIi)%*%Xi+
                                                           sum(P1IPII%*%qTPIX)-sum(P1IPII*t(qTPIX)))/sigma2
      }
    }
  } else {
    nQTL <- 0
    IOC <- IOM <- matrix(0, ncd+nx+1, ncd+nx+1)
  }

  for(i in 1:ncd){
    Di <- D.matrix[, i]
    PDD <- t(invector)%*%PI%*%((Di%x%t(rep(1, ncd)))*D.matrix)
    PTD <- t(invector)%*%TPI%*%Di/sigma2
    IOC0 <- c(rep(0, nQTL), PDD, PTD, t(Di)%*%t(PI)%*%X)
    IOC[i+nQTL, ] <- IOC0
  }

  PTX0 <- t(iqvector)%*%t(TPI)%*%X/sigma2
  PTX <- rbind(c(ind/(2*sigma2), PTX0), cbind(t(PTX0), t(X)%*%X))
  IOC[(ncd+nQTL+1):(ncd+nQTL+nx+1), ] <- cbind(t(IOC[1:(ncd+nQTL), (ncd+nQTL+1):(ncd+nQTL+nx+1)]), PTX)
  IOC <- IOC/sigma2

  for(i in 1:ncd){
    Di <- D.matrix[, i]
    for(j in i:ncd){
      Dj <- D.matrix[, j]
      IOM[nQTL+j, nQTL+i] <- IOM[nQTL+i, nQTL+j] <- (t(invector)%*%(T.matrix*TPI)%*%(Di*Dj)+sum((TPI%*%Dj)%*%t(TPI%*%Di))-
                                                       sum((TPI%*%Dj)*(TPI%*%Di)))/(sigma2)^2
    }
    IOM[nQTL+ncd+1, nQTL+i] <- IOM[nQTL+i, nQTL+ncd+1] <- (t(invector)%*%(S.matrix*TPI)%*%Di+sum((TPI%*%Di)%*%t(SPII))-
                                                             sum((TPI%*%Di)*SPII))/(sigma2)^2

    for(j in 1:nx){
      Xi <- X[, j]
      qTPIX <- t(iqvector)%*%t(TPI)*Xi
      IOM[nQTL+ncd+1+j, nQTL+i] <- IOM[nQTL+i, nQTL+ncd+1+j] <- (t(Xi)%*%(T.matrix*TPI)%*%Di+sum((TPI%*%Di)%*%qTPIX)-
                                                                   sum((TPI%*%Di)*t(qTPIX)))/(sigma2)^2
    }
  }

  IOM[nQTL+ncd+1, nQTL+ncd+1] <- (t(invector)%*%(S.matrix*S.matrix*PI)%*%iqvector+sum(SPII%*%t(SPII))-
                                    sum(SPII*SPII))/(sigma2)^2

  for(i in 1:nx){
    Xi <- X[, i]
    qTPIX <- t(iqvector)%*%t(TPI)*Xi
    IOM[nQTL+ncd+1, nQTL+ncd+1+i] <- IOM[nQTL+ncd+1+i, nQTL+ncd+1] <- (t(iqvector)%*%t(S.matrix*TPI)%*%Xi+
                                                                         sum(SPII%*%qTPIX)-sum(SPII*t(qTPIX)))/(sigma2)^2
    for(j in i:nx){
      Xj <- X[, j]
      qTPIXj <- t(iqvector)%*%t(TPI)*Xj
      IOM[nQTL+ncd+1+j, nQTL+ncd+1+i] <- IOM[nQTL+ncd+1+i, nQTL+ncd+1+j] <- (t(iqvector)%*%t(T.matrix*TPI)%*%(Xi*Xj)+
                                                                               sum(t(qTPIXj)%*%qTPIX)-sum(qTPIXj*qTPIX))/(sigma2)^2
    }
  }

  IOI <- IOC-IOM
  avc.matrix <- solve(IOI)
  EMvar <- diag(avc.matrix)
  if(sum(EMvar < 0) > 0){
    EMvar[EMvar < 0] <- 0
    warning("The result contains parameters whose approximate variance cannot be calculated.")
  }

  if(var.pos){
    nvar <- c(paste("QTL", 1:nQTL, sep = ""), colnames(D.matrix), "residual.var", colnames(X))
  } else {
    nvar <- c(colnames(D.matrix), "residual.var", colnames(X))
  }
  names(EMvar) <- colnames(avc.matrix) <- row.names(avc.matrix) <- nvar

  result <- EM
  result[[length(EM)+1]] <- avc.matrix
  result[[length(EM)+2]] <- EMvar
  names(result) <- c(names(EM), "avc.matrix", "EMvar")

  return(result)
}