R/runMod1G.R

In espenmei/Bio1G: Fits biometric models in OpenMx

#' Run a biometric model based on 1 generation of anchestors
#' @param vrs string vector with name of variables
#' @param Alist list of OpenMx components for additive genetic model
#' @param Tlist list of OpenMx components for twin/pregnancy environmental model
#' @param Elist list of OpenMx components for unique environmental model
#' @param D_M dataframe with mother identifiers
#' @param D_F dataframe with father identifiers
#' @param D_P dataframe with pregnancy identifiers and mz dummy
#' @param D_C dataframe with pregnancy identifiers, mz dummy, covariates and variables to be analyzed
#' @export
#' @examples
#' library(OpenMx)
#' 
#' data(D_5var) # Load example dataframe
#' 
#' nv = 5
#' vrs = paste0("Y", 1:nv)
#'
#' D_C = D_5var[, c(vrs, "Pid", "mz", "female")]
#' D_P = D_5var %>%
#' group_by(Pid) %>%
#' summarise(Mid = unique(Mid),
#'           Fid = unique(Fid),
#'           mz = unique(mz)) %>%
#' ungroup() %>%
#' data.frame()
#' D_M = data.frame(Mid = unique(D_P$Mid))
#' D_F = data.frame(Fid = unique(D_P$Fid))
#'
#'# Specify an "independent pathway" model
#'
#' # Additive genetic model
#' lCA = mxMatrix("Full", nv, 1, T, 0.5, labF("lCA", nv, 1), name = "lCA")
#' lSA = mxMatrix("Diag", nv, nv, T, 0.5, labD("lSA", nv), name = "lSA")
#' VA = mxAlgebra(lCA%*%t(lCA) + lSA%*%t(lSA), name = "VA")
#' # Twin specific environmental model
#' lCT = mxMatrix("Full", nv, 1, T, 0.5, labF("lCT", nv, 1), name = "lCT")
#' lST = mxMatrix("Diag", nv, nv, T, 0.5, labD("lST", nv), name = "lST")
#' VT = mxAlgebra(lCT%*%t(lCT) + lST%*%t(lST), name = "VT")
#' # Unique environmental model
#' lCE = mxMatrix("Full", nv, 1, T, 0.5, labF("lCE", nv, 1), name = "lCE")
#' lSE = mxMatrix("Diag", nv, nv, T, 0.5, labD("lSE", nv), name = "lSE")
#' VE = mxAlgebra(lCE%*%t(lCE) + lSE%*%t(lSE), name = "VE")
#' # Mean model
#' B = mxMatrix("Full", nv, 2, T, 0, labF("b", nv, 2), name = "B")
#' X = mxMatrix("Full", 2, 1, F, 1, c("x1", "data.female"), name = "X")
#' My = mxAlgebra(B%*%X, name = "My")
#'
#' Apars = list(lCA, lSA, VA)
#' Tpars = list(lCT, lST, VT)
#' Epars = list(lCE, lSE, VE)
#' Mpars = list(B, X, My)
#' fit = runMod1G(vrs, Apars, Tpars, Epars, Mpars, D_M, D_F, D_P, D_C)
#' summary(fit)
runMod1G = function(vrs = NULL,
                    Alist = NULL, Tlist = NULL, Elist = NULL, Mlist = NULL,
                    D_M = NULL, D_F = NULL, D_P = NULL, D_C = NULL) {
  
  if(is.null(vrs)) {
    stop("Variable names must be supplied")
  }
  
  nv = length(vrs)
  
  if(is.null(D_M) | is.null(D_F) | is.null(D_P) | is.null(D_C)) {
    stop("A dataframe must be supplied for each level")
  }
  
  if(is.null(Alist) & is.null(Tlist) & is.null(Elist)) {
    stop("Some kind of model must be specified")
  }
  
  if(is.null(Alist)) {
    Alist = mxMatrix("Zero", nv, nv, name = "AV")
  }
  
  if(is.null(Tlist)) {
    Tlist = mxMatrix("Zero", nv, nv, name = "TV")
  }
  
  if(is.null(Elist)) {
    Elist = mxMatrix("Zero", nv, nv, name = "EV")
  }
  
  if(is.null(Mlist)) {
    Mlist = mxMatrix("Zero", nv, nv, name = "My")
  }
  
  fit = Mod1G(nv, vrs,
              Alist, Tlist, Elist, Mlist,
              D_M, D_F, D_P, D_C)
  
  return(fit)
  
}