grmsem.fit {grmsem}

grmsem model fitting function

Description

This function fits a grmsem model.

Usage

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
)

Arguments

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.

Details

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.

Value

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)

Examples

#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)