Description Usage Arguments Value Reasons for Function Errors References See Also Examples
This function generates a correlated system of M
equations representing a system of repeated measures at
M
time points. The equations may contain 1) ordinal (r ≥q 2 categories), continuous (normal, nonnormal, and mixture distributions),
count (regular and zeroinflated, Poisson and Negative Binomial) independent variables X;
2) continuous error terms E; 3) a discrete time variable Time; and 4) random effects U. The assumptions are that
1) there are at least 2 equations, 2) the independent variables, random effect terms, and error terms are uncorrelated,
3) each equation has an error term, 4) all error terms have a continuous nonmixture distribution or all have a continuous mixture
distribution, 5) all random effects are continuous, and 6) growth is linear (with respect to time). The random effects may be a
random intercept, a random slope for time, or a random slope for any of the X variables. Continuous variables are simulated
using either Fleishman's thirdorder (method = "Fleishman"
, doi: 10.1007/BF02293811)
or Headrick's fifthorder (method = "Polynomial"
, doi: 10.1016/S01679473(02)000725) power method transformation (PMT).
Simulation occurs at the componentlevel for continuous mixture distributions. The target correlation matrix is specified in terms of
correlations with components of continuous mixture variables. These components are transformed into
the desired mixture variables using random multinomial variables based on the mixing probabilities. The X terms can be the same across equations (i.e., modeling sex or height) or may be timevarying covariates. The
equations may contain different numbers of X terms (i.e., a covariate could be missing for a given equation).
The outcomes Y are generated using a hierarchical linear models (HLM) approach, which allows the data to be structured in at least two levels.
Level1 is the repeated measure (time or condition) and other subjectlevel variables. Level1 is nested
within Level2, which describes the average of the outcome (the intercept) and growth (slope for time) as a function
of grouplevel variables. The first level captures the withinsubject variation, while the second level describes the betweensubjects variability.
Using a HLM provides a way to determine if: a) subjects differ at a specific time point with respect to the dependent
variable, b) growth rates differ across conditions, or c) growth rates differ across subjects. Random effects
describe deviation at the subjectlevel from the average (fixed) effect described by the slope coefficients (betas). See the
The Hierarchical Linear Models Approach for a System of Correlated Equations with Multiple Variable Types vignette for a
description of the HLM model. The user can specify subjectlevel X terms, and each subjectlevel X term is crossed with
all grouplevel X terms. The equations may also contain interactions between X variables. Interactions specified in
int.var
between two grouplevel covariates are themselves considered grouplevel covariates and will be crossed with subjectlevel
covariates. Interactions between two subjectlevel covariates are considered subjectlevel covariates and will be crossed with
grouplevel covariates. Since Time
is a subjectlevel variable, each grouplevel term is crossed with Time
unless otherwise specified.
Random effects may be added for the intercept, time slope, or effects of any of the covariates. The type of random intercept and
time slope (i.e., nonmixture or mixture) is specified in rand.int
and rand.tsl
. This type may vary by equation.
The random effects for independent variables are specified in rand.var
and may also contain a combination of nonmixture and
mixture continuous distributions. If the parameter lists are of length M + 1
, the random effects are the same variables across equations and the
correlation for the effects corr.u
is a matrix. If the parameter lists are of length 2 * M
, the random effects are different variables
across equations and the correlation for the effects corr.u
is a list.
The independent variables, interactions, Time
effect, random effects, and error terms are summed
together to produce the outcomes Y. The beta coefficients may be the same or differ across equations. The user specifies the betas for
the independent variables in betas
, for the interactions between two grouplevel or two subjectlevel covariates in betas.int
,
for the groupsubject level interactions in betas.subj
, and for the Time
interactions in betas.tint
.
Setting a coefficient to 0 will eliminate that term. The user also provides the correlations 1) between E terms; 2) between X
variables within each outcome Y_p, p = 1, ..., M
, and between outcome pairs; and 3) between U variables within each
outcome Y_p, p = 1, ..., M
, and between outcome pairs. The order of the independent variables in corr.x
must be
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
). The order of the random effects in corr.u
must be 1st random intercept, 2nd random
time slope, 3rd continuous nonmixture random effects, and 4th components of continuous mixture random effects.
The variables are generated from multivariate normal variables with intermediate correlations calculated using
intercorr2
, which employs correlation method 2. See SimCorrMix
for a
description of the correlation method and the techniques used to generate each variable type. The order of the variables returned is
1st covariates X (as specified in corr.x
), 2nd groupgroup or subjectsubject interactions (ordered as in int.var
),
3rd subjectgroup interactions (1st by subjectlevel variable as specified in subj.var
, 2nd by covariate as specified in corr.x
),
and 4th time interactions (either as specified in corr.x
with grouplevel covariates or in tint.var
).
This function contains no parameter checks in order to decrease simulation time. That should be done first using
checkpar
. Summaries of the system can be obtained using summary_sys
. The
Correlated Systems of Statistical Equations with Multiple Variable Types demonstrates examples.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17  corrsys2(n = 10000, M = NULL, Time = NULL, method = c("Fleishman",
"Polynomial"), error_type = c("non_mix", "mix"), means = list(),
vars = list(), skews = list(), skurts = list(), fifths = list(),
sixths = list(), Six = list(), mix_pis = list(), mix_mus = list(),
mix_sigmas = list(), mix_skews = list(), mix_skurts = list(),
mix_fifths = list(), mix_sixths = list(), mix_Six = list(),
marginal = list(), support = list(), lam = list(), p_zip = list(),
pois_eps = list(), size = list(), prob = list(), mu = list(),
p_zinb = list(), nb_eps = list(), corr.x = list(), corr.e = NULL,
same.var = NULL, subj.var = NULL, int.var = NULL, tint.var = NULL,
betas.0 = NULL, betas = list(), betas.subj = list(),
betas.int = list(), betas.t = NULL, betas.tint = list(),
rand.int = c("none", "non_mix", "mix"), rand.tsl = c("none", "non_mix",
"mix"), rand.var = NULL, corr.u = list(), seed = 1234,
use.nearPD = TRUE, eigmin = 0, adjgrad = FALSE, B1 = NULL,
tau = 0.5, tol = 0.1, steps = 100, msteps = 10, errorloop = FALSE,
epsilon = 0.001, maxit = 1000, quiet = FALSE)

n 
the sample size (i.e. the length of each simulated variable; default = 10000) 
M 
the number of dependent variables Y (outcomes); equivalently, the number of equations in the system 
Time 
a vector of values to use for time; each subject receives the same time value; if 
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 if there are random effects, a list of length 
vars 
a list of same length and order as 
skews 
if no random effects, a list of length if there are random effects, a list of length 
skurts 
a list of same length and order as 
fifths 
a list of same length and order as 
sixths 
a list of same length and order as 
Six 
a list of length
keep 
mix_pis 
list of length

mix_mus 
list of same length and order as

mix_sigmas 
list of same length and order as

mix_skews 
list of same length and order as

mix_skurts 
list of same length and order as

mix_fifths 
list of same length and order as

mix_sixths 
list of same length and order as

mix_Six 
a list of same length and order as pth component of qth component of 
marginal 
a list of length 
support 
a list of length 
lam 
list of length 
p_zip 
a list of vectors of probabilities of structural zeros (not including zeros from the Poisson distribution) for the
zeroinflated Poisson variables (see 
pois_eps 
list of length 
size 
list of length 
prob 
list of length 
mu 
list of length 
p_zinb 
a vector of probabilities of structural zeros (not including zeros from the NB distribution) for the zeroinflated NB variables
(see 
nb_eps 
list of length 
corr.x 
list of length 
corr.e 
correlation matrix for continuous nonmixture or components of mixture error terms 
same.var 
either a vector or a matrix; if a vector, 
subj.var 
matrix where 1st column is outcome index ( 
int.var 
matrix where 1st column is outcome index ( 
tint.var 
matrix where 1st column is outcome index ( 
betas.0 
vector of length 
betas 
list of length 
betas.subj 
list of length 
betas.int 
list of length 
betas.t 
vector of length 
betas.tint 
list of length 
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 
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 
rand.var 
matrix where 1st column is outcome index ( 
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 correlations are specified in terms of components of mixture variables (if present);
order is 1st random intercept (if 
seed 
the seed value for random number generation (default = 1234) 
use.nearPD 
TRUE to convert the overall intermediate correlation matrix formed by the X (for all outcomes and independent
variables), E, or the random effects to the nearest positive definite matrix with 
eigmin 
minimum replacement eigenvalue if overall intermediate correlation matrix is not positivedefinite (default = 0) 
adjgrad 
TRUE to use 
B1 
the initial matrix for algorithm; if NULL, uses a scaled initial matrix with diagonal elements 
tau 
parameter used to calculate theta (default = 0.5) 
tol 
maximum error for Frobenius norm distance between new matrix and original matrix (default = 0.1) 
steps 
maximum number of steps for k (default = 100) 
msteps 
maximum number of steps for m (default = 10) 
errorloop 
if TRUE, uses 
epsilon 
the maximum acceptable error between the final and target correlation matrices (default = 0.001)
in the calculation of ordinal intermediate correlations with 
maxit 
the maximum number of iterations to use (default = 1000) in the calculation of ordinal
intermediate correlations with 
quiet 
if FALSE prints messages, if TRUE suppresses messages 
A list with the following components:
Y
matrix with n
rows and M
columns of outcomes
X
list of length M
containing X_{ord(pj)}, X_{cont(pj)}, X_{comp(pj)}, X_{pois(pj)}, X_{nb(pj)}
X_all
list of length M
containing X_{ord(pj)}, X_{cont(pj)}, X_{mix(pj)}, X_{pois(pj)}, X_{nb(pj)}, X interactions as indicated by
int.var
, subjectgroup level term interactions as indicated by subj.var
, Time_p, and Time interactions as indicated by
tint.var
; order is 1st covariates X (as specified in corr.x
), 2nd groupgroup or subjectsubject interactions (ordered as in int.var
),
3rd subjectgroup interactions (1st by subjectlevel variable as specified in subj.var
, 2nd by covariate as specified in
corr.x
), and 4th time interactions (either as specified in corr.x
with grouplevel covariates or in tint.var
)
E
matrix with n
rows containing continuous nonmixture or components of continuous mixture error terms
E_mix
matrix with n
rows containing continuous mixture error terms
Sigma_X0
matrix of intermediate correlations calculated by intercorr2
Sigma_X
matrix of intermediate correlations after nearPD
or adj_grad
function has been used;
applied to generate the normal variables transformed to get the desired distributions
Error_Time
the time in minutes required to use the error loop
Time
the total simulation time in minutes
niter
a matrix of the number of iterations used in the error loop
If continuous variables are produced:
constants
a list of maximum length 2 * M
, the 1st M
components are data.frames of the constants for the
X_{cont(pj)}, X_comp(pj) and E_p, the 2nd M
components are for random effects (if present),
SixCorr
a list of maximum length 2 * M
, the 1st M
components are lists of sixth cumulant correction
values used to obtain valid pdf's for the X_{cont(pj)}, X_comp(pj), and E_p, the 2nd M
components are for random effects (if present),
valid.pdf
a list of maximum length 2 * M
of vectors where the ith element is "TRUE" if the constants for the ith
continuous variable generate a valid pdf, else "FALSE"; the 1st M
components are for the X_{cont(pj)},
X_comp(pj), and E_p, the 2nd M
components are for random effects (if present)
If random effects are produced:
U
a list of length M
containing matrices of continuous nonmixture and components of mixture random effects,
U_all
a list of length M
containing matrices of continuous nonmixture and mixture random effects,
V
a list of length M
containing matrices of design matrices for random effects,
rmeans2
and rvars2
the means and variances of the nonmixture and components reordered in accordance with the
random intercept and time slope types (input for summary_sys
)
1) The most likely cause for function errors is that the parameter inputs are mispecified. Using checkpar
prior to simulation can help decrease these errors.
2) Another reason for error is that no solutions to fleish
or
poly
converged when using find_constants
. If this happens,
the simulation will stop. It may help to first use find_constants
for each continuous variable to
determine if a sixth cumulant correction value is needed. If the standardized cumulants are obtained from calc_theory
, the user
may need to use rounded values as inputs (i.e. skews = round(skews, 8)
). For example, in order to ensure that skew is exactly 0 for symmetric distributions.
3) The kurtosis for a continuous variable may be outside the region of possible values. There is an associated lower kurtosis boundary for
associated with a given skew (for Fleishman's method) or skew and fifth and sixth cumulants (for Headrick's method). Use
calc_lower_skurt
to determine the boundary for a given set of cumulants.
Barbiero A & Ferrari PA (2015). Simulation of correlated Poisson variables. Applied Stochastic Models in Business and Industry, 31:66980. doi: 10.1002/asmb.2072.
Barbiero A & Ferrari PA (2015). GenOrd: Simulation of Discrete Random Variables with Given
Correlation Matrix and Marginal Distributions. R package version 1.4.0.
https://CRAN.Rproject.org/package=GenOrd
Davenport JW, Bezder JC, & Hathaway RJ (1988). Parameter Estimation for Finite Mixture Distributions. Computers & Mathematics with Applications, 15(10):81928.
Demirtas H (2006). A method for multivariate ordinal data generation given marginal distributions and correlations. Journal of Statistical
Computation and Simulation, 76(11):10171025.
doi: 10.1080/10629360600569246.
Demirtas H (2014). Joint Generation of Binary and Nonnormal Continuous Data. Biometrics & Biostatistics, S12.
Demirtas H, Hedeker D, & Mermelstein RJ (2012). Simulation of massive public health data by power polynomials. Statistics in Medicine, 31(27):33373346. doi: 10.1002/sim.5362.
Everitt BS (1996). An Introduction to Finite Mixture Distributions. Statistical Methods in Medical Research, 5(2):107127. doi: 10.1177/096228029600500202.
Ferrari PA & Barbiero A (2012). Simulating ordinal data. Multivariate Behavioral Research, 47(4): 566589. doi: 10.1080/00273171.2012.692630.
Fialkowski AC (2017). SimMultiCorrData: Simulation of Correlated Data with Multiple Variable Types. R package version 0.2.1. https://CRAN.Rproject.org/package=SimMultiCorrData.
Fialkowski AC (2018). SimCorrMix: Simulation of Correlated Data of Multiple Variable Types including Continuous and Count Mixture Distributions. R package version 0.1.0. https://github.com/AFialkowski/SimCorrMix
Fleishman AI (1978). A Method for Simulating Nonnormal Distributions. Psychometrika, 43:521532. doi: 10.1007/BF02293811.
Headrick TC (2002). Fast Fifthorder Polynomial Transforms for Generating Univariate and Multivariate Nonnormal Distributions. Computational Statistics & Data Analysis, 40(4):685711. doi: 10.1016/S01679473(02)000725. (ScienceDirect)
Headrick TC (2004). On Polynomial Transformations for Simulating Multivariate Nonnormal Distributions. Journal of Modern Applied Statistical Methods, 3(1):6571. doi: 10.22237/jmasm/1083370080.
Headrick TC, Kowalchuk RK (2007). The Power Method Transformation: Its Probability Density Function, Distribution Function, and Its Further Use for Fitting Data. Journal of Statistical Computation and Simulation, 77:229249. doi: 10.1080/10629360600605065.
Headrick TC, Sawilowsky SS (1999). Simulating Correlated Nonnormal Distributions: Extending the Fleishman Power Method. Psychometrika, 64:2535. doi: 10.1007/BF02294317.
Headrick TC, Sheng Y, & Hodis FA (2007). Numerical Computing and Graphics for the Power Method Transformation Using
Mathematica. Journal of Statistical Software, 19(3):1  17.
doi: 10.18637/jss.v019.i03.
Higham N (2002). Computing the nearest correlation matrix  a problem from finance; IMA Journal of Numerical Analysis 22:329343.
Ismail N & Zamani H (2013). Estimation of Claim Count Data Using Negative Binomial, Generalized Poisson, ZeroInflated Negative Binomial and ZeroInflated Generalized Poisson Regression Models. Casualty Actuarial Society EForum 41(20):128.
Kincaid C (2005). Guidelines for Selecting the Covariance Structure in Mixed Model Analysis. Computational Statistics and Data Analysis, 198(30):18.
Lambert D (1992). ZeroInflated Poisson Regression, with an Application to Defects in Manufacturing. Technometrics 34(1):114.
Lininger M, Spybrook J, & Cheatham CC (2015). Hierarchical Linear Model: Thinking Outside the Traditional RepeatedMeasures AnalysisofVariance Box. Journal of Athletic Training, 50(4):438441. doi: 10.4085/1062605049.5.09.
McCulloch CE, Searle SR, Neuhaus JM (2008). Generalized, Linear, and Mixed Models (2nd ed.). Wiley Series in Probability and Statistics. Hoboken, New Jersey: John Wiley & Sons, Inc.
Olsson U, Drasgow F, & Dorans NJ (1982). The Polyserial Correlation Coefficient. Psychometrika, 47(3):33747. doi: 10.1007/BF02294164.
Pearson RK (2011). Exploring Data in Engineering, the Sciences, and Medicine. In. New York: Oxford University Press.
Schork NJ, Allison DB, & Thiel B (1996). Mixture Distributions in Human Genetics Research. Statistical Methods in Medical Research, 5:155178. doi: 10.1177/096228029600500204.
Vale CD & Maurelli VA (1983). Simulating Multivariate Nonnormal Distributions. Psychometrika, 48:465471. doi: 10.1007/BF02293687.
Van Der Leeden R (1998). Multilevel Analysis of Repeated Measures Data. Quality & Quantity, 32(1):1529.
Yee TW (2017). VGAM: Vector Generalized Linear and Additive Models.
https://CRAN.Rproject.org/package=VGAM.
Zhang X, Mallick H, & Yi N (2016). ZeroInflated Negative Binomial Regression for Differential Abundance Testing in Microbiome Studies. Journal of Bioinformatics and Genomics 2(2):19. doi: 10.18454/jbg.2016.2.2.1.
find_constants
, intercorr2
,
checkpar
, summary_sys
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140  M < 3
B < calc_theory("Beta", c(4, 1.5))
skews < lapply(seq_len(M), function(x) B[3])
skurts < lapply(seq_len(M), function(x) B[4])
fifths < lapply(seq_len(M), function(x) B[5])
sixths < lapply(seq_len(M), function(x) B[6])
Six < lapply(seq_len(M), function(x) list(0.03))
corr.e < matrix(c(1, 0.4, 0.4^2, 0.4, 1, 0.4, 0.4^2, 0.4, 1), M, M,
byrow = TRUE)
means < lapply(seq_len(M), function(x) B[1])
vars < lapply(seq_len(M), function(x) B[2]^2)
marginal < list(0.3, 0.4, 0.5)
support < lapply(seq_len(M), function(x) list(0:1))
corr.x < list(list(matrix(1, 1, 1), matrix(0.4, 1, 1), matrix(0.4, 1, 1)),
list(matrix(0.4, 1, 1), matrix(1, 1, 1), matrix(0.4, 1, 1)),
list(matrix(0.4, 1, 1), matrix(0.4, 1, 1), matrix(1, 1, 1)))
betas < list(0.5)
betas.t < 1
betas.tint < list(0.25)
Sys2 < corrsys2(10000, M, Time = 1:M, "Polynomial", "non_mix", means, vars,
skews, skurts, fifths, sixths, Six, marginal = marginal, support = support,
corr.x = corr.x, corr.e = corr.e, betas = betas, betas.t = betas.t,
betas.tint = betas.tint, quiet = TRUE)
## Not run:
seed < 276
n < 10000
M < 3
Time < 1:M
# Error terms have a beta(4, 1.5) distribution with an AR(1, p = 0.4)
# correlation structure
B < calc_theory("Beta", c(4, 1.5))
skews < lapply(seq_len(M), function(x) B[3])
skurts < lapply(seq_len(M), function(x) B[4])
fifths < lapply(seq_len(M), function(x) B[5])
sixths < lapply(seq_len(M), function(x) B[6])
Six < lapply(seq_len(M), function(x) list(0.03))
error_type < "non_mix"
corr.e < matrix(c(1, 0.4, 0.4^2, 0.4, 1, 0.4, 0.4^2, 0.4, 1), M, M,
byrow = TRUE)
# 1 continuous mixture of Normal(2, 1) and Normal(2, 1) for each Y
mix_pis < lapply(seq_len(M), function(x) list(c(0.4, 0.6)))
mix_mus < lapply(seq_len(M), function(x) list(c(2, 2)))
mix_sigmas < lapply(seq_len(M), function(x) list(c(1, 1)))
mix_skews < lapply(seq_len(M), function(x) list(c(0, 0)))
mix_skurts < lapply(seq_len(M), function(x) list(c(0, 0)))
mix_fifths < lapply(seq_len(M), function(x) list(c(0, 0)))
mix_sixths < lapply(seq_len(M), function(x) list(c(0, 0)))
mix_Six < list()
Nstcum < calc_mixmoments(mix_pis[[1]][[1]], mix_mus[[1]][[1]],
mix_sigmas[[1]][[1]], mix_skews[[1]][[1]], mix_skurts[[1]][[1]],
mix_fifths[[1]][[1]], mix_sixths[[1]][[1]])
means < lapply(seq_len(M), function(x) c(Nstcum[1], B[1]))
vars < lapply(seq_len(M), function(x) c(Nstcum[2]^2, B[2]^2))
# 1 binary variable for each Y
marginal < lapply(seq_len(M), function(x) list(0.4))
support < list(NULL, list(c(0, 1)), NULL)
# 1 Poisson variable for each Y
lam < list(1, 5, 10)
# Y2 and Y3 are zeroinflated Poisson variables
p_zip < list(NULL, 0.05, 0.1)
# 1 NB variable for each Y
size < list(10, 15, 20)
prob < list(0.3, 0.4, 0.5)
# either prob or mu is required (not both)
mu < mapply(function(x, y) x * (1  y)/y, size, prob, SIMPLIFY = FALSE)
# Y2 and Y3 are zeroinflated NB variables
p_zinb < list(NULL, 0.05, 0.1)
pois_eps < nb_eps < list()
# The 2nd (the normal mixture) variable is the same across Y
same.var < 2
# Create the correlation matrix in terms of the components of the normal
# mixture
K < 5
corr.x < list()
corr.x[[1]] < list(matrix(0.1, K, K), matrix(0.2, K, K), matrix(0.3, K, K))
diag(corr.x[[1]][[1]]) < 1
# set correlation between components to 0
corr.x[[1]][[1]][2:3, 2:3] < diag(2)
# set correlations with the same variable equal across outcomes
corr.x[[1]][[2]][, same.var] < corr.x[[1]][[3]][, same.var] <
corr.x[[1]][[1]][, same.var]
corr.x[[2]] < list(t(corr.x[[1]][[2]]), matrix(0.35, K, K),
matrix(0.4, K, K))
diag(corr.x[[2]][[2]]) < 1
corr.x[[2]][[2]][2:3, 2:3] < diag(2)
corr.x[[2]][[2]][, same.var] < corr.x[[2]][[3]][, same.var] <
t(corr.x[[1]][[2]][same.var, ])
corr.x[[2]][[3]][same.var, ] < corr.x[[1]][[3]][same.var, ]
corr.x[[2]][[2]][same.var, ] < t(corr.x[[2]][[2]][, same.var])
corr.x[[3]] < list(t(corr.x[[1]][[3]]), t(corr.x[[2]][[3]]),
matrix(0.5, K, K))
diag(corr.x[[3]][[3]]) < 1
corr.x[[3]][[3]][2:3, 2:3] < diag(2)
corr.x[[3]][[3]][, same.var] < t(corr.x[[1]][[3]][same.var, ])
corr.x[[3]][[3]][same.var, ] < t(corr.x[[3]][[3]][, same.var])
# The 2nd and 3rd variables of each Y are subjectlevel variables
subj.var < matrix(c(1, 2, 1, 3, 2, 2, 2, 3, 3, 2, 3, 3), 6, 2, byrow = TRUE)
int.var < tint.var < NULL
betas.0 < 0
betas < list(seq(0.5, 0.5 + (K  2) * 0.25, 0.25))
betas.subj < list(seq(0.5, 0.5 + (K  2) * 0.1, 0.1))
betas.int < list()
betas.t < 1
betas.tint < list(c(0.25, 0.5))
method < "Polynomial"
# Check parameter inputs
checkpar(M, method, error_type, means, vars, skews, skurts, fifths, sixths,
Six, mix_pis, mix_mus, mix_sigmas, mix_skews, mix_skurts, mix_fifths,
mix_sixths, mix_Six, marginal, support, lam, p_zip, pois_eps, size, prob,
mu, p_zinb, nb_eps, corr.x, corr.yx = list(), corr.e, same.var, subj.var,
int.var, tint.var, betas.0, betas, betas.subj, betas.int, betas.t,
betas.tint)
# Simulated system using correlation method 2
N < corrsys2(n, M, Time, method, error_type, means, vars, skews, skurts,
fifths, sixths, Six, mix_pis, mix_mus, mix_sigmas, mix_skews, mix_skurts,
mix_fifths, mix_sixths, mix_Six, marginal, support, lam, p_zip, pois_eps,
size, prob, mu, p_zinb, nb_eps, corr.x, corr.e, same.var, subj.var,
int.var, tint.var, betas.0, betas, betas.subj, betas.int, betas.t,
betas.tint, seed = seed, use.nearPD = FALSE)
# Summarize the results
S < summary_sys(N$Y, N$E, E_mix = NULL, N$X, N$X_all, M, method, means,
vars, skews, skurts, fifths, sixths, mix_pis, mix_mus, mix_sigmas,
mix_skews, mix_skurts, mix_fifths, mix_sixths, marginal, support, lam,
p_zip, size, prob, mu, p_zinb, corr.x, corr.e)
## End(Not run)

Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.