## Packages ----
(dplyr)
(reshape2)
(nlme)
(mvtnorm)
(e1071)
(knitr)
(car)
(lme4)
(Matrix)
(magic)
(MASS)
(MCMCpack)
## Functions ----
<- (n, k, P, , , = , pi_par = , trt = 1, missing_type, ){
p <- P
n_pattern <- ()
covar_x <- ((n), (n, 1, 0.3))
(covar_x) <- ("x",1:k)
xmat <- ((1,n),covar_x)
y_mat <- (0, n, p)
( == "MVN")
eps <- rmvnorm(n, = (0,p), = )
( == "MVT")
# eps <- rmvt(n, delta = rep(0,p), sigma = sigma, df = 3)
eps <- rmvt(n, delta = (0,p), = /5*3, = 5)
( == "MVGamma"){
u <- rmvnorm(n, = (0,p), = ())
var_eps <- ()
eps <- (0, n, p)
(j 1:p){
shape_para <- 2
scale_para <- (var_eps[j]/shape_para)
eps,j] <- ((u,j]), shape = shape_para, = scale_para) - shape_para*scale_para
}
}
y_mean <- (, 1, (x) (xmat%*%x))
y_mat <- y_mean + eps
mean_y <- (y_mat)
cov_eps <- (eps)
cov_y <- (y_mat)
db <- (y_mat)
(db) <- ("y",1:p)
db <- (db, covar_x)
db$num <- (1:n)
db$id <- (trt, "-", 1:n)
(missing_type == "MCAR"){
db$ <- (1:p, size = n, = , prob = )
}
(missing_type %in% ("MAR","MNAR")){
(pi_par[3] != 0 & missing_type == "MAR"){
("This is not a MAR setting")
}
logit_inv <- (x){
(x) / (1 + (x))
}
.pattern <- (1, n)
(i_missing p:2){
.score <- ((1, db,(i_missing - 1,i_missing)])) %*% pi_par
.pi <- logit_inv((.score))
.pattern <- ( (n = n, size = 1, prob = .pi) == 1, i_missing, .pattern)
}
db$ <- .pattern
}
db_comb <- db
db_comb$trt <- trt
(i 1:(db_comb)){
(j 2:p){
(db_comb$[i] <= j & db_comb$[i] > 1){
db_comb[i, j] <-
}
}
}
db_long <- reshape2::melt(db_comb, id.vars = ("id","pattern", "num", ("x", 1:k), "trt"),
variable.name = ("time") , value.name = "aval")
db_long <- db_long %>% group_by(id) %>%
mutate(
= ()) %>%
ungroup()
((db_comb = db_comb, db_long = db_long,
mean_y = mean_y, cov_eps = cov_eps, cov_y = cov_y))
}
sim_one_extreme <- (n, k, P, , , = , pi_par = , trt = 1, missing_type, , outlier = "none"){
p <- P
n_pattern <- ()
covar_x <- ((n), (n, 1, 0.3))
(covar_x) <- ("x",1:k)
xmat <- ((1,n),covar_x)
y_mat <- (0, n, p)
( == "MVN")
eps <- rmvnorm(n, = (0,p), = )
( == "MVT")
# eps <- rmvt(n, delta = rep(0,p), sigma = sigma, df = 3)
eps <- rmvt(n, delta = (0,p), = /5*3, = 5)
( == "MVGamma"){
u <- rmvnorm(n, = (0,p), = ())
var_eps <- ()
eps <- (0, n, p)
(j 1:p){
shape_para <- 2
scale_para <- (var_eps[j]/shape_para)
eps,j] <- ((u,j]), shape = shape_para, = scale_para) - shape_para*scale_para
}
}
y_mean <- (, 1, (x) (xmat%*%x))
y_mat <- y_mean + eps
mean_y <- (y_mat)
cov_eps <- (eps)
cov_y <- (y_mat)
db <- (y_mat)
(db) <- ("y",1:p)
db <- (db, covar_x)
db$num <- (1:n)
db$id <- (trt, "-", 1:n)
(missing_type == "MCAR"){
db$ <- (1:p, size = n, = , prob = )
}
(missing_type %in% ("MAR","MNAR")){
(pi_par[3] != 0 & missing_type == "MAR"){
("This is not a MAR setting")
}
logit_inv <- (x){
(x) / (1 + (x))
}
.pattern <- (1, n)
(i_missing p:2){
.score <- ((1, db,(i_missing - 1,i_missing)])) %*% pi_par
.pi <- logit_inv((.score))
.pattern <- ( (n = n, size = 1, prob = .pi) == 1, i_missing, .pattern)
}
db$ <- .pattern
}
db_comb <- db
db_comb$trt <- trt
(i 1:(db_comb)){
(j 2:p){
(db_comb$[i] <= j & db_comb$[i] > 1){
db_comb[i, j] <-
}
}
}
(outlier != "none"){
ind_max <- (db_comb,p] == (db_comb,p], na.rm = ))
(outlier == "Type 1"){
db_comb[ind_max, p] <- 3*db_comb[ind_max, p]
}
(outlier == "Type 2"){
db_comb[ind_max, 2:p] <- 3*db_comb[ind_max, 2:p]
}
}
db_long <- reshape2::melt(db_comb, id.vars = ("id","pattern", "num", ("x", 1:k), "trt"),
variable.name = ("time") , value.name = "aval")
db_long <- db_long %>% group_by(id) %>%
mutate(
= (),
trt = trt) %>%
ungroup()
((db_comb = db_comb, db_long = db_long))
}
<- (db_comb, db_long, M){
n1 <- ((db_comb$trt == 1))
n2 <- ((db_comb$trt == 2))
p <- ((db_long$))
obs_pi <- 1 - (((db_comb,p] & db_comb$trt == 1))/n1,
((db_comb,p] & db_comb$trt == 2))/n2)
names_dep_var <- (db_comb, (1:p)])
names_covar <- (db_comb, ((p+1):(p+k))])
db_comb_ctl <- db_comb[which(db_comb$trt == 1),]
db_comb_trt <- db_comb[which(db_comb$trt == 2),]
# names_reg_var <- c(names_covar, names_dep_var[1])
# For control group, conduct sequential for obs; conduct ~ X+Y1 for imputed
dep_var <- ()
(j 1:p){
(j == 1){
dep_var[j]] <- names_covar
}
{
dep_var[j]] <- (dep_var[j-1]], names_dep_var[j-1])
}
}
name_all_var <- (db_comb)
imp_fn <- (){
fit_time_ctl <- ()
fit_time_trt <- ()
beta_time_ctl <- ()
beta_time_trt <- ()
beta_cov_ctl <- ()
beta_cov_trt <- ()
sigma_time_ctl <- (0)
sigma_time_trt <- (0)
db_imp_ctl <- db_comb_ctl
db_imp_trt <- db_comb_trt
(j 1:p){
# control
fit_time_ctl[j]] <- fit_model((dep_var[j]],
name_all_var[j]),
= db_imp_ctl)
beta_time_ctl[j]] <- fit_time_ctl[j]]$
beta_cov_ctl[j]] <- (fit_time_ctl[j]])
sigma_time_ctl[j] <- (fit_time_ctl[j]])$^2
# treatment
fit_time_trt[j]] <- fit_model((dep_var[j]],
name_all_var[j]),
= db_imp_trt)
beta_time_trt[j]] <- fit_time_trt[j]]$
beta_cov_trt[j]] <- (fit_time_trt[j]])
sigma_time_trt[j] <- (fit_time_trt[j]])$^2
(j > 1){
sigma_ini <- sigma_time_ctl[j]
df_fit <- (fit_time_ctl[j]])$[1]
sigma_imp_ctl <- (df_fit-2)*sigma_ini/(1, = df_fit)
imp_mean_ctl <- (fit_time_ctl[j]], newdata = db_imp_ctl)[which(db_imp_ctl$ <= j & db_imp_ctl$ > 1)]
imp_sd_ctl <- ((sigma_imp_ctl), ((db_imp_ctl$ <= j & db_imp_ctl$ > 1)))
db_imp_ctl[which(db_imp_ctl$ <= j & db_imp_ctl$ > 1),j] <- (((db_imp_ctl$ <= j & db_imp_ctl$ > 1)), = imp_mean_ctl, = imp_sd_ctl)
sigma_ini <- sigma_time_trt[j]
df_fit <- (fit_time_trt[j]])$[1]
sigma_imp_trt <- (df_fit-2)*sigma_ini/(1, = df_fit)
imp_mean_trt <- (fit_time_trt[j]], newdata = db_imp_trt)[which(db_imp_trt$ <= j & db_imp_trt$ > 1)]
imp_sd_trt <- ((sigma_imp_trt), ((db_imp_trt$ <= j & db_imp_trt$ > 1)))
db_imp_trt[which(db_imp_trt$ <= j & db_imp_trt$ > 1),j] <- (((db_imp_trt$ <= j & db_imp_trt$ > 1)), = imp_mean_trt, = imp_sd_trt)
}
}
beta_imp_ctl <- ()
beta_imp_trt <- ()
sigma_imp_ctl <- (0)
control_term1 <- (0)
trt_term1 <- (0)
control_term2 <- (0, k, p)
trt_term2 <- (0, k, p)
(j 1:p){
beta_imp_ctl[j]] <- (rmvnorm(1, = beta_time_ctl[j]], = beta_cov_ctl[j]]))
beta_imp_trt[j]] <- (rmvnorm(1, = beta_time_trt[j]], = beta_cov_trt[j]]))
# sigma_ini <- sigma_time_ctl[j]
# df_fit <- summary(fit_time_ctl[[j]])$df[1]
# sigma_imp_ctl[j] <- (df_fit-2)*sigma_ini/rchisq(1, df = df_fit)
(j == 1){
control_term1[j] <- beta_imp_ctl[j]][1]
control_term2,j] <- beta_imp_ctl[j]][2:(k+1)]
trt_term1[j] <- beta_imp_trt[j]][1]
trt_term2,j] <- beta_imp_trt[j]][2:(k+1)]
}
{
sum_term1 <- 0
sum_term2 <- (0,k)
sum_term3 <- 0
sum_term4 <- (0,k)
(i in (k+2):(k+j)){
sum_term1 <- beta_imp_ctl[j]][i]*control_term1[i-k-1] + sum_term1
sum_term2 <- beta_imp_ctl[j]][i]*control_term2,i-k-1] + sum_term2
sum_term3 <- beta_imp_trt[j]][i]*trt_term1[i-k-1] + sum_term3
sum_term4 <- beta_imp_trt[j]][i]*trt_term2,i-k-1] + sum_term4
}
control_term1[j] <- beta_imp_ctl[j]][1] + sum_term1
control_term2,j] <- beta_imp_ctl[j]][2:(k+1)] + sum_term2
trt_term1[j] <- beta_imp_trt[j]][1] + sum_term3
trt_term2,j] <- beta_imp_trt[j]][2:(k+1)] + sum_term4
}
}
sigma_coef <- (p)
(j 1:p){
(j < p){
(ind in (j+1):p){
sum_term <- 0
(i j:(ind-1)){
sum_term <- beta_time_ctl[ind]][i+k+1]*sigma_coef[j,i] + sum_term
}
sigma_coef[j,ind] <- sum_term
}
}
}
component_mat <- (sigma_coef, 2, (x) x^2*sigma_time_ctl)
sigma_total <- (component_mat, 2, )
sigma_ctl <- (sigma_total)
(j1 1:(p-1)){
(j2 in (j1+1):p){
(ind 1:(j2-1)){
sigma_ctl[j1, j2] <- beta_time_ctl[j2]][ind+1+k]*sigma_ctl[j1, ind] + sigma_ctl[j1, j2]
}
sigma_ctl[j2, j1] <- sigma_ctl[j1, j2]
}
}
imp_value_ctl <- db_comb_ctl,p]
imp_value_trt <- db_comb_trt,p]
sigma_pattern <- (0)
mean_ctl <- (0, n1)
mean_trt <- (0, n2)
sigma_mi_ctl <- rwish(p+1, sigma_ctl)/(p+1)
(j 2:p){
sigma_mo <- sigma_mi_ctl[j:p, 1:(j-1)]%*%(sigma_mi_ctl[1:(j-1), 1:(j-1)])
sigma_pattern[j-1] <- (sigma_mi_ctl[j:p, j:p] - sigma_mo%*%sigma_mi_ctl[1:(j-1),j:p])[p-j+1, p-j+1]
(((db_comb_ctl$ == j)) > 0){
pattern_mat_ctl <- db_comb_ctl[which(db_comb_ctl$ == j), ]
mu_miss_ctl <- control_term1[j:p] + (pattern_mat_ctl,(p+1):(p+k)], 1, (x) ((control_term2,j:p], k, (p-j)+1), 2, (y) (x*y)))
mu_obs_ctl <- control_term1[1:(j-1)] + (pattern_mat_ctl,(p+1):(p+k)], 1, (x) ((control_term2,1:(j-1)], k, j-1), 2, (y) (x*y)))
mean_pattern_ctl <- mu_miss_ctl + (((pattern_mat_ctl,1:(j-1)] - (mu_obs_ctl)), (pattern_mat_ctl), j-1), 1, (x) sigma_mo%*%x)
mean_pattern_ctl <- (mean_pattern_ctl, p-j+1, (pattern_mat_ctl))
mean_ctl[which(db_comb_ctl$ == j)] <- mean_pattern_ctl[p-j+1,]
imp_value_ctl[which(db_comb_ctl$ == j)] <- (mean_pattern_ctl[p-j+1,], (x) (1, x, (sigma_pattern[j-1])))
}
(((db_comb_trt$ == j)) > 0){
pattern_mat_trt <- db_comb_trt[which(db_comb_trt$ == j), ]
mu_miss_trt <- control_term1[j:p] + (pattern_mat_trt,(p+1):(p+k)], 1, (x) ((control_term2,j:p], k, (p-j)+1), 2, (y) (x*y)))
mu_obs_trt <- trt_term1[1:(j-1)] + (pattern_mat_trt,(p+1):(p+k)], 1, (x) ((trt_term2,1:(j-1)], k, j-1), 2, (y) (x*y)))
mean_pattern_trt <- mu_miss_trt + (((pattern_mat_trt,1:(j-1)] - (mu_obs_trt)), (pattern_mat_trt), j-1), 1, (x) sigma_mo%*%x)
mean_pattern_trt <- (mean_pattern_trt, p-j+1, (pattern_mat_trt))
mean_trt[which(db_comb_trt$ == j)] <- mean_pattern_trt[p-j+1,]
imp_value_trt[which(db_comb_trt$ == j)] <- (mean_pattern_trt[p-j+1,], (x) (1, x, (sigma_pattern[j-1])))
}
}
imp_value <- (imp_value_ctl, imp_value_trt)
(imp_value)
}
imp_value <- (0, n1+n2, M)
(m 1:M){
imp_value,m] <- imp_fn()
}
# change from baseline
base_value <- db_comb,1]
chg_imp <- (imp_value, 2, (x) x - base_value)
# chg_imp_ctl <- chg_imp[1:n1,]
# chg_imp_trt <- chg_imp[(n1+1):(n1+n2),]
(chg_imp = (chg_imp))
}
<- (db_comb, db_long, chg_imp, M, fit_model){
n1 <- ((db_comb$trt == 1))
n2 <- ((db_comb$trt == 2))
p <- ((db_long$))
obs_pi <- 1 - (((db_comb,p] & db_comb$trt == 1))/n1,
((db_comb,p] & db_comb$trt == 2))/n2)
names_dep_var <- (db_comb, 1:p])
names_covar <- (db_comb, (p+1):(p+k)])
db_comb_ctl <- db_comb[which(db_comb$trt == 1),]
db_comb_trt <- db_comb[which(db_comb$trt == 2),]
chg_imp_ctl <- chg_imp,1:n1]
chg_imp_trt <- chg_imp,(n1+1):(n1+n2)]
## point estimates
# (1) mean estimator
# (i) simple average
theta1 <- (chg_imp_ctl)
theta2 <- (chg_imp_trt)
theta_diff <- theta2 - theta1
mean_est <- (theta1, theta2, theta_diff)
# Rubin's estimate
wm_ctl <- ((chg_imp_ctl, 1, )/n1)
bm_ctl <- ((chg_imp_ctl, 1, ))
var_rubin_ctl <- wm_ctl + (1+1/M)*bm_ctl
wm_trt <- ((chg_imp_trt, 1, )/n2)
bm_trt <- ((chg_imp_trt, 1, ))
var_rubin_trt <- wm_trt + (1+1/M)*bm_trt
wm_diff <- wm_ctl + wm_trt
bm_diff <- ((chg_imp_trt, 1, ) - (chg_imp_ctl, 1, ))
var_rubin_diff <- wm_diff + (1+1/M)*bm_diff
var_mean_rubin <- (var_rubin_ctl, var_rubin_trt, var_rubin_diff)
# (ii) robust regression for each group
db_reg_ctl <- db_comb_ctl,(1,(p+1):(p+k))]
db_reg_trt <- db_comb_trt,(1,(p+1):(p+k))]
adj_mean <- (db_comb,(p+1):(p+k)])
coef1_mat <- ((1, adj_mean), 1, k+1)
coef2_mat <- ((1, adj_mean), 1, k+1)
var_xbar <- (db_comb,(p+1):(p+k)], 2, )/(n1+n2)
rr_fn <- (y){
z1 <- y[1:n1]
z2 <- y(n1+1):(n1+n2)]
db_reg_ctl$wt <- z1
db_reg_trt$wt <- z2
# using Huber weight
fit_ctl <- fit_model(wt ~ x1 + x2, = db_reg_ctl)
# fit_ctl <- lm(wt ~ x1 + x2, data = db_reg_ctl)
# using bisquare weight
# fit_ctl <- rlm(wt ~ x1 + x2, data = db_reg_ctl,
# psi = psi.bisquare)
beta_ctl <- (fit_ctl)
fit_trt <- fit_model(wt ~ x1 + x2, = db_reg_trt)
# fit_trt <- lm(wt ~ x1 + x2, data = db_reg_trt)
beta_trt <- (fit_trt)
theta1 <- beta_ctl[1] + (beta_ctl[2:(k+1)]*adj_mean)
theta2 <- beta_trt[1] + (beta_trt[2:(k+1)]*adj_mean)
theta_diff <- theta2 - theta1
mean_est <- (theta1, theta2, theta_diff)
cov_beta_ctl <- (fit_ctl)
cov_beta_trt <- (fit_trt)
var_theta1 <- (coef1_mat%*%cov_beta_ctl%*%(coef1_mat)) +
(beta_ctl[2:(k+1)]^2*var_xbar)
var_theta2 <- (coef2_mat%*%cov_beta_trt%*%(coef2_mat)) +
(beta_trt[2:(k+1)]^2*var_xbar)
var_diff <- (coef1_mat%*%cov_beta_ctl%*%(coef1_mat)) +
(coef2_mat%*%cov_beta_trt%*%(coef2_mat)) +
((beta_trt[2:(k+1)] - beta_ctl[2:(k+1)])^2*var_xbar)
var_est <- (var_theta1, var_theta2, var_diff)
((mean_est, var_est))
}
mi_reg_mat <- (chg_imp, 1, rr_fn)
reg_est <- (mi_reg_mat)[1:3]
# Rubin's estimate
wm_ctl <- (mi_reg_mat)[4]
bm_ctl <- (mi_reg_mat[1,])
var_rubin_ctl <- wm_ctl + (1+1/M)*bm_ctl
wm_trt <- (mi_reg_mat)[5]
bm_trt <- (mi_reg_mat[2,])
var_rubin_trt <- wm_trt + (1+1/M)*bm_trt
wm_diff <- (mi_reg_mat)[6]
bm_diff <- (mi_reg_mat[3,])
var_rubin_diff <- wm_diff + (1+1/M)*bm_diff
var_reg_rubin <- (var_rubin_ctl, var_rubin_trt, var_rubin_diff)
# (iii) robust regression for the whole data
db_reg <- db_comb,("x1", "x2", "trt")]
coef1_mat <- ((1, adj_mean, 0), 1, k+2)
coef2_mat <- ((1, adj_mean, 1), 1, k+2)
coef3_mat <- coef2_mat - coef1_mat
var_xbar <- (db_comb,(p+1):(p+k)], 2, )/(n1+n2)
rr2_fn <- (y){
z <- y
db_reg$wt <- z
# (1)
# using Huber weight
fit_rr <- fit_model(wt ~ x1 + x2 + (trt), = db_reg)
# fit_rr <- lm(wt ~ x1 + x2 + factor(trt), data = db_reg)
# using bisquare weight
# fit_ctl <- rlm(wt ~ x1 + x2, data = db_reg_ctl,
# psi = psi.bisquare)
beta_rr <- (fit_rr)
theta1 <- beta_rr[1] + (beta_rr[2:(k+1)]*adj_mean)
theta2 <- beta_rr[1] + (beta_rr[2:(k+1)]*adj_mean) + beta_rr[k+2]
theta_diff <- theta2 - theta1
mean_est <- (theta1, theta2, theta_diff)
cov_beta_rr <- (fit_rr)
var_theta1 <- (coef1_mat%*%cov_beta_rr%*%(coef1_mat)) +
(beta_rr[2:(k+1)]^2*var_xbar)
var_theta2 <- (coef2_mat%*%cov_beta_rr%*%(coef2_mat)) +
(beta_rr[2:(k+1)]^2*var_xbar)
var_diff <- (coef3_mat%*%cov_beta_rr%*%(coef3_mat))
var_est <- (var_theta1, var_theta2, var_diff)
((mean_est, var_est))
}
mi_reg2_mat <- (chg_imp, 1, rr2_fn)
reg2_est <- (mi_reg2_mat)[1:3]
# Rubin's estimate
wm_ctl <- (mi_reg2_mat)[4]
bm_ctl <- (mi_reg2_mat[1,])
var_rubin_ctl <- wm_ctl + (1+1/M)*bm_ctl
wm_trt <- (mi_reg2_mat)[5]
bm_trt <- (mi_reg2_mat[2,])
var_rubin_trt <- wm_trt + (1+1/M)*bm_trt
wm_diff <- (mi_reg2_mat)[6]
bm_diff <- (mi_reg2_mat[3,])
var_rubin_diff <- wm_diff + (1+1/M)*bm_diff
var_reg2_rubin <- (var_rubin_ctl, var_rubin_trt, var_rubin_diff)
((mean_est = mean_est,
var_mean_rubin = var_mean_rubin,
reg_est = reg_est,
var_reg_rubin = var_reg_rubin,
reg2_est = reg2_est,
var_reg2_rubin = var_reg2_rubin))
}
rubin_mean <- (db_comb, db_long, chg_imp, fit_model){
n1 <- ((db_comb$trt == 1))
n2 <- ((db_comb$trt == 2))
p <- ((db_long$))
obs_pi <- 1 - (((db_comb,p] & db_comb$trt == 1))/n1,
((db_comb,p] & db_comb$trt == 2))/n2)
names_dep_var <- (db_comb, (1:p)])
names_covar <- (db_comb, ((p+1):(p+k))])
db_comb_ctl <- db_comb[which(db_comb$trt == 1),]
db_comb_trt <- db_comb[which(db_comb$trt == 2),]
chg_imp_ctl <- chg_imp,1:n1]
chg_imp_trt <- chg_imp,(n1+1):(n1+n2)]
## point estimates
# (1) mean estimator
# (i) simple average
theta1 <- (chg_imp_ctl)
theta2 <- (chg_imp_trt)
theta_diff <- theta2 - theta1
mean_est <- (theta1, theta2, theta_diff)
# (ii) robust regression for each group
db_reg_ctl <- db_comb_ctl,(p+1):(p+k)]
db_reg_trt <- db_comb_trt,(p+1):(p+k)]
adj_mean <- (db_comb,(p+1):(p+k)])
rr_fn <- (y){
z1 <- y[1:n1]
z2 <- y(n1+1):(n1+n2)]
db_reg_ctl$wt <- z1
db_reg_trt$wt <- z2
# using Huber weight
fit_ctl <- fit_model(wt ~ x1 + x2, = db_reg_ctl)
# fit_ctl <- lm(wt ~ x1 + x2, data = db_reg_ctl)
# using bisquare weight
# fit_ctl <- rlm(wt ~ x1 + x2, data = db_reg_ctl,
# psi = psi.bisquare)
beta_ctl <- (fit_ctl)
fit_trt <- fit_model(wt ~ x1 + x2, = db_reg_trt)
# fit_trt <- lm(wt ~ x1 + x2, data = db_reg_trt)
beta_trt <- (fit_trt)
theta1 <- beta_ctl[1] + (beta_ctl[2:(k+1)]*adj_mean)
theta2 <- beta_trt[1] + (beta_trt[2:(k+1)]*adj_mean)
theta_diff <- theta2 - theta1
mean_est <- (theta1, theta2, theta_diff)
(mean_est)
}
mi_reg_mat <- (chg_imp, 1, rr_fn)
reg_est <- (mi_reg_mat)
# (iii) robust regression for the whole data
db_reg <- db_comb,("x1", "x2", "trt")]
coef1_mat <- ((1, adj_mean, 0), 1, k+2)
coef2_mat <- ((1, adj_mean, 1), 1, k+2)
coef3_mat <- coef2_mat - coef1_mat
var_xbar <- (db_comb,(p+1):(p+k)], 2, )/(n1+n2)
rr2_fn <- (y){
z <- y
db_reg$wt <- z
# (1)
# using Huber weight
fit_rr <- fit_model(wt ~ x1 + x2 + (trt), = db_reg)
# fit_rr <- lm(wt ~ x1 + x2 + factor(trt), data = db_reg)
# using bisquare weight
# fit_ctl <- rlm(wt ~ x1 + x2, data = db_reg_ctl,
# psi = psi.bisquare)
beta_rr <- (fit_rr)
theta1 <- beta_rr[1] + (beta_rr[2:(k+1)]*adj_mean)
theta2 <- beta_rr[1] + (beta_rr[2:(k+1)]*adj_mean) + beta_rr[k+2]
theta_diff <- theta2 - theta1
mean_est <- (theta1, theta2, theta_diff)
(mean_est)
}
mi_reg2_mat <- (chg_imp, 1, rr2_fn)
reg2_est <- (mi_reg2_mat)
((mean_est = mean_est,
reg_est = reg_est,
reg2_est = reg2_est))
}
nonpara_mi_fn <- (db_comb, db_long, M, B){
db_comb_ctl <- db_comb[which(db_comb$trt == 1),]
db_comb_trt <- db_comb[which(db_comb$trt == 2),]
mean_boot <- (0, 3, B)
reg_boot <- (0, 3, B)
reg2_boot <- (0, 3, B)
(b 1:B){
(b)
id_boot <- ((db_comb$id), = )
data_boot <- (id_new = 1:(id_boot), id = id_boot)
data_boot <- (data_boot, db_comb, all.x = )
data_boot <- (data_boot,-(1:2)], data_boot,1:2])
(data_boot) <- ((data_boot)[1:((data_boot)-2)],"id_old", "id")
data_long <- (id_new = 1:(id_boot), id = id_boot)
data_long <- (data_long, db_long, all.x = )
(data_long) <- ("id_old", "id", (data_long)[3:(data_long)])
chg_imp <- (data_boot, data_long, M)
mi_res <- rubin_mean(data_boot, data_long, chg_imp, fit_model)
mean_boot,b] <- mi_res$mean_est
reg_boot,b] <- mi_res$reg_est
reg2_boot,b] <- mi_res$reg2_est
}
var_mean_boot <- (mean_boot, 1, )
var_reg_boot <- (reg_boot, 1, )
var_reg2_boot <- (reg2_boot, 1, )
((var_mean_boot = var_mean_boot,
var_reg_boot = var_reg_boot,
var_reg2_boot = var_reg2_boot))
}
main<-(seed, case){
#the main fucntion to run for one simulation
(seed)
(missing_type == "MCAR"){
tmp1 <- (n = N, k = k, P = p, = Sigma, = beta_ctl, = Pi_ctl, trt = 1, missing_type = missing_type, = )
tmp2 <- (n = N, k = k, P = p, = Sigma, = beta_trt, = Pi_trt, trt = 2, missing_type = missing_type, = )
}
(missing_type %in% ("MAR", "MNAR")){
(case == "Case0"){
tmp1 <- (n = N, k = k, P = p, = Sigma, = beta_ctl, pi_par = phi_ctl, trt = 1, missing_type = missing_type, = )
tmp2 <- (n = N, k = k, P = p, = Sigma, = beta_trt, pi_par = phi_trt, trt = 2, missing_type = missing_type, = )
}
(case == "Case1"){
tmp1 <- sim_one_extreme(n = N, k = k, P = p, = Sigma, = beta_ctl, pi_par = phi_ctl, trt = 1, missing_type = missing_type, = )
tmp2 <- sim_one_extreme(n = N, k = k, P = p, = Sigma, = beta_trt, pi_par = phi_trt, trt = 2, missing_type = missing_type, = )
}
(case == "Case2"){
tmp1 <- (n = N, k = k, P = p, = Sigma, = beta_ctl, pi_par = phi_ctl, trt = 1, missing_type = missing_type, = )
tmp2 <- sim_one_extreme(n = N, k = k, P = p, = Sigma, = beta_trt, pi_par = phi_trt, trt = 2, missing_type = missing_type, = )
}
(case == "Case3"){
tmp1 <- sim_one_extreme(n = N, k = k, P = p, = Sigma, = beta_ctl, pi_par = phi_ctl, trt = 1, missing_type = missing_type, = )
tmp2 <- (n = N, k = k, P = p, = Sigma, = beta_trt, pi_par = phi_trt, trt = 2, missing_type = missing_type, = )
}
(case == "Case4"){
tmp1 <- sim_one_extreme(n = N, k = k, P = p, = Sigma, = beta_ctl, pi_par = phi_ctl, trt = 1, missing_type = missing_type, = , outlier = "Type 2")
tmp2 <- sim_one_extreme(n = N, k = k, P = p, = Sigma, = beta_trt, pi_par = phi_trt, trt = 2, missing_type = missing_type, = , outlier = "Type 2")
}
(case == "Case5"){
tmp1 <- (n = N, k = k, P = p, = Sigma, = beta_ctl, pi_par = phi_ctl, trt = 1, missing_type = missing_type, = )
tmp2 <- sim_one_extreme(n = N, k = k, P = p, = Sigma, = beta_trt, pi_par = phi_trt, trt = 2, missing_type = missing_type, = , outlier = "Type 2")
}
(case == "Case6"){
tmp1 <- sim_one_extreme(n = N, k = k, P = p, = Sigma, = beta_ctl, pi_par = phi_ctl, trt = 1, missing_type = missing_type, = , outlier = "Type 2")
tmp2 <- (n = N, k = k, P = p, = Sigma, = beta_trt, pi_par = phi_trt, trt = 2, missing_type = missing_type, = )
}
}
db_comb <- (tmp1[1]], tmp2[1]]) # to fit lm
db_long <- (tmp1[2]], tmp2[2]]) # to fit mmrm
## multiple imputation
chg_imp <- (db_comb, db_long, M)
mi_res <- (db_comb, db_long, chg_imp, M, fit_model)
minp_res <- nonpara_mi_fn(db_comb, db_long, M, B)
mean_mi <- mi_res$mean_est
reg_mi <- mi_res$reg_est
reg2_mi <- mi_res$reg2_est
var_rubin_mean <- mi_res$var_mean_rubin
var_rubin_reg <- mi_res$var_reg_rubin
var_rubin_reg2 <- mi_res$var_reg2_rubin
var_boot_mean <- minp_res$var_mean_boot
var_boot_reg <- minp_res$var_reg_boot
var_boot_reg2 <- minp_res$var_reg2_boot
((mean_mi = mean_mi,
reg_mi = reg_mi,
reg2_mi = reg2_mi,
var_rubin_mean = var_rubin_mean,
var_rubin_reg = var_rubin_reg,
var_rubin_reg2 = var_rubin_reg2,
var_boot_mean = var_boot_mean,
var_boot_reg = var_boot_reg,
var_boot_reg2 = var_boot_reg2))
}
## Data ----
N <- 100
k <- 2 # dimension of covariates (omit intercept)
p <- 5 # number of visits
mu_beta_ctl <- (0, 1, 2, 3, 4)
# mu_beta_trt <- c(0, 1.3, 2.8, 4, 5.5)
mu_beta_trt <- (0, 1.3, 2.3, 3.5, 4.8)
(123)
beta_ctl <- ((k+1, mu_beta_ctl[1], 1), (k+1,mu_beta_ctl[2], 1),
(k+1, mu_beta_ctl[3], 1), (k+1,mu_beta_ctl[4], 1),
(k+1, mu_beta_ctl[5], 1))
# beta_trt <- rbind(rnorm(k+1, mu_beta_trt[1], 1), rnorm(k+1,mu_beta_trt[2], 1),
# rnorm(k+1, mu_beta_trt[3], 1), rnorm(k+1,mu_beta_trt[4], 1),
# rnorm(k+1, mu_beta_trt[5], 1))
# beta_trt[1,] <- beta_ctl[1,]
# # Under H0
beta_trt <- beta_ctl
<- (2.0, 1.8, 2.0, 2.1, 2.2)
corr <- (
(1, 0.6, 0.3, 0.2, 0.1,
0.6, 1, 0.7, 0.5, 0.2,
0.3, 0.7, 1, 0.6, 0.4,
0.2, 0.5, 0.6, 1, 0.5,
0.1, 0.2, 0.4, 0.5, 1), 5, 5)
Sigma <- () %*% corr %*% ()
## Missing types
# (1) MCAR
# missing_type = "MCAR"
# (i) case 1
# Pi_ctl <- c(80, rep(5,4)) / 100
# Pi_trt <- c(80, rep(5,4)) / 100
# (ii) case 2
# Pi_ctl <- c(80, rep(5,4)) / 100
# Pi_trt <- c(70, 9, rep(7,3)) / 100
# (iii) case 3
# Pi_ctl <- c(70, 9, rep(7,3)) / 100
# Pi_trt <- c(80, rep(5,4)) / 100
# (iv) case 4
# Pi_ctl <- c(60, rep(10,4)) / 100
# Pi_trt <- c(60, rep(10,4)) / 100
# (2) MAR
missing_type = "MAR"
## (i) ctl > trt, approx 0.8, 0.7
# phi_ctl <- c(-3.2, 0.2, 0)
# phi_trt <- c(-3.5, 0.2, 0)
## (ii) ctl = trt, approx 0.8
# phi_ctl <- c(-3.5, 0.2, 0)
# phi_trt <- c(-3.6, 0.2, 0)
# ## (iii) trt > ctl, approxiamte 0.7, 0.8
# phi_ctl <- c(-2.8, 0.2, 0)
# phi_trt <- c(-4.0, 0.2, 0)
## Under H0, ctl = trt, approx 0.8
phi_ctl <- (-3.5, 0.2, 0)
phi_trt <- (-3.5, 0.2, 0)
## Distribution
= "MVN"
# distribution = "MVT"
# distribution = "MVGamma"
## True value ----
( == "MVN"){
# true_mean <- c(5.526485, 6.604056, 1.077571) # MVN
true_mean <- (5.526485, 5.526485, 0)
}
( == "MVT"){
true_mean <- (5.526554, 6.607928, 1.081374) # MVT
# true_mean <- c(5.526554, 5.526554, 0)
}
( == "MVGamma"){
# true_mean <- c(5.526218, 6.603291, 1.077073) # MVGamma
true_mean <- (5.526218, 5.526218, 0)
}
## Analysis model ----
# fit_model <- lm
fit_model <- MASS::rlm # RR
# fit_model <- Rfit::rfit # rank regression
## Simulation results ----
M <- 10
B <- 100
# Get task id for each simulation job (not run in Rstudio Serve)
task_id <- (("SGE_TASK_ID"))
res <- main(seed = task_id, case = "Case6")
# Save all the objects into i.Rdata
(res, = (task_id, ".Rdata"))
