R/data_read_selection_long.R
In AnGabrio/missingHE: Missing Outcome Data in Health Economic Evaluation
Defines functions
data_read_selection_long
Documented in
data_read_selection_long
#' A function to read and re-arrange the data in different ways
#'
#' This internal function imports the data and outputs only those variables that are needed to run the model
#' according to the information provided by the user.
#' @param data A data frame in which to find variables supplied in \code{model.eff} and \code{model.cost}. Among these,
#' effectiveness, cost and treatment indicator (only two arms) variables must always be provided and named 'u', 'c' and 't' respectively.
#' @param model.eff A formula expression in conventional \code{R} linear modelling syntax. The response must be a health economics
#' effectiveness outcome ('u') whose name must correspond to that used in \code{data}, and
#' any covariates are given on the right-hand side. If there are no covariates, specify \code{1} on the right hand side.
#' By default, covariates are placed on the "location" parameter of the distribution through a linear model.
#' Random effects can also be specified for each model parameter.
#' @param model.cost A formula expression in conventional \code{R} linear modelling syntax. The response must be a health economics
#' cost outcome ('c') whose name must correspond to that used in \code{data}, and any covariates are given on the right-hand side.
#' If there are no covariates, specify \code{1} on the right hand side. By default, covariates are placed on the "location"
#' parameter of the distribution through a linear model. Random effects can also be specified for each model parameter.
#' @param model.mu A formula expression in conventional \code{R} linear modelling syntax. The response must be indicated with the
#' term 'mu'(missing effects) and any covariates used to estimate the probability of missing effects are given on the right-hand side.
#' If there are no covariates, specify \code{1} on the right hand side. By default, covariates are placed on the "probability" parameter for the missing effects through a logistic-linear model.
#' Random effects can also be specified for each model parameter.
#' @param model.mc A formula expression in conventional \code{R} linear modelling syntax. The response must be indicated with the
#' term 'mc'(missing costs) and any covariates used to estimate the probability of missing costs should be given on the right-hand side.
#' If there are no covariates, specify \code{1} on the right hand side. By default, covariates are placed on the "probability" parameter for the missing costs through a logistic-linear model.
#' Random effects can also be specified for each model parameter.
#' @param type Type of missingness mechanism assumed. Choices are Missing At Random (MAR) and Missing Not At Random (MNAR).
#' @param center Logical. If \code{center} is \code{TRUE} all the covariates in the model are centered.
#' @keywords read data
#' @importFrom stats na.omit sd as.formula model.matrix model.frame model.response terms aggregate reshape
#' @export
#' @examples
#' #Internal function only
#' #no examples
#' #
#' #
data_read_selection_long <- function(data, model.eff, model.cost, model.mu, model.mc, type, center) {
if(is.data.frame(data) == FALSE) {
stop("object data must be provided as data frame")
}
if(any(names(data) == "u") == TRUE & any(names(data) == "c") == TRUE) {
u <- as.name("u")
c <- as.name("c")
}
data.wide <- reshape(data, v.names = c("u", "c"), timevar = c("time"),
direction = "wide", idvar = "id", sep = ".")
data.wide_u <- data.wide[, grepl('u\\.|id', colnames(data.wide))]
wide_u_index <- grepl('u\\.|id', colnames(data.wide_u))
data.wide_u <- data.wide_u[, wide_u_index]
data.wide_u_obs_index <- as.data.frame(which(is.na(data.wide_u) == FALSE, arr.ind = TRUE))
data.wide_u_obs_index <- merge(data.wide_u_obs_index, aggregate(col ~ row, data.wide_u_obs_index, max))
data.wide_u_obs_index <- unique(data.wide_u_obs_index)
names(data.wide_u_obs_index) <- c("id", "drop_u")
data.wide_u_drop <- merge(data.wide, data.wide_u_obs_index, by.x = "id", all.x = TRUE)
data.wide_u_drop$drop_u <- data.wide_u_drop$drop_u - 1
data.wide_c <- data.wide[, grepl('c\\.|id', colnames(data.wide))]
wide_c_index <- grepl('c\\.|id', colnames(data.wide_c))
data.wide_c <- data.wide_c[, wide_c_index]
data.wide_c_obs_index <- as.data.frame(which(is.na(data.wide_c) == FALSE, arr.ind = TRUE))
data.wide_c_obs_index <- merge(data.wide_c_obs_index, aggregate(col ~ row, data.wide_c_obs_index, max))
data.wide_c_obs_index <- unique(data.wide_c_obs_index)
names(data.wide_c_obs_index) <- c("id", "drop_c")
data.wide_uc_drop <- merge(data.wide_u_drop, data.wide_c_obs_index, by.x = "id", all.x = TRUE)
data.wide_uc_drop$drop_c <- data.wide_uc_drop$drop_c - 1
data.wide_uc_drop <- data.wide_uc_drop[, which(names(data.wide_uc_drop) %in% c("id","drop_u","drop_c"))]
data.wide_drop <- merge(data.wide, data.wide_uc_drop, by = "id", all = TRUE)
u_index_names <- colnames(data.wide_drop)[grepl('u\\.', colnames(data.wide_drop))]
c_index_names <- colnames(data.wide_drop)[grepl('c\\.', colnames(data.wide_drop))]
data.long <- reshape(data.wide_drop, varying = c(u_index_names,c_index_names),
direction = "long", idvar = "id", sep = ".")
time_max <- max(data.long$time)
count_mis_u <- rep(NA, length(data.long$drop_u))
count_mis_u <- ifelse(is.na(data.long$u) == FALSE, 1, count_mis_u)
count_mis_u <- ifelse(data.long$drop_u == 1 & is.na(data.long$u) == TRUE, 3, count_mis_u)
count_mis_u <- ifelse(data.long$drop_u <= data.long$time & is.na(data.long$u) == TRUE, 3, count_mis_u)
count_mis_u <- ifelse(data.long$drop_u > data.long$time & is.na(data.long$u) == TRUE, 2, count_mis_u)
count_mis_c <- rep(NA, length(data.long$drop_c))
count_mis_c <- ifelse(is.na(data.long$c) == FALSE, 1, count_mis_c)
count_mis_c <- ifelse(data.long$drop_c == 1 & is.na(data.long$c) == TRUE, 3, count_mis_c)
count_mis_c <- ifelse(data.long$drop_c <= data.long$time & is.na(data.long$c) == TRUE, 3, count_mis_c)
count_mis_c <- ifelse(data.long$drop_c > data.long$time & is.na(data.long$c) == TRUE, 2, count_mis_c)
data <- data.long
my_subset <- function() {
drop_u <- drop_c <- NULL
cov_matrix <- subset(data, select = -c(u, c, drop_u, drop_c, time, t))[data$time == 1,]
cov_matrix
}
cov_matrix <- my_subset()
cov_matrix <- cov_matrix[!unlist(vapply(cov_matrix, anyNA, logical(1)))]
is.formula <- function (x) { inherits(x, "formula") }
if(is.formula(model.eff) == FALSE | is.formula(model.cost) == FALSE) {
stop("model.eff and/or model.cost must be formula objects")
}
if(is.logical(center) == FALSE) { stop("center must be either TRUE or FALSE") }
fixed_u <- nobars_(model.eff)
fixed_c <- nobars_(model.cost)
random_u <- fb(model.eff)
random_c <- fb(model.cost)
fname_re_u_coeff <- as.formula(paste("u", "0", sep=" ~ "))
fname_re_c_coeff <- as.formula(paste("c", "0", sep=" ~ "))
fname_re_mu_coeff <- as.formula(paste("mu", "0", sep=" ~ "))
fname_re_mc_coeff <- as.formula(paste("mc", "0", sep=" ~ "))
clusn_u <- clusn_c <- NULL
clusn_mu <- clusn_mc <- NULL
if(!is.null(random_u) & length(random_u) > 1 | !is.null(random_c) & length(random_c) > 1) {
stop("random effects can be included in the formula only through a single expression within brackets")
}
if(all(names(model.frame(fixed_u, data = data)) %in% c("u", names(cov_matrix))) == FALSE |
all(names(model.frame(fixed_c, data = data)) %in% c("c", "u", names(cov_matrix))) == FALSE) {
stop("partially-observed covariates cannot be included in the fixed effects model")
}
if(all(names(model.frame(fixed_u, data = data)) %in% names(data)) == FALSE |
all(names(model.frame(fixed_c, data = data)) %in% names(data)) == FALSE) {
stop("you must provide names in the formula that correspond to those in the data")
}
if("u" %in% labels(terms(fixed_u)) | "c" %in% labels(terms(fixed_c))) {
stop("please remove 'u' from the right hand side of model.eff and/or 'c' from the right hand side of model.cost")
}
if(names(model.frame(fixed_u, data = data)[1]) != "u") {
stop("you must set 'u' as the response in the formula model.eff")
}
if("c" %in% names(model.frame(fixed_u, data = data))) {
stop("dependence allowed only through the cost model; please remove 'c' from model.eff")
}
if(names(model.frame(fixed_c, data = data)[1]) != "c") {
stop("you must set 'c' as the response in the formula model.cost")
}
if("u" %in% labels(terms(fixed_c))) {
if(length(grep(":u", labels(terms(fixed_c)))) != 0 | length(grep("u:", labels(terms(fixed_c)))) != 0) {
stop("no interaction effects for 'u' is allowed")
}
}
if("t" %in% names(model.frame(fixed_c, data = data)) | "t" %in% names(model.frame(fixed_u, data = data))) {
stop("treatment indicator must be provided only in the data. Please remove 't' from 'model.eff' and/or 'model.cost'")
}
index_mis_u <- which(is.na(data$u[data$time == 1]))
index_mis_c <- which(is.na(data$c[data$time == 1]))
data2 <- data
data2$mu <- data2$drop_u
data2$mc <- data2$drop_c
data.t1 <- data[data$time == 1, ]
data.t1$u[is.na(data.t1$u)] <- -999999
data.t1$c[is.na(data.t1$c)] <- -999999
mf_u_fixed <- model.frame(formula = fixed_u, data = data.t1)
mf_c_fixed <- model.frame(formula = fixed_c, data = data.t1)
terms <- NULL
x_u_fixed <- model.matrix(attr(mf_u_fixed, "terms"), data = mf_u_fixed)
x_c_fixed <- model.matrix(attr(mf_c_fixed, "terms"), data = mf_c_fixed)
if("u" %in% names(mf_c_fixed)){
mf_c_fixed$u[index_mis_u] <- NA
}
name_re_u_coeff <- NULL
name_re_c_coeff <- NULL
if(!is.null(random_u)){
name_re_u_coeff <- sub("\\|.*", "", random_u)
if(grepl("0 + 1", name_re_u_coeff, fixed = TRUE) == TRUE) { stop("Either remove or add the random intercept")}
name_clus_u <- sub('.*\\|', '', random_u)
if(lengths(strsplit(name_clus_u, " ")) > 2) { stop("a single clustering variable must selected for each formula") }
name_clus_u <- gsub(" ", "", name_clus_u, fixed = TRUE)
if(!name_clus_u %in% names(cov_matrix)) { stop("the clustering variable must be among the variables in the dataset") }
if(strsplit(name_re_u_coeff, "")[[1]][1] == 0) {
no_random_int_u <- TRUE} else {no_random_int_u <- FALSE }
if(no_random_int_u == TRUE) {
name_re_u_coeff <- sub("[0]", "", name_re_u_coeff)
name_re_u_coeff <- sub("[+]", "", name_re_u_coeff)
}
if(name_re_u_coeff == "" | name_re_u_coeff == " ") { stop("please state for which variables the random effects are assumed") }
fname_re_u_coeff <- as.formula(paste("u", name_re_u_coeff, sep = " ~ "))
if(all(names(model.frame(fname_re_u_coeff, data = data.t1)) %in% c("0", "1", names(model.frame(fixed_u, data = data.t1)))) == FALSE) {
stop("only covariates defined as fixed effects can be included in the random effects model")
}
if("u" %in% labels(terms(fname_re_u_coeff))) {
stop("please remove 'u' from the random effects expression of model.eff")
}
if("c" %in% labels(terms(fname_re_u_coeff))) {
stop("dependence allowed only through the cost model; please remove 'c' from model.eff")
}
mf_u_random <- model.frame(formula = fname_re_u_coeff, data = data.t1)
x_u_random <- model.matrix(attr(mf_u_random, "terms"), data = mf_u_random)
if(no_random_int_u == TRUE) {
x_u_random <- as.matrix(x_u_random[, !colnames(x_u_random) == "(Intercept)"])
if(is.null(colnames(x_u_random)) == TRUE & dim(x_u_random)[2] == 1) {
colnames(x_u_random) <- gsub(" ", "", name_re_u_coeff)
}
}
clus_u <- data[, name_clus_u]
if(!is.factor(clus_u)) { stop("clustering variables must be defined as factors") }
clusn_u <- as.numeric(clus_u)
if(!all(diff(sort(unique(clusn_u))) == 1) | !min(clusn_u) == 1) {
stop("ordered levels of clustering variables must not have gaps and must start from 1")
}
}
if(!is.null(random_c)){
name_re_c_coeff <- sub("\\|.*", "", random_c)
if(grepl("0 + 1", name_re_c_coeff, fixed = TRUE) == TRUE) { stop("Either remove or add the random intercept")}
name_clus_c <- sub('.*\\|', '', random_c)
if(lengths(strsplit(name_clus_c, " ")) > 2) { stop("a single clustering variable must selected for each formula") }
name_clus_c <- gsub(" ", "", name_clus_c, fixed = TRUE)
if(!name_clus_c %in% names(cov_matrix)) { stop("the clustering variable must be among the variables in the dataset") }
if(strsplit(name_re_c_coeff, "")[[1]][1] == 0) {
no_random_int_c <- TRUE} else {no_random_int_c <- FALSE }
if(no_random_int_c == TRUE) {
name_re_c_coeff <- sub("[0]", "", name_re_c_coeff)
name_re_c_coeff <- sub("[+]", "", name_re_c_coeff)
}
if(name_re_c_coeff == "" | name_re_c_coeff == " ") { stop("please state for which variables the random effects are assumed") }
if(gsub(" ", "", name_re_c_coeff) == "u" & no_random_int_c == FALSE) {name_re_c_coeff <- "1 + u" }
fname_re_c_coeff <- as.formula(paste("c", name_re_c_coeff, sep = " ~ "))
if(all(names(model.frame(fname_re_c_coeff, data = data.t1)) %in% c("0", "1", names(model.frame(fixed_c, data = data.t1)))) == FALSE) {
stop("only covariates defined as fixed effects can be included in the random effects model")
}
if("c" %in% labels(terms(fname_re_c_coeff))) {
stop("please remove 'c' from the random effects expression of model.cost")
}
if("u" %in% labels(terms(fname_re_c_coeff))) {
if(length(grep(":u", labels(terms(fname_re_c_coeff)))) != 0 | length(grep("u:", labels(terms(fname_re_c_coeff)))) != 0) {
stop("no interaction effects for 'u' is allowed")
}
}
mf_c_random <- model.frame(formula = fname_re_c_coeff, data = data.t1)
x_c_random <- model.matrix(attr(mf_c_random, "terms"), data = mf_c_random)
if("u" %in% labels(terms(fname_re_c_coeff)) & length(labels(terms(fname_re_c_coeff))) == 1) {
x_c_random <- subset(x_c_random, select = -c(u))
}
if(no_random_int_c == TRUE) {
x_c_random <- as.matrix(x_c_random[, !colnames(x_c_random) == "(Intercept)"])
if(is.null(colnames(x_c_random)) == TRUE & dim(x_c_random)[2] == 1) {
colnames(x_c_random) <- gsub(" ", "", name_re_c_coeff)
}
}
clus_c <- data[, name_clus_c]
if(!is.factor(clus_c)) { stop("clustering variables must be defined as factors") }
clusn_c <- as.numeric(clus_c)
if(!all(diff(sort(unique(clusn_c))) == 1) | !min(clusn_c) == 1) {
stop("ordered levels of clustering variables must not have gaps and must start from 1")
}
}
N1 <- N2 <- c()
N1 <- length(data.t1$u[data.t1$t == 1])
N2 <- length(data.t1$u[data.t1$t == 2])
N <- c(N1, N2)
n.time <- length(unique(data$time))
eff1.mat <- cost1.mat <- matrix(NA, nrow = N1, ncol = n.time)
eff2.mat <- cost2.mat <- matrix(NA, nrow = N2, ncol = n.time)
m_eff <- count_mis_u
m_cost <- count_mis_c
time1.mat <- matrix(NA, nrow = N1, ncol = n.time)
time2.mat <- matrix(NA, nrow = N2, ncol = n.time)
m_eff1.mat <- m_cost1.mat <- matrix(NA, nrow = N1, ncol = n.time)
m_eff2.mat <- m_cost2.mat <- matrix(NA, nrow = N2, ncol = n.time)
for(i in 1:n.time) {
eff1.mat[, i] <- data$u[data$time == i & data$t == 1]
eff2.mat[, i] <- data$u[data$time == i & data$t == 2]
cost1.mat[, i] <- data$c[data$time == i & data$t == 1]
cost2.mat[, i] <- data$c[data$time == i & data$t == 2]
time1.mat[, i] <- data$time[data$time == i & data$t == 1]
time2.mat[, i] <- data$time[data$time == i & data$t == 2]
m_eff1.mat[, i] <- count_mis_u[data$time == i & data$t == 1]
m_eff2.mat[, i] <- count_mis_u[data$time == i & data$t == 2]
m_cost1.mat[, i] <- count_mis_c[data$time == i & data$t == 1]
m_cost2.mat[, i] <- count_mis_c[data$time == i & data$t == 2]
}
effects <- list(eff1.mat, eff2.mat)
costs <- list(cost1.mat, cost2.mat)
time <- list(time1.mat, time2.mat)
m_eff <- list(m_eff1.mat, m_eff2.mat)
m_cost <- list(m_cost1.mat, m_cost2.mat)
N1_cc <- N2_cc <- N1_mis <- N2_mis <- matrix(NA, nrow = 2, ncol = n.time)
eff1_cc <- eff2_cc <- cost1_cc <- cost2_cc <- vector(mode = "list", length = n.time)
for(i in 1:n.time){
N1_cc[1, i] <- length(as.vector(apply(eff1.mat, 2, na.omit)[[i]]))
N1_cc[2, i] <- length(as.vector(apply(cost1.mat, 2, na.omit)[[i]]))
N2_cc[1, i] <- length(as.vector(apply(eff2.mat, 2, na.omit)[[i]]))
N2_cc[2, i] <- length(as.vector(apply(cost2.mat, 2, na.omit)[[i]]))
eff1_cc[[i]] <- as.vector(apply(eff1.mat, 2, na.omit)[[i]])
cost1_cc[[i]] <- as.vector(apply(cost1.mat, 2, na.omit)[[i]])
eff2_cc[[i]] <- as.vector(apply(eff2.mat, 2, na.omit)[[i]])
cost2_cc[[i]] <- as.vector(apply(cost2.mat, 2, na.omit)[[i]])
}
N_cc <- list(N1_cc, N2_cc)
N1_mis <- N1 - N1_cc
N2_mis <- N2 - N2_cc
N_mis <- list(N1_mis, N2_mis)
eff_cc <- list(eff1_cc, eff2_cc)
cost_cc <- list(cost1_cc, cost2_cc)
cov1_u_fixed <- as.data.frame(x_u_fixed[data.t1$t == 1, ])
names(cov1_u_fixed) <- colnames(x_u_fixed)
cov2_u_fixed <- as.data.frame(x_u_fixed[data.t1$t == 2, ])
names(cov2_u_fixed) <- colnames(x_u_fixed)
cov_u_fixed <- list(cov1_u_fixed, cov2_u_fixed)
x_c_hold_fixed <- x_c_fixed
if("u" %in% colnames(x_c_hold_fixed)) {
x_c_fixed <- subset(x_c_hold_fixed, select = -c(u))
}
cov1_c_fixed <- as.data.frame(x_c_fixed[data.t1$t == 1, ])
names(cov1_c_fixed) <- colnames(x_c_fixed)
cov2_c_fixed <- as.data.frame(x_c_fixed[data.t1$t == 2, ])
names(cov2_c_fixed) <- colnames(x_c_fixed)
cov_c_fixed <- list(cov1_c_fixed, cov2_c_fixed)
covu_fixed <- list(cov1_u_fixed, cov2_u_fixed)
mean_cov_u_fixed <- list(apply(as.matrix(cov1_u_fixed), 2, mean), apply(as.matrix(cov2_u_fixed), 2, mean))
covc_fixed <- list(cov1_c_fixed, cov2_c_fixed)
mean_cov_c_fixed <- list(apply(as.matrix(cov1_c_fixed), 2, mean), apply(as.matrix(cov2_c_fixed), 2, mean))
cov1_u_center_fixed <- as.data.frame(scale(cov1_u_fixed, scale = FALSE))
cov2_u_center_fixed <- as.data.frame(scale(cov2_u_fixed, scale = FALSE))
cov1_u_center_fixed[, 1] <- rep(1, nrow(cov1_u_fixed))
cov2_u_center_fixed[, 1] <- rep(1, nrow(cov2_u_fixed))
cov_u_center_fixed <- list(cov1_u_center_fixed, cov2_u_center_fixed)
mean_cov_u_center_fixed <- list(apply(as.matrix(cov1_u_center_fixed), 2, mean), apply(as.matrix(cov2_u_center_fixed), 2, mean))
cov1_c_center_fixed <- as.data.frame(scale(cov1_c_fixed, scale = FALSE))
cov2_c_center_fixed <- as.data.frame(scale(cov2_c_fixed, scale = FALSE))
cov1_c_center_fixed[, 1] <- rep(1, nrow(cov1_c_fixed))
cov2_c_center_fixed[, 1] <- rep(1, nrow(cov2_c_fixed))
cov_c_center_fixed <- list(cov1_c_center_fixed, cov2_c_center_fixed)
mean_cov_c_center_fixed <- list(apply(as.matrix(cov1_c_center_fixed), 2, mean), apply(as.matrix(cov2_c_center_fixed), 2, mean))
if(center == TRUE) {
cov_u_fixed <- cov_u_center_fixed
cov_c_fixed <- cov_c_center_fixed
mean_cov_u_fixed <- mean_cov_u_center_fixed
mean_cov_c_fixed <- mean_cov_c_center_fixed
}
if(!is.null(random_u)){
cov1_u_random <- as.data.frame(x_u_random[data.t1$t == 1, ])
names(cov1_u_random) <- colnames(x_u_random)
cov2_u_random <- as.data.frame(x_u_random[data.t1$t == 2, ])
names(cov2_u_random) <- colnames(x_u_random)
cov_u_random <- list(cov1_u_random, cov2_u_random)
covu_random <- list(cov1_u_random, cov2_u_random)
mean_cov_u_random <- list(apply(as.matrix(cov1_u_random), 2, mean), apply(as.matrix(cov2_u_random), 2, mean))
cov1_u_center_random <- as.data.frame(scale(cov1_u_random, scale = FALSE))
cov2_u_center_random <- as.data.frame(scale(cov2_u_random, scale = FALSE))
if(no_random_int_u == FALSE) {
cov1_u_center_random[, 1] <- rep(1, nrow(cov1_u_random))
cov2_u_center_random[, 1] <- rep(1, nrow(cov2_u_random))
}
cov_u_center_random <- list(cov1_u_center_random, cov2_u_center_random)
mean_cov_u_center_random <- list(apply(as.matrix(cov1_u_center_random), 2, mean), apply(as.matrix(cov2_u_center_random), 2, mean))
if(center == TRUE) {
cov_u_random <- cov_u_center_random
mean_cov_u_random <- mean_cov_u_center_random
}
clusn_u1 <- clusn_u[data.t1$t == 1]
clusn_u1 <- factor(clusn_u1, levels = unique(clusn_u1))
clusn_u2 <- clusn_u[data.t1$t == 2]
clusn_u2 <- factor(clusn_u2, levels = unique(clusn_u2))
}
if(!is.null(random_c)){
x_c_hold_random <- x_c_random
if("u" %in% colnames(x_c_hold_random)) {
x_c_random <- subset(x_c_hold_random, select = -c(u))
}
cov1_c_random <- as.data.frame(x_c_random[data.t1$t == 1, ])
names(cov1_c_random) <- colnames(x_c_random)
cov2_c_random <- as.data.frame(x_c_random[data.t1$t == 2, ])
names(cov2_c_random) <- colnames(x_c_random)
cov_c_random <- list(cov1_c_random, cov2_c_random)
covc_random <- list(cov1_c_random, cov2_c_random)
mean_cov_c_random <- list(apply(as.matrix(cov1_c_random), 2, mean), apply(as.matrix(cov2_c_random), 2, mean))
cov1_c_center_random <- as.data.frame(scale(cov1_c_random, scale = FALSE))
cov2_c_center_random <- as.data.frame(scale(cov2_c_random, scale = FALSE))
if(no_random_int_c == FALSE) {
cov1_c_center_random[, 1] <- rep(1, nrow(cov1_c_random))
cov2_c_center_random[, 1] <- rep(1, nrow(cov2_c_random))
}
cov_c_center_random <- list(cov1_c_center_random, cov2_c_center_random)
mean_cov_c_center_random <- list(apply(as.matrix(cov1_c_center_random), 2, mean), apply(as.matrix(cov2_c_center_random), 2, mean))
if(center == TRUE) {
cov_c_random <- cov_c_center_random
mean_cov_c_random <- mean_cov_c_center_random
}
clusn_c1 <- clusn_c[data.t1$t == 1]
clusn_c1 <- factor(clusn_c1, levels = unique(clusn_c1))
clusn_c2 <- clusn_c[data.t1$t == 2]
clusn_c2 <- factor(clusn_c2, levels = unique(clusn_c2))
}
data2.t1 <- data[data$time == 1, ]
data2.t1$u[is.na(data2.t1$u) == TRUE] <- -999999
data2.t1$c[is.na(data2.t1$c) == TRUE] <- -999999
data2.t1$mu <- c(m_eff1.mat[,1], m_eff2.mat[,1])
data2.t1$mc <- c(m_cost1.mat[,1], m_cost2.mat[,1])
if(!is.formula(model.mu) | !is.formula(model.mc)) {
stop("model.mu and/or model.mc must be formula objects")
}
fixed_mu <- nobars_(model.mu)
fixed_mc <- nobars_(model.mc)
random_mu <- fb(model.mu)
random_mc <- fb(model.mc)
if(!is.null(random_mu) & length(random_mu) > 1 | !is.null(random_mc) & length(random_mc) > 1) {
stop("random effects can be included in the formula only through a single expression within brackets")
}
if(all(names(model.frame(fixed_mu, data = data2.t1)) %in% c("mu", "u", names(cov_matrix))) == FALSE |
all(names(model.frame(fixed_mc, data = data2.t1)) %in% c("mc", "c", names(cov_matrix))) == FALSE) {
stop("partially-observed covariates cannot be included in the model")
}
if(all(names(model.frame(fixed_mu, data = data2.t1)) %in% names(data2.t1)) == FALSE |
all(names(model.frame(fixed_mc, data = data2.t1)) %in% names(data2.t1)) == FALSE) {
stop("you must provide names in the formula that correspond to those in the data")
}
if(names(model.frame(fixed_mu, data = data2.t1)[1]) != "mu") {
stop("you must set 'mu' as the response in the formula model.mu")
}
if(names(model.frame(fixed_mc, data = data2.t1)[1]) != "mc") {
stop("you must set 'mc' as the response in the formula model.mc")
}
if("t" %in% names(model.frame(fixed_mc, data = data2.t1)) | "t" %in% names(model.frame(fixed_mu, data = data2.t1))) {
stop("treatment indicator must be provided only in the data. Please remove 't' from 'model.mu' and/or 'model.mc'")
}
if("c" %in% names(model.frame(fixed_mu, data = data2.t1)) | "u" %in% names(model.frame(fixed_mc, data = data2.t1))) {
stop("please remove 'u' from model.mc and/or remove 'c' from model.mu")
}
if("u" %in% labels(terms(fixed_mu))) {
if(length(grep(":u", labels(terms(fixed_mu)))) != 0 | length(grep("u:", labels(terms(fixed_mu)))) != 0) {
stop("no interaction effects for 'u' is allowed")
}
}
if("c" %in% labels(terms(fixed_mc))) {
if(length(grep(":c", labels(terms(fixed_mc)))) != 0 | length(grep("c:", labels(terms(fixed_mc)))) != 0) {
stop("no interaction effects for 'c' is allowed")
}
}
name_re_mu_coeff <- NULL
name_re_mc_coeff <- NULL
if(!is.null(random_mu)){
name_re_mu_coeff <- sub("\\|.*", "", random_mu)
if(grepl("0 + 1", name_re_mu_coeff, fixed = TRUE) == TRUE) { stop("Either remove or add the random intercept")}
name_clus_mu <- sub('.*\\|', '', random_mu)
if(lengths(strsplit(name_clus_mu, " ")) > 2) {stop("a single clustering variable must be selected for each formula") }
name_clus_mu <- gsub(" ", "", name_clus_mu, fixed = TRUE)
if(!name_clus_mu %in% names(cov_matrix)) { stop("the clustering variable must be among the variables in the dataset") }
if(strsplit(name_re_mu_coeff, "")[[1]][1] == 0) {
no_random_int_mu <- TRUE} else {no_random_int_mu <- FALSE }
if(no_random_int_mu == TRUE) {
name_re_mu_coeff <- sub("[0]", "", name_re_mu_coeff)
name_re_mu_coeff <- sub("[+]", "", name_re_mu_coeff)
}
if(name_re_mu_coeff == "" | name_re_mu_coeff == " ") { stop("please state for which variables the random effects are assumed") }
if(gsub(" ", "", name_re_mu_coeff) == "u" & no_random_int_mu == FALSE) {name_re_mu_coeff <- "1 + u" }
fname_re_mu_coeff <- as.formula(paste("mu", name_re_mu_coeff, sep = " ~ "))
if(all(names(model.frame(fname_re_mu_coeff, data = data2)) %in% c("0","1", names(model.frame(fixed_mu, data = data2)))) == FALSE) {
stop("only covariates inlcued as fixed effects can be included in the random effects model")
}
if("mu" %in% labels(terms(fname_re_mu_coeff))) {
stop("please remove 'mu' from the right hand side of model.mu")
}
if("u" %in% labels(terms(fname_re_mu_coeff))) {
if(length(grep(":u", labels(terms(fname_re_mu_coeff)))) != 0 | length(grep("u:", labels(terms(fname_re_mu_coeff)))) != 0) {
stop("no interaction effects for 'u' are allowed")
}
}
clus_mu <- data[, name_clus_mu]
if(!is.factor(clus_mu)) { stop("clustering variables must be defined as factors") }
clusn_mu <- as.numeric(clus_mu)
if(!all(diff(sort(unique(clusn_mu))) == 1) | !min(clusn_mu) == 1) {
stop("ordered levels of clustering variables must not have gaps and must start from 1")
}
}
if(!is.null(random_mc)){
name_re_mc_coeff <- sub("\\|.*", "", random_mc)
if(grepl("0 + 1", name_re_mc_coeff, fixed = TRUE) == TRUE) { stop("Either remove or add the random intercept")}
name_clus_mc <- sub('.*\\|', '', random_mc)
if(lengths(strsplit(name_clus_mc, " ")) > 2) {stop("a single clustering variable must be selected for each formula") }
name_clus_mc <- gsub(" ", "", name_clus_mc, fixed = TRUE)
if(!name_clus_mc %in% names(cov_matrix)) { stop("the clustering variable must be among the variables in the dataset") }
if(strsplit(name_re_mc_coeff, "")[[1]][1] == 0) {
no_random_int_mc <- TRUE} else {no_random_int_mc <- FALSE }
if(no_random_int_mc == TRUE) {
name_re_mc_coeff <- sub("[0]", "", name_re_mc_coeff)
name_re_mc_coeff <- sub("[+]", "", name_re_mc_coeff)
}
if(name_re_mc_coeff == "" | name_re_mc_coeff == " ") { stop("please state for which variables the random effects are assumed") }
if(gsub(" ", "", name_re_mc_coeff) == "c" & no_random_int_mc == FALSE) {name_re_mc_coeff <- "1 + c" }
fname_re_mc_coeff <- as.formula(paste("mc", name_re_mc_coeff, sep=" ~ "))
if(all(names(model.frame(fname_re_mc_coeff, data = data2)) %in% c("0","1", names(model.frame(fixed_mc, data = data2)))) == FALSE) {
stop("only covariates inlcued as fixed effects can be included in the random effects model")
}
if("mc" %in% labels(terms(fname_re_mc_coeff))) {
stop("please remove 'mc' from the right hand side of model.mc")
}
if("c" %in% labels(terms(fname_re_mc_coeff))) {
if(length(grep(":c", labels(terms(fname_re_mc_coeff)))) != 0 | length(grep("c:", labels(terms(fname_re_mc_coeff)))) != 0) {
stop("no interaction effects for 'c' is allowed")
}
}
clus_mc <- data[, name_clus_mc]
if(!is.factor(clus_mc)) { stop("clustering variables must be defined as factors") }
clusn_mc <- as.numeric(clus_mc)
if(!all(diff(sort(unique(clusn_mc))) == 1) | !min(clusn_mc) == 1) {
stop("ordered levels of clustering variables must not have gaps and must start from 1")
}
}
mf_mu_fixed <- model.frame(formula = fixed_mu, data = data2.t1)
mf_mc_fixed <- model.frame(formula = fixed_mc, data = data2.t1)
z_u_fixed <- model.matrix(attr(mf_mu_fixed, "terms"), data = mf_mu_fixed)
z_c_fixed <- model.matrix(attr(mf_mc_fixed, "terms"), data = mf_mc_fixed)
z_u_hold_fixed <- z_u_fixed
if("u" %in% colnames(z_u_hold_fixed)) {
z_u_fixed <- subset(z_u_hold_fixed, select = -c(u))
}
z_c_hold_fixed <- z_c_fixed
if("c" %in% colnames(z_c_hold_fixed)) {
z_c_fixed <- subset(z_c_hold_fixed, select = -c(c))
}
covz1_u_fixed <- as.data.frame(z_u_fixed[data2.t1$t == 1, ])
names(covz1_u_fixed) <- colnames(z_u_fixed)
covz2_u_fixed <- as.data.frame(z_u_fixed[data2.t1$t == 2, ])
names(covz2_u_fixed) <- colnames(z_u_fixed)
covz_u_fixed <- list(covz1_u_fixed, covz2_u_fixed)
covzu_fixed <- list(covz1_u_fixed, covz2_u_fixed)
mean_covz_u_fixed <- list(apply(as.matrix(covz1_u_fixed), 2, mean), apply(as.matrix(covz2_u_fixed), 2, mean))
names(covz_u_fixed) <- names(mean_covz_u_fixed) <- c("Control", "Intervention")
covz1_c_fixed <- as.data.frame(z_c_fixed[data2.t1$t == 1, ])
names(covz1_c_fixed) <- colnames(z_c_fixed)
covz2_c_fixed <- as.data.frame(z_c_fixed[data2.t1$t == 2, ])
names(covz2_c_fixed) <- colnames(z_c_fixed)
covz_c_fixed <- list(covz1_c_fixed, covz2_c_fixed)
covzc_fixed <- list(covz1_c_fixed, covz2_c_fixed)
mean_covz_c_fixed <- list(apply(as.matrix(covz1_c_fixed), 2, mean), apply(as.matrix(covz2_c_fixed), 2, mean))
covz1_u_center_fixed <- as.data.frame(scale(covz1_u_fixed, scale = FALSE))
covz2_u_center_fixed <- as.data.frame(scale(covz2_u_fixed, scale = FALSE))
covz1_u_center_fixed[, 1] <- rep(1, nrow(covz1_u_fixed))
covz2_u_center_fixed[, 1] <- rep(1, nrow(covz2_u_fixed))
covz_u_center_fixed <- list(covz1_u_center_fixed, covz2_u_center_fixed)
mean_covz_u_center_fixed <- list(apply(as.matrix(covz1_u_center_fixed), 2, mean), apply(as.matrix(covz2_u_center_fixed), 2, mean))
covz1_c_center_fixed <- as.data.frame(scale(covz1_c_fixed, scale = FALSE))
covz2_c_center_fixed <- as.data.frame(scale(covz2_c_fixed, scale = FALSE))
covz1_c_center_fixed[, 1] <- rep(1, nrow(covz1_c_fixed))
covz2_c_center_fixed[, 1] <- rep(1, nrow(covz2_c_fixed))
covz_c_center_fixed <- list(covz1_c_center_fixed, covz2_c_center_fixed)
mean_covz_c_center_fixed <- list(apply(as.matrix(covz1_c_center_fixed), 2, mean), apply(as.matrix(covz2_c_center_fixed), 2, mean))
if(center == TRUE) {
covz_u_fixed <- covz_u_center_fixed
covz_c_fixed <- covz_c_center_fixed
mean_covz_u_fixed <- mean_covz_u_center_fixed
mean_covz_c_fixed <- mean_covz_c_center_fixed
}
if(!is.null(random_mu)){
mf_mu_random <- model.frame(formula = fname_re_mu_coeff, data = data2.t1)
z_u_random <- model.matrix(attr(mf_mu_random, "terms"), data = mf_mu_random)
if(no_random_int_mu == TRUE) {
z_u_random <- as.matrix(z_u_random[, !colnames(z_u_random) == "(Intercept)"])
if(dim(z_u_random)[2] == 1) { colnames(z_u_random) <- colnames(model.matrix(attr(mf_mu_random, "terms"), data = mf_mu_random))[2] }
}
z_u_hold_random <- z_u_random
if("u" %in% colnames(z_u_hold_random)) {
z_u_random <- subset(z_u_hold_random, select = -c(u))
}
covz1_u_random <- as.data.frame(z_u_random[data2.t1$t == 1, ])
names(covz1_u_random) <- colnames(z_u_random)
covz2_u_random <- as.data.frame(z_u_random[data2.t1$t == 2, ])
names(covz2_u_random) <- colnames(z_u_random)
covz_u_random <- list(covz1_u_random, covz2_u_random)
covzu_random <- list(covz1_u_random, covz2_u_random)
mean_covz_u_random <- list(apply(as.matrix(covz1_u_random), 2, mean), apply(as.matrix(covz2_u_random), 2, mean))
names(covz_u_random) <- names(mean_covz_u_random) <- c("Control", "Intervention")
covz1_u_center_random <- as.data.frame(scale(covz1_u_random, scale = FALSE))
covz2_u_center_random <- as.data.frame(scale(covz2_u_random, scale = FALSE))
if(no_random_int_mu == FALSE) {
covz1_u_center_random[, 1] <- rep(1, nrow(covz1_u_random))
covz2_u_center_random[, 1] <- rep(1, nrow(covz2_u_random))
}
covz_u_center_random <- list(covz1_u_center_random, covz2_u_center_random)
mean_covz_u_center_random <- list(apply(as.matrix(covz1_u_center_random), 2, mean), apply(as.matrix(covz2_u_center_random), 2, mean))
if(center == TRUE) {
covz_u_random <- covz_u_center_random
mean_covz_u_random <- mean_covz_u_center_random
}
clusn_mu1 <- clusn_mu[data2.t1$t == 1]
clusn_mu1 <- factor(clusn_mu1, levels = unique(clusn_mu1))
clusn_mu2 <- clusn_mu[data2.t1$t == 2]
clusn_mu2 <- factor(clusn_mu2, levels = unique(clusn_mu2))
}
if(!is.null(random_mc)){
mf_mc_random <- model.frame(formula = fname_re_mc_coeff, data = data2.t1)
z_c_random <- model.matrix(attr(mf_mc_random, "terms"), data = mf_mc_random)
if(no_random_int_mc == TRUE) {
z_c_random <- as.matrix(z_c_random[, !colnames(z_c_random) == "(Intercept)"])
if(dim(z_c_random)[2] == 1) { colnames(z_c_random) <- colnames(model.matrix(attr(mf_mc_random, "terms"), data = mf_mc_random))[2] }
}
z_c_hold_random <- z_c_random
if("c" %in% colnames(z_c_hold_random)) {
z_c_random <- subset(z_c_hold_random, select = -c(c))
}
covz1_c_random <- as.data.frame(z_c_random[data2.t1$t == 1, ])
names(covz1_c_random) <- colnames(z_c_random)
covz2_c_random <- as.data.frame(z_c_random[data2.t1$t == 2, ])
names(covz2_c_random) <- colnames(z_c_random)
covz_c_random <- list(covz1_c_random, covz2_c_random)
covzc_random <- list(covz1_c_random, covz2_c_random)
mean_covz_c_random <- list(apply(as.matrix(covz1_c_random), 2, mean), apply(as.matrix(covz2_c_random), 2, mean))
names(covz_c_random) <- names(mean_covz_c_random) <- c("Control", "Intervention")
covz1_c_center_random <- as.data.frame(scale(covz1_c_random, scale = FALSE))
covz2_c_center_random <- as.data.frame(scale(covz2_c_random, scale = FALSE))
if(no_random_int_mc == FALSE) {
covz1_c_center_random[, 1] <- rep(1, nrow(covz1_c_random))
covz2_c_center_random[, 1] <- rep(1, nrow(covz2_c_random))
}
covz_c_center_random <- list(covz1_c_center_random, covz2_c_center_random)
mean_covz_c_center_random <- list(apply(as.matrix(covz1_c_center_random), 2, mean), apply(as.matrix(covz2_c_center_random), 2, mean))
if(center == TRUE) {
covz_c_random <- covz_c_center_random
mean_covz_c_random <- mean_covz_c_center_random
}
clusn_mc1 <- clusn_mc[data2.t1$t == 1]
clusn_mc1 <- factor(clusn_mc1, levels = unique(clusn_mc1))
clusn_mc2 <- clusn_mc[data2.t1$t == 2]
clusn_mc2 <- factor(clusn_mc2, levels = unique(clusn_mc2))
}
names(covz_u_fixed) <- names(mean_covz_u_fixed) <- names(covz_c_fixed) <- names(mean_covz_c_fixed) <- c("Control", "Intervention")
names(cov_u_fixed) <- names(cov_c_fixed) <- names(mean_cov_u_fixed) <- names(mean_cov_c_fixed) <- c("Control", "Intervention")
if(!is.null(random_c)) {
names(cov_c_random) <- names(mean_cov_c_random) <- c("Control", "Intervention")
clusn_c <- list("Control" = clusn_c1, "Intervention" = clusn_c2)
} else {cov_c_random <- mean_cov_c_random <- NULL}
if(!is.null(random_mc)) {
names(covz_c_random) <- names(mean_covz_c_random) <- c("Control", "Intervention")
clusn_mc <- list("Control" = clusn_mc1, "Intervention" = clusn_mc2)
} else {covz_c_random <- mean_covz_c_random <- NULL}
if(!is.null(random_u)) {
names(cov_u_random) <- names(mean_cov_u_random) <- c("Control", "Intervention")
clusn_u <- list("Control" = clusn_u1, "Intervention" = clusn_u2)
} else {cov_u_random <- mean_cov_u_random <- NULL}
if(!is.null(random_mu)) {
names(covz_u_random) <- names(mean_covz_u_random) <- c("Control", "Intervention")
clusn_mu <- list("Control" = clusn_mu1, "Intervention" = clusn_mu2)
} else {covz_u_random <- mean_covz_u_random <- NULL}
names(m_eff) <- names(m_cost) <- names(time) <- c("Control", "Intervention")
names(effects) <- names(costs) <- names(eff_cc) <- names(cost_cc) <- c("Control", "Intervention")
data_raw <- list("raw_effects" = effects, "raw_costs" = costs, "raw_effects_cc" = eff_cc, "raw_costs_cc" = cost_cc, "arm_lengths" = N,
"arm_lengths_cc" = N_cc, "arm_missing_data" = N_mis, "missing_effects" = m_eff, "missing_costs" = m_cost, "time" = time,
"covariates_effects_fixed" = cov_u_fixed, "covariates_costs_fixed" = cov_c_fixed, "mean_cov_effects_fixed" = mean_cov_u_fixed, "mean_cov_costs_fixed" = mean_cov_c_fixed,
"covariates_missing_effects_fixed" = covz_u_fixed, "mean_cov_missing_effects_fixed" = mean_covz_u_fixed, "covariates_missing_costs_fixed" = covz_c_fixed,
"mean_cov_missing_costs_fixed" = mean_covz_c_fixed, "covariates_effects_random" = cov_u_random, "covariates_costs_random" = cov_c_random, "mean_cov_effects_random" = mean_cov_u_random, "mean_cov_costs_random" = mean_cov_c_random,
"covariates_missing_effects_random" = covz_u_random, "mean_cov_missing_effects_random" = mean_covz_u_random, "covariates_missing_costs_random" = covz_c_random,
"mean_cov_missing_costs_random" = mean_covz_c_random, "clus_u" = clusn_u, "clus_c" = clusn_c, "clus_mu" = clusn_mu, "clus_mc" = clusn_mc, "data_ind" = data2.t1, "data_long" = data)
model_formula <- list("mf_model.u_fixed" = fixed_u, "mf_model.c_fixed" = fixed_c, "mf_model.mu_fixed" = fixed_mu, "mf_model.mc_fixed" = fixed_mc,
"mf_model.u_random" = fname_re_u_coeff, "mf_model.c_random" = fname_re_c_coeff, "mf_model.mu_random" = fname_re_mu_coeff, "mf_model.mc_random" = fname_re_mc_coeff)
data_list <- list("data_raw" = data_raw, "model_formula" = model_formula)
return(data_list)
}
AnGabrio/missingHE documentation built on March 22, 2023, 12:55 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions data_read_selection_long
Documented in data_read_selection_long
#' A function to read and re-arrange the data in different ways
#'
#' This internal function imports the data and outputs only those variables that are needed to run the model
#' according to the information provided by the user.
#' @param data A data frame in which to find variables supplied in \code{model.eff} and \code{model.cost}. Among these,
#' effectiveness, cost and treatment indicator (only two arms) variables must always be provided and named 'u', 'c' and 't' respectively.
#' @param model.eff A formula expression in conventional \code{R} linear modelling syntax. The response must be a health economics
#' effectiveness outcome ('u') whose name must correspond to that used in \code{data}, and
#' any covariates are given on the right-hand side. If there are no covariates, specify \code{1} on the right hand side.
#' By default, covariates are placed on the "location" parameter of the distribution through a linear model.
#' Random effects can also be specified for each model parameter.
#' @param model.cost A formula expression in conventional \code{R} linear modelling syntax. The response must be a health economics
#' cost outcome ('c') whose name must correspond to that used in \code{data}, and any covariates are given on the right-hand side.
#' If there are no covariates, specify \code{1} on the right hand side. By default, covariates are placed on the "location"
#' parameter of the distribution through a linear model. Random effects can also be specified for each model parameter.
#' @param model.mu A formula expression in conventional \code{R} linear modelling syntax. The response must be indicated with the
#' term 'mu'(missing effects) and any covariates used to estimate the probability of missing effects are given on the right-hand side.
#' If there are no covariates, specify \code{1} on the right hand side. By default, covariates are placed on the "probability" parameter for the missing effects through a logistic-linear model.
#' Random effects can also be specified for each model parameter.
#' @param model.mc A formula expression in conventional \code{R} linear modelling syntax. The response must be indicated with the
#' term 'mc'(missing costs) and any covariates used to estimate the probability of missing costs should be given on the right-hand side.
#' If there are no covariates, specify \code{1} on the right hand side. By default, covariates are placed on the "probability" parameter for the missing costs through a logistic-linear model.
#' Random effects can also be specified for each model parameter.
#' @param type Type of missingness mechanism assumed. Choices are Missing At Random (MAR) and Missing Not At Random (MNAR).
#' @param center Logical. If \code{center} is \code{TRUE} all the covariates in the model are centered.
#' @keywords read data
#' @importFrom stats na.omit sd as.formula model.matrix model.frame model.response terms aggregate reshape
#' @export
#' @examples
#' #Internal function only
#' #no examples
#' #
#' #
data_read_selection_long <- function(data, model.eff, model.cost, model.mu, model.mc, type, center) {
if(is.data.frame(data) == FALSE) {
stop("object data must be provided as data frame")
}
if(any(names(data) == "u") == TRUE & any(names(data) == "c") == TRUE) {
u <- as.name("u")
c <- as.name("c")
}
data.wide <- reshape(data, v.names = c("u", "c"), timevar = c("time"),
direction = "wide", idvar = "id", sep = ".")
data.wide_u <- data.wide[, grepl('u\\.|id', colnames(data.wide))]
wide_u_index <- grepl('u\\.|id', colnames(data.wide_u))
data.wide_u <- data.wide_u[, wide_u_index]
data.wide_u_obs_index <- as.data.frame(which(is.na(data.wide_u) == FALSE, arr.ind = TRUE))
data.wide_u_obs_index <- merge(data.wide_u_obs_index, aggregate(col ~ row, data.wide_u_obs_index, max))
data.wide_u_obs_index <- unique(data.wide_u_obs_index)
names(data.wide_u_obs_index) <- c("id", "drop_u")
data.wide_u_drop <- merge(data.wide, data.wide_u_obs_index, by.x = "id", all.x = TRUE)
data.wide_u_drop$drop_u <- data.wide_u_drop$drop_u - 1
data.wide_c <- data.wide[, grepl('c\\.|id', colnames(data.wide))]
wide_c_index <- grepl('c\\.|id', colnames(data.wide_c))
data.wide_c <- data.wide_c[, wide_c_index]
data.wide_c_obs_index <- as.data.frame(which(is.na(data.wide_c) == FALSE, arr.ind = TRUE))
data.wide_c_obs_index <- merge(data.wide_c_obs_index, aggregate(col ~ row, data.wide_c_obs_index, max))
data.wide_c_obs_index <- unique(data.wide_c_obs_index)
names(data.wide_c_obs_index) <- c("id", "drop_c")
data.wide_uc_drop <- merge(data.wide_u_drop, data.wide_c_obs_index, by.x = "id", all.x = TRUE)
data.wide_uc_drop$drop_c <- data.wide_uc_drop$drop_c - 1
data.wide_uc_drop <- data.wide_uc_drop[, which(names(data.wide_uc_drop) %in% c("id","drop_u","drop_c"))]
data.wide_drop <- merge(data.wide, data.wide_uc_drop, by = "id", all = TRUE)
u_index_names <- colnames(data.wide_drop)[grepl('u\\.', colnames(data.wide_drop))]
c_index_names <- colnames(data.wide_drop)[grepl('c\\.', colnames(data.wide_drop))]
data.long <- reshape(data.wide_drop, varying = c(u_index_names,c_index_names),
direction = "long", idvar = "id", sep = ".")
time_max <- max(data.long$time)
count_mis_u <- rep(NA, length(data.long$drop_u))
count_mis_u <- ifelse(is.na(data.long$u) == FALSE, 1, count_mis_u)
count_mis_u <- ifelse(data.long$drop_u == 1 & is.na(data.long$u) == TRUE, 3, count_mis_u)
count_mis_u <- ifelse(data.long$drop_u <= data.long$time & is.na(data.long$u) == TRUE, 3, count_mis_u)
count_mis_u <- ifelse(data.long$drop_u > data.long$time & is.na(data.long$u) == TRUE, 2, count_mis_u)
count_mis_c <- rep(NA, length(data.long$drop_c))
count_mis_c <- ifelse(is.na(data.long$c) == FALSE, 1, count_mis_c)
count_mis_c <- ifelse(data.long$drop_c == 1 & is.na(data.long$c) == TRUE, 3, count_mis_c)
count_mis_c <- ifelse(data.long$drop_c <= data.long$time & is.na(data.long$c) == TRUE, 3, count_mis_c)
count_mis_c <- ifelse(data.long$drop_c > data.long$time & is.na(data.long$c) == TRUE, 2, count_mis_c)
data <- data.long
my_subset <- function() {
drop_u <- drop_c <- NULL
cov_matrix <- subset(data, select = -c(u, c, drop_u, drop_c, time, t))[data$time == 1,]
cov_matrix
}
cov_matrix <- my_subset()
cov_matrix <- cov_matrix[!unlist(vapply(cov_matrix, anyNA, logical(1)))]
is.formula <- function (x) { inherits(x, "formula") }
if(is.formula(model.eff) == FALSE | is.formula(model.cost) == FALSE) {
stop("model.eff and/or model.cost must be formula objects")
}
if(is.logical(center) == FALSE) { stop("center must be either TRUE or FALSE") }
fixed_u <- nobars_(model.eff)
fixed_c <- nobars_(model.cost)
random_u <- fb(model.eff)
random_c <- fb(model.cost)
fname_re_u_coeff <- as.formula(paste("u", "0", sep=" ~ "))
fname_re_c_coeff <- as.formula(paste("c", "0", sep=" ~ "))
fname_re_mu_coeff <- as.formula(paste("mu", "0", sep=" ~ "))
fname_re_mc_coeff <- as.formula(paste("mc", "0", sep=" ~ "))
clusn_u <- clusn_c <- NULL
clusn_mu <- clusn_mc <- NULL
if(!is.null(random_u) & length(random_u) > 1 | !is.null(random_c) & length(random_c) > 1) {
stop("random effects can be included in the formula only through a single expression within brackets")
}
if(all(names(model.frame(fixed_u, data = data)) %in% c("u", names(cov_matrix))) == FALSE |
all(names(model.frame(fixed_c, data = data)) %in% c("c", "u", names(cov_matrix))) == FALSE) {
stop("partially-observed covariates cannot be included in the fixed effects model")
}
if(all(names(model.frame(fixed_u, data = data)) %in% names(data)) == FALSE |
all(names(model.frame(fixed_c, data = data)) %in% names(data)) == FALSE) {
stop("you must provide names in the formula that correspond to those in the data")
}
if("u" %in% labels(terms(fixed_u)) | "c" %in% labels(terms(fixed_c))) {
stop("please remove 'u' from the right hand side of model.eff and/or 'c' from the right hand side of model.cost")
}
if(names(model.frame(fixed_u, data = data)[1]) != "u") {
stop("you must set 'u' as the response in the formula model.eff")
}
if("c" %in% names(model.frame(fixed_u, data = data))) {
stop("dependence allowed only through the cost model; please remove 'c' from model.eff")
}
if(names(model.frame(fixed_c, data = data)[1]) != "c") {
stop("you must set 'c' as the response in the formula model.cost")
}
if("u" %in% labels(terms(fixed_c))) {
if(length(grep(":u", labels(terms(fixed_c)))) != 0 | length(grep("u:", labels(terms(fixed_c)))) != 0) {
stop("no interaction effects for 'u' is allowed")
}
}
if("t" %in% names(model.frame(fixed_c, data = data)) | "t" %in% names(model.frame(fixed_u, data = data))) {
stop("treatment indicator must be provided only in the data. Please remove 't' from 'model.eff' and/or 'model.cost'")
}
index_mis_u <- which(is.na(data$u[data$time == 1]))
index_mis_c <- which(is.na(data$c[data$time == 1]))
data2 <- data
data2$mu <- data2$drop_u
data2$mc <- data2$drop_c
data.t1 <- data[data$time == 1, ]
data.t1$u[is.na(data.t1$u)] <- -999999
data.t1$c[is.na(data.t1$c)] <- -999999
mf_u_fixed <- model.frame(formula = fixed_u, data = data.t1)
mf_c_fixed <- model.frame(formula = fixed_c, data = data.t1)
terms <- NULL
x_u_fixed <- model.matrix(attr(mf_u_fixed, "terms"), data = mf_u_fixed)
x_c_fixed <- model.matrix(attr(mf_c_fixed, "terms"), data = mf_c_fixed)
if("u" %in% names(mf_c_fixed)){
mf_c_fixed$u[index_mis_u] <- NA
}
name_re_u_coeff <- NULL
name_re_c_coeff <- NULL
if(!is.null(random_u)){
name_re_u_coeff <- sub("\\|.*", "", random_u)
if(grepl("0 + 1", name_re_u_coeff, fixed = TRUE) == TRUE) { stop("Either remove or add the random intercept")}
name_clus_u <- sub('.*\\|', '', random_u)
if(lengths(strsplit(name_clus_u, " ")) > 2) { stop("a single clustering variable must selected for each formula") }
name_clus_u <- gsub(" ", "", name_clus_u, fixed = TRUE)
if(!name_clus_u %in% names(cov_matrix)) { stop("the clustering variable must be among the variables in the dataset") }
if(strsplit(name_re_u_coeff, "")[[1]][1] == 0) {
no_random_int_u <- TRUE} else {no_random_int_u <- FALSE }
if(no_random_int_u == TRUE) {
name_re_u_coeff <- sub("[0]", "", name_re_u_coeff)
name_re_u_coeff <- sub("[+]", "", name_re_u_coeff)
}
if(name_re_u_coeff == "" | name_re_u_coeff == " ") { stop("please state for which variables the random effects are assumed") }
fname_re_u_coeff <- as.formula(paste("u", name_re_u_coeff, sep = " ~ "))
if(all(names(model.frame(fname_re_u_coeff, data = data.t1)) %in% c("0", "1", names(model.frame(fixed_u, data = data.t1)))) == FALSE) {
stop("only covariates defined as fixed effects can be included in the random effects model")
}
if("u" %in% labels(terms(fname_re_u_coeff))) {
stop("please remove 'u' from the random effects expression of model.eff")
}
if("c" %in% labels(terms(fname_re_u_coeff))) {
stop("dependence allowed only through the cost model; please remove 'c' from model.eff")
}
mf_u_random <- model.frame(formula = fname_re_u_coeff, data = data.t1)
x_u_random <- model.matrix(attr(mf_u_random, "terms"), data = mf_u_random)
if(no_random_int_u == TRUE) {
x_u_random <- as.matrix(x_u_random[, !colnames(x_u_random) == "(Intercept)"])
if(is.null(colnames(x_u_random)) == TRUE & dim(x_u_random)[2] == 1) {
colnames(x_u_random) <- gsub(" ", "", name_re_u_coeff)
}
}
clus_u <- data[, name_clus_u]
if(!is.factor(clus_u)) { stop("clustering variables must be defined as factors") }
clusn_u <- as.numeric(clus_u)
if(!all(diff(sort(unique(clusn_u))) == 1) | !min(clusn_u) == 1) {
stop("ordered levels of clustering variables must not have gaps and must start from 1")
}
}
if(!is.null(random_c)){
name_re_c_coeff <- sub("\\|.*", "", random_c)
if(grepl("0 + 1", name_re_c_coeff, fixed = TRUE) == TRUE) { stop("Either remove or add the random intercept")}
name_clus_c <- sub('.*\\|', '', random_c)
if(lengths(strsplit(name_clus_c, " ")) > 2) { stop("a single clustering variable must selected for each formula") }
name_clus_c <- gsub(" ", "", name_clus_c, fixed = TRUE)
if(!name_clus_c %in% names(cov_matrix)) { stop("the clustering variable must be among the variables in the dataset") }
if(strsplit(name_re_c_coeff, "")[[1]][1] == 0) {
no_random_int_c <- TRUE} else {no_random_int_c <- FALSE }
if(no_random_int_c == TRUE) {
name_re_c_coeff <- sub("[0]", "", name_re_c_coeff)
name_re_c_coeff <- sub("[+]", "", name_re_c_coeff)
}
if(name_re_c_coeff == "" | name_re_c_coeff == " ") { stop("please state for which variables the random effects are assumed") }
if(gsub(" ", "", name_re_c_coeff) == "u" & no_random_int_c == FALSE) {name_re_c_coeff <- "1 + u" }
fname_re_c_coeff <- as.formula(paste("c", name_re_c_coeff, sep = " ~ "))
if(all(names(model.frame(fname_re_c_coeff, data = data.t1)) %in% c("0", "1", names(model.frame(fixed_c, data = data.t1)))) == FALSE) {
stop("only covariates defined as fixed effects can be included in the random effects model")
}
if("c" %in% labels(terms(fname_re_c_coeff))) {
stop("please remove 'c' from the random effects expression of model.cost")
}
if("u" %in% labels(terms(fname_re_c_coeff))) {
if(length(grep(":u", labels(terms(fname_re_c_coeff)))) != 0 | length(grep("u:", labels(terms(fname_re_c_coeff)))) != 0) {
stop("no interaction effects for 'u' is allowed")
}
}
mf_c_random <- model.frame(formula = fname_re_c_coeff, data = data.t1)
x_c_random <- model.matrix(attr(mf_c_random, "terms"), data = mf_c_random)
if("u" %in% labels(terms(fname_re_c_coeff)) & length(labels(terms(fname_re_c_coeff))) == 1) {
x_c_random <- subset(x_c_random, select = -c(u))
}
if(no_random_int_c == TRUE) {
x_c_random <- as.matrix(x_c_random[, !colnames(x_c_random) == "(Intercept)"])
if(is.null(colnames(x_c_random)) == TRUE & dim(x_c_random)[2] == 1) {
colnames(x_c_random) <- gsub(" ", "", name_re_c_coeff)
}
}
clus_c <- data[, name_clus_c]
if(!is.factor(clus_c)) { stop("clustering variables must be defined as factors") }
clusn_c <- as.numeric(clus_c)
if(!all(diff(sort(unique(clusn_c))) == 1) | !min(clusn_c) == 1) {
stop("ordered levels of clustering variables must not have gaps and must start from 1")
}
}
N1 <- N2 <- c()
N1 <- length(data.t1$u[data.t1$t == 1])
N2 <- length(data.t1$u[data.t1$t == 2])
N <- c(N1, N2)
n.time <- length(unique(data$time))
eff1.mat <- cost1.mat <- matrix(NA, nrow = N1, ncol = n.time)
eff2.mat <- cost2.mat <- matrix(NA, nrow = N2, ncol = n.time)
m_eff <- count_mis_u
m_cost <- count_mis_c
time1.mat <- matrix(NA, nrow = N1, ncol = n.time)
time2.mat <- matrix(NA, nrow = N2, ncol = n.time)
m_eff1.mat <- m_cost1.mat <- matrix(NA, nrow = N1, ncol = n.time)
m_eff2.mat <- m_cost2.mat <- matrix(NA, nrow = N2, ncol = n.time)
for(i in 1:n.time) {
eff1.mat[, i] <- data$u[data$time == i & data$t == 1]
eff2.mat[, i] <- data$u[data$time == i & data$t == 2]
cost1.mat[, i] <- data$c[data$time == i & data$t == 1]
cost2.mat[, i] <- data$c[data$time == i & data$t == 2]
time1.mat[, i] <- data$time[data$time == i & data$t == 1]
time2.mat[, i] <- data$time[data$time == i & data$t == 2]
m_eff1.mat[, i] <- count_mis_u[data$time == i & data$t == 1]
m_eff2.mat[, i] <- count_mis_u[data$time == i & data$t == 2]
m_cost1.mat[, i] <- count_mis_c[data$time == i & data$t == 1]
m_cost2.mat[, i] <- count_mis_c[data$time == i & data$t == 2]
}
effects <- list(eff1.mat, eff2.mat)
costs <- list(cost1.mat, cost2.mat)
time <- list(time1.mat, time2.mat)
m_eff <- list(m_eff1.mat, m_eff2.mat)
m_cost <- list(m_cost1.mat, m_cost2.mat)
N1_cc <- N2_cc <- N1_mis <- N2_mis <- matrix(NA, nrow = 2, ncol = n.time)
eff1_cc <- eff2_cc <- cost1_cc <- cost2_cc <- vector(mode = "list", length = n.time)
for(i in 1:n.time){
N1_cc[1, i] <- length(as.vector(apply(eff1.mat, 2, na.omit)[[i]]))
N1_cc[2, i] <- length(as.vector(apply(cost1.mat, 2, na.omit)[[i]]))
N2_cc[1, i] <- length(as.vector(apply(eff2.mat, 2, na.omit)[[i]]))
N2_cc[2, i] <- length(as.vector(apply(cost2.mat, 2, na.omit)[[i]]))
eff1_cc[[i]] <- as.vector(apply(eff1.mat, 2, na.omit)[[i]])
cost1_cc[[i]] <- as.vector(apply(cost1.mat, 2, na.omit)[[i]])
eff2_cc[[i]] <- as.vector(apply(eff2.mat, 2, na.omit)[[i]])
cost2_cc[[i]] <- as.vector(apply(cost2.mat, 2, na.omit)[[i]])
}
N_cc <- list(N1_cc, N2_cc)
N1_mis <- N1 - N1_cc
N2_mis <- N2 - N2_cc
N_mis <- list(N1_mis, N2_mis)
eff_cc <- list(eff1_cc, eff2_cc)
cost_cc <- list(cost1_cc, cost2_cc)
cov1_u_fixed <- as.data.frame(x_u_fixed[data.t1$t == 1, ])
names(cov1_u_fixed) <- colnames(x_u_fixed)
cov2_u_fixed <- as.data.frame(x_u_fixed[data.t1$t == 2, ])
names(cov2_u_fixed) <- colnames(x_u_fixed)
cov_u_fixed <- list(cov1_u_fixed, cov2_u_fixed)
x_c_hold_fixed <- x_c_fixed
if("u" %in% colnames(x_c_hold_fixed)) {
x_c_fixed <- subset(x_c_hold_fixed, select = -c(u))
}
cov1_c_fixed <- as.data.frame(x_c_fixed[data.t1$t == 1, ])
names(cov1_c_fixed) <- colnames(x_c_fixed)
cov2_c_fixed <- as.data.frame(x_c_fixed[data.t1$t == 2, ])
names(cov2_c_fixed) <- colnames(x_c_fixed)
cov_c_fixed <- list(cov1_c_fixed, cov2_c_fixed)
covu_fixed <- list(cov1_u_fixed, cov2_u_fixed)
mean_cov_u_fixed <- list(apply(as.matrix(cov1_u_fixed), 2, mean), apply(as.matrix(cov2_u_fixed), 2, mean))
covc_fixed <- list(cov1_c_fixed, cov2_c_fixed)
mean_cov_c_fixed <- list(apply(as.matrix(cov1_c_fixed), 2, mean), apply(as.matrix(cov2_c_fixed), 2, mean))
cov1_u_center_fixed <- as.data.frame(scale(cov1_u_fixed, scale = FALSE))
cov2_u_center_fixed <- as.data.frame(scale(cov2_u_fixed, scale = FALSE))
cov1_u_center_fixed[, 1] <- rep(1, nrow(cov1_u_fixed))
cov2_u_center_fixed[, 1] <- rep(1, nrow(cov2_u_fixed))
cov_u_center_fixed <- list(cov1_u_center_fixed, cov2_u_center_fixed)
mean_cov_u_center_fixed <- list(apply(as.matrix(cov1_u_center_fixed), 2, mean), apply(as.matrix(cov2_u_center_fixed), 2, mean))
cov1_c_center_fixed <- as.data.frame(scale(cov1_c_fixed, scale = FALSE))
cov2_c_center_fixed <- as.data.frame(scale(cov2_c_fixed, scale = FALSE))
cov1_c_center_fixed[, 1] <- rep(1, nrow(cov1_c_fixed))
cov2_c_center_fixed[, 1] <- rep(1, nrow(cov2_c_fixed))
cov_c_center_fixed <- list(cov1_c_center_fixed, cov2_c_center_fixed)
mean_cov_c_center_fixed <- list(apply(as.matrix(cov1_c_center_fixed), 2, mean), apply(as.matrix(cov2_c_center_fixed), 2, mean))
if(center == TRUE) {
cov_u_fixed <- cov_u_center_fixed
cov_c_fixed <- cov_c_center_fixed
mean_cov_u_fixed <- mean_cov_u_center_fixed
mean_cov_c_fixed <- mean_cov_c_center_fixed
}
if(!is.null(random_u)){
cov1_u_random <- as.data.frame(x_u_random[data.t1$t == 1, ])
names(cov1_u_random) <- colnames(x_u_random)
cov2_u_random <- as.data.frame(x_u_random[data.t1$t == 2, ])
names(cov2_u_random) <- colnames(x_u_random)
cov_u_random <- list(cov1_u_random, cov2_u_random)
covu_random <- list(cov1_u_random, cov2_u_random)
mean_cov_u_random <- list(apply(as.matrix(cov1_u_random), 2, mean), apply(as.matrix(cov2_u_random), 2, mean))
cov1_u_center_random <- as.data.frame(scale(cov1_u_random, scale = FALSE))
cov2_u_center_random <- as.data.frame(scale(cov2_u_random, scale = FALSE))
if(no_random_int_u == FALSE) {
cov1_u_center_random[, 1] <- rep(1, nrow(cov1_u_random))
cov2_u_center_random[, 1] <- rep(1, nrow(cov2_u_random))
}
cov_u_center_random <- list(cov1_u_center_random, cov2_u_center_random)
mean_cov_u_center_random <- list(apply(as.matrix(cov1_u_center_random), 2, mean), apply(as.matrix(cov2_u_center_random), 2, mean))
if(center == TRUE) {
cov_u_random <- cov_u_center_random
mean_cov_u_random <- mean_cov_u_center_random
}
clusn_u1 <- clusn_u[data.t1$t == 1]
clusn_u1 <- factor(clusn_u1, levels = unique(clusn_u1))
clusn_u2 <- clusn_u[data.t1$t == 2]
clusn_u2 <- factor(clusn_u2, levels = unique(clusn_u2))
}
if(!is.null(random_c)){
x_c_hold_random <- x_c_random
if("u" %in% colnames(x_c_hold_random)) {
x_c_random <- subset(x_c_hold_random, select = -c(u))
}
cov1_c_random <- as.data.frame(x_c_random[data.t1$t == 1, ])
names(cov1_c_random) <- colnames(x_c_random)
cov2_c_random <- as.data.frame(x_c_random[data.t1$t == 2, ])
names(cov2_c_random) <- colnames(x_c_random)
cov_c_random <- list(cov1_c_random, cov2_c_random)
covc_random <- list(cov1_c_random, cov2_c_random)
mean_cov_c_random <- list(apply(as.matrix(cov1_c_random), 2, mean), apply(as.matrix(cov2_c_random), 2, mean))
cov1_c_center_random <- as.data.frame(scale(cov1_c_random, scale = FALSE))
cov2_c_center_random <- as.data.frame(scale(cov2_c_random, scale = FALSE))
if(no_random_int_c == FALSE) {
cov1_c_center_random[, 1] <- rep(1, nrow(cov1_c_random))
cov2_c_center_random[, 1] <- rep(1, nrow(cov2_c_random))
}
cov_c_center_random <- list(cov1_c_center_random, cov2_c_center_random)
mean_cov_c_center_random <- list(apply(as.matrix(cov1_c_center_random), 2, mean), apply(as.matrix(cov2_c_center_random), 2, mean))
if(center == TRUE) {
cov_c_random <- cov_c_center_random
mean_cov_c_random <- mean_cov_c_center_random
}
clusn_c1 <- clusn_c[data.t1$t == 1]
clusn_c1 <- factor(clusn_c1, levels = unique(clusn_c1))
clusn_c2 <- clusn_c[data.t1$t == 2]
clusn_c2 <- factor(clusn_c2, levels = unique(clusn_c2))
}
data2.t1 <- data[data$time == 1, ]
data2.t1$u[is.na(data2.t1$u) == TRUE] <- -999999
data2.t1$c[is.na(data2.t1$c) == TRUE] <- -999999
data2.t1$mu <- c(m_eff1.mat[,1], m_eff2.mat[,1])
data2.t1$mc <- c(m_cost1.mat[,1], m_cost2.mat[,1])
if(!is.formula(model.mu) | !is.formula(model.mc)) {
stop("model.mu and/or model.mc must be formula objects")
}
fixed_mu <- nobars_(model.mu)
fixed_mc <- nobars_(model.mc)
random_mu <- fb(model.mu)
random_mc <- fb(model.mc)
if(!is.null(random_mu) & length(random_mu) > 1 | !is.null(random_mc) & length(random_mc) > 1) {
stop("random effects can be included in the formula only through a single expression within brackets")
}
if(all(names(model.frame(fixed_mu, data = data2.t1)) %in% c("mu", "u", names(cov_matrix))) == FALSE |
all(names(model.frame(fixed_mc, data = data2.t1)) %in% c("mc", "c", names(cov_matrix))) == FALSE) {
stop("partially-observed covariates cannot be included in the model")
}
if(all(names(model.frame(fixed_mu, data = data2.t1)) %in% names(data2.t1)) == FALSE |
all(names(model.frame(fixed_mc, data = data2.t1)) %in% names(data2.t1)) == FALSE) {
stop("you must provide names in the formula that correspond to those in the data")
}
if(names(model.frame(fixed_mu, data = data2.t1)[1]) != "mu") {
stop("you must set 'mu' as the response in the formula model.mu")
}
if(names(model.frame(fixed_mc, data = data2.t1)[1]) != "mc") {
stop("you must set 'mc' as the response in the formula model.mc")
}
if("t" %in% names(model.frame(fixed_mc, data = data2.t1)) | "t" %in% names(model.frame(fixed_mu, data = data2.t1))) {
stop("treatment indicator must be provided only in the data. Please remove 't' from 'model.mu' and/or 'model.mc'")
}
if("c" %in% names(model.frame(fixed_mu, data = data2.t1)) | "u" %in% names(model.frame(fixed_mc, data = data2.t1))) {
stop("please remove 'u' from model.mc and/or remove 'c' from model.mu")
}
if("u" %in% labels(terms(fixed_mu))) {
if(length(grep(":u", labels(terms(fixed_mu)))) != 0 | length(grep("u:", labels(terms(fixed_mu)))) != 0) {
stop("no interaction effects for 'u' is allowed")
}
}
if("c" %in% labels(terms(fixed_mc))) {
if(length(grep(":c", labels(terms(fixed_mc)))) != 0 | length(grep("c:", labels(terms(fixed_mc)))) != 0) {
stop("no interaction effects for 'c' is allowed")
}
}
name_re_mu_coeff <- NULL
name_re_mc_coeff <- NULL
if(!is.null(random_mu)){
name_re_mu_coeff <- sub("\\|.*", "", random_mu)
if(grepl("0 + 1", name_re_mu_coeff, fixed = TRUE) == TRUE) { stop("Either remove or add the random intercept")}
name_clus_mu <- sub('.*\\|', '', random_mu)
if(lengths(strsplit(name_clus_mu, " ")) > 2) {stop("a single clustering variable must be selected for each formula") }
name_clus_mu <- gsub(" ", "", name_clus_mu, fixed = TRUE)
if(!name_clus_mu %in% names(cov_matrix)) { stop("the clustering variable must be among the variables in the dataset") }
if(strsplit(name_re_mu_coeff, "")[[1]][1] == 0) {
no_random_int_mu <- TRUE} else {no_random_int_mu <- FALSE }
if(no_random_int_mu == TRUE) {
name_re_mu_coeff <- sub("[0]", "", name_re_mu_coeff)
name_re_mu_coeff <- sub("[+]", "", name_re_mu_coeff)
}
if(name_re_mu_coeff == "" | name_re_mu_coeff == " ") { stop("please state for which variables the random effects are assumed") }
if(gsub(" ", "", name_re_mu_coeff) == "u" & no_random_int_mu == FALSE) {name_re_mu_coeff <- "1 + u" }
fname_re_mu_coeff <- as.formula(paste("mu", name_re_mu_coeff, sep = " ~ "))
if(all(names(model.frame(fname_re_mu_coeff, data = data2)) %in% c("0","1", names(model.frame(fixed_mu, data = data2)))) == FALSE) {
stop("only covariates inlcued as fixed effects can be included in the random effects model")
}
if("mu" %in% labels(terms(fname_re_mu_coeff))) {
stop("please remove 'mu' from the right hand side of model.mu")
}
if("u" %in% labels(terms(fname_re_mu_coeff))) {
if(length(grep(":u", labels(terms(fname_re_mu_coeff)))) != 0 | length(grep("u:", labels(terms(fname_re_mu_coeff)))) != 0) {
stop("no interaction effects for 'u' are allowed")
}
}
clus_mu <- data[, name_clus_mu]
if(!is.factor(clus_mu)) { stop("clustering variables must be defined as factors") }
clusn_mu <- as.numeric(clus_mu)
if(!all(diff(sort(unique(clusn_mu))) == 1) | !min(clusn_mu) == 1) {
stop("ordered levels of clustering variables must not have gaps and must start from 1")
}
}
if(!is.null(random_mc)){
name_re_mc_coeff <- sub("\\|.*", "", random_mc)
if(grepl("0 + 1", name_re_mc_coeff, fixed = TRUE) == TRUE) { stop("Either remove or add the random intercept")}
name_clus_mc <- sub('.*\\|', '', random_mc)
if(lengths(strsplit(name_clus_mc, " ")) > 2) {stop("a single clustering variable must be selected for each formula") }
name_clus_mc <- gsub(" ", "", name_clus_mc, fixed = TRUE)
if(!name_clus_mc %in% names(cov_matrix)) { stop("the clustering variable must be among the variables in the dataset") }
if(strsplit(name_re_mc_coeff, "")[[1]][1] == 0) {
no_random_int_mc <- TRUE} else {no_random_int_mc <- FALSE }
if(no_random_int_mc == TRUE) {
name_re_mc_coeff <- sub("[0]", "", name_re_mc_coeff)
name_re_mc_coeff <- sub("[+]", "", name_re_mc_coeff)
}
if(name_re_mc_coeff == "" | name_re_mc_coeff == " ") { stop("please state for which variables the random effects are assumed") }
if(gsub(" ", "", name_re_mc_coeff) == "c" & no_random_int_mc == FALSE) {name_re_mc_coeff <- "1 + c" }
fname_re_mc_coeff <- as.formula(paste("mc", name_re_mc_coeff, sep=" ~ "))
if(all(names(model.frame(fname_re_mc_coeff, data = data2)) %in% c("0","1", names(model.frame(fixed_mc, data = data2)))) == FALSE) {
stop("only covariates inlcued as fixed effects can be included in the random effects model")
}
if("mc" %in% labels(terms(fname_re_mc_coeff))) {
stop("please remove 'mc' from the right hand side of model.mc")
}
if("c" %in% labels(terms(fname_re_mc_coeff))) {
if(length(grep(":c", labels(terms(fname_re_mc_coeff)))) != 0 | length(grep("c:", labels(terms(fname_re_mc_coeff)))) != 0) {
stop("no interaction effects for 'c' is allowed")
}
}
clus_mc <- data[, name_clus_mc]
if(!is.factor(clus_mc)) { stop("clustering variables must be defined as factors") }
clusn_mc <- as.numeric(clus_mc)
if(!all(diff(sort(unique(clusn_mc))) == 1) | !min(clusn_mc) == 1) {
stop("ordered levels of clustering variables must not have gaps and must start from 1")
}
}
mf_mu_fixed <- model.frame(formula = fixed_mu, data = data2.t1)
mf_mc_fixed <- model.frame(formula = fixed_mc, data = data2.t1)
z_u_fixed <- model.matrix(attr(mf_mu_fixed, "terms"), data = mf_mu_fixed)
z_c_fixed <- model.matrix(attr(mf_mc_fixed, "terms"), data = mf_mc_fixed)
z_u_hold_fixed <- z_u_fixed
if("u" %in% colnames(z_u_hold_fixed)) {
z_u_fixed <- subset(z_u_hold_fixed, select = -c(u))
}
z_c_hold_fixed <- z_c_fixed
if("c" %in% colnames(z_c_hold_fixed)) {
z_c_fixed <- subset(z_c_hold_fixed, select = -c(c))
}
covz1_u_fixed <- as.data.frame(z_u_fixed[data2.t1$t == 1, ])
names(covz1_u_fixed) <- colnames(z_u_fixed)
covz2_u_fixed <- as.data.frame(z_u_fixed[data2.t1$t == 2, ])
names(covz2_u_fixed) <- colnames(z_u_fixed)
covz_u_fixed <- list(covz1_u_fixed, covz2_u_fixed)
covzu_fixed <- list(covz1_u_fixed, covz2_u_fixed)
mean_covz_u_fixed <- list(apply(as.matrix(covz1_u_fixed), 2, mean), apply(as.matrix(covz2_u_fixed), 2, mean))
names(covz_u_fixed) <- names(mean_covz_u_fixed) <- c("Control", "Intervention")
covz1_c_fixed <- as.data.frame(z_c_fixed[data2.t1$t == 1, ])
names(covz1_c_fixed) <- colnames(z_c_fixed)
covz2_c_fixed <- as.data.frame(z_c_fixed[data2.t1$t == 2, ])
names(covz2_c_fixed) <- colnames(z_c_fixed)
covz_c_fixed <- list(covz1_c_fixed, covz2_c_fixed)
covzc_fixed <- list(covz1_c_fixed, covz2_c_fixed)
mean_covz_c_fixed <- list(apply(as.matrix(covz1_c_fixed), 2, mean), apply(as.matrix(covz2_c_fixed), 2, mean))
covz1_u_center_fixed <- as.data.frame(scale(covz1_u_fixed, scale = FALSE))
covz2_u_center_fixed <- as.data.frame(scale(covz2_u_fixed, scale = FALSE))
covz1_u_center_fixed[, 1] <- rep(1, nrow(covz1_u_fixed))
covz2_u_center_fixed[, 1] <- rep(1, nrow(covz2_u_fixed))
covz_u_center_fixed <- list(covz1_u_center_fixed, covz2_u_center_fixed)
mean_covz_u_center_fixed <- list(apply(as.matrix(covz1_u_center_fixed), 2, mean), apply(as.matrix(covz2_u_center_fixed), 2, mean))
covz1_c_center_fixed <- as.data.frame(scale(covz1_c_fixed, scale = FALSE))
covz2_c_center_fixed <- as.data.frame(scale(covz2_c_fixed, scale = FALSE))
covz1_c_center_fixed[, 1] <- rep(1, nrow(covz1_c_fixed))
covz2_c_center_fixed[, 1] <- rep(1, nrow(covz2_c_fixed))
covz_c_center_fixed <- list(covz1_c_center_fixed, covz2_c_center_fixed)
mean_covz_c_center_fixed <- list(apply(as.matrix(covz1_c_center_fixed), 2, mean), apply(as.matrix(covz2_c_center_fixed), 2, mean))
if(center == TRUE) {
covz_u_fixed <- covz_u_center_fixed
covz_c_fixed <- covz_c_center_fixed
mean_covz_u_fixed <- mean_covz_u_center_fixed
mean_covz_c_fixed <- mean_covz_c_center_fixed
}
if(!is.null(random_mu)){
mf_mu_random <- model.frame(formula = fname_re_mu_coeff, data = data2.t1)
z_u_random <- model.matrix(attr(mf_mu_random, "terms"), data = mf_mu_random)
if(no_random_int_mu == TRUE) {
z_u_random <- as.matrix(z_u_random[, !colnames(z_u_random) == "(Intercept)"])
if(dim(z_u_random)[2] == 1) { colnames(z_u_random) <- colnames(model.matrix(attr(mf_mu_random, "terms"), data = mf_mu_random))[2] }
}
z_u_hold_random <- z_u_random
if("u" %in% colnames(z_u_hold_random)) {
z_u_random <- subset(z_u_hold_random, select = -c(u))
}
covz1_u_random <- as.data.frame(z_u_random[data2.t1$t == 1, ])
names(covz1_u_random) <- colnames(z_u_random)
covz2_u_random <- as.data.frame(z_u_random[data2.t1$t == 2, ])
names(covz2_u_random) <- colnames(z_u_random)
covz_u_random <- list(covz1_u_random, covz2_u_random)
covzu_random <- list(covz1_u_random, covz2_u_random)
mean_covz_u_random <- list(apply(as.matrix(covz1_u_random), 2, mean), apply(as.matrix(covz2_u_random), 2, mean))
names(covz_u_random) <- names(mean_covz_u_random) <- c("Control", "Intervention")
covz1_u_center_random <- as.data.frame(scale(covz1_u_random, scale = FALSE))
covz2_u_center_random <- as.data.frame(scale(covz2_u_random, scale = FALSE))
if(no_random_int_mu == FALSE) {
covz1_u_center_random[, 1] <- rep(1, nrow(covz1_u_random))
covz2_u_center_random[, 1] <- rep(1, nrow(covz2_u_random))
}
covz_u_center_random <- list(covz1_u_center_random, covz2_u_center_random)
mean_covz_u_center_random <- list(apply(as.matrix(covz1_u_center_random), 2, mean), apply(as.matrix(covz2_u_center_random), 2, mean))
if(center == TRUE) {
covz_u_random <- covz_u_center_random
mean_covz_u_random <- mean_covz_u_center_random
}
clusn_mu1 <- clusn_mu[data2.t1$t == 1]
clusn_mu1 <- factor(clusn_mu1, levels = unique(clusn_mu1))
clusn_mu2 <- clusn_mu[data2.t1$t == 2]
clusn_mu2 <- factor(clusn_mu2, levels = unique(clusn_mu2))
}
if(!is.null(random_mc)){
mf_mc_random <- model.frame(formula = fname_re_mc_coeff, data = data2.t1)
z_c_random <- model.matrix(attr(mf_mc_random, "terms"), data = mf_mc_random)
if(no_random_int_mc == TRUE) {
z_c_random <- as.matrix(z_c_random[, !colnames(z_c_random) == "(Intercept)"])
if(dim(z_c_random)[2] == 1) { colnames(z_c_random) <- colnames(model.matrix(attr(mf_mc_random, "terms"), data = mf_mc_random))[2] }
}
z_c_hold_random <- z_c_random
if("c" %in% colnames(z_c_hold_random)) {
z_c_random <- subset(z_c_hold_random, select = -c(c))
}
covz1_c_random <- as.data.frame(z_c_random[data2.t1$t == 1, ])
names(covz1_c_random) <- colnames(z_c_random)
covz2_c_random <- as.data.frame(z_c_random[data2.t1$t == 2, ])
names(covz2_c_random) <- colnames(z_c_random)
covz_c_random <- list(covz1_c_random, covz2_c_random)
covzc_random <- list(covz1_c_random, covz2_c_random)
mean_covz_c_random <- list(apply(as.matrix(covz1_c_random), 2, mean), apply(as.matrix(covz2_c_random), 2, mean))
names(covz_c_random) <- names(mean_covz_c_random) <- c("Control", "Intervention")
covz1_c_center_random <- as.data.frame(scale(covz1_c_random, scale = FALSE))
covz2_c_center_random <- as.data.frame(scale(covz2_c_random, scale = FALSE))
if(no_random_int_mc == FALSE) {
covz1_c_center_random[, 1] <- rep(1, nrow(covz1_c_random))
covz2_c_center_random[, 1] <- rep(1, nrow(covz2_c_random))
}
covz_c_center_random <- list(covz1_c_center_random, covz2_c_center_random)
mean_covz_c_center_random <- list(apply(as.matrix(covz1_c_center_random), 2, mean), apply(as.matrix(covz2_c_center_random), 2, mean))
if(center == TRUE) {
covz_c_random <- covz_c_center_random
mean_covz_c_random <- mean_covz_c_center_random
}
clusn_mc1 <- clusn_mc[data2.t1$t == 1]
clusn_mc1 <- factor(clusn_mc1, levels = unique(clusn_mc1))
clusn_mc2 <- clusn_mc[data2.t1$t == 2]
clusn_mc2 <- factor(clusn_mc2, levels = unique(clusn_mc2))
}
names(covz_u_fixed) <- names(mean_covz_u_fixed) <- names(covz_c_fixed) <- names(mean_covz_c_fixed) <- c("Control", "Intervention")
names(cov_u_fixed) <- names(cov_c_fixed) <- names(mean_cov_u_fixed) <- names(mean_cov_c_fixed) <- c("Control", "Intervention")
if(!is.null(random_c)) {
names(cov_c_random) <- names(mean_cov_c_random) <- c("Control", "Intervention")
clusn_c <- list("Control" = clusn_c1, "Intervention" = clusn_c2)
} else {cov_c_random <- mean_cov_c_random <- NULL}
if(!is.null(random_mc)) {
names(covz_c_random) <- names(mean_covz_c_random) <- c("Control", "Intervention")
clusn_mc <- list("Control" = clusn_mc1, "Intervention" = clusn_mc2)
} else {covz_c_random <- mean_covz_c_random <- NULL}
if(!is.null(random_u)) {
names(cov_u_random) <- names(mean_cov_u_random) <- c("Control", "Intervention")
clusn_u <- list("Control" = clusn_u1, "Intervention" = clusn_u2)
} else {cov_u_random <- mean_cov_u_random <- NULL}
if(!is.null(random_mu)) {
names(covz_u_random) <- names(mean_covz_u_random) <- c("Control", "Intervention")
clusn_mu <- list("Control" = clusn_mu1, "Intervention" = clusn_mu2)
} else {covz_u_random <- mean_covz_u_random <- NULL}
names(m_eff) <- names(m_cost) <- names(time) <- c("Control", "Intervention")
names(effects) <- names(costs) <- names(eff_cc) <- names(cost_cc) <- c("Control", "Intervention")
data_raw <- list("raw_effects" = effects, "raw_costs" = costs, "raw_effects_cc" = eff_cc, "raw_costs_cc" = cost_cc, "arm_lengths" = N,
"arm_lengths_cc" = N_cc, "arm_missing_data" = N_mis, "missing_effects" = m_eff, "missing_costs" = m_cost, "time" = time,
"covariates_effects_fixed" = cov_u_fixed, "covariates_costs_fixed" = cov_c_fixed, "mean_cov_effects_fixed" = mean_cov_u_fixed, "mean_cov_costs_fixed" = mean_cov_c_fixed,
"covariates_missing_effects_fixed" = covz_u_fixed, "mean_cov_missing_effects_fixed" = mean_covz_u_fixed, "covariates_missing_costs_fixed" = covz_c_fixed,
"mean_cov_missing_costs_fixed" = mean_covz_c_fixed, "covariates_effects_random" = cov_u_random, "covariates_costs_random" = cov_c_random, "mean_cov_effects_random" = mean_cov_u_random, "mean_cov_costs_random" = mean_cov_c_random,
"covariates_missing_effects_random" = covz_u_random, "mean_cov_missing_effects_random" = mean_covz_u_random, "covariates_missing_costs_random" = covz_c_random,
"mean_cov_missing_costs_random" = mean_covz_c_random, "clus_u" = clusn_u, "clus_c" = clusn_c, "clus_mu" = clusn_mu, "clus_mc" = clusn_mc, "data_ind" = data2.t1, "data_long" = data)
model_formula <- list("mf_model.u_fixed" = fixed_u, "mf_model.c_fixed" = fixed_c, "mf_model.mu_fixed" = fixed_mu, "mf_model.mc_fixed" = fixed_mc,
"mf_model.u_random" = fname_re_u_coeff, "mf_model.c_random" = fname_re_c_coeff, "mf_model.mu_random" = fname_re_mu_coeff, "mf_model.mc_random" = fname_re_mc_coeff)
data_list <- list("data_raw" = data_raw, "model_formula" = model_formula)
return(data_list)
}
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.