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
)phenotype file as R dataframe, even for single vectors (columns: k phenotypes, rows: ni individuals in same order as G). No default.
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.
vector of free parameters for genetic factor loadings (free:1, not-free:0). Default NULL, all parameters are estimated.
vector of free parameters for residual factor loadings (free:1, not-free:0). Default NULL, all parameters are estimated.
vector of starting values for genetic factor loadings. Default NULL, all starting values are generated.
vector of starting values for residual factor loadings. Default NULL, all starting values are generated.
estimation of the loglikelihood using optim BFGS (TRUE/FALSE). Default FALSE.
estimation of standard errors by recalculating the Hessian matrix. Default FALSE.
number of cores for multi-threaded calculations (numeric). Default 1.
grmsem model selection. Options: "Cholesky", "IP", "IPC", "DS". Default "Cholesky".
listwise complete observations across all phenotypes (all NA are excluded). Default FALSE.
print output of the model.fit function including estimates (printest.txt) that can be used as starting values. Default FALSE.
cluster type. Options: "PSOCK", "FORK". Default "PSOCK".
optimisation function from stats or optimParallel libraries. Options: "optim", "optimParallel". Default "optim".
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.fit returns a grmsem.fit list object consisting of:
list of input parameters
matrix of the model specification
list output of the maximum likelihood estimation, if LogL TRUE
dataframe of fitted grmsem model with estimated parameters and SEs, if estSE TRUE
variance/covariance matrix
number of phenotypes
total number of observations across all phenotypes
number of observations per phenotype
number of individuals with at least one phenotype
type of grmsem model
vector of phenotype names
model.in list of input parameters:
a - genetic, e - residual parameters
parameter label
starting values
free parameters
model.fit list output of the maximum likelihood estimation:
output via optim()
estimated parameters: factor loadings for `Cholesky`, `IP` and `IPC` models, but variance components for `DS`
loglikelihood
optim() calls
optim() convergence
optim() message
model.out data.frame of fitted grmsem model:
parameter label
estimated parameters
gradient
SE
Z (Wald)
p (Wald)
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.
# NOT RUN {
#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)
# }
# NOT RUN {
out <- grmsem.fit(ph.small, G.small, LogL = FALSE, estSE = FALSE)
# }
# NOT RUN {
#Run a Cholesky model:
# }
# NOT RUN {
out <- grmsem.fit(ph.small, G.small, LogL = TRUE, estSE = TRUE)
# }
Run the code above in your browser using DataLab