-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy path5_make_tables.R
128 lines (114 loc) · 5.45 KB
/
5_make_tables.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
rm(list = ls())
# un-comment if running individually
setwd("~/Documents/GitHub/ACEimpute/Manuscript_Simulations/Misspecified_Imputation_Model/")
# load packages
library(tidyverse)
library(latex2exp)
library(xtable)
# create table to display simulation results
create_table <- function(sim_results, true_params) {
bias_df <- sim_results %>%
dplyr::select(-dplyr::starts_with("se_")) %>%
tidyr::gather("param_method", "est", -1) %>%
dplyr::mutate(data = ifelse(test = grepl(pattern = "_fd_", x = param_method),
yes = "Full Data",
no = ifelse(test = grepl(pattern = "_cc_", x = param_method),
yes = "Complete Case",
no = ifelse(test = grepl(pattern = "_cmi_mi_", x = param_method),
yes = "CMI-MI",
no = "CMI"))),
method = paste0(toupper(sub(".*_", "", param_method)), " (", data, ")"),
param = sub("_.*", "", param_method),
truth = unlist(true_params[param]),
bias = est - truth)
see_df <- sim_results %>%
dplyr::select(sim, starts_with("se_")) %>%
dplyr::mutate(sim = 1:dplyr::n()) %>%
tidyr::gather("param_method", "se_est", dplyr::starts_with("se_")) %>%
dplyr::mutate(data = ifelse(test = grepl(pattern = "_fd_", x = param_method),
yes = "Full Data",
no = ifelse(test = grepl(pattern = "_cc_", x = param_method),
yes = "Complete Case",
no = ifelse(test = grepl(pattern = "_cmi_mi_", x = param_method),
yes = "CMI-MI",
no = "CMI"))),
param_method = sub("se_", "", param_method))
mse_df <- bias_df %>%
left_join(see_df, by = c("sim", "data", "param_method")) %>%
mutate(mse = bias^2 + se_est^2,
cover = I(est + qnorm(0.025)*se_est < truth & est + qnorm(0.975)*se_est > truth))
mse_df <- mse_df %>%
dplyr::filter(!(param == "sigma2" & est < 0)) %>%
dplyr::mutate(
data = factor(x = data,
levels = c("CMI", "CMI-MI", "Complete Case", "Full Data"),
labels = c("CMI", "CMI-MI", "Complete Case", "Full Data")),
Method = factor(x = method,
levels = c("ACE (CMI)", "REML (CMI)", "REML (CMI-MI)", "REML (Complete Case)", "REML (Full Data)"),
labels = c("ACE", "CMI", "CMI-MI", "CCA", "Oracle")),
Parameter = factor(param,
levels = c("alpha", "beta1", "sigma2"),
labels = c("alpha", "beta1", "sigma2")))
summary_df = mse_df %>%
group_by(Parameter, Method) %>%
summarize(Bias = mean(bias),
SEE = mean(se_est, na.rm = T),
ESE = sd(est),
MSE = mean(bias^2),
CPr = mean(cover))
return(summary_df)
}
# true parameter values
truth = list(alpha = 1, beta1 = 1, sigma2 = 1)
# create light censoring results table
results_25perc = read.csv("sim_results/full_data_reml_estimates.csv") %>%
select(-censored)
results_25perc = read.csv("sim_results/ace_estimates_25perc.csv") %>%
select(-censored) %>%
left_join(results_25perc, ., by = "sim")
results_25perc = read.csv("sim_results/cmi_mi_estimates_25perc.csv") %>%
select(-censored) %>%
left_join(results_25perc, ., by = "sim")
table_25 = create_table(sim_results = results_25perc, true_params = truth)
# create medium censoring results table
results_50perc = read.csv("sim_results/full_data_reml_estimates.csv") %>%
select(-censored)
results_50perc = read.csv("sim_results/ace_estimates_50perc.csv") %>%
select(-censored) %>%
left_join(results_50perc, ., by = "sim")
results_50perc = read.csv("sim_results/cmi_mi_estimates_50perc.csv") %>%
select(-censored) %>%
left_join(results_50perc, ., by = "sim")
table_50 = create_table(sim_results = results_50perc, true_params = truth)
# create heavy censoring results table
results_75perc = read.csv("sim_results/full_data_reml_estimates.csv") %>%
select(-censored)
results_75perc = read.csv("sim_results/ace_estimates_75perc.csv") %>%
select(-censored) %>%
left_join(results_75perc, ., by = "sim")
results_75perc = read.csv("sim_results/cmi_mi_estimates_75perc.csv") %>%
select(-censored) %>%
left_join(results_75perc, ., by = "sim")
table_75 = create_table(sim_results = results_75perc, true_params = truth)
sink("tables/medium_and_heavy_xtable.txt")
# combine medium and heavy censoring results
rbind(table_50, table_75) %>%
cbind(data.frame(Censoring = c("Medium", rep("", 8),
"Heavy", rep("", 8))), .) %>%
mutate(Parameter = rep(c("alpha", rep("", 2),
"beta", rep("", 2),
"sigma2", rep("", 2)), 2)) %>%
# print xtable for LaTeX use
xtable(x = ., digits = 3) %>%
print(include.rownames = F)
sink()
sink("tables/light_xtable.txt")
table_25 %>%
cbind(data.frame(Censoring = c("Light", rep("", 8))), .) %>%
mutate(Parameter = c("alpha", rep("", 2),
"beta", rep("", 2),
"sigma2", rep("", 2))) %>%
# print xtable for LaTeX use
xtable(x = ., digits = 3) %>%
print(include.rownames = F)
sink()