SimRepeat (version 0.1.0)

summary_sys: Summary of Correlated Systems of Variables

Description

This function summarizes the results of nonnormsys, corrsys, or corrsys2. The inputs are either the simulated variables or inputs for those functions. See their documentation for more information. If only selected descriptions are desired, keep the non-relevant parameter inputs at their defaults. For example, if only a description of the error terms are desired, error_type = "non_mix", and method = "Polynomial", specify E, M, method, means, vars, skews, skurts, fifths, sixths, corr.e.

Usage

summary_sys(Y = NULL, E = NULL, E_mix = NULL, X = list(),
  X_all = list(), M = NULL, method = c("Fleishman", "Polynomial"),
  means = list(), vars = list(), skews = list(), skurts = list(),
  fifths = list(), sixths = list(), mix_pis = list(), mix_mus = list(),
  mix_sigmas = list(), mix_skews = list(), mix_skurts = list(),
  mix_fifths = list(), mix_sixths = list(), marginal = list(),
  support = list(), lam = list(), p_zip = list(), size = list(),
  prob = list(), mu = list(), p_zinb = list(), corr.x = list(),
  corr.e = NULL, U = list(), U_all = list(), rand.int = c("none",
  "non_mix", "mix"), rand.tsl = c("none", "non_mix", "mix"),
  corr.u = list(), rmeans2 = list(), rvars2 = list())

Arguments

Y

the matrix of outcomes simulated with corrsys or corrsys2

E

the matrix of continuous non-mixture or components of mixture error terms

E_mix

the matrix of continuous mixture error terms

X

a list of length M where X[[p]] = cbind(X_cat(pj), X_cont(pj), X_comp(pj), X_pois(pj), X_nb(pj)); keep X[[p]] = NULL if \(Y_p\) has no independent variables

X_all

a list of length M where X_all[[p]] contains all independent variables, interactions, and time for \(Y_p\); keep X_all[[p]] = NULL if \(Y_p\) has no independent variables

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 third-order polynomial transformation and "Polynomial" uses Headrick's fifth-order transformation

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 non-mixture (\(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 non-mixture 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 non-mixture 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 non-mixture (\(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 non-mixture 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 non-mixture 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 non-mixture 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 non-mixture 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 non-mixture distributions; not necessary for 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 j-th 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 j-th 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)}\)

marginal

a list of length M, with the p-th 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 j-th ordinal variable for outcome \(Y_p\); if the variable can take r values, the vector will contain r - 1 probabilities (the r-th 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 p-th 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 j-th ordinal variable for outcome \(Y_p\); if not provided, the default for r categories is 1, ..., r

lam

list of length M, p-th 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 zero-inflated 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 zero-inflated 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 zero-inflated Poisson distribution; if p_zip is in (-(exp(lam) - 1)^(-1), 0), \(Y_{pois}\) has a zero-deflated Poisson distribution and p_zip is not a probability; if p_zip = -(exp(lam) - 1)^(-1), \(Y_{pois}\) has a positive-Poisson distribution (see VGAM::dpospois); order is 1st regular Poisson and 2nd zero-inflated 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

size

list of length M, p-th component a vector of size parameters for the Negative Binomial variables for outcome \(Y_p\) (see stats::nbinom); order is 1st regular NB and 2nd zero-inflated 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, p-th component a vector of success probabilities for the Negative Binomial variables for outcome \(Y_p\) (see stats::nbinom); order is 1st regular NB and 2nd zero-inflated 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, p-th component a vector of mean values for the Negative Binomial variables for outcome \(Y_p\) (see stats::nbinom); order is 1st regular NB and 2nd zero-inflated 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 zero-inflated 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 zero-inflated 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 zero-deflated NB distribution and p_zinb is not a probability; if p_zinb = -prob^size/(1 - prob^size), \(Y_{nb}\) has a positive-NB distribution (see VGAM::dposnegbin); order is 1st regular NB and 2nd zero-inflated 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

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 non-mixture (same order as in skews), 3rd components of continuous mixture (same order as in mix_pis), 4th regular Poisson, 5th zero-inflated Poisson (same order as in lam), 6th regular NB, and 7th zero-inflated 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 non-symmetric 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.e

correlation matrix for continuous non-mixture or components of mixture error terms

U

a list of length M of continuous non-mixture and components of mixture random effects

U_all

a list of length M of continuous non-mixture and mixture random effects

rand.int

"none" (default) if no random intercept term for all outcomes, "non_mix" if all random intercepts have a continuous non-mixture 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 non-mixture 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 non-mixture distribution, "mix" if all random time slopes have a continuous mixture distribution; also can be a vector of length M as in rand.int

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 non-symmetric 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 non-mixture distributions, 4th other random slopes with mixture distributions

rmeans2

a list returned from corrsys or corrsys2 which has the non-mixture and component means ordered according to types of random intercept and time slope

rvars2

a list returned like rmeans

Value

A list with the following components:

cont_sum_y a data.frame summarizing the simulated distributions of the \(Y_p\),

cont_sum_e a data.frame summarizing the simulated distributions of the non-mixture or components of mixture \(E_p\),

target_sum_e a data.frame summarizing the target distributions of the non-mixture or components of mixture \(E_p\),

mix_sum_e a data.frame summarizing the simulated distributions of the mixture \(E_p\),

target_mix_e a data.frame summarizing the target distributions of the mixture \(E_p\),

rho.y correlation matrix of dimension M x M for \(Y_p\)

rho.e correlation matrix for the non-mixture or components of mixture \(E_p\)

rho.emix correlation matrix for the mixture \(E_p\)

rho.ye matrix with correlations between \(Y_p\) (rows) and the non-mixture or components of mixture \(E_p\) (columns)

rho.yemix matrix with correlations between \(Y_p\) (rows) and the mixture \(E_p\) (columns)

sum_xall a data.frame summarizing X_all without the Time variable,

rho.yx a list of length M, where rho.yx[[p]] is matrix of correlations between \(Y\) (rows) and X[[p]] = \({X_ord(pj), X_cont(pj), X_comp(pj), X_pois(pj), X_nb(pj)}\) (columns)

rho.yxall a list of length M, where rho.yx[[p]] is matrix of correlations between \(Y\) (rows) and X_all[[p]] (columns) not including Time

rho.x a list of length M of lists of length M where rho.x[[p]][[q]] = cor(cbind(X[[p]], X[[q]])) if p!= q or rho.x[[p]][[q]] = cor(X[[p]])) if p = q, where X[[p]] = \({X_ord(pj), X_cont(pj), X_comp(pj), X_pois(pj), X_nb(pj)}\)

rho.xall a list of length M of lists of length M where rho.xall[[p]][[q]] = cor(cbind(X_all[[p]], X_all[[q]])) if p!= q or rho.xall[[p]][[q]] = cor(X_all[[p]])) if p = q, not including Time

maxerr a list of length M containing a vector of length M with the maximum correlation errors between outcomes, maxerr[[p]]][q] = abs(max(corr.x[[p]][[q]] - rho.x[[p]][[q]]))

Additional components vary based on the type of simulated variables:

If ordinal variables are produced: ord_sum_x a list where ord_sum_x[[j]] is a data.frame summarizing \(X_{ord(pj)}\) for all p = 1, ..., M

If continuous variables are produced: cont_sum_x a data.frame summarizing the simulated distributions of the \(X_{cont(pj)}\) and \(X_comp(pj)\),

target_sum_x a data.frame summarizing the target distributions of the \(X_{cont(pj)}\) and \(X_comp(pj)\),

mix_sum_x a data.frame summarizing the simulated distributions of the \(X_{mix(pj)}\),

target_mix_x a data.frame summarizing the target distributions of the \(X_{mix(pj)}\)

If Poisson variables are produced: pois_sum_x a data.frame summarizing the simulated distributions of the \(X_{pois(pj)}\)

If Negative Binomial variables are produced: nb_sum_x a data.frame summarizing the simulated distributions of the \(X_{nb(pj)}\)

If random effects are produced: cont_sum_u a data.frame summarizing the simulated distributions of the \(U_{cont(pj)}\) and \(U_{comp(pj)}\),

target_sum_u a data.frame summarizing the target distributions of the \(U_{cont(pj)}\) and \(U_{comp(pj)}\),

sum_uall a data.frame summarizing the simulated distributions of U_all,

mix_sum_u a data.frame summarizing the simulated distributions of the \(U_{mix(pj)}\),

target_mix_u a data.frame summarizing the target distributions of the \(U_{mix(pj)}\),

rho.u list of length M, each component a list of length M; rho.u[[p]][[q]] = cor(cbind(U[[p]], U[[q]])) if p != q or rho.u[[p]][[q]] = cor(U[[p]])) if p = q

rho.uall list of length M, each component a list of length M; rho.uall[[p]][[q]] = cor(cbind(U_all[[p]], U_all[[q]])) if p != q or rho.uall[[p]][[q]] = cor(U_all[[p]])) if p = q

maxerr_u list of length M containing a vector of length M with the maximum correlation errors for \(U\) between outcomes maxerr_u[[p]]][q] = abs(max(corr.u[[p]][[q]] - rho.u[[p]][[q]]))

References

See references for SimRepeat.

See Also

nonnormsys, corrsys, corrsys2

Examples

# NOT RUN {
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)
Sys1 <- corrsys(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)
Sum1 <- summary_sys(Sys1$Y, Sys1$E, E_mix = NULL, Sys1$X, Sys1$X_all, M,
  "Polynomial", means, vars, skews, skurts, fifths, sixths,
  marginal = marginal, support = support, corr.x = corr.x, corr.e = corr.e)

# }
# 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 zero-inflated 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 zero-inflated NB variables
p_zinb <- list(NULL, 0.05, 0.1)

# 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 subject-level 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 = list(),
  size, prob, mu, p_zinb, nb_eps = list(), 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 1
N <- corrsys(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, size,
  prob, mu, p_zinb, 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)
S$sum_xall
S$maxerr
# }
# NOT RUN {
# }