Last data update: 2014.03.03

 R: Create MxExpectationGREML Object
mxExpectationGREMLR Documentation

Create MxExpectationGREML Object

Description

This function creates a new MxExpectationGREML object.

Usage

mxExpectationGREML(V, yvars=character(0), Xvars=list(), addOnes=TRUE, blockByPheno=TRUE, 
                  staggerZeroes=TRUE, dataset.is.yX=FALSE, casesToDropFromV=integer(0))

Arguments

V

Character string; the name of the MxAlgebra or MxMatrix to serve as the 'V' matrix (the model-expected covariance matrix). Internally, the 'V' matrix is assumed to be symmetric, and its elements above the main diagonal are ignored.

yvars, Xvars, addOnes, blockByPheno, staggerZeroes

Passed to mxGREMLDataHandler().

dataset.is.yX

Logical; defaults to FALSE. If TRUE, then the first column of the raw dataset is taken as-is to be the 'y' phenotype vector, and the remaining columns are taken as-is to be the 'X' matrix of covariates. In this case, mxGREMLDataHandler() is never internally called at runtime, and all other arguments besides V and casesToDropFromV are ignored.

casesToDropFromV

Integer vector. Its elements are the numbers of the rows and columns of covariance matrix 'V' to be dropped at runtime, usually because they correspond to rows of 'y' or 'X' that contained missing observations. By default, no cases are dropped from 'V.' Ignored unless dataset.is.yX=TRUE.

Details

"GREML" stands for "genomic-relatedness-matrix restricted maximum-likelihood." In the strictest sense of the term, it refers to genetic variance-component estimation from matrices of subjects' pairwise degree of genetic relatedness, as calculated from genome-wide marker data. It is from this original motivation that some of the terminology originates, such as calling 'y' the "phenotype" vector. However, OpenMx's implementation of GREML is applicable for analyses from any subject-matter domain, and in which the following assumptions are reasonable:

  1. Conditional on 'X' (the covariates), the phenotype vector (response variable) 'y' is a single realization from a multivariate-normal distribution having (in general) a dense covariance matrix, 'V.'

  2. The parameters of the covariance matrix, such as variance components, are of primary interest.

  3. The random effects are normally distributed.

  4. Weighted least-squares regression, using the inverse of 'V' as a weight matrix, is an adequate model for the phenotypic means. Note that the regression coefficients are not actually free parameters to be numerically optimized.

Computationally, the chief distinguishing feature of an OpenMx GREML analysis is that the phenotype vector, 'y,' is a single realization of a random vector that, in general, cannot be partitioned into independent subvectors. For this reason, definition variables are not compatible (and should be unnecessary with) GREML expectation. GREML expectation can still be used if the covariance matrix is sparse, but as of this writing, OpenMx does not take advantage of the sparseness to improved performance. Because of the limitations of restricted maximum likelihood, GREML expectation is presently incompatible with ordinal variables.

Value

Returns a new object of class MxExpectationGREML.

References

One of the first uses of the acronym "GREML":
Benjamin DJ, Cesarini D, van der Loos MJHM, Dawes CT, Koellinger PD, et al. (2012) The genetic architecture of economic and political preferences. Proceedings of the National Academy of Sciences 109: 8026-8031. doi: 10.1073/pnas.1120666109

The OpenMx User's guide can be found at http://openmx.psyc.virginia.edu/documentation.

See Also

See MxExpectationGREML for the S4 class created by mxExpectationGREML(). More information about the OpenMx package may be found here.

Examples

dat <- cbind(rnorm(100),rep(1,100))
colnames(dat) <- c("y","x")

ge <- mxExpectationGREML(V="V",yvars="y",Xvars=list("X"),addOnes=FALSE)
gff <- mxFitFunctionGREML(dV=c(ve="I"))
plan <- mxComputeSequence(freeSet=c("Ve"),steps=list(
  mxComputeNewtonRaphson(fitfunction="fitfunction"),
  mxComputeOnce('fitfunction',
    c('fit','gradient','hessian','ihessian')),
  mxComputeStandardError(),
  mxComputeReportDeriv(),
  mxComputeReportExpectation()
))

testmod <- mxModel(
  "GREMLtest",
  mxData(observed = dat, type="raw"),
  mxMatrix(type = "Full", nrow = 1, ncol=1, free=TRUE,
    values = 1, labels = "ve", lbound = 0.0001, name = "Ve"),
  mxMatrix("Iden",nrow=100,name="I",condenseSlots=TRUE),
  mxAlgebra(I %x% Ve,name="V"),
  ge,
  gff,
  plan
)
str(testmod)

Results


R version 3.3.1 (2016-06-21) -- "Bug in Your Hair"
Copyright (C) 2016 The R Foundation for Statistical Computing
Platform: x86_64-pc-linux-gnu (64-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> library(OpenMx)
Loading required package: digest
Loading required package: MASS
Loading required package: Matrix
Loading required package: Rcpp
Loading required package: parallel

Attaching package: 'OpenMx'

The following objects are masked from 'package:Matrix':

    %&%, expm

> png(filename="/home/ddbj/snapshot/RGM3/R_CC/result/OpenMx/mxExpectationGREML.Rd_%03d_medium.png", width=480, height=480)
> ### Name: mxExpectationGREML
> ### Title: Create MxExpectationGREML Object
> ### Aliases: mxExpectationGREML
> 
> ### ** Examples
> 
> dat <- cbind(rnorm(100),rep(1,100))
> colnames(dat) <- c("y","x")
> 
> ge <- mxExpectationGREML(V="V",yvars="y",Xvars=list("X"),addOnes=FALSE)
> gff <- mxFitFunctionGREML(dV=c(ve="I"))
> plan <- mxComputeSequence(freeSet=c("Ve"),steps=list(
+   mxComputeNewtonRaphson(fitfunction="fitfunction"),
+   mxComputeOnce('fitfunction',
+     c('fit','gradient','hessian','ihessian')),
+   mxComputeStandardError(),
+   mxComputeReportDeriv(),
+   mxComputeReportExpectation()
+ ))
> 
> testmod <- mxModel(
+   "GREMLtest",
+   mxData(observed = dat, type="raw"),
+   mxMatrix(type = "Full", nrow = 1, ncol=1, free=TRUE,
+     values = 1, labels = "ve", lbound = 0.0001, name = "Ve"),
+   mxMatrix("Iden",nrow=100,name="I",condenseSlots=TRUE),
+   mxAlgebra(I %x% Ve,name="V"),
+   ge,
+   gff,
+   plan
+ )
> str(testmod)
Formal class 'MxModel' [package "OpenMx"] with 21 slots
  ..@ name             : chr "GREMLtest"
  ..@ matrices         :List of 2
  .. ..$ Ve:Formal class 'FullMatrix' [package "OpenMx"] with 13 slots
  .. .. .. ..@ name           : chr "Ve"
  .. .. .. ..@ values         : num [1, 1] 1
  .. .. .. ..@ labels         : chr [1, 1] "ve"
  .. .. .. ..@ free           : logi [1, 1] TRUE
  .. .. .. ..@ lbound         : num [1, 1] 1e-04
  .. .. .. ..@ ubound         : num [1, 1] NA
  .. .. .. ..@ .squareBrackets: num[0 , 0 ] 
  .. .. .. ..@ .persist       : logi TRUE
  .. .. .. ..@ .condenseSlots : logi FALSE
  .. .. .. ..@ display        : chr(0) 
  .. .. .. ..@ dependencies   : int(0) 
  .. .. .. ..@ joinModel      : chr NA
  .. .. .. ..@ joinKey        : chr NA
  .. ..$ I :Formal class 'IdenMatrix' [package "OpenMx"] with 13 slots
  .. .. .. ..@ name           : chr "I"
  .. .. .. ..@ values         : num [1:100, 1:100] 1 0 0 0 0 0 0 0 0 0 ...
  .. .. .. ..@ labels         : chr [1, 1] NA
  .. .. .. ..@ free           : logi [1, 1] FALSE
  .. .. .. ..@ lbound         : num [1, 1] NA
  .. .. .. ..@ ubound         : num [1, 1] NA
  .. .. .. ..@ .squareBrackets: num[0 , 0 ] 
  .. .. .. ..@ .persist       : logi TRUE
  .. .. .. ..@ .condenseSlots : logi TRUE
  .. .. .. ..@ display        : chr(0) 
  .. .. .. ..@ dependencies   : int(0) 
  .. .. .. ..@ joinModel      : chr NA
  .. .. .. ..@ joinKey        : chr NA
  ..@ algebras         :List of 1
  .. ..$ V:Formal class 'MxAlgebra' [package "OpenMx"] with 7 slots
  .. .. .. ..@ formula  : language I %x% Ve
  .. .. .. ..@ name     : chr "V"
  .. .. .. ..@ fixed    : logi FALSE
  .. .. .. ..@ .dimnames: NULL
  .. .. .. ..@ result   : num[0 , 0 ] 
  .. .. .. ..@ joinModel: chr NA
  .. .. .. ..@ joinKey  : chr NA
  ..@ constraints      : list()
  ..@ intervals        : list()
  ..@ latentVars       : chr(0) 
  ..@ manifestVars     : chr(0) 
  ..@ data             :Formal class 'MxDataStatic' [package "OpenMx"] with 17 slots
  .. .. ..@ observed            : num [1:100, 1:2] 0.433 -1.698 -0.694 -1.097 0.298 ...
  .. .. .. ..- attr(*, "dimnames")=List of 2
  .. .. .. .. ..$ : NULL
  .. .. .. .. ..$ : chr [1:2] "y" "x"
  .. .. ..@ means               : num [1, 1] NA
  .. .. ..@ type                : chr "raw"
  .. .. ..@ numObs              : num 100
  .. .. ..@ acov                : num [1, 1] NA
  .. .. ..@ fullWeight          : num [1, 1] NA
  .. .. ..@ thresholds          : num [1, 1] NA
  .. .. ..@ thresholdColumns    : int(0) 
  .. .. ..@ thresholdLevels     : int(0) 
  .. .. ..@ indexVector         : int(0) 
  .. .. ..@ identicalDefVars    : int(0) 
  .. .. ..@ identicalMissingness: int(0) 
  .. .. ..@ identicalRows       : int(0) 
  .. .. ..@ .isSorted           : logi FALSE
  .. .. ..@ .needSort           : logi TRUE
  .. .. ..@ primaryKey          : chr NA
  .. .. ..@ name                : chr "data"
  ..@ submodels        : list()
  ..@ expectation      :Formal class 'MxExpectationGREML' [package "OpenMx"] with 22 slots
  .. .. ..@ V            : chr "V"
  .. .. ..@ yvars        : chr "y"
  .. .. ..@ Xvars        :List of 1
  .. .. .. ..$ : chr "X"
  .. .. ..@ addOnes      : logi FALSE
  .. .. ..@ blockByPheno : logi TRUE
  .. .. ..@ staggerZeroes: logi TRUE
  .. .. ..@ dataset.is.yX: logi FALSE
  .. .. ..@ X            : num [1, 1] NA
  .. .. ..@ y            : NULL
  .. .. ..@ yXcolnames   : chr(0) 
  .. .. ..@ casesToDrop  : int(0) 
  .. .. ..@ b            : num[0 , 0 ] 
  .. .. ..@ bcov         : num[0 , 0 ] 
  .. .. ..@ numFixEff    : int(0) 
  .. .. ..@ dims         : chr "foo"
  .. .. ..@ numStats     : num(0) 
  .. .. ..@ dataColumns  : num(0) 
  .. .. ..@ name         : chr "expectation"
  .. .. ..@ data         : int NA
  .. .. ..@ .runDims     : chr(0) 
  .. .. ..@ output       : list()
  .. .. ..@ debug        : list()
  ..@ fitfunction      :Formal class 'MxFitFunctionGREML' [package "OpenMx"] with 14 slots
  .. .. ..@ dV            : Named chr "I"
  .. .. .. ..- attr(*, "names")= chr "ve"
  .. .. ..@ dVnames       : chr "ve"
  .. .. ..@ MLfit         : num 0
  .. .. ..@ numObs        : int 0
  .. .. ..@ aug           : chr(0) 
  .. .. ..@ augGrad       : chr(0) 
  .. .. ..@ augHess       : chr(0) 
  .. .. ..@ info          : list()
  .. .. ..@ dependencies  : int(0) 
  .. .. ..@ expectation   : int(0) 
  .. .. ..@ vector        : logi FALSE
  .. .. ..@ rowDiagnostics: logi(0) 
  .. .. ..@ result        : num[0 , 0 ] 
  .. .. ..@ name          : chr "fitfunction"
  ..@ compute          :Formal class 'MxComputeSequence' [package "OpenMx"] with 8 slots
  .. .. ..@ independent: logi FALSE
  .. .. ..@ steps      :List of 5
  .. .. .. ..$ :Formal class 'MxComputeNewtonRaphson' [package "OpenMx"] with 10 slots
  .. .. .. .. .. ..@ fitfunction: chr "fitfunction"
  .. .. .. .. .. ..@ maxIter    : int 100
  .. .. .. .. .. ..@ tolerance  : num 1e-12
  .. .. .. .. .. ..@ verbose    : int 0
  .. .. .. .. .. ..@ id         : int(0) 
  .. .. .. .. .. ..@ freeSet    : chr NA
  .. .. .. .. .. ..@ output     : list()
  .. .. .. .. .. ..@ debug      : list()
  .. .. .. .. .. ..@ .persist   : logi TRUE
  .. .. .. .. .. ..@ name       : chr "compute"
  .. .. .. ..$ :Formal class 'MxComputeOnce' [package "OpenMx"] with 11 slots
  .. .. .. .. .. ..@ from       : chr "fitfunction"
  .. .. .. .. .. ..@ what       : chr [1:4] "fit" "gradient" "hessian" "ihessian"
  .. .. .. .. .. ..@ how        : NULL
  .. .. .. .. .. ..@ verbose    : int 0
  .. .. .. .. .. ..@ .is.bestfit: logi FALSE
  .. .. .. .. .. ..@ id         : int(0) 
  .. .. .. .. .. ..@ freeSet    : chr NA
  .. .. .. .. .. ..@ output     : list()
  .. .. .. .. .. ..@ debug      : list()
  .. .. .. .. .. ..@ .persist   : logi TRUE
  .. .. .. .. .. ..@ name       : chr "compute"
  .. .. .. ..$ :Formal class 'MxComputeStandardError' [package "OpenMx"] with 6 slots
  .. .. .. .. .. ..@ id      : int(0) 
  .. .. .. .. .. ..@ freeSet : chr NA
  .. .. .. .. .. ..@ output  : list()
  .. .. .. .. .. ..@ debug   : list()
  .. .. .. .. .. ..@ .persist: logi TRUE
  .. .. .. .. .. ..@ name    : chr "compute"
  .. .. .. ..$ :Formal class 'MxComputeReportDeriv' [package "OpenMx"] with 6 slots
  .. .. .. .. .. ..@ id      : int(0) 
  .. .. .. .. .. ..@ freeSet : chr NA
  .. .. .. .. .. ..@ output  : list()
  .. .. .. .. .. ..@ debug   : list()
  .. .. .. .. .. ..@ .persist: logi TRUE
  .. .. .. .. .. ..@ name    : chr "compute"
  .. .. .. ..$ :Formal class 'MxComputeReportExpectation' [package "OpenMx"] with 6 slots
  .. .. .. .. .. ..@ id      : int(0) 
  .. .. .. .. .. ..@ freeSet : chr NA
  .. .. .. .. .. ..@ output  : list()
  .. .. .. .. .. ..@ debug   : list()
  .. .. .. .. .. ..@ .persist: logi TRUE
  .. .. .. .. .. ..@ name    : chr "compute"
  .. .. ..@ id         : int(0) 
  .. .. ..@ freeSet    : chr "Ve"
  .. .. ..@ output     : list()
  .. .. ..@ debug      : list()
  .. .. ..@ .persist   : logi TRUE
  .. .. ..@ name       : chr "compute"
  ..@ independent      : logi FALSE
  ..@ options          : list()
  ..@ output           : list()
  ..@ runstate         : list()
  ..@ .newobjects      : logi FALSE
  ..@ .resetdata       : logi FALSE
  ..@ .wasRun          : logi FALSE
  ..@ .modifiedSinceRun: logi TRUE
  ..@ .version         :Classes 'package_version', 'numeric_version'  hidden list of 1
  .. ..$ : int [1:3] 2 6 7
> 
> 
> 
> 
> 
> dev.off()
null device 
          1 
>


Created & Maintained by Osamu Ogasawara (osamu.ogasawara@gmail.com) and IMSBIO