Nothing
#' @title Summary of Correlated Systems of Variables
#'
#' @description This function summarizes the results of \code{\link[SimRepeat]{nonnormsys}}, \code{\link[SimRepeat]{corrsys}}, or
#' \code{\link[SimRepeat]{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, \code{error_type = "non_mix"}, and
#' \code{method = "Polynomial"}, specify \code{E, M, method, means, vars, skews, skurts, fifths, sixths, corr.e}.
#'
#' @param Y the matrix of outcomes simulated with \code{corrsys} or \code{corrsys2}
#' @param E the matrix of continuous non-mixture or components of mixture error terms
#' @param E_mix the matrix of continuous mixture error terms
#' @param X a list of length \code{M} where \code{X[[p]] = cbind(X_cat(pj), X_cont(pj), X_comp(pj), X_pois(pj), X_nb(pj))}; keep
#' \code{X[[p]] = NULL} if \eqn{Y_p} has no independent variables
#' @param X_all a list of length \code{M} where \code{X_all[[p]]} contains all independent variables, interactions, and time for \eqn{Y_p}; keep
#' \code{X_all[[p]] = NULL} if \eqn{Y_p} has no independent variables
#' @param M the number of dependent variables \eqn{Y} (outcomes); equivalently, the number of equations in the system
#' @param 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
#' @param means if no random effects, a list of length \code{M} where \code{means[[p]]} contains a vector of means for the continuous independent variables
#' in equation p with non-mixture (\eqn{X_{cont}}) or mixture (\eqn{X_{mix}}) distributions and for the error terms (\eqn{E});
#' order in vector is \eqn{X_{cont}, X_{mix}, E}
#'
#' if there are random effects, a list of length \code{M + 1} if the effects are the same across equations or \code{2 * M} if they differ;
#' where \code{means[M + 1]} or \code{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 \eqn{U_0} (if \code{rand.int != "none"}), 2nd random time slope \eqn{U_1} (if \code{rand.tsl != "none"}),
#' 3rd other random slopes with non-mixture distributions \eqn{U_{cont}}, 4th other random slopes with mixture distributions \eqn{U_{mix}}
#' @param vars a list of same length and order as \code{means} containing vectors of variances for the continuous variables, error terms, and any random effects
#' @param skews if no random effects, a list of length \code{M} where \code{skews[[p]]} contains a vector of skew values for the continuous independent variables
#' in equation p with non-mixture (\eqn{X_{cont}}) distributions and for \eqn{E} if \code{error_type = "non_mix"}; order in vector is \eqn{X_{cont}, E}
#'
#' if there are random effects, a list of length \code{M + 1} if the effects are the same across equations or \code{2 * M} if they differ;
#' where \code{skews[M + 1]} or \code{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 \eqn{U_0} (if \code{rand.int = "non_mix"}), 2nd random time slope \eqn{U_1} (if \code{rand.tsl = "non_mix"}),
#' 3rd other random slopes with non-mixture distributions \eqn{U_{cont}}
#' @param skurts a list of same length and order as \code{skews} containing vectors of standardized kurtoses (kurtosis - 3) for the continuous variables,
#' error terms, and any random effects with non-mixture distributions
#' @param fifths a list of same length and order as \code{skews} containing vectors of standardized fifth cumulants for the continuous variables,
#' error terms, and any random effects with non-mixture distributions; not necessary for \code{method = "Fleishman"}
#' @param sixths a list of same length and order as \code{skews} containing vectors of standardized sixth cumulants for the continuous variables,
#' error terms, and any random effects with non-mixture distributions; not necessary for \code{method = "Fleishman"}
#' @param mix_pis list of length \code{M}, \code{M + 1} or \code{2 * M}, where \code{mix_pis[1:M]} are for \eqn{X_{cont}, E} (if \code{error_type = "mix"}) and
#' \code{mix_pis[M + 1]} or \code{mix_pis[(M + 1):(2 * M)]} are for mixture \eqn{U}; use \code{mix_pis[[p]] = NULL} if equation p has no continuous mixture terms
#' if \code{error_type = "non_mix"} and there are only random effects (i.e., \code{length(corr.x) = 0}), use \code{mix_pis[1:M] = NULL} so that
#' \code{mix_pis[M + 1]} or \code{mix_pis[(M + 1):(2 * M)]} describes the mixture \eqn{U};
#'
#' \code{mix_pis[[p]][[j]]} is a vector of mixing probabilities of the component distributions for \eqn{X_{mix(pj)}}, the j-th mixture covariate for outcome \eqn{Y_p};
#' the last vector in \code{mix_pis[[p]]} is for \eqn{E_p} (if \code{error_type = "mix"}); components should be ordered as in \code{corr.x}
#'
#' \code{mix_pis[[M + p]][[j]]} is a vector of mixing probabilities of the component distributions for \eqn{U_{(pj)}}, the j-th random effect with a mixture distribution
#' for outcome \eqn{Y_p}; order is 1st random intercept (if \code{rand.int = "mix"}), 2nd random time slope (if \code{rand.tsl = "mix"}),
#' 3rd other random slopes with mixture distributions; components should be ordered as in \code{corr.u}
#' @param mix_mus list of same length and order as \code{mix_pis};
#'
#' \code{mix_mus[[p]][[j]]} is a vector of means of the component distributions for \eqn{X_{mix(pj)}},
#' the last vector in \code{mix_mus[[p]]} is for \eqn{E_p} (if \code{error_type = "mix"})
#'
#' \code{mix_mus[[p]][[j]]} is a vector of means of the component distributions for \eqn{U_{mix(pj)}}
#' @param mix_sigmas list of same length and order as \code{mix_pis};
#'
#' \code{mix_sigmas[[p]][[j]]} is a vector of standard deviations of the component distributions for \eqn{X_{mix(pj)}},
#' the last vector in \code{mix_sigmas[[p]]} is for \eqn{E_p} (if \code{error_type = "mix"})
#'
#' \code{mix_sigmas[[p]][[j]]} is a vector of standard deviations of the component distributions for \eqn{U_{mix(pj)}}
#' @param mix_skews list of same length and order as \code{mix_pis};
#'
#' \code{mix_skews[[p]][[j]]} is a vector of skew values of the component distributions for \eqn{X_{mix(pj)}},
#' the last vector in \code{mix_skews[[p]]} is for \eqn{E_p} (if \code{error_type = "mix"})
#'
#' \code{mix_skews[[p]][[j]]} is a vector of skew values of the component distributions for \eqn{U_{mix(pj)}}
#' @param mix_skurts list of same length and order as \code{mix_pis};
#'
#' \code{mix_skurts[[p]][[j]]} is a vector of standardized kurtoses of the component distributions for \eqn{X_{mix(pj)}},
#' the last vector in \code{mix_skurts[[p]]} is for \eqn{E_p} (if \code{error_type = "mix"})
#'
#' \code{mix_skurts[[p]][[j]]} is a vector of standardized kurtoses of the component distributions for \eqn{U_{mix(pj)}}
#' @param mix_fifths list of same length and order as \code{mix_pis}; not necessary for \code{method = "Fleishman"};
#'
#' \code{mix_fifths[[p]][[j]]} is a vector of standardized fifth cumulants of the component distributions for \eqn{X_{mix(pj)}},
#' the last vector in \code{mix_fifths[[p]]} is for \eqn{E_p} (if \code{error_type = "mix"})
#'
#' \code{mix_fifths[[p]][[j]]} is a vector of standardized fifth cumulants of the component distributions for \eqn{U_{mix(pj)}}
#' @param mix_sixths list of same length and order as \code{mix_pis}; not necessary for \code{method = "Fleishman"};
#'
#' \code{mix_sixths[[p]][[j]]} is a vector of standardized sixth cumulants of the component distributions for \eqn{X_{mix(pj)}},
#' the last vector in \code{mix_sixths[[p]]} is for \eqn{E_p} (if \code{error_type = "mix"})
#'
#' \code{mix_sixths[[p]][[j]]} is a vector of standardized sixth cumulants of the component distributions for \eqn{U_{mix(pj)}}
#' @param marginal a list of length \code{M}, with the p-th component a list of cumulative probabilities for the ordinal variables
#' associated with outcome \eqn{Y_p} (use \code{marginal[[p]] = NULL} if outcome \eqn{Y_p} has no ordinal variables);
#' \code{marginal[[p]][[j]]} is a vector of the cumulative probabilities defining the marginal distribution of \eqn{X_{ord(pj)}},
#' the j-th ordinal variable for outcome \eqn{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;
#' \code{length(marginal[[p]])} can differ across outcomes; the order should be the same as in \code{corr.x}
#' @param support a list of length \code{M}, with the p-th component a list of support values for the ordinal variables associated
#' with outcome \eqn{Y_p}; use \code{support[[p]] = NULL} if outcome \eqn{Y_p} has no ordinal variables;
#' \code{support[[p]][[j]]} is a vector of the support values defining the marginal distribution of \eqn{X_{ord(pj)}},
#' the j-th ordinal variable for outcome \eqn{Y_p}; if not provided, the default for r categories is 1, ..., r
#' @param lam list of length \code{M}, p-th component a vector of lambda (means > 0) values for Poisson variables for outcome \eqn{Y_p}
#' (see \code{stats::dpois}); order is 1st regular Poisson and 2nd zero-inflated Poisson; use \code{lam[[p]] = NULL} if outcome \eqn{Y_p} has no Poisson variables;
#' \code{length(lam[[p]])} can differ across outcomes; the order should be the same as in \code{corr.x}
#' @param 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 \code{VGAM::dzipois}); if \code{p_zip} = 0, \eqn{Y_{pois}} has a regular Poisson
#' distribution; if \code{p_zip} is in (0, 1), \eqn{Y_{pois}} has a zero-inflated Poisson distribution;
#' if \code{p_zip} is in \code{(-(exp(lam) - 1)^(-1), 0)}, \eqn{Y_{pois}} has a zero-deflated Poisson distribution and \code{p_zip}
#' is not a probability; if \code{p_zip = -(exp(lam) - 1)^(-1)}, \eqn{Y_{pois}} has a positive-Poisson distribution
#' (see \code{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 \code{M}, all Poisson variables in equation p
#' given \code{p_zip[p]}; otherwise, missing values are set to 0 and ordered 1st
#' @param size list of length \code{M}, p-th component a vector of size parameters for the Negative Binomial variables for outcome \eqn{Y_p}
#' (see \code{stats::nbinom}); order is 1st regular NB and 2nd zero-inflated NB; use \code{size[[p]] = NULL} if outcome \eqn{Y_p} has no Negative Binomial variables;
#' \code{length(size[[p]])} can differ across outcomes; the order should be the same as in \code{corr.x}
#' @param prob list of length \code{M}, p-th component a vector of success probabilities for the Negative Binomial variables for outcome \eqn{Y_p}
#' (see \code{stats::nbinom}); order is 1st regular NB and 2nd zero-inflated NB; use \code{prob[[p]] = NULL} if outcome \eqn{Y_p} has no Negative Binomial variables;
#' \code{length(prob[[p]])} can differ across outcomes; the order should be the same as in \code{corr.x}
#' @param mu list of length \code{M}, p-th component a vector of mean values for the Negative Binomial variables for outcome \eqn{Y_p}
#' (see \code{stats::nbinom}); order is 1st regular NB and 2nd zero-inflated NB; use \code{mu[[p]] = NULL} if outcome \eqn{Y_p} has no Negative Binomial variables;
#' \code{length(mu[[p]])} can differ across outcomes; the order should be the same as in \code{corr.x}; for zero-inflated NB variables,
#' this refers to the mean of the NB distribution (see \code{VGAM::dzinegbin})
#' (*Note: either \code{prob} or \code{mu} should be supplied for all Negative Binomial variables, not a mixture)
#' @param p_zinb a vector of probabilities of structural zeros (not including zeros from the NB distribution) for the zero-inflated NB variables
#' (see \code{VGAM::dzinegbin}); if \code{p_zinb} = 0, \eqn{Y_{nb}} has a regular NB distribution;
#' if \code{p_zinb} is in \code{(-prob^size/(1 - prob^size),} \code{0)}, \eqn{Y_{nb}} has a zero-deflated NB distribution and \code{p_zinb}
#' is not a probability; if \code{p_zinb = -prob^size/(1 - prob^size)}, \eqn{Y_{nb}} has a positive-NB distribution (see
#' \code{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 \code{M}, all NB variables in equation p
#' given \code{p_zinb[p]}; otherwise, missing values are set to 0 and ordered 1st
#' @param corr.x list of length \code{M}, each component a list of length \code{M}; \code{corr.x[[p]][[q]]} is matrix of correlations
#' for independent variables in equations p (\eqn{X_{(pj)}} for outcome \eqn{Y_p}) and q (\eqn{X_{(qj)}} for outcome \eqn{Y_q});
#' order: 1st ordinal (same order as in \code{marginal}), 2nd continuous non-mixture (same order as in \code{skews}),
#' 3rd components of continuous mixture (same order as in \code{mix_pis}), 4th regular Poisson, 5th zero-inflated Poisson (same order as in \code{lam}),
#' 6th regular NB, and 7th zero-inflated NB (same order as in \code{size});
#' if p = q, \code{corr.x[[p]][[q]]} is a correlation matrix with \code{nrow(corr.x[[p]][[q]])} = # \eqn{X_{(pj)}} for outcome \eqn{Y_p};
#' if p != q, \code{corr.x[[p]][[q]]} is a non-symmetric matrix of correlations where rows correspond to covariates for \eqn{Y_p}
#' so that \code{nrow(corr.x[[p]][[q]])} = # \eqn{X_{(pj)}} for outcome \eqn{Y_p} and
#' columns correspond to covariates for \eqn{Y_q} so that \code{ncol(corr.x[[p]][[q]])} = # \eqn{X_{(qj)}} for outcome \eqn{Y_q};
#' use \code{corr.x[[p]][[q]] = NULL} if equation q has no \eqn{X_{(qj)}}; use \code{corr.x[[p]] = NULL} if equation p has no \eqn{X_{(pj)}}
#' @param corr.e correlation matrix for continuous non-mixture or components of mixture error terms
#' @param U a list of length \code{M} of continuous non-mixture and components of mixture random effects
#' @param U_all a list of length \code{M} of continuous non-mixture and mixture random effects
#' @param 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 \code{M} containing a combination (i.e., \code{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)
#' @param 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 \code{M} as in \code{rand.int}
#' @param corr.u if the random effects are the same variables across equations, a matrix of correlations for \eqn{U};
#' if the random effects are different variables across equations, a list of length \code{M}, each component a list of length \code{M};
#' \code{corr.u[[p]][[q]]} is matrix of correlations for random effects in equations p (\eqn{U_{(pj)}} for outcome \eqn{Y_p}) and
#' q (\eqn{U_{(qj)}} for outcome \eqn{Y_q});
#' if p = q, \code{corr.u[[p]][[q]]} is a correlation matrix with \code{nrow(corr.u[[p]][[q]])} = # \eqn{U_{(pj)}} for outcome \eqn{Y_p};
#' if p != q, \code{corr.u[[p]][[q]]} is a non-symmetric matrix of correlations where rows correspond to \eqn{U_{(pj)}} for \eqn{Y_p}
#' so that \code{nrow(corr.u[[p]][[q]])} = # \eqn{U_{(pj)}} for outcome \eqn{Y_p} and
#' columns correspond to \eqn{U_{(qj)}} for \eqn{Y_q} so that \code{ncol(corr.u[[p]][[q]])} = # \eqn{U_{(qj)}} for outcome \eqn{Y_q};
#' the number of random effects for \eqn{Y_p} is taken from \code{nrow(corr.u[[p]][[1]])} so that if there should be random effects,
#' there must be entries for \code{corr.u};
#' use \code{corr.u[[p]][[q]] = NULL} if equation q has no \eqn{U_{(qj)}}; use \code{corr.u[[p]] = NULL} if equation p has no \eqn{U_{(pj)}};
#'
#' correlations are specified in terms of components of mixture variables (if present);
#' order is 1st random intercept (if \code{rand.int != "none"}), 2nd random time slope (if \code{rand.tsl != "none"}),
#' 3rd other random slopes with non-mixture distributions, 4th other random slopes with mixture distributions
#' @param rmeans2 a list returned from \code{corrsys} or \code{corrsys2} which has the non-mixture and component means ordered according to
#' types of random intercept and time slope
#' @param rvars2 a list returned like \code{rmeans}
#' @import SimMultiCorrData
#' @import SimCorrMix
#' @importFrom stats cor dbeta dbinom dchisq density dexp df dgamma dlnorm dlogis dmultinom dnbinom dnorm dpois dt dunif dweibull ecdf
#' median pbeta pbinom pchisq pexp pf pgamma plnorm plogis pnbinom pnorm ppois pt punif pweibull qbeta qbinom qchisq qexp qf qgamma
#' qlnorm qlogis qnbinom qnorm qpois qt quantile qunif qweibull rbeta rbinom rchisq rexp rf rgamma rlnorm rlogis rmultinom rnbinom
#' rnorm rpois rt runif rweibull sd uniroot var
#' @import utils
#' @importFrom Matrix nearPD
#' @importFrom VGAM qzipois qzinegbin
#' @export
#' @keywords summary
#' @seealso \code{\link[SimRepeat]{nonnormsys}}, \code{\link[SimRepeat]{corrsys}}, \code{\link[SimRepeat]{corrsys2}}
#'
#' @return A list with the following components:
#'
#' @return \code{cont_sum_y} a data.frame summarizing the simulated distributions of the \eqn{Y_p},
#'
#' @return \code{cont_sum_e} a data.frame summarizing the simulated distributions of the non-mixture or components of mixture \eqn{E_p},
#'
#' @return \code{target_sum_e} a data.frame summarizing the target distributions of the non-mixture or components of mixture \eqn{E_p},
#'
#' @return \code{mix_sum_e} a data.frame summarizing the simulated distributions of the mixture \eqn{E_p},
#'
#' @return \code{target_mix_e} a data.frame summarizing the target distributions of the mixture \eqn{E_p},
#'
#' @return \code{rho.y} correlation matrix of dimension \code{M x M} for \eqn{Y_p}
#'
#' @return \code{rho.e} correlation matrix for the non-mixture or components of mixture \eqn{E_p}
#'
#' @return \code{rho.emix} correlation matrix for the mixture \eqn{E_p}
#'
#' @return \code{rho.ye} matrix with correlations between \eqn{Y_p} (rows) and the non-mixture or components of mixture \eqn{E_p} (columns)
#'
#' @return \code{rho.yemix} matrix with correlations between \eqn{Y_p} (rows) and the mixture \eqn{E_p} (columns)
#'
#' @return \code{sum_xall} a data.frame summarizing \code{X_all} without the Time variable,
#'
#' @return \code{rho.yx} a list of length \code{M}, where \code{rho.yx[[p]]} is matrix of correlations
#' between \eqn{Y} (rows) and \code{X[[p]]} = \eqn{{X_ord(pj), X_cont(pj), X_comp(pj), X_pois(pj), X_nb(pj)}} (columns)
#'
#' @return \code{rho.yxall} a list of length \code{M}, where \code{rho.yx[[p]]} is matrix of correlations
#' between \eqn{Y} (rows) and \code{X_all[[p]]} (columns) not including \code{Time}
#'
#' @return \code{rho.x} a list of length \code{M} of lists of length \code{M} where
#' \code{rho.x[[p]][[q]] = cor(cbind(X[[p]], X[[q]]))} if p!= q or
#' \code{rho.x[[p]][[q]] = cor(X[[p]]))} if p = q, where \code{X[[p]]} = \eqn{{X_ord(pj), X_cont(pj), X_comp(pj), X_pois(pj), X_nb(pj)}}
#'
#' @return \code{rho.xall} a list of length \code{M} of lists of length \code{M} where
#' \code{rho.xall[[p]][[q]] = cor(cbind(X_all[[p]], X_all[[q]]))} if p!= q or
#' \code{rho.xall[[p]][[q]] = cor(X_all[[p]]))} if p = q, not including \code{Time}
#'
#' @return \code{maxerr} a list of length \code{M} containing a vector of length \code{M} with the maximum correlation errors between outcomes,
#' \code{maxerr[[p]]][q] = abs(max(corr.x[[p]][[q]] - rho.x[[p]][[q]]))}
#'
#' @return Additional components vary based on the type of simulated variables:
#'
#' @return If \bold{ordinal variables} are produced:
#' \code{ord_sum_x} a list where \code{ord_sum_x[[j]]} is a data.frame summarizing \eqn{X_{ord(pj)}} for all p = 1, ..., \code{M}
#'
#' @return If \bold{continuous variables} are produced:
#' \code{cont_sum_x} a data.frame summarizing the simulated distributions of the \eqn{X_{cont(pj)}} and \eqn{X_comp(pj)},
#'
#' \code{target_sum_x} a data.frame summarizing the target distributions of the \eqn{X_{cont(pj)}} and \eqn{X_comp(pj)},
#'
#' \code{mix_sum_x} a data.frame summarizing the simulated distributions of the \eqn{X_{mix(pj)}},
#'
#' \code{target_mix_x} a data.frame summarizing the target distributions of the \eqn{X_{mix(pj)}}
#'
#' @return If \bold{Poisson variables} are produced:
#' \code{pois_sum_x} a data.frame summarizing the simulated distributions of the \eqn{X_{pois(pj)}}
#'
#' @return If \bold{Negative Binomial variables} are produced:
#' \code{nb_sum_x} a data.frame summarizing the simulated distributions of the \eqn{X_{nb(pj)}}
#'
#' @return If \bold{random effects} are produced:
#' \code{cont_sum_u} a data.frame summarizing the simulated distributions of the \eqn{U_{cont(pj)}} and \eqn{U_{comp(pj)}},
#'
#' \code{target_sum_u} a data.frame summarizing the target distributions of the \eqn{U_{cont(pj)}} and \eqn{U_{comp(pj)}},
#'
#' \code{sum_uall} a data.frame summarizing the simulated distributions of \code{U_all},
#'
#' \code{mix_sum_u} a data.frame summarizing the simulated distributions of the \eqn{U_{mix(pj)}},
#'
#' \code{target_mix_u} a data.frame summarizing the target distributions of the \eqn{U_{mix(pj)}},
#'
#' \code{rho.u} list of length \code{M}, each component a list of length \code{M};
#' \code{rho.u[[p]][[q]] = cor(cbind(U[[p]], U[[q]]))} if p != q or \code{rho.u[[p]][[q]] = cor(U[[p]]))} if p = q
#'
#' \code{rho.uall} list of length \code{M}, each component a list of length \code{M};
#' \code{rho.uall[[p]][[q]] = cor(cbind(U_all[[p]], U_all[[q]]))} if p != q or \code{rho.uall[[p]][[q]] = cor(U_all[[p]]))} if p = q
#'
#' \code{maxerr_u} list of length \code{M} containing a vector of length \code{M} with the maximum correlation errors for \eqn{U} between outcomes
#' \code{maxerr_u[[p]]][q] = abs(max(corr.u[[p]][[q]] - rho.u[[p]][[q]]))}
#'
#' @examples
#' 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)
#'
#' \dontrun{
#' 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
#' }
#'
#' @references
#' See references for \code{\link[SimRepeat]{SimRepeat}}.
#'
summary_sys <- function(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()) {
error_type <- ifelse(is.null(E_mix) & is.null(E), "none",
ifelse(is.null(E_mix), "non_mix", "mix"))
if (length(X_all) == 0 & length(X) > 0) X_all <- X
K.cat <- rep(0, M)
K.pois <- rep(0, M)
K.nb <- rep(0, M)
if (length(marginal) > 0) K.cat <- lengths(marginal)
if (length(lam) > 0) K.pois <- lengths(lam)
if (length(size) > 0) K.nb <- lengths(size)
if (max(K.pois) > 0) {
p_zip0 <- p_zip
p_zip <- list()
if (class(p_zip0) == "numeric" & length(p_zip0) == 1) {
for (i in 1:M) {
if (K.pois[i] == 0) p_zip <- append(p_zip, list(NULL)) else
p_zip[[i]] <- rep(p_zip0, K.pois[i])
}
} else if (class(p_zip0) == "numeric" & length(p_zip0) == M) {
for (i in 1:M) {
if (K.pois[i] == 0) p_zip <- append(p_zip, list(NULL)) else
p_zip[[i]] <- rep(p_zip0[i], K.pois[i])
}
} else if (class(p_zip0) == "list" & length(p_zip0) == M) {
for (i in 1:M) {
if (K.pois[i] == 0) {
p_zip <- append(p_zip, list(NULL))
next
}
if (length(p_zip0[[i]]) == K.pois[i])
p_zip[[i]] <- p_zip0[[i]] else
p_zip[[i]] <- c(rep(0, K.pois[i] - length(p_zip0[[i]])),
p_zip0[[i]])
}
} else {
p_zip <- lapply(K.pois, function(x) if (x == 0) NULL else rep(0, x))
}
}
if (max(K.nb) > 0) {
if (length(prob) > 0)
mu <- mapply(function(x, y) x * (1 - y)/y, size, prob, SIMPLIFY = FALSE)
p_zinb0 <- p_zinb
p_zinb <- list()
if (class(p_zinb0) == "numeric" & length(p_zinb0) == 1) {
for (i in 1:M) {
if (K.nb[i] == 0) p_zinb <- append(p_zinb, list(NULL)) else
p_zinb[[i]] <- rep(p_zinb0, K.nb[i])
}
} else if (class(p_zinb0) == "numeric" & length(p_zinb0) == M) {
for (i in 1:M) {
if (K.nb[i] == 0) p_zinb <- append(p_zinb, list(NULL)) else
p_zinb[[i]] <- rep(p_zinb0[i], K.nb[i])
}
} else if (class(p_zinb0) == "list" & length(p_zinb0) == M) {
for (i in 1:M) {
if (K.nb[i] == 0) {
p_zinb <- append(p_zinb, list(NULL))
next
}
if (length(p_zinb0[[i]]) == K.nb[i])
p_zinb[[i]] <- p_zinb0[[i]] else
p_zinb[[i]] <- c(rep(0, K.nb[i] - length(p_zinb0[[i]])),
p_zinb0[[i]])
}
} else {
p_zinb <- lapply(K.nb, function(x) if (x == 0) NULL else rep(0, x))
}
}
K.r <- rep(0, M)
if (class(corr.u) == "list" & length(corr.u) > 0)
K.r <- sapply(mapply('[', corr.u, lengths(corr.u), SIMPLIFY = FALSE),
function(x) if (is.null(x)) 0 else nrow(x[[1]]))
if (class(corr.u) == "matrix") K.r <- rep(ncol(corr.u), M)
rmix_pis <- list()
if (length(means) %in% c(2 * M, M + 1)) {
means <- means[1:M]
vars <- vars[1:M]
}
if (length(skews) == (M + 1)) {
rskews <- skews[M + 1]
rskurts <- skurts[M + 1]
if (method == "Polynomial") {
rfifths <- fifths[M + 1]
rsixths <- sixths[M + 1]
}
}
if (length(skews) == (2 * M)) {
rskews <- skews[(M + 1):(2 * M)]
rskurts <- skurts[(M + 1):(2 * M)]
if (method == "Polynomial") {
rfifths <- fifths[(M + 1):(2 * M)]
rsixths <- sixths[(M + 1):(2 * M)]
}
}
if (length(skews) %in% c(2 * M, M + 1)) {
skews <- skews[1:M]
skurts <- skurts[1:M]
if (method == "Polynomial") {
fifths <- fifths[1:M]
sixths <- sixths[1:M]
}
}
if (length(mix_pis) == (M + 1)) {
rmix_pis <- mix_pis[M + 1]
rmix_mus <- mix_mus[M + 1]
rmix_sigmas <- mix_sigmas[M + 1]
rmix_skews <- mix_skews[M + 1]
rmix_skurts <- mix_skurts[M + 1]
if (method == "Polynomial") {
rmix_fifths <- mix_fifths[M + 1]
rmix_sixths <- mix_sixths[M + 1]
}
}
if (length(mix_pis) == (2 * M)) {
rmix_pis <- mix_pis[(M + 1):(2 * M)]
rmix_mus <- mix_mus[(M + 1):(2 * M)]
rmix_sigmas <- mix_sigmas[(M + 1):(2 * M)]
rmix_skews <- mix_skews[(M + 1):(2 * M)]
rmix_skurts <- mix_skurts[(M + 1):(2 * M)]
if (method == "Polynomial") {
rmix_fifths <- mix_fifths[(M + 1):(2 * M)]
rmix_sixths <- mix_sixths[(M + 1):(2 * M)]
}
}
if (length(mix_pis) %in% c(2 * M, M + 1)) {
mix_pis <- mix_pis[1:M]
mix_mus <- mix_mus[1:M]
mix_sigmas <- mix_sigmas[1:M]
mix_skews <- mix_skews[1:M]
mix_skurts <- mix_skurts[1:M]
if (method == "Polynomial") {
mix_fifths <- mix_fifths[1:M]
mix_sixths <- mix_sixths[1:M]
}
}
K.mix <- rep(0, M)
K.comp <- rep(0, M)
K.error <- rep(0, M)
if (length(mix_pis) > 0) {
K.mix <- lengths(mix_pis)
K.comp <- sapply(lapply(mix_pis, unlist), length)
k.comp <- c(0, cumsum(K.comp))
}
K.cont <- lengths(vars) - K.mix
k.cont <- c(0, cumsum(K.cont))
k.mix <- c(0, cumsum(K.mix))
if (error_type == "mix") {
K.error <- sapply(mix_pis, function(x) lengths(x[length(x)]))
k.error <- c(0, cumsum(K.error))
K.x <- K.cont + K.comp - K.error + K.cat + K.pois + K.nb
K.mix2 <- K.mix - 1
} else {
K.x <- K.cont + K.comp - 1 + K.cat + K.pois + K.nb
K.mix2 <- K.mix
}
if (length(marginal) > 0) {
if (length(support) == 0) {
for (m in 1:M) {
support <- append(support, list(NULL))
if (K.cat[m] > 0) support[[m]] <- lapply(marginal[[m]],
function(x) 1:(length(x) + 1))
}
} else {
for (m in 1:M) {
if (K.cat[m] > 0 & is.null(support[[m]])) {
support[[m]] <- lapply(marginal[[m]], function(x) 1:(length(x) + 1))
} else if (K.cat[m] > 0) {
for (i in 1:K.cat[m]) {
if (class(support[[m]]) != "list")
support[[m]] <- append(support[[m]], list(NULL))
if (length(support[[m]][[i]]) != (length(marginal[[m]][[i]]) + 1))
support[[m]][[i]] <- 1:(length(marginal[[m]][[i]]) + 1)
}
} else support <- append(support, list(NULL))
}
}
}
if (!is.null(Y)) {
rho.y <- cor(Y)
mom <- apply(Y, 2, calc_moments)
medians <- apply(Y, 2, median)
mins <- apply(Y, 2, min)
maxs <- apply(Y, 2, max)
cont_sum_y <- as.data.frame(cbind(1:M, rep(nrow(Y), M), mom[1, ], mom[2, ],
medians, mins, maxs, mom[3, ], mom[4, ], mom[5, ], mom[6, ]))
colnames(cont_sum_y) <- c("Outcome", "N", "Mean", "SD", "Median",
"Min", "Max", "Skew", "Skurtosis", "Fifth", "Sixth")
result <- list(cont_sum_y = cont_sum_y, rho.y = rho.y)
}
if (error_type != "none") {
mom <- apply(E, 2, calc_moments)
medians <- apply(E, 2, median)
mins <- apply(E, 2, min)
maxs <- apply(E, 2, max)
if (error_type == "non_mix") {
e.means <- mapply('[[', means, lengths(means))
e.vars <- mapply('[[', vars, lengths(vars))
cont_sum_e <- as.data.frame(cbind(1:M, rep(nrow(E), M), mom[1, ],
mom[2, ], medians, mins, maxs, mom[3, ], mom[4, ], mom[5, ], mom[6, ]))
colnames(cont_sum_e) <- c("Outcome", "N", "Mean", "SD", "Median",
"Min", "Max", "Skew", "Skurtosis", "Fifth", "Sixth")
if (method == "Fleishman") {
target_sum_e <- as.data.frame(cbind(1:M, e.means, sqrt(e.vars),
mapply('[[', skews, lengths(skews)),
mapply('[[', skurts, lengths(skurts))))
colnames(target_sum_e) <- c("Outcome", "Mean", "SD", "Skew",
"Skurtosis")
} else {
target_sum_e <- as.data.frame(cbind(1:M, e.means, sqrt(e.vars),
mapply('[[', skews, lengths(skews)),
mapply('[[', skurts, lengths(skurts)),
mapply('[[', fifths, lengths(fifths)),
mapply('[[', sixths, lengths(sixths))))
colnames(target_sum_e) <- c("Outcome", "Mean", "SD", "Skew",
"Skurtosis", "Fifth", "Sixth")
}
}
if (error_type == "mix") {
e.means <- unlist(mapply('[', mix_mus, lengths(mix_mus)))
e.vars <- (unlist(mapply('[', mix_sigmas, lengths(mix_sigmas))))^2
cont_sum_e <- as.data.frame(cbind(rep(1:M, K.error),
unlist(lapply(K.error, seq)), rep(nrow(E), ncol(mom)), mom[1, ],
mom[2, ], medians, mins, maxs, mom[3, ], mom[4, ], mom[5, ], mom[6, ]))
colnames(cont_sum_e) <- c("Outcome", "Component", "N", "Mean", "SD",
"Median", "Min", "Max", "Skew", "Skurtosis", "Fifth", "Sixth")
rho.emix <- cor(E_mix)
mom <- apply(E_mix, 2, calc_moments)
medians <- apply(E_mix, 2, median)
mins <- apply(E_mix, 2, min)
maxs <- apply(E_mix, 2, max)
mix_sum_e <- as.data.frame(cbind(1:M, rep(nrow(E_mix), M), mom[1, ],
mom[2, ], medians, mins, maxs, mom[3, ], mom[4, ], mom[5, ], mom[6, ]))
colnames(mix_sum_e) <- c("Outcome", "N", "Mean", "SD", "Median",
"Min", "Max", "Skew", "Skurtosis", "Fifth", "Sixth")
target_mix_e <- NULL
if (method == "Fleishman") {
target_sum_e <- as.data.frame(cbind(rep(1:M, K.error),
unlist(lapply(K.error, seq)), e.means, sqrt(e.vars),
unlist(mapply('[', mix_skews, lengths(mix_skews))),
unlist(mapply('[', mix_skurts, lengths(mix_skurts)))))
colnames(target_sum_e) <- c("Outcome", "Component", "Mean", "SD",
"Skew", "Skurtosis")
for (i in 1:M) {
target_mix_e <- rbind(target_mix_e,
calc_mixmoments(mix_pis[[i]][[length(mix_pis[[i]])]],
mix_mus[[i]][[length(mix_pis[[i]])]],
mix_sigmas[[i]][[length(mix_pis[[i]])]],
mix_skews[[i]][[length(mix_pis[[i]])]],
mix_skurts[[i]][[length(mix_pis[[i]])]]))
}
target_mix_e <- as.data.frame(cbind(1:M, target_mix_e))
colnames(target_mix_e) <- c("Outcome", "Mean", "SD", "Skew",
"Skurtosis")
} else {
target_sum_e <- as.data.frame(cbind(rep(1:M, K.error),
unlist(lapply(K.error, seq)), e.means, sqrt(e.vars),
unlist(mapply('[', mix_skews, lengths(mix_skews))),
unlist(mapply('[', mix_skurts, lengths(mix_skurts))),
unlist(mapply('[', mix_fifths, lengths(mix_fifths))),
unlist(mapply('[', mix_sixths, lengths(mix_sixths)))))
colnames(target_sum_e) <- c("Outcome", "Component", "Mean", "SD",
"Skew", "Skurtosis", "Fifth", "Sixth")
for (i in 1:M) {
target_mix_e <- rbind(target_mix_e,
calc_mixmoments(mix_pis[[i]][[length(mix_pis[[i]])]],
mix_mus[[i]][[length(mix_pis[[i]])]],
mix_sigmas[[i]][[length(mix_pis[[i]])]],
mix_skews[[i]][[length(mix_pis[[i]])]],
mix_skurts[[i]][[length(mix_pis[[i]])]],
mix_fifths[[i]][[length(mix_pis[[i]])]],
mix_sixths[[i]][[length(mix_pis[[i]])]]))
}
target_mix_e <- as.data.frame(cbind(1:M, target_mix_e))
colnames(target_mix_e) <- c("Outcome", "Mean", "SD", "Skew",
"Skurtosis", "Fifth", "Sixth")
}
rownames(target_mix_e) <- rownames(mix_sum_e)
}
rownames(target_sum_e) <- rownames(cont_sum_e)
rho.e <- cor(E)
result <- append(result, list(cont_sum_e = cont_sum_e,
target_sum_e = target_sum_e, rho.e = rho.e))
if (error_type == "mix") {
result <- append(result, list(mix_sum_e = mix_sum_e,
target_mix_e = target_mix_e, rho.emix = rho.emix))
}
if (!is.null(Y)) {
rho.ye <- cor(cbind(Y, E))[1:M, (M + 1):(M + ncol(E)), drop = FALSE]
result <- append(result, list(rho.ye = rho.ye))
if (error_type == "mix") {
rho.yemix <- cor(cbind(Y, E_mix))[1:M, (M + 1):(M + ncol(E_mix)),
drop = FALSE]
result <- append(result, list(rho.yemix = rho.yemix))
}
}
}
if (max(K.x) > 0) {
n <- max(unlist(lapply(X, nrow)))
for (i in 1:M) {
if (!is.null(X_all[[i]])) {
if (var(X_all[[i]][, ncol(X_all[[i]])]) == 0)
X_all[[i]] <- X_all[[i]][, -ncol(X_all[[i]]), drop = FALSE]
}
}
if (length(marginal) > 0) {
ord_sum_x <- list()
for (j in 1:max(K.cat)) {
ord_sum_x <- append(ord_sum_x, list(NULL))
for (i in 1:M) {
if (K.cat[i] == 0 | K.cat[i] < j) next
csum <- cumsum(table(X[[i]][, j]))/n
ord_sum_x[[j]] <- rbind(ord_sum_x[[j]],
cbind(rep(i, length(support[[i]][[j]])), support[[i]][[j]],
append(marginal[[i]][[j]], 1),
c(csum, rep(0, length(support[[i]][[j]]) - length(csum)))))
}
rownames(ord_sum_x[[j]]) <- 1:nrow(ord_sum_x[[j]])
ord_sum_x[[j]] <- as.data.frame(ord_sum_x[[j]])
colnames(ord_sum_x[[j]]) <- c("Outcome", "Support", "Target",
"Simulated")
names(ord_sum_x)[j] <- paste("O", j, sep = "")
}
result <- append(result, list(ord_sum_x = ord_sum_x))
}
sum_xall <- list()
for (i in 1:M) {
sum_xall <- append(sum_xall, list(NULL))
if (K.x[i] > 0) {
mom <- apply(X_all[[i]], 2, calc_moments)
medians <- apply(X_all[[i]], 2, median)
mins <- apply(X_all[[i]], 2, min)
maxs <- apply(X_all[[i]], 2, max)
sum_xall[[i]] <- cbind(rep(i, ncol(X_all[[i]])),
1:ncol(X_all[[i]]), rep(n, ncol(X_all[[i]])), mom[1, ], mom[2, ],
medians, mins, maxs, mom[3, ], mom[4, ], mom[5, ], mom[6, ])
rownames(sum_xall[[i]]) <- paste(colnames(X_all[[i]]), ":", i, sep = "")
}
}
sum_xall <- as.data.frame(do.call(rbind, sum_xall))
colnames(sum_xall) <- c("Outcome", "X", "N", "Mean", "SD", "Median",
"Min", "Max", "Skew", "Skurtosis", "Fifth", "Sixth")
if (max(K.x - K.cat - K.pois - K.nb) > 0) {
cont_sum_x <- list()
target_sum_x <- list()
mix_sum_x <- list()
target_mix_x <- list()
Y_comp <- list()
for (i in 1:M) {
cont_sum_x <- append(cont_sum_x, list(NULL))
target_sum_x <- append(target_sum_x, list(NULL))
mix_sum_x <- append(mix_sum_x, list(NULL))
target_mix_x <- append(target_mix_x, list(NULL))
Y_comp <- append(Y_comp, list(NULL))
if ((K.x[i] - K.cat[i] - K.pois[i] - K.nb[i]) > 0) {
Y_comp[[i]] <- X[[i]][, (K.cat[i] + 1):(ncol(X[[i]]) - K.pois[i] -
K.nb[i]), drop = FALSE]
mom <- apply(Y_comp[[i]], 2, calc_moments)
medians <- apply(Y_comp[[i]], 2, median)
mins <- apply(Y_comp[[i]], 2, min)
maxs <- apply(Y_comp[[i]], 2, max)
cont_sum_x[[i]] <- cbind(rep(i, ncol(Y_comp[[i]])),
1:ncol(Y_comp[[i]]), rep(n, ncol(Y_comp[[i]])), mom[1, ], mom[2, ],
medians, mins, maxs, mom[3, ], mom[4, ], mom[5, ], mom[6, ])
rownames(cont_sum_x[[i]]) <- paste("cont", rep(i, ncol(Y_comp[[i]])),
"_", 1:ncol(Y_comp[[i]]), sep = "")
ind <- ncol(Y_comp[[i]]) - K.comp[i] + K.error[i]
if (method == "Fleishman") {
if (ind > 0) {
target_sum_x[[i]] <- cbind(means[[i]][1:ind],
sqrt(vars[[i]][1:ind]), skews[[i]][1:ind],
skurts[[i]][1:ind])
}
if (K.mix2[i] > 0) {
target_sum_x[[i]] <- rbind(target_sum_x[[i]],
cbind(unlist(mix_mus[[i]][1:K.mix2[i]]),
unlist(mix_sigmas[[i]][1:K.mix2[i]]),
unlist(mix_skews[[i]][1:K.mix2[i]]),
unlist(mix_skurts[[i]][1:K.mix2[i]])))
}
} else {
if (ind > 0) {
target_sum_x[[i]] <- cbind(means[[i]][1:ind],
sqrt(vars[[i]][1:ind]), skews[[i]][1:ind], skurts[[i]][1:ind],
fifths[[i]][1:ind], sixths[[i]][1:ind])
}
if (K.mix2[i] > 0) {
target_sum_x[[i]] <- rbind(target_sum_x[[i]],
cbind(unlist(mix_mus[[i]][1:K.mix2[i]]),
unlist(mix_sigmas[[i]][1:K.mix2[i]]),
unlist(mix_skews[[i]][1:K.mix2[i]]),
unlist(mix_skurts[[i]][1:K.mix2[i]]),
unlist(mix_fifths[[i]][1:K.mix2[i]]),
unlist(mix_sixths[[i]][1:K.mix2[i]])))
}
}
target_sum_x[[i]] <- cbind(rep(i, ncol(Y_comp[[i]])),
1:ncol(Y_comp[[i]]), target_sum_x[[i]])
rownames(target_sum_x[[i]]) <- rownames(cont_sum_x[[i]])
if (K.mix2[i] > 0) {
mom <- apply(X_all[[i]][, (K.cat[i] + ind + 1):(K.cat[i] + ind +
K.mix2[i]), drop = FALSE], 2, calc_moments)
medians <- apply(X_all[[i]][, (K.cat[i] + ind + 1):(K.cat[i] + ind +
K.mix2[i]), drop = FALSE], 2, median)
mins <- apply(X_all[[i]][, (K.cat[i] + ind + 1):(K.cat[i] + ind +
K.mix2[i]), drop = FALSE], 2, min)
maxs <- apply(X_all[[i]][, (K.cat[i] + ind + 1):(K.cat[i] + ind +
K.mix2[i]), drop = FALSE], 2, max)
mix_sum_x[[i]] <- cbind(rep(i, K.mix2[i]), 1:K.mix2[i],
rep(n, K.mix2[i]), mom[1, ], mom[2, ],
medians, mins, maxs, mom[3, ], mom[4, ], mom[5, ], mom[6, ])
rownames(mix_sum_x[[i]]) <- paste("mix", rep(i, K.mix2[i]), "_",
1:K.mix2[i], sep = "")
if (method == "Fleishman") {
for (j in 1:K.mix2[i]) {
target_mix_x[[i]] <- rbind(target_mix_x[[i]],
calc_mixmoments(mix_pis[[i]][[j]], mix_mus[[i]][[j]],
mix_sigmas[[i]][[j]], mix_skews[[i]][[j]],
mix_skurts[[i]][[j]]))
}
} else {
for (j in 1:K.mix2[i]) {
target_mix_x[[i]] <- rbind(target_mix_x[[i]],
calc_mixmoments(mix_pis[[i]][[j]], mix_mus[[i]][[j]],
mix_sigmas[[i]][[j]], mix_skews[[i]][[j]],
mix_skurts[[i]][[j]], mix_fifths[[i]][[j]],
mix_sixths[[i]][[j]]))
}
}
target_mix_x[[i]] <- cbind(rep(i, K.mix2[i]), 1:K.mix2[i],
target_mix_x[[i]])
rownames(target_mix_x[[i]]) <- rownames(mix_sum_x[[i]])
}
}
}
cont_sum_x <- as.data.frame(do.call(rbind, cont_sum_x))
colnames(cont_sum_x) <- c("Outcome", "X", "N", "Mean", "SD", "Median",
"Min", "Max", "Skew", "Skurtosis", "Fifth", "Sixth")
target_sum_x <- as.data.frame(do.call(rbind, target_sum_x))
if (method == "Fleishman") {
colnames(target_sum_x) <- c("Outcome", "X", "Mean", "SD", "Skew",
"Skurtosis")
} else {
colnames(target_sum_x) <- c("Outcome", "X", "Mean", "SD", "Skew",
"Skurtosis", "Fifth", "Sixth")
}
result <- append(result, list(cont_sum_x = cont_sum_x,
target_sum_x = target_sum_x, sum_xall = sum_xall))
if (max(K.mix2) > 0) {
target_mix_x <- as.data.frame(do.call(rbind, target_mix_x))
colnames(target_mix_x) <- colnames(target_sum_x)
mix_sum_x <- as.data.frame(do.call(rbind, mix_sum_x))
colnames(mix_sum_x) <- colnames(cont_sum_x)
result <- append(result, list(mix_sum_x = mix_sum_x,
target_mix_x = target_mix_x))
}
}
if (length(lam) > 0) {
pois_sum_x <- list()
for (i in 1:M) {
if (K.pois[i] == 0) {
pois_sum_x <- append(pois_sum_x, list(NULL))
} else {
mom <- apply(X[[i]][, (ncol(X[[i]]) - K.pois[i] - K.nb[i] +
1):(ncol(X[[i]]) - K.nb[i]), drop = FALSE], 2, calc_moments)
medians <- apply(X[[i]][, (ncol(X[[i]]) - K.pois[i] - K.nb[i] +
1):(ncol(X[[i]]) - K.nb[i]), drop = FALSE], 2, median)
mins <- apply(X[[i]][, (ncol(X[[i]]) - K.pois[i] - K.nb[i] +
1):(ncol(X[[i]]) - K.nb[i]), drop = FALSE], 2, min)
maxs <- apply(X[[i]][, (ncol(X[[i]]) - K.pois[i] - K.nb[i] +
1):(ncol(X[[i]]) - K.nb[i]), drop = FALSE], 2, max)
p_0 <- apply(X[[i]][, (ncol(X[[i]]) - K.pois[i] - K.nb[i] +
1):(ncol(X[[i]]) - K.nb[i]), drop = FALSE], 2,
function(x) sum(x == 0)/n)
pois_sum_x[[i]] <- cbind(rep(i, K.pois[i]), 1:K.pois[i],
rep(n, K.pois[i]), p_0,
mapply(function(x, y) dzipois(0, x, y), lam[[i]], p_zip[[i]]),
mom[1, ], (1 - p_zip[[i]]) * lam[[i]], mom[2, ]^2,
lam[[i]] + (lam[[i]]^2) * p_zip[[i]]/(1 - p_zip[[i]]),
medians, mins, maxs, mom[3, ], mom[4, ])
rownames(pois_sum_x[[i]]) <- paste("pois", rep(i, K.pois[i]), "_",
1:K.pois[i], sep = "")
}
}
pois_sum_x <- as.data.frame(do.call(rbind, pois_sum_x))
colnames(pois_sum_x) <- c("Outcome", "X", "N", "P0", "Exp_P0", "Mean",
"Exp_Mean", "Var", "Exp_Var", "Median", "Min", "Max", "Skew",
"Skurtosis")
result <- append(result, list(pois_sum_x = pois_sum_x))
}
if (length(size) > 0) {
nb_sum_x <- list()
for (i in 1:M) {
if (K.nb[i] == 0) {
nb_sum_x <- append(nb_sum_x, list(NULL))
} else {
prob <- size[[i]]/(mu[[i]] + size[[i]])
p_0 <- apply(X[[i]][, (ncol(X[[i]]) - K.nb[i] + 1):ncol(X[[i]]),
drop = FALSE], 2, function(x) sum(x == 0)/n)
mom <- apply(X[[i]][, (ncol(X[[i]]) - K.nb[i] + 1):ncol(X[[i]]),
drop = FALSE], 2, calc_moments)
medians <- apply(X[[i]][, (ncol(X[[i]]) - K.nb[i] + 1):ncol(X[[i]]),
drop = FALSE], 2, median)
mins <- apply(X[[i]][, (ncol(X[[i]]) - K.nb[i] + 1):ncol(X[[i]]),
drop = FALSE], 2, min)
maxs <- apply(X[[i]][, (ncol(X[[i]]) - K.nb[i] + 1):ncol(X[[i]]),
drop = FALSE], 2, max)
nb_sum_x[[i]] <- cbind(rep(i, K.nb[i]), 1:K.nb[i], rep(n, K.nb[i]),
p_0, mapply(function(x, y, z) dzinegbin(0, x, munb = y, pstr0 = z),
size[[i]], mu[[i]], p_zinb[[i]]),
prob, mom[1, ], (1 - p_zinb[[i]]) * mu[[i]], mom[2, ]^2,
(1 - p_zinb[[i]]) * mu[[i]] * (1 + mu[[i]] *
(p_zinb[[i]] + 1/size[[i]])), medians, mins, maxs, mom[3, ],
mom[4, ])
rownames(nb_sum_x[[i]]) <- paste("nb", rep(i, K.nb[i]), "_",
1:K.nb[i], sep = "")
}
}
nb_sum_x <- as.data.frame(do.call(rbind, nb_sum_x))
colnames(nb_sum_x) <- c("Outcome", "X", "N", "P0", "Exp_P0", "Prob",
"Mean", "Exp_Mean", "Var", "Exp_Var", "Median", "Min", "Max", "Skew",
"Skurtosis")
result <- append(result, list(nb_sum_x = nb_sum_x))
}
rho.x <- list()
rho.xall <- list()
for (p in 1:M) {
if (K.x[p] == 0) {
rho.x <- append(rho.x, list(NULL))
rho.xall <- append(rho.xall, list(NULL))
next
}
rho.x[[p]] <- list()
rho.xall[[p]] <- list()
for (q in 1:M) {
if (K.x[q] == 0) {
rho.x[[p]] <- append(rho.x[[p]], list(NULL))
rho.xall[[p]] <- append(rho.xall[[p]], list(NULL))
next
}
if (q > p) {
rho.x[[p]][[q]] <- cor(cbind(X[[p]], X[[q]]))
rho.xall[[p]][[q]] <- cor(cbind(X_all[[p]], X_all[[q]]))
} else if (q == p) {
rho.x[[p]][[q]] <- cor(X[[p]])
rho.xall[[p]][[q]] <- cor(X_all[[p]])
} else {
rho.x[[p]][[q]] <- t(rho.x[[q]][[p]])
rho.xall[[p]][[q]] <- t(rho.xall[[q]][[p]])
}
}
}
result <- append(result, list(rho.x = rho.x, rho.xall = rho.xall))
if (!is.null(Y)) {
rho.yx <- list()
rho.yxall <- list()
for (p in 1:M) {
if (K.x[p] == 0) {
rho.yx <- append(rho.yx, list(NULL))
rho.yxall <- append(rho.yxall, list(NULL))
next
}
rho.yx[[p]] <- cor(cbind(Y, X[[p]]))[1:M, (M + 1):(M + ncol(X[[p]])),
drop = FALSE]
rho.yxall[[p]] <-
cor(cbind(Y, X_all[[p]]))[1:M, (M + 1):(M + ncol(X_all[[p]])),
drop = FALSE]
}
result <- append(result, list(rho.yx = rho.yx, rho.yxall = rho.yxall))
}
}
if (length(corr.x) > 0) {
emax2 <- list()
for (p in 1:M) {
if (K.x[p] == 0) {
emax2 <- append(emax2, list(NULL))
next
}
emax2[[p]] <- numeric(M)
for (q in 1:M) {
if (K.x[q] == 0) {
emax2[[p]][q] <- NA
next
}
if (q > p) {
emax2[[p]][q] <- max(abs(rho.x[[p]][[q]][1:K.x[p],
(K.x[p] + 1):(K.x[p] + K.x[q])] - corr.x[[p]][[q]]))
} else if (q == p) {
emax2[[p]][q] <- max(abs(rho.x[[p]][[q]] - corr.x[[p]][[q]]))
} else {
emax2[[p]][q] <- emax2[[q]][p]
}
}
}
result <- append(result, list(maxerr = emax2))
}
if (length(rmeans2) > 0) {
if (length(rand.int) != M) rand.int <- rep(rand.int[1], M)
if (length(rand.tsl) != M) rand.tsl <- rep(rand.tsl[1], M)
M0 <- ifelse(class(corr.u) == "list", M, 1)
K.rmix <- rep(0, M0)
K.rcomp <- rep(0, M0)
if (length(rmix_pis) > 0) {
K.rmix <- lengths(rmix_pis)
K.rcmix <- unlist(lapply(rmix_pis, lengths))
K.rcomp <- sapply(lapply(rmix_pis, unlist), length)
k.rcomp <- c(0, cumsum(K.rcomp))
}
K.rcont <- sapply(rvars2,
function(x) if (is.null(x[[1]])) 0 else length(x)) - K.rcomp
k.rcont <- c(0, cumsum(K.rcont))
k.rmix <- c(0, cumsum(K.rmix))
rskews2 <- list()
rskurts2 <- list()
if (method == "Polynomial") {
rfifths2 <- list()
rsixths2 <- list()
}
for (i in 1:M0) {
rskews2[[i]] <- append(rskews2, list(NULL))
rskurts2[[i]] <- append(rskurts2, list(NULL))
if (method == "Polynomial") {
rfifths2[[i]] <- append(rfifths2, list(NULL))
rsixths2[[i]] <- append(rsixths2, list(NULL))
}
if (K.r[i] == 0) next
if (K.rmix[i] > 0) ind <- c(0, cumsum(lengths(rmix_pis[[i]])))
if ((rand.int[i] == "none" & rand.tsl[i] == "mix") |
(rand.int[i] == "mix" & rand.tsl[i] == "none")) {
if (K.rcont[i] > 0) {
rskews2[[i]] <- c(rmix_skews[[i]][[1]], rskews[[i]],
unlist(rmix_skews[[i]][-1]))
rskurts2[[i]] <- c(rmix_skurts[[i]][[1]], rskurts[[i]],
unlist(rmix_skurts[[i]][-1]))
if (method == "Polynomial") {
rfifths2[[i]] <- c(rmix_fifths[[i]][[1]], rfifths[[i]],
unlist(rmix_fifths[[i]][-1]))
rsixths2[[i]] <- c(rmix_sixths[[i]][[1]], rsixths[[i]],
unlist(rmix_sixths[[i]][-1]))
}
} else {
rskews2[[i]] <- unlist(rmix_skews[[i]])
rskurts2[[i]] <- unlist(rmix_skurts[[i]])
if (method == "Polynomial") {
rfifths2[[i]] <- unlist(rmix_fifths[[i]])
rsixths2[[i]] <- unlist(rmix_sixths[[i]])
}
}
} else if (rand.int[i] == "non_mix" & rand.tsl[i] == "mix") {
if (K.rcont[i] > 1) {
rskews2[[i]] <- c(rskews[[i]][1], rmix_skews[[i]][[1]],
rskews[[i]][-1], unlist(rmix_skews[[i]][-1]))
rskurts2[[i]] <- c(rskurts[[i]][1], rmix_skurts[[i]][[1]],
rskurts[[i]][-1], unlist(rmix_skurts[[i]][-1]))
if (method == "Polynomial") {
rfifths2[[i]] <- c(rfifths[[i]][1], rmix_fifths[[i]][[1]],
rfifths[[i]][-1], unlist(rmix_fifths[[i]][-1]))
rsixths2[[i]] <- c(rsixths[[i]][1], rmix_sixths[[i]][[1]],
rsixths[[i]][-1], unlist(rmix_sixths[[i]][-1]))
}
} else {
rskews2[[i]] <- c(rskews[[i]][1], unlist(rmix_skews[[i]]))
rskurts2[[i]] <- c(rskurts[[i]][1], unlist(rmix_skurts[[i]]))
if (method == "Polynomial") {
rfifths2[[i]] <- c(rfifths[[i]][1], unlist(rmix_fifths[[i]]))
rsixths2[[i]] <- c(rsixths[[i]][1], unlist(rmix_sixths[[i]]))
}
}
} else if (rand.int[i] == "mix" & rand.tsl[i] == "non_mix") {
rskews2[[i]] <- c(rmix_skews[[i]][[1]], rskews[[i]],
unlist(rmix_skews[[i]][-1]))
rskurts2[[i]] <- c(rmix_skurts[[i]][[1]], rskurts[[i]],
unlist(rmix_skurts[[i]][-1]))
if (method == "Polynomial") {
rfifths2[[i]] <- c(rmix_fifths[[i]][[1]], rfifths[[i]],
unlist(rmix_fifths[[i]][-1]))
rsixths2[[i]] <- c(rmix_sixths[[i]][[1]], rsixths[[i]],
unlist(rmix_sixths[[i]][-1]))
}
} else if (rand.int[i] == "mix" & rand.tsl[i] == "mix") {
if (K.rcont[i] > 0) {
rskews2[[i]] <- c(unlist(rmix_skews[[i]][1:2]),
rskews[[i]], unlist(rmix_skews[[i]][-c(1:2)]))
rskurts2[[i]] <- c(unlist(rmix_skurts[[i]][1:2]),
rskurts[[i]], unlist(rmix_skurts[[i]][-c(1:2)]))
if (method == "Polynomial") {
rfifths2[[i]] <- c(unlist(rmix_fifths[[i]][1:2]),
rfifths[[i]], unlist(rmix_fifths[[i]][-c(1:2)]))
rsixths2[[i]] <- c(unlist(rmix_sixths[[i]][1:2]),
rsixths[[i]], unlist(rmix_sixths[[i]][-c(1:2)]))
}
} else {
rskews2[[i]] <- unlist(rmix_skews[[i]])
rskurts2[[i]] <- unlist(rmix_skurts[[i]])
if (method == "Polynomial") {
rfifths2[[i]] <- unlist(rmix_fifths[[i]])
rsixths2[[i]] <- unlist(rmix_sixths[[i]])
}
}
} else {
if (K.rcont[i] > 0) {
rskews2[[i]] <- rskews[[i]]
rskurts2[[i]] <- rskurts[[i]]
if (method == "Polynomial") {
rfifths2[[i]] <- rfifths[[i]]
rsixths2[[i]] <- rsixths[[i]]
}
}
if (K.rmix[i] > 0) {
rskews2[[i]] <- c(rskews2[[i]], unlist(rmix_skews[[i]]))
rskurts2[[i]] <- c(rskurts2[[i]], unlist(rmix_skurts[[i]]))
if (method == "Polynomial") {
rfifths2[[i]] <- c(rfifths2[[i]], unlist(rmix_fifths[[i]]))
rsixths2[[i]] <- c(rsixths2[[i]], unlist(rmix_sixths[[i]]))
}
}
}
}
if (max(K.rmix) > 0) {
U_mix <- list()
for (i in 1:M0) {
if (K.rmix[i] == 0) {
U_mix <- append(U_mix, list(NULL))
next
}
if ((rand.int[i] == "none" & rand.tsl[i] == "mix") |
(rand.int[i] == "mix" & rand.tsl[i] == "none") |
(rand.int[i] == "mix" & rand.tsl[i] == "non_mix")) {
U_mix[[i]] <- U_all[[i]][, 1, drop = FALSE]
if (K.rmix[i] > 1)
U_mix[[i]] <- cbind(U_mix[[i]],
U_all[[i]][, (ncol(U_all[[i]]) - K.rmix[i] +
2):ncol(U_all[[i]]), drop = FALSE])
} else if (rand.int[i] == "non_mix" & rand.tsl[i] == "mix") {
U_mix[[i]] <- U_all[[i]][, 2, drop = FALSE]
if (K.rmix[i] > 1)
U_mix[[i]] <- cbind(U_mix[[i]],
U_all[[i]][, (ncol(U_all[[i]]) - K.rmix[i] +
2):ncol(U_all[[i]])])
} else if (rand.int[i] == "mix" & rand.tsl[i] == "mix") {
U_mix[[i]] <- U_all[[i]][, 1:2]
if (K.rmix[i] > 2)
U_mix[[i]] <- cbind(U_mix[[i]],
U_all[[i]][, (ncol(U_all[[i]]) - K.rmix[i] +
3):ncol(U_all[[i]]), drop = FALSE])
} else {
U_mix[[i]] <- U_all[[i]][, (ncol(U_all[[i]]) - K.rmix[i] +
1):ncol(U_all[[i]]), drop = FALSE]
}
}
}
rho.u <- list()
for (i in 1:M0) {
if (K.r[i] == 0) {
rho.u <- append(rho.u, list(NULL))
next
}
rho.u[[i]] <- list()
for (j in 1:M0) {
if (K.r[j] == 0) {
rho.u[[i]] <- append(rho.u[[i]], list(NULL))
next
}
if (j > i) {
rhou <- cor(cbind(U[[i]], U[[j]]))
rho.u[[i]][[j]] <- rhou[1:K.r[i], (K.r[i] + 1):ncol(rhou),
drop = FALSE]
} else if (j == i) {
rho.u[[i]][[j]] <- cor(U[[i]])
} else {
rho.u[[i]][[j]] <- t(rho.u[[j]][[i]])
}
}
}
cont_sum_u <- list()
target_sum_u <- list()
mix_sum_u <- list()
target_mix_u <- list()
sum_uall <- list()
for (i in 1:M0) {
sum_uall <- append(sum_uall, list(NULL))
cont_sum_u <- append(cont_sum_u, list(NULL))
target_sum_u <- append(target_sum_u, list(NULL))
mix_sum_u <- append(mix_sum_u, list(NULL))
target_mix_u <- append(target_mix_u, list(NULL))
if (K.r[i] > 0) {
mom <- apply(U_all[[i]], 2, calc_moments)
medians <- apply(U_all[[i]], 2, median)
mins <- apply(U_all[[i]], 2, min)
maxs <- apply(U_all[[i]], 2, max)
sum_uall[[i]] <- cbind(rep(i, ncol(U_all[[i]])), 1:ncol(U_all[[i]]),
rep(nrow(U_all[[i]]), ncol(U_all[[i]])), mom[1, ], mom[2, ], medians,
mins, maxs, mom[3, ], mom[4, ], mom[5, ], mom[6, ])
rownames(sum_uall[[i]]) <- colnames(U_all[[i]])
mom <- apply(U[[i]], 2, calc_moments)
medians <- apply(U[[i]], 2, median)
mins <- apply(U[[i]], 2, min)
maxs <- apply(U[[i]], 2, max)
cont_sum_u[[i]] <- cbind(rep(i, K.r[i]), 1:K.r[i],
rep(nrow(U[[i]]), ncol(U[[i]])), mom[1, ], mom[2, ], medians, mins,
maxs, mom[3, ], mom[4, ], mom[5, ], mom[6, ])
rownames(cont_sum_u[[i]]) <- paste("cont", rep(i, ncol(U[[i]])),
"_", 1:ncol(U[[i]]), sep = "")
if (method == "Fleishman") {
target_sum_u[[i]] <- cbind(rmeans2[[i]], sqrt(rvars2[[i]]),
rskews2[[i]], rskurts2[[i]])
} else {
target_sum_u[[i]] <- cbind(rmeans2[[i]], sqrt(rvars2[[i]]),
rskews2[[i]], rskurts2[[i]], rfifths2[[i]], rsixths2[[i]])
}
target_sum_u[[i]] <- cbind(rep(i, K.r[i]), 1:K.r[i], target_sum_u[[i]])
rownames(target_sum_u[[i]]) <- rownames(cont_sum_u[[i]])
if (K.rmix[i] > 0) {
mom <- apply(U_mix[[i]], 2, calc_moments)
medians <- apply(U_mix[[i]], 2, median)
mins <- apply(U_mix[[i]], 2, min)
maxs <- apply(U_mix[[i]], 2, max)
mix_sum_u[[i]] <- cbind(rep(i, K.rmix[i]), 1:K.rmix[i],
rep(nrow(U_mix[[i]]), ncol(U_mix[[i]])), mom[1, ], mom[2, ],
medians, mins, maxs, mom[3, ], mom[4, ], mom[5, ], mom[6, ])
rownames(mix_sum_u[[i]]) <- paste("mix", rep(i, K.rmix[i]),
"_", 1:K.rmix[i], sep = "")
if (method == "Fleishman") {
for (j in 1:K.rmix[i]) {
target_mix_u[[i]] <- rbind(target_mix_u[[i]],
calc_mixmoments(rmix_pis[[i]][[j]], rmix_mus[[i]][[j]],
rmix_sigmas[[i]][[j]], rmix_skews[[i]][[j]],
rmix_skurts[[i]][[j]]))
}
} else {
for (j in 1:K.rmix[i]) {
target_mix_u[[i]] <- rbind(target_mix_u[[i]],
calc_mixmoments(rmix_pis[[i]][[j]], rmix_mus[[i]][[j]],
rmix_sigmas[[i]][[j]], rmix_skews[[i]][[j]],
rmix_skurts[[i]][[j]], rmix_fifths[[i]][[j]],
rmix_sixths[[i]][[j]]))
}
}
target_mix_u[[i]] <- cbind(rep(i, K.rmix[i]), 1:K.rmix[i],
target_mix_u[[i]])
rownames(target_mix_u[[i]]) <- rownames(mix_sum_u[[i]])
}
}
}
cont_sum_u <- as.data.frame(do.call(rbind, cont_sum_u))
colnames(cont_sum_u) <- c("Outcome", "U", "N", "Mean", "SD", "Median",
"Min", "Max", "Skew", "Skurtosis", "Fifth", "Sixth")
sum_uall <- as.data.frame(do.call(rbind, sum_uall))
colnames(sum_uall) <- c("Outcome", "U", "N", "Mean", "SD", "Median",
"Min", "Max", "Skew", "Skurtosis", "Fifth", "Sixth")
target_sum_u <- as.data.frame(do.call(rbind, target_sum_u))
if (method == "Fleishman") {
colnames(target_sum_u) <- c("Outcome", "U", "Mean", "SD", "Skew",
"Skurtosis")
} else {
colnames(target_sum_u) <- c("Outcome", "U", "Mean", "SD", "Skew",
"Skurtosis", "Fifth", "Sixth")
}
result <- append(result, list(target_sum_u = target_sum_u,
cont_sum_u = cont_sum_u, sum_uall = sum_uall, rho.u = rho.u))
if (max(K.rmix) > 0) {
rho.uall <- list()
for (i in 1:M0) {
if (K.rmix[i] == 0) {
rho.uall <- append(rho.uall, list(NULL))
next
}
rho.uall[[i]] <- list()
for (j in 1:M0) {
if (K.rmix[j] == 0) {
rho.uall[[i]] <- append(rho.uall[[i]], list(NULL))
next
}
if (j >= i) {
rhou <- cor(cbind(U_all[[i]], U_all[[j]]))
if (j == i) rho.uall[[i]][[j]] <-
rhou[1:ncol(U_all[[i]]), 1:ncol(U_all[[i]]), drop = FALSE]
if (j > i) rho.uall[[i]][[j]] <-
rhou[1:ncol(U_all[[i]]), (ncol(U_all[[i]]) + 1):ncol(rhou),
drop = FALSE]
}
if (j < i) rho.uall[[i]][[j]] <- t(rho.uall[[j]][[i]])
}
}
target_mix_u <- as.data.frame(do.call(rbind, target_mix_u))
colnames(target_mix_u) <- colnames(target_sum_u)
mix_sum_u <- as.data.frame(do.call(rbind, mix_sum_u))
colnames(mix_sum_u) <- colnames(cont_sum_u)
result <- append(result, list(rho.uall = rho.uall,
target_mix_u = target_mix_u, mix_sum_u = mix_sum_u))
}
}
if (class(corr.u) == "matrix") {
emax3 <- max(abs(rho.u[[1]][[1]] - corr.u))
result <- append(result, list(maxerr_u = emax3))
}
if (class(corr.u) == "list" & length(corr.u) > 0) {
emax3 <- list()
for (i in 1:M) {
if (K.r[i] == 0) {
emax3 <- append(emax3, list(NULL))
next
}
emax3[[i]] <- numeric(M)
for (j in 1:M) {
if (K.r[j] == 0) {
emax3[[i]][j] <- NA
next
}
if (j > i) {
emax3[[i]][j] <- max(abs(rho.u[[i]][[j]] - corr.u[[i]][[j]]))
} else if (j == i) {
emax3[[i]][j] <- max(abs(rho.u[[i]][[j]] - corr.u[[i]][[j]]))
} else {
emax3[[i]][j] <- emax3[[j]][i]
}
}
}
result <- append(result, list(maxerr_u = emax3))
}
result
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.