grmsem.fit: grmsem model fitting function
In grmsem: Genetic-Relationship-Matrix Structural Equation Modelling (GRMSEM)

Description Usage Arguments Details Value Examples

This function fits a grmsem model.

grmsem.fit(
  ph,
  G,
  A.v.free = NULL,
  E.v.free = NULL,
  A.v.startval = NULL,
  E.v.startval = NULL,
  LogL = FALSE,
  estSE = FALSE,
  cores = 1,
  model = "Cholesky",
  compl.ph = FALSE,
  printest = FALSE,
  cluster = "PSOCK",
  optim = "optim",
  verbose = FALSE
)

`ph`	phenotype file as R dataframe, even for single vectors (columns: k phenotypes, rows: ni individuals in same order as G). No default.
`G`	GRM matrix as provided by the grm.input or grm.bin.input.R function. Use the same order of individuals as in ph. No default.
`A.v.free`	vector of free parameters for genetic factor loadings (free:1, not-free:0). Default NULL, all parameters are estimated.
`E.v.free`	vector of free parameters for residual factor loadings (free:1, not-free:0). Default NULL, all parameters are estimated.
`A.v.startval`	vector of starting values for genetic factor loadings. Default NULL, all starting values are generated.
`E.v.startval`	vector of starting values for residual factor loadings. Default NULL, all starting values are generated.
`LogL`	estimation of the loglikelihood using optim BFGS (TRUE/FALSE). Default FALSE.
`estSE`	estimation of standard errors by recalculating the Hessian matrix. Default FALSE.
`cores`	number of cores for multi-threaded calculations (numeric). Default 1.
`model`	grmsem model selection. Options: "Cholesky", "IP", "IPC", "DS". Default "Cholesky".
`compl.ph`	listwise complete observations across all phenotypes (all NA are excluded). Default FALSE.
`printest`	print output of the model.fit function including estimates (printest.txt) that can be used as starting values. Default FALSE.
`cluster`	cluster type. Options: "PSOCK", "FORK". Default "PSOCK".
`optim`	optimisation function from stats or optimParallel libraries. Options: "optim", "optimParallel". Default "optim".
`verbose`	additional model fit information: (i) phenotype vector, (ii) n of GRM and corresponding I matrix when data are missing, (iii) Hessian if `estSE` `TRUE`. Default FALSE.

grmsem models estimate genetic (A) and residual (E) variance/covariance of quantitative traits (AE model), where E in GRM-based methods can capture both, shared and unique residual influences. The estimation of parameters and their SEs is performed with the function grmsem.fit. Specifically, the loglikelihood is estimated with stats::optim and the BFGS (Broyden-Fletcher-Goldfarb-Shannon) approach and the variance/covariance matrix of estimated parameters with numDeriv::hessian. The statistical significance of estimated parameters is assessed using a Wald test, assuming multivariate normality.

grmsem.fit allows fitting different models describing the underlying multivariate genetic architecture of quantitative traits, as captured by a genetic-relationship-matrix (GRM), using structural equation modelling (SEM) techniques and a maximum likelihood approach. The user can fit multiple predefined model structures to the data. A Cholesky decomposition, Independent Pathway, and hybrid Independent Pathway/Cholesky models can be fitted by setting the model option to Cholesky, IP or IPC, respectively. In addition, the Cholesky model can be re-parametrised as Direct Symmetric model, estimating genetic and residual covariances directly, using the model option DS. Each model can be adapted by the user by setting free parameters (A.v.free and E.v.free options) and starting values (A.v.startval and E.v.startval options).

Input parameters are returned as model.in list object. Output from the maximum likelihood estimation is also given as list model.fit and the fitted grmsem model with estimated parameters and SEs is returned as dataframe model.out. The returned grmsem.fit object can be used to estimate genetic and residual covariance and correlations (grmsem.var function), bivariate heritabilities (grmsem.biher function), and factorial co-heritabilities and co-environmentalities (grmsem.fcoher function). All estimated parameters of the fitted grmsem model can also be standardised using the function grmsem.stpar.

Listwise complete observations can be selected with the option compl.ph=TRUE. Otherwise, grmsem.fit fits, like GREML, all available data to the model with the default option compl.ph FALSE. Using the option LogL=FALSE, the user can check the model input parameters and formula without a maximum likelihood estimation. Using the option estSE=FALSE, the user can carry out a maximum likelihood estimation without the estimation of SEs for estimated parameters that require calculating the Hessian. Note that grmsem.fit should preferably be run in parallel, by setting the cores option to the required number of cores.

When grmsem.fit is called with LogL TRUE, the user will see also the iterations of the stats::optim loglikelihood estimation, which are not included in the exported grmsem.fit list object.

grmsem.fit returns a grmsem.fit list object consisting of:

`model.in`	list of input parameters
`formula`	matrix of the model specification
`model.fit`	list output of the maximum likelihood estimation, if `LogL` `TRUE`
`model.out`	dataframe of fitted grmsem model with estimated parameters and SEs, if `estSE` `TRUE`
`VCOV`	variance/covariance matrix
`k`	number of phenotypes
`n`	total number of observations across all phenotypes
`n.obs`	number of observations per phenotype
`n.ind`	number of individuals with at least one phenotype
`model`	type of grmsem model
`ph.nms`	vector of phenotype names

model.in list of input parameters:

`part`	a - genetic, e - residual parameters
`label`	parameter label
`value`	starting values
`freepar`	free parameters

model.fit list output of the maximum likelihood estimation:

`optimisation`	output via optim()
`estimates`	estimated parameters: factor loadings for 'Cholesky', 'IP' and 'IPC' models, but variance components for 'DS'
`LL`	loglikelihood
`calls`	optim() calls
`convergence`	optim() convergence
`message`	optim() message

model.out data.frame of fitted grmsem model:

`label`	parameter label
`estimates`	estimated parameters
`gradient`	gradient
`se`	SE
`Z`	Z (Wald)
`p`	p (Wald)

#Set up a Cholesky model: Model formula and total number of parameters
#ph.small is a trivariate phenotype file for 100 individuals in the same order as G.small
#nrow = 100, ncol = 3
#(runtime should be less than one minute)
out <- grmsem.fit(ph.small, G.small, LogL = FALSE, estSE = FALSE)

#Run a Cholesky model:
out <- grmsem.fit(ph.small, G.small, LogL = TRUE, estSE = TRUE)