M 
the number of dependent variables Y (outcomes); equivalently, the number of equations in the system

method 
the PMT method used to generate all continuous variables, including independent variables (covariates), error terms, and random effects;
"Fleishman" uses Fleishman's thirdorder polynomial transformation and "Polynomial" uses Headrick's fifthorder transformation

error_type 
"non_mix" if all error terms have continuous nonmixture distributions, "mix" if all error terms have continuous mixture distributions

means 
if no random effects, a list of length M where means[[p]] contains a vector of means for the continuous independent variables
in equation p with nonmixture (X_{cont}) or mixture (X_{mix}) distributions and for the error terms (E);
order in vector is X_{cont}, X_{mix}, E
if there are random effects, a list of length M + 1 if the effects are the same across equations or 2 * M if they differ;
where means[M + 1] or means[(M + 1):(2 * M)] are vectors of means for all random effects with continuous nonmixture or mixture distributions;
order in vector is 1st random intercept U_0 (if rand.int != "none" ), 2nd random time slope U_1 (if rand.tsl != "none" ),
3rd other random slopes with nonmixture distributions U_{cont}, 4th other random slopes with mixture distributions U_{mix}

vars 
a list of same length and order as means containing vectors of variances for the continuous variables, error terms, and any random effects

skews 
if no random effects, a list of length M where skews[[p]] contains a vector of skew values for the continuous independent variables
in equation p with nonmixture (X_{cont}) distributions and for E if error_type = "non_mix" ; order in vector is X_{cont}, E
if there are random effects, a list of length M + 1 if the effects are the same across equations or 2 * M if they differ;
where skews[M + 1] or skews[(M + 1):(2 * M)] are vectors of skew values for all random effects with continuous nonmixture distributions;
order in vector is 1st random intercept U_0 (if rand.int = "non_mix" ), 2nd random time slope U_1 (if rand.tsl = "non_mix" ),
3rd other random slopes with nonmixture distributions U_{cont}

skurts 
a list of same length and order as skews containing vectors of standardized kurtoses (kurtosis  3) for the continuous variables,
error terms, and any random effects with nonmixture distributions

fifths 
a list of same length and order as skews containing vectors of standardized fifth cumulants for the continuous variables,
error terms, and any random effects with nonmixture distributions; not necessary for method = "Fleishman"

sixths 
a list of same length and order as skews containing vectors of standardized sixth cumulants for the continuous variables,
error terms, and any random effects with nonmixture distributions; not necessary for method = "Fleishman"

Six 
a list of length M , M + 1 , or 2 * M , where Six[1:M] are for X_{cont}, E (if error_type = "non_mix" ) and
Six[M + 1] or Six[(M + 1):(2 * M)] are for nonmixture U;
if error_type = "mix" and there are only random effects (i.e., length(corr.x) = 0 ), use Six[1:M] = rep(list(NULL), M) so that
Six[M + 1] or Six[(M + 1):(2 * M)] describes the nonmixture U;
Six[[p]][[j]] is a vector of sixth cumulant correction values to aid in finding a valid PDF for X_{cont(pj)}, the
jth continuous nonmixture covariate for outcome Y_p; the last vector in Six[[p]] is for E_p (if error_type = "non_mix" );
use Six[[p]][[j]] = NULL if no correction desired for X_{cont(pj)};
use Six[[p]] = NULL if no correction desired for any continuous nonmixture covariate or error term in equation p
Six[[M + p]][[j]] is a vector of sixth cumulant correction values to aid in finding a valid PDF for U_{(pj)}, the
jth nonmixture random effect for outcome Y_p; use Six[[M + p]][[j]] = NULL if no correction desired for U_{(pj)};
use Six[[M + p]] = NULL if no correction desired for any continuous nonmixture random effect in equation p
keep Six = list() if no corrections desired for all equations or if method = "Fleishman"

mix_pis 
list of length M , M + 1 or 2 * M , where mix_pis[1:M] are for X_{cont}, E (if error_type = "mix" ) and
mix_pis[M + 1] or mix_pis[(M + 1):(2 * M)] are for mixture U; use mix_pis[[p]] = NULL if equation p has no continuous mixture terms
if error_type = "non_mix" and there are only random effects (i.e., length(corr.x) = 0 ), use mix_pis[1:M] = NULL so that
mix_pis[M + 1] or mix_pis[(M + 1):(2 * M)] describes the mixture U;
mix_pis[[p]][[j]] is a vector of mixing probabilities of the component distributions for X_{mix(pj)}, the jth mixture covariate for outcome Y_p;
the last vector in mix_pis[[p]] is for E_p (if error_type = "mix" ); components should be ordered as in corr.x
mix_pis[[M + p]][[j]] is a vector of mixing probabilities of the component distributions for U_{(pj)}, the jth random effect with a mixture distribution
for outcome Y_p; order is 1st random intercept (if rand.int = "mix" ), 2nd random time slope (if rand.tsl = "mix" ),
3rd other random slopes with mixture distributions; components should be ordered as in corr.u

mix_mus 
list of same length and order as mix_pis ;
mix_mus[[p]][[j]] is a vector of means of the component distributions for X_{mix(pj)},
the last vector in mix_mus[[p]] is for E_p (if error_type = "mix" )
mix_mus[[p]][[j]] is a vector of means of the component distributions for U_{mix(pj)}

mix_sigmas 
list of same length and order as mix_pis ;
mix_sigmas[[p]][[j]] is a vector of standard deviations of the component distributions for X_{mix(pj)},
the last vector in mix_sigmas[[p]] is for E_p (if error_type = "mix" )
mix_sigmas[[p]][[j]] is a vector of standard deviations of the component distributions for U_{mix(pj)}

mix_skews 
list of same length and order as mix_pis ;
mix_skews[[p]][[j]] is a vector of skew values of the component distributions for X_{mix(pj)},
the last vector in mix_skews[[p]] is for E_p (if error_type = "mix" )
mix_skews[[p]][[j]] is a vector of skew values of the component distributions for U_{mix(pj)}

mix_skurts 
list of same length and order as mix_pis ;
mix_skurts[[p]][[j]] is a vector of standardized kurtoses of the component distributions for X_{mix(pj)},
the last vector in mix_skurts[[p]] is for E_p (if error_type = "mix" )
mix_skurts[[p]][[j]] is a vector of standardized kurtoses of the component distributions for U_{mix(pj)}

mix_fifths 
list of same length and order as mix_pis ; not necessary for method = "Fleishman" ;
mix_fifths[[p]][[j]] is a vector of standardized fifth cumulants of the component distributions for X_{mix(pj)},
the last vector in mix_fifths[[p]] is for E_p (if error_type = "mix" )
mix_fifths[[p]][[j]] is a vector of standardized fifth cumulants of the component distributions for U_{mix(pj)}

mix_sixths 
list of same length and order as mix_pis ; not necessary for method = "Fleishman" ;
mix_sixths[[p]][[j]] is a vector of standardized sixth cumulants of the component distributions for X_{mix(pj)},
the last vector in mix_sixths[[p]] is for E_p (if error_type = "mix" )
mix_sixths[[p]][[j]] is a vector of standardized sixth cumulants of the component distributions for U_{mix(pj)}

mix_Six 
a list of same length and order as mix_pis ;
keep mix_Six = list() if no corrections desired for all equations or if method = "Fleishman"
pth component of mix_Six[1:M] is a list of length equal to the total number of component distributions for the X_{mix(p)}
and E_p (if error_type = "mix" );
mix_Six[[p]][[j]] is a vector of sixth cumulant corrections for the jth component distribution (i.e., if there are 2
continuous mixture independent variables for Y_p, where X_{mix(p1)} has 2 components and X_{mix(p2)} has 3
components, then length(mix_Six[[p]]) = 5 and mix_Six[[p]][[3]] would correspond to the 1st component of
X_{mix(p2)}); use mix_Six[[p]][[j]] = NULL if no correction desired for that component;
use mix_Six[[p]] = NULL if no correction desired for any component of X_{mix(p)} and E_p
qth component of mix_Six[M + 1] or mix_Six[(M + 1):(2 * M)] is a list of length equal to the total number of component distributions for
the U_{mix(q)}; mix_Six[[q]][[j]] is a vector of sixth cumulant corrections for the jth component distribution; use
mix_Six[[q]][[j]] = NULL if no correction desired for that component;
use mix_Six[[q]] = NULL if no correction desired for any component of U_{mix(q)}

marginal 
a list of length M , with the pth component a list of cumulative probabilities for the ordinal variables
associated with outcome Y_p (use marginal[[p]] = NULL if outcome Y_p has no ordinal variables);
marginal[[p]][[j]] is a vector of the cumulative probabilities defining the marginal distribution of X_{ord(pj)},
the jth ordinal variable for outcome Y_p; if the variable can take r values, the vector will contain r  1 probabilities
(the rth is assumed to be 1); for binary variables, the probability is the probability of the 1st category, which has the smaller support value;
length(marginal[[p]]) can differ across outcomes; the order should be the same as in corr.x

support 
a list of length M , with the pth component a list of support values for the ordinal variables associated
with outcome Y_p; use support[[p]] = NULL if outcome Y_p has no ordinal variables;
support[[p]][[j]] is a vector of the support values defining the marginal distribution of X_{ord(pj)},
the jth ordinal variable for outcome Y_p; if not provided, the default for r categories is 1, ..., r

lam 
list of length M , pth component a vector of lambda (means > 0) values for Poisson variables for outcome Y_p
(see stats::dpois ); order is 1st regular Poisson and 2nd zeroinflated Poisson; use lam[[p]] = NULL if outcome Y_p has no Poisson variables;
length(lam[[p]]) can differ across outcomes; the order should be the same as in corr.x

p_zip 
a list of vectors of probabilities of structural zeros (not including zeros from the Poisson distribution) for the
zeroinflated Poisson variables (see VGAM::dzipois ); if p_zip = 0, Y_{pois} has a regular Poisson
distribution; if p_zip is in (0, 1), Y_{pois} has a zeroinflated Poisson distribution;
if p_zip is in ((exp(lam)  1)^(1), 0) , Y_{pois} has a zerodeflated Poisson distribution and p_zip
is not a probability; if p_zip = (exp(lam)  1)^(1) , Y_{pois} has a positivePoisson distribution
(see VGAM::dpospois ); order is 1st regular Poisson and 2nd zeroinflated Poisson;
if a single number, all Poisson variables given this value; if a vector of length M , all Poisson variables in equation p
given p_zip[p] ; otherwise, missing values are set to 0 and ordered 1st

pois_eps 
list of length M , pth component a vector of length lam[[p]] containing cumulative probability truncation values
used to calculate intermediate correlations involving Poisson variables; order is 1st regular Poisson and 2nd zeroinflated Poisson;
if a single number, all Poisson variables given this value; if a vector of length M , all Poisson variables in equation p
given pois_eps[p] ; otherwise, missing values are set to 0.0001 and ordered 1st

size 
list of length M , pth component a vector of size parameters for the Negative Binomial variables for outcome Y_p
(see stats::dnbinom ); order is 1st regular NB and 2nd zeroinflated NB; use size[[p]] = NULL if outcome Y_p has no Negative Binomial variables;
length(size[[p]]) can differ across outcomes; the order should be the same as in corr.x

prob 
list of length M , pth component a vector of success probabilities for the Negative Binomial variables for outcome Y_p
(see stats::dnbinom ); order is 1st regular NB and 2nd zeroinflated NB; use prob[[p]] = NULL if outcome Y_p has no Negative Binomial variables;
length(prob[[p]]) can differ across outcomes; the order should be the same as in corr.x

mu 
list of length M , pth component a vector of mean values for the Negative Binomial variables for outcome Y_p
(see stats::dnbinom ); order is 1st regular NB and 2nd zeroinflated NB; use mu[[p]] = NULL if outcome Y_p has no Negative Binomial variables;
length(mu[[p]]) can differ across outcomes; the order should be the same as in corr.x ; for zeroinflated NB variables,
this refers to the mean of the NB distribution (see VGAM::dzinegbin )
(*Note: either prob or mu should be supplied for all Negative Binomial variables, not a mixture)

p_zinb 
a vector of probabilities of structural zeros (not including zeros from the NB distribution) for the zeroinflated NB variables
(see VGAM::dzinegbin ); if p_zinb = 0, Y_{nb} has a regular NB distribution;
if p_zinb is in (prob^size/(1  prob^size), 0) , Y_{nb} has a zerodeflated NB distribution and p_zinb
is not a probability; if p_zinb = prob^size/(1  prob^size) , Y_{nb} has a positiveNB distribution (see
VGAM::dposnegbin ); order is 1st regular NB and 2nd zeroinflated NB;
if a single number, all NB variables given this value; if a vector of length M , all NB variables in equation p
given p_zinb[p] ; otherwise, missing values are set to 0 and ordered 1st

nb_eps 
list of length M , pth component a vector of length size[[p]] containing cumulative probability truncation values
used to calculate intermediate correlations involving Negative Binomial variables; order is 1st regular NB and 2nd zeroinflated NB;
if a single number, all NB variables given this value; if a vector of length M , all NB variables in equation p
given nb_eps[p] ; otherwise, missing values are set to 0.0001 and ordered 1st

corr.x 
list of length M , each component a list of length M ; corr.x[[p]][[q]] is matrix of correlations
for independent variables in equations p (X_{(pj)} for outcome Y_p) and q (X_{(qj)} for outcome Y_q);
order: 1st ordinal (same order as in marginal ), 2nd continuous nonmixture (same order as in skews ),
3rd components of continuous mixture (same order as in mix_pis ), 4th regular Poisson, 5th zeroinflated Poisson (same order as in lam ),
6th regular NB, and 7th zeroinflated NB (same order as in size );
if p = q, corr.x[[p]][[q]] is a correlation matrix with nrow(corr.x[[p]][[q]]) = # X_{(pj)} for outcome Y_p;
if p != q, corr.x[[p]][[q]] is a nonsymmetric matrix of correlations where rows correspond to covariates for Y_p
so that nrow(corr.x[[p]][[q]]) = # X_{(pj)} for outcome Y_p and
columns correspond to covariates for Y_q so that ncol(corr.x[[p]][[q]]) = # X_{(qj)} for outcome Y_q;
use corr.x[[p]][[q]] = NULL if equation q has no X_{(qj)}; use corr.x[[p]] = NULL if equation p has no X_{(pj)}

corr.yx 
input for nonnormsys only; a list of length M , where the pth component is a 1 row matrix of correlations between Y_p and X_{(pj)};
if there are mixture variables and the betas are desired in terms of these (and not the components), then corr.yx
should be specified in terms of correlations between outcomes and nonmixture or mixture variables, and the number of columns of the matrices
of corr.yx should not match the dimensions of the matrices in corr.x ; if the betas are desired in terms of
the components, then corr.yx should be specified in terms of correlations between outcomes and nonmixture or components of
mixture variables, and the number of columns of the matrices of corr.yx should match the dimensions of the matrices in corr.x ;
use corr.yx[[p]] = NULL if equation p has no X_{(pj)}

corr.e 
correlation matrix for continuous nonmixture or components of mixture error terms

same.var 
either a vector or a matrix; if a vector, same.var includes column numbers of corr.x[[1]][[1]]
corresponding to independent variables that should be identical across equations; these terms must have the same indices for all
p = 1, ..., M ; i.e., if the 1st ordinal variable represents sex which should be the same for each equation, then
same.var[1] = 1 since ordinal variables are 1st in corr.x[[1]][[1]] and sex is the 1st ordinal variable, and
the 1st term for all other outcomes must also be sex;
if a matrix, columns 1 and 2 are outcome p and column index in corr.x[[p]][[p]] for 1st instance of variable,
columns 3 and 4 are outcome q and column index in corr.x[[q]][[q]] for subsequent instances of variable; i.e., if
1st term for all outcomes is sex and M = 3 , then same.var = matrix(c(1, 1, 2, 1, 1, 1, 3, 1), 2, 4, byrow = TRUE) ; the
independent variable index corresponds to ordinal, continuous nonmixture, component of continuous mixture, Poisson, or
NB variable

subj.var 
matrix where 1st column is outcome index (p = 1, ..., M ), 2nd column is independent variable index
corresponding to covariate which is a a subjectlevel term (not including time), including timevarying covariates;
the independent variable index corresponds to ordinal, continuous nonmixture, continuous mixture (not mixture component), Poisson, or
NB variable; assumes all other variables are grouplevel terms; these subjectlevel terms are used to form interactions with the group
level terms

int.var 
matrix where 1st column is outcome index (p = 1, ..., M ), 2nd and 3rd columns are indices
corresponding to independent variables to form interactions between; this includes all interactions that are not accounted for by
a subjectgroup level interaction (as indicated by subj.var ) or by a timecovariate interaction (as indicated by
tint.var ); ex: 1, 2, 3 indicates that for outcome 1, the 2nd and 3rd independent variables form an interaction term;
the independent variable index corresponds to ordinal, continuous nonmixture, continuous mixture (not mixture component), Poisson, or
NB variable

tint.var 
matrix where 1st column is outcome index (p = 1, ..., M ), 2nd column is index of
independent variable to form interaction with time; if tint.var = NULL or no X_{(pj)} are indicated for outcome Y_p,
this includes all grouplevel variables (variables not indicated as subjectlevel variables in subj.var ), else includes only
terms indicated by 2nd column (i.e., in order to include subjectlevel variables);
ex: 1, 1 indicates that for outcome 1, the 1st independent variable has an interaction with time;
the independent variable index corresponds to ordinal, continuous nonmixture, continuous mixture (not mixture component), Poisson, or
NB variable

betas.0 
vector of length M containing intercepts, if NULL all set equal to 0; if length 1, all intercepts set to
betas.0

betas 
list of length M , pth component a vector of coefficients for outcome Y_p, including group and subjectlevel terms;
order is order of variables in corr.x[[p]][[p]] ; if betas = list() , all set to 0 so that all Y only have intercept
and/or interaction terms plus error terms; if all outcomes have the same betas, use list of length 1; if Y_p only has intercept
and/or interaction terms, set betas[[p]] = NULL ; if there are continuous mixture variables, beta is for mixture variable
(not for components)

betas.subj 
list of length M , pth component a vector of coefficients for interaction terms between grouplevel terms and
subjectlevel terms given in subj.var ; order is the same order as given in subj.var ; if all outcomes have the same betas,
use list of length 1; if Y_p only has grouplevel terms, set betas.subj[[p]] = NULL ;
if there are continuous mixture variables, beta is for mixture variable (not for components)

betas.int 
list of length M , pth component a vector of coefficients for interaction terms indicated in int.var ;
order is the same order as given in int.var ; if all outcomes have the same betas, use list of length 1;
if Y_p has none, set betas.int[[p]] = NULL ;
if there are continuous mixture variables, beta is for mixture variable (not for components)

betas.t 
vector of length M of coefficients for time terms, if NULL all set equal to 1;
if length 1, all intercepts set to betas.t

betas.tint 
list of length M , pth component a vector of coefficients for interaction terms indicated in tint.var ;
order is the same order as given in tint.var ; if all outcomes have the same betas, use list of length 1;
if Y_p has none, set betas.tint[[p]] = NULL ;
if there are continuous mixture variables, beta is for mixture variable (not for components)

rand.int 
"none" (default) if no random intercept term for all outcomes, "non_mix" if all random intercepts have a continuous
nonmixture distribution, "mix" if all random intercepts have a continuous mixture distribution;
also can be a vector of length M containing a combination (i.e., c("non_mix", "mix", "none") if the 1st has a nonmixture
distribution, the 2nd has a mixture distribution, and 3rd outcome has no random intercept)

rand.tsl 
"none" (default) if no random slope for time for all outcomes, "non_mix" if all random time slopes have a
continuous nonmixture distribution, "mix" if all random time slopes have a continuous mixture distribution; also can
be a vector of length M as in rand.int

rand.var 
matrix where 1st column is outcome index (p = 1, ..., M ), 2nd column is independent variable index
corresponding to covariate to assign random effect to (not including the random intercept or time slope if present);
the independent variable index corresponds to ordinal, continuous nonmixture, continuous mixture (not mixture component), Poisson, or
NB variable; order is 1st continuous nonmixture and 2nd continuous mixture random effects; note that the order of the rows corresponds
to the order of the random effects in corr.u not the order of the
independent variable so that a continuous mixture covariate with a nonmixture random effect would be ordered before a
continuous nonmixture covariate with a mixture random effect (the 2nd column of rand.var indicates the specific covariate)

corr.u 
if the random effects are the same variables across equations, a matrix of correlations for U;
if the random effects are different variables across equations, a list of length M , each component a list of length M ;
corr.u[[p]][[q]] is matrix of correlations for random effects in equations p (U_{(pj)} for outcome Y_p) and
q (U_{(qj)} for outcome Y_q);
if p = q, corr.u[[p]][[q]] is a correlation matrix with nrow(corr.u[[p]][[q]]) = # U_{(pj)} for outcome Y_p;
if p != q, corr.u[[p]][[q]] is a nonsymmetric matrix of correlations where rows correspond to U_{(pj)} for Y_p
so that nrow(corr.u[[p]][[q]]) = # U_{(pj)} for outcome Y_p and
columns correspond to U_{(qj)} for Y_q so that ncol(corr.u[[p]][[q]]) = # U_{(qj)} for outcome Y_q;
the number of random effects for Y_p is taken from nrow(corr.u[[p]][[1]]) so that if there should be random effects,
there must be entries for corr.u ;
use corr.u[[p]][[q]] = NULL if equation q has no U_{(qj)}; use corr.u[[p]] = NULL if equation p has no U_{(pj)};
correlations are specified in terms of components of mixture variables (if present);
order is 1st random intercept (if rand.int != "none" ), 2nd random time slope (if rand.tsl != "none" ),
3rd other random slopes with nonmixture distributions, 4th other random slopes with mixture distributions

quiet 
if FALSE prints messages, if TRUE suppresses messages
