# 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 {
# }
Run the code above in your browser using DataLab