Nothing
#' Simulate data based on a model and parameter distributions
#'
#' @param design a design dataset. See example
#' @param model A function with the first argument the simulation design, i.e. a dataset with the columns ... The second argument to this function is a dataset with parameters for every individual. This can be supplied by the user, or generated by this sim_data if theta and omega_mat are supplied.
#' @param theta vector of fixed effect parameters
#' @param omega_mat vector of between subject random effects, specified as lower triangle
#' @param par_names A character vector linking the parameters in the model to the variables in the dataset. See example.
#' @param par_values parameter values
#' @param draw_iiv draw between subject random effects?
#' @param error see example
#' @param n number of simulations to perform
#' @return a vector of simulated dependent variables (for us in the VPC plotting function)
#' @family aggregate functions
#' @seealso \code{\link{vpc}}
#' @export
#' @details
#' This function generates the simulated dependent values for use in the VPC plotting function.
sim_data <- function (design = cbind(id = c(1,1,1), idv = c(0,1,2)),
model = function(x) { return(x$alpha + x$beta) },
theta,
omega_mat,
par_names,
par_values = NULL,
draw_iiv = "mvrnorm",
error = list(proportional = 0, additive = 0, exponential = 0),
n=100) {
if (is.null(par_values)) {
param <- draw_params_mvr( # draw parameter values. can also be just population values, or specified manually ()
ids = 1:n,
n_sim = n,
theta,
omega_mat = triangle_to_full(omega_mat),
par_names = par_names)
} else {
param <- par_values
}
sim_des <- do.call("rbind", replicate(n, design, simplify = FALSE))
sim_des$sim <- rep(1:n, each=nrow(design[,1]))
sim_des$join <- paste(sim_des$sim, sim_des$id, sep="_")
param$join <- paste(param$sim, param$id, sep="_")
tmp <- dplyr::as_tibble(merge(sim_des, param,
by.x="join", by.y="join"))
tmp_pred <- cbind(data.frame(design), matrix(rep(theta, each=nrow(design[,1])), ncol=length(theta)))
colnames(tmp_pred)[length(tmp_pred)-length(par_names)+1:3] <- par_names
tmp$dv <- add_noise(model(tmp), ruv = error)
tmp$pred <- rep(model(tmp_pred), n)
colnames(tmp) <- gsub("\\.x", "", colnames(tmp))
tmp %>%
dplyr::arrange_("sim", "id", "time")
}
sim_data_tte <- function (fit, t_cens = NULL, n = 100) {
fit$coefficients <- as.list(fit$coefficients)
dat <- data.frame(model.matrix(fit))
for (i in seq(fit$coefficients)) { fit$coefficients[[i]] <- as.numeric (fit$coefficients[[i]]) }
fact <- as.matrix(attr(fit$terms, "factors"))
parm <- t(fact) %*% as.numeric(fit$coefficients)
tmp.single <- data.frame (
par = exp(as.numeric(fit$coefficients[1]) + as.matrix(dat[,rownames(parm)]) %*% parm),
dv = 1
)
tmp <- do.call("rbind", replicate(n, tmp.single, simplify = FALSE))
tmp$sim <- rep(1:n, each=nrow(tmp.single[,1]))
if (!fit$dist %in% c("exponential", "weibull")) {
cat (paste("Simulation of ", fit$dist, "distribution not yet implemented, sorry."))
return()
}
if (fit$dist == "exponential") {
tmp$t = rweibull(nrow(dat[,1]) * n, shape = 1, scale = tmp$par)
# or using: tmp$t = rexp(length(design$id), 1/tmp$par)
}
if (fit$dist == "weibull") {
# annoyinly, the survreg and rweibull mix up the shape/scale parameter names and also take the inverse!!!
tmp$t = rweibull(nrow(dat[,1]) * n, shape = 1/fit$scale, scale = tmp$par)
}
if (sum(tmp$t > t_cens) > 0) {
tmp[tmp$t > t_cens,]$t <- t_cens
}
out <- c()
for (i in 1:n) {
km_fit <- compute_kaplan(tmp[tmp$sim == i,])
idx_new <- idx + length(unique(tmp[tmp$sim == i,]$t))-1
out <- rbind(out, cbind(i, km_fit$time, km_fit$surv))
}
colnames(out) <- c("sim", "time", "dv")
return(dplyr::as_tibble(data.frame(out)))
}
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.