In this R Markdown file we conduct an econometric analysis of kaizen evidence to estimate 10 incremental effects and two scale parameters (i.e., one for coma comparisons and one for kaizen tasks) of the US EQ-5D-Y-3L data. This software is free for private, non-commercial use only. Please contact the authors for commercial usage.
We estimate a hybrid model where we combine paired comparisons and kaizen tasks data, for the complete and per block data. We then use the bootstrap method to estimate the standard errors. Finally, we analyze the relationship between different sample sizes and the standard errors per block.
#setwd("C:/Users/maksat/Box Sync/USF PC/EuroQol/Kaisen task/paper1/Wave 1/markdown")
setwd("~/markdown")
# clear memory
rm(list = ls())
# use libraries
library(reshape2) # for command dcast and melt
library(dplyr) # for command %>%
library(cfdecomp) # for cluster.resample
library(doParallel) # for doParallel
library(foreach) # for doParallel
library(iterators) # for doParallel
library(parallel) # for doParallel
library(knitr) # for kableBelow, we present data notations and definitions of new variables created.
Subject, i = 1 to N; task, t = 1 to T;
T1 is the number of coma comparison tasks and paired comparison tasks including warmup
T2 is the number of kaizen tasks including warmup
‘set’ is an indicator of set type:
0 paired or coma comparison
1 first kaizen response (1 improvement)
2 second kaizen response (2 improvements)
3 third kaizen response (3 improvements)
‘MO’ is Mobility, Levels 1 to 3
‘SC’ is Self-Care, Levels 1 to 3
‘UA’ is Usual Activities, Levels 1 to 3
‘PD’ is Pain/Discomfort, Levels 1 to 3
‘AD’ is Anxiety Depression, Levels 1 to 3
‘COMA’ is an indicator for coma, 1 if coma; 0 otherwise
‘gamma1’ is the log-scale factor of paired comparisons
‘gamma2’ is the log-scale factor of kaizen tasks
# specify the setting
N = 631 ## number of survey record
Z = 230 ## number of variables per record
Z1 = 131 ## the number of respondent-specific variables
Z2 = 28 ## the number of task-specific variables
Z3 = 71 ## the number of temporal variables (e.g., page times)
T1 = 17 ## number of paired comparisons per respondent
T2 = 11 ## number of kaizen tasks per respondent
spec = c(N, Z, Z1, Z2, Z3, T1, T2) ## summary of setting specification# log-likelihood function of logit with two gammas
llf_g <- function(param0, X) {
param1 <- matrix(param0[1:11],nrow=1)
gamma1 <- as.matrix(
exp(param0[12])*as.matrix(as.numeric(X$set==0),ncol=1)+
exp(param0[13])*as.matrix(as.numeric(X$set!=0),ncol=1),ncol=1)
suffixes <- c(".A1", ".A2", ".A3", ".A4", ".A5", ".A6")
v <- matrix(0, nrow = nrow(X), ncol = length(suffixes)) # value of objects
for (i in seq_along(suffixes)) {
col_indices <- grep(paste0(suffixes[i], "$"), names(X))
subset_X <- X[, col_indices]
v[, i] <- as.matrix(subset_X) %*% t(param1)
v[, i] <- exp(gamma1[,1]*(1-v[, i]))
v[v==exp(gamma1[,1])] <- 0
}
return(colSums(log( v[,1] / (v %*% matrix(data = 1, nrow = 6, ncol = 1))) ))
}
# block estimation function
doBlock <- function(i) {
param0 <- matrix(param0[i, 1:ncol(param0)],nrow=1)
colnames(param0) <- c("MO1","MO2","SC1","SC2","UA1","UA2","PD1","PD2","AD1","AD2",
"COMA","gamma1","gamma2")
start_values <- c(unlist(param0[1, 1:ncol(param0)]))
data_subset <- boot[boot$subjectindex == i, ]
model <- maxLik::maxLik(llf_g, X = data_subset, start = start_values,
method = method_ml, fixed = "COMA",
control = list(iterlim = iteration_ml,
tol = tol_ml, printLevel = 2))
return(matrix(model$estimate[1:13],nrow=1))
}
# block and sample size estimation function
doBlock_SS <- function(i) {
param0 <- matrix(param0[i, 1:ncol(param0)],nrow=1)
colnames(param0) <- c("MO1","MO2","SC1","SC2","UA1","UA2","PD1","PD2","AD1","AD2",
"COMA","gamma1","gamma2")
start_values <- c(unlist(param0[1, 1:ncol(param0)]))
model_size <- list() # Estimates by sample size
result_size <- list()
for (n in c(150,125,100,75,50,40,30,20)) {
data_subset <- boot[(boot$subjectindex == i) & (boot$newID <= n), ]
model_size[[n]] <- maxLik::maxLik(llf_g, X = data_subset, start = start_values,
method = method_ml, fixed = "COMA",
control = list(iterlim = iteration_ml,
tol = tol_ml, printLevel = 2))
param0[1, 1:13] <- matrix(model_size[[n]]$estimate[1:13],nrow=1)
result_size[[n]] <- param0
}
return(do.call(cbind, result_size))
}# load the respondent data
resp_i = read.csv("resp_i_231010.csv")
resp_i <- resp_i[resp_i$incomplete==0,] ## Only completes
# load the task data
task_it = read.csv("task_it_231012.csv")In this section we recode the data to create the variables defined above
# Format from wide to long
id_var <- c("id","kz","tnum2","set")
task_ita = data.frame(melt(task_it, id.vars = id_var, value.name = "A"))
# Create attribute variables MO1 to AD2, LLS and COMA
prefixes <- c("MO", "SC", "UA", "PD", "AD")
for (i in c(1:5)) {
task_ita[paste0(prefixes[[i]], 1)] <- as.numeric(as.numeric(
substr(task_ita$A, i, i)) > 1 )*as.numeric(task_ita$A!="44444")
task_ita[paste0(prefixes[[i]], 2)] <- as.numeric(as.numeric(
substr(task_ita$A, i, i)) > 2 )*as.numeric(task_ita$A!="44444")
}
task_ita$COMA <- as.numeric(task_ita$A=="44444") #create COMA
# Format from long to wide
task_ita <- task_ita[,!(names(task_ita) %in% c("A"))]
attr_var <- c("MO1","MO2","SC1","SC2","UA1","UA2","PD1","PD2","AD1","AD2","COMA")
logit_it <- reshape(task_ita, direction = "wide", idvar = id_var,
timevar = "variable", v.names = attr_var)
logit_it[is.na(logit_it)] <- 0
# Append the block indicator
logit_it <- merge(logit_it, subset(resp_i, select = c(id, subjectindex)),
by="id", all.x = T)
# Select only the coma comparison and kaizen responses
logit_it <- logit_it[(logit_it$set==0 & logit_it$tnum2<=6)|(logit_it$set>0),] In this section, we will estimate 10 incremental effects and two scale parameters by incorporating the coma comparisons with the evidence from the kaizen tasks.
seed <- 231012 # random number generator
tol_ml <- 1e-8 # tolerance criterion for maximum likelihood
method_ml <- "BFGS" # maximization estimation technique for maximum likelihood
iteration_ml <- 150 # number of iteration for maximum likelihoodUsing the conditional logit (CL) model, we will estimate the preference weights β and ancillary parameters (e.g., σ) on a coma scale using the preference evidence from the kaizen tasks.
# CC and KZ Hybrid (two gammas)
g2_var <- c("MO1","MO2","SC1","SC2","UA1","UA2","PD1","PD2","AD1","AD2",
"COMA","gamma1","gamma2")
param0 <- cbind(t(as.matrix(rep(0.1,10), nrow=1)),
t(as.matrix(c(1,0.5,1.5),nrow=1)))
colnames(param0) <- g2_var
g2_model <- (maxLik::maxLik(llf_g, X=logit_it,
start = c(unlist(param0[1,1:ncol(param0)])),
method = method_ml, fixed = "COMA",
control = list(iterlim=iteration_ml,
tol=tol_ml, printLevel=0)))
summary(g2_model) # BFGS results## --------------------------------------------
## Maximum Likelihood estimation
## BFGS maximization, 108 iterations
## Return code 0: successful convergence
## Log-Likelihood: -28092.48
## 12 free parameters
## Estimates:
## Estimate Std. error t value Pr(> t)
## MO1 0.095968 0.006087 15.767 < 2e-16 ***
## MO2 0.152068 0.007057 21.550 < 2e-16 ***
## SC1 0.042086 0.006071 6.932 4.14e-12 ***
## SC2 0.090359 0.006183 14.613 < 2e-16 ***
## UA1 0.059681 0.005927 10.069 < 2e-16 ***
## UA2 0.127845 0.006976 18.327 < 2e-16 ***
## PD1 0.269529 0.009363 28.786 < 2e-16 ***
## PD2 0.331517 0.011849 27.980 < 2e-16 ***
## AD1 0.026954 0.006189 4.355 1.33e-05 ***
## AD2 0.139972 0.006670 20.984 < 2e-16 ***
## COMA 1.000000 0.000000 NA NA
## gamma1 0.510108 0.131691 3.874 0.000107 ***
## gamma2 1.542827 0.038874 39.688 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## --------------------------------------------Next, we analyze the effect of each block on the results.
# Assess the effects of block on the results
param0 <- t(g2_model$estimate[1:13])
colnames(param0) <- g2_var
start_values <- c(unlist(param0[1, 1:ncol(param0)]))
model_list <- list() # Estimates by block
# Iterate through subject indices
for (i in 1:4) {
data_subset <- logit_it[logit_it$subjectindex == i, ]
model_list[[i]] <- maxLik::maxLik(llf_g, X = data_subset, start = start_values,
method = method_ml, fixed = "COMA",
control = list(iterlim = iteration_ml,
tol = tol_ml, printLevel = 0))
}summary(model_list[[1]]) # Block 1## --------------------------------------------
## Maximum Likelihood estimation
## BFGS maximization, 98 iterations
## Return code 0: successful convergence
## Log-Likelihood: -7167.535
## 12 free parameters
## Estimates:
## Estimate Std. error t value Pr(> t)
## MO1 0.10370 0.01188 8.727 < 2e-16 ***
## MO2 0.15348 0.01217 12.615 < 2e-16 ***
## SC1 0.04244 0.01285 3.302 0.000961 ***
## SC2 0.07662 0.01245 6.157 7.42e-10 ***
## UA1 0.07381 0.01228 6.011 1.85e-09 ***
## UA2 0.10191 0.01433 7.109 1.17e-12 ***
## PD1 0.30156 0.01775 16.986 < 2e-16 ***
## PD2 0.29953 0.01751 17.102 < 2e-16 ***
## AD1 0.03381 0.01285 2.630 0.008533 **
## AD2 0.13621 0.01285 10.603 < 2e-16 ***
## COMA 1.00000 0.00000 NA NA
## gamma1 1.06633 0.24295 4.389 1.14e-05 ***
## gamma2 1.48868 0.06422 23.183 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## --------------------------------------------summary(model_list[[2]]) # Block 2## --------------------------------------------
## Maximum Likelihood estimation
## BFGS maximization, 90 iterations
## Return code 0: successful convergence
## Log-Likelihood: -7124.736
## 12 free parameters
## Estimates:
## Estimate Std. error t value Pr(> t)
## MO1 0.14169 0.01440 9.837 < 2e-16 ***
## MO2 0.10904 0.01469 7.425 1.13e-13 ***
## SC1 0.03541 0.01394 2.541 0.0111 *
## SC2 0.10695 0.01266 8.446 < 2e-16 ***
## UA1 0.05933 0.01301 4.560 5.13e-06 ***
## UA2 0.14178 0.01482 9.569 < 2e-16 ***
## PD1 0.29750 0.02003 14.855 < 2e-16 ***
## PD2 0.31644 0.02286 13.845 < 2e-16 ***
## AD1 0.03253 0.01326 2.453 0.0142 *
## AD2 0.11206 0.01294 8.659 < 2e-16 ***
## COMA 1.00000 0.00000 NA NA
## gamma1 0.51835 0.50749 1.021 0.3071
## gamma2 1.47169 0.07643 19.255 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## --------------------------------------------summary(model_list[[3]]) # Block 3## --------------------------------------------
## Maximum Likelihood estimation
## BFGS maximization, 89 iterations
## Return code 0: successful convergence
## Log-Likelihood: -6980.211
## 12 free parameters
## Estimates:
## Estimate Std. error t value Pr(> t)
## MO1 0.07947 0.01334 5.956 2.58e-09 ***
## MO2 0.14987 0.01638 9.149 < 2e-16 ***
## SC1 0.01760 0.01357 1.297 0.195
## SC2 0.08015 0.01323 6.060 1.36e-09 ***
## UA1 0.04815 0.01183 4.069 4.72e-05 ***
## UA2 0.10256 0.01515 6.772 1.27e-11 ***
## PD1 0.23520 0.02224 10.575 < 2e-16 ***
## PD2 0.34112 0.03111 10.966 < 2e-16 ***
## AD1 0.02037 0.01333 1.527 0.127
## AD2 0.12570 0.01478 8.508 < 2e-16 ***
## COMA 1.00000 0.00000 NA NA
## gamma1 0.46376 0.31583 1.468 0.142
## gamma2 1.53532 0.09843 15.599 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## --------------------------------------------summary(model_list[[4]]) # Block 4## --------------------------------------------
## Maximum Likelihood estimation
## BFGS maximization, 99 iterations
## Return code 0: successful convergence
## Log-Likelihood: -6757.411
## 12 free parameters
## Estimates:
## Estimate Std. error t value Pr(> t)
## MO1 0.08181 0.01121 7.298 2.93e-13 ***
## MO2 0.18273 0.01604 11.390 < 2e-16 ***
## SC1 0.06559 0.01132 5.796 6.79e-09 ***
## SC2 0.09707 0.01259 7.710 1.26e-14 ***
## UA1 0.07469 0.01172 6.375 1.83e-10 ***
## UA2 0.15342 0.01454 10.549 < 2e-16 ***
## PD1 0.28349 0.01718 16.500 < 2e-16 ***
## PD2 0.37891 0.02463 15.383 < 2e-16 ***
## AD1 0.03955 0.01189 3.327 0.000878 ***
## AD2 0.18583 0.01359 13.675 < 2e-16 ***
## COMA 1.00000 0.00000 NA NA
## gamma1 0.53134 0.17285 3.074 0.002112 **
## gamma2 1.60380 0.07697 20.838 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## --------------------------------------------We now start the parallel computing to compute the standard errors (SE) using bootstrap method. In our bootstrap method, we will re-sample with replacement and compute the SE for 10 preference weights and 2 scaling parameters.
we first conduct the bootstrap method to the complete data.
reps = 1000
r <- matrix(0,reps,13) # results matrix
param0 <- t(g2_model$estimate[1:13]) # starting values
colnames(param0) <- g2_var
colnames(r) <- g2_var
start_values <- c(unlist(param0[1, 1:ncol(param0)]))
start_time <- Sys.time() # track computational time
for (i in 1:reps) { # About 8 hours
boot <- rbind(cluster.resample(logit_it[logit_it$subjectindex == 1, ],
cluster.name="id", size = NA, newID = F),
cluster.resample(logit_it[logit_it$subjectindex == 2, ],
cluster.name="id", size = NA, newID = F),
cluster.resample(logit_it[logit_it$subjectindex == 3, ],
cluster.name="id", size = NA, newID = F),
cluster.resample(logit_it[logit_it$subjectindex == 4, ],
cluster.name="id", size = NA, newID = F))
final_model <- (maxLik::maxLik(llf_g, X=boot,
start = start_values, method = method_ml, fixed = "COMA",
control = list(iterlim=iteration_ml, tol=tol_ml, printLevel=0)))
r[i,1:13] <- final_model$estimate[1:13]
end_time <- Sys.time()
cat(paste(i, ifelse(i %% 10 == 0,paste(" ",round((end_time - start_time),
digits=2), "\n")," ")))
}boot231016a <-as.data.frame(r)
boot_se = apply(boot231016a, 2, function(x) sd(x))
results_final = round(cbind(g2_model[["estimate"]],boot_se), digits = 3)
colnames(results_final) = c("Estimates", "Bootstrap SE")
kable(results_final, align = "c", caption = "**Summary Statistics of Bootstrap
Method for Final Data**", row.names = TRUE, escape = FALSE, centering = T)| Estimates | Bootstrap SE | |
|---|---|---|
| MO1 | 0.096 | 0.010 |
| MO2 | 0.152 | 0.012 |
| SC1 | 0.042 | 0.009 |
| SC2 | 0.090 | 0.010 |
| UA1 | 0.060 | 0.010 |
| UA2 | 0.128 | 0.010 |
| PD1 | 0.270 | 0.015 |
| PD2 | 0.332 | 0.020 |
| AD1 | 0.027 | 0.011 |
| AD2 | 0.140 | 0.013 |
| COMA | 1.000 | 0.000 |
| gamma1 | 0.510 | 0.131 |
| gamma2 | 1.543 | 0.063 |
we now conduct the bootstrap method for each block.
param0 <- rbind(t(model_list[[1]]$estimate[1:13]),
t(model_list[[2]]$estimate[1:13]),
t(model_list[[3]]$estimate[1:13]),
t(model_list[[4]]$estimate[1:13])) # starting values
colnames(param0) = g2_var
# setup parallel processing
# detectCores()
cl <- makeCluster(detectCores() - 1, type="PSOCK")
registerDoParallel(cl)
r <- matrix(0,reps,52) # results matrix
start_time <- Sys.time() # track computational time
for (i in 1:reps) { # About 3 hours
boot <- rbind(cluster.resample(logit_it[logit_it$subjectindex == 1, ],
cluster.name="id", size = NA, newID = F),
cluster.resample(logit_it[logit_it$subjectindex == 2, ],
cluster.name="id", size = NA, newID = F),
cluster.resample(logit_it[logit_it$subjectindex == 3, ],
cluster.name="id", size = NA, newID = F),
cluster.resample(logit_it[logit_it$subjectindex == 4, ],
cluster.name="id", size = NA, newID = F))
clusterExport(cl=cl, c('boot','llf_g','method_ml',
'iteration_ml','tol_ml','param0'))
r[i,1:52] <- t(matrix(unlist(parLapply(cl, 1:4, doBlock)),ncol=1))
end_time <- Sys.time()
cat(paste(i, ifelse(i %% 10 == 0,paste(" ",round((end_time - start_time),
digits=2), "\n")," ")))
}
stopCluster(cl)boot231016b <-as.data.frame(r) # save the results
boot_se = apply(boot231016b, 2, function(x) sd(x))
boot_se = t(matrix(boot_se, ncol = 13, byrow = T))
results_block = round(cbind(model_list[[1]][["estimate"]], boot_se[,1],
model_list[[2]][["estimate"]], boot_se[,2],
model_list[[3]][["estimate"]], boot_se[,3],
model_list[[4]][["estimate"]], boot_se[,4]), digits = 3)
colnames(results_block) = c("Est1", "SE1", "Est2", "SE2", "Est3", "SE3",
"Est4", "SE4")
kable(results_block, align = "c", caption = "**Summary Statistics of Bootstrap Method
per Block**", row.names = TRUE, escape = FALSE, centering = T)| Est1 | SE1 | Est2 | SE2 | Est3 | SE3 | Est4 | SE4 | |
|---|---|---|---|---|---|---|---|---|
| MO1 | 0.104 | 0.019 | 0.142 | 0.023 | 0.079 | 0.019 | 0.082 | 0.018 |
| MO2 | 0.153 | 0.020 | 0.109 | 0.024 | 0.150 | 0.025 | 0.183 | 0.025 |
| SC1 | 0.042 | 0.016 | 0.035 | 0.019 | 0.018 | 0.020 | 0.066 | 0.017 |
| SC2 | 0.077 | 0.019 | 0.107 | 0.022 | 0.080 | 0.022 | 0.097 | 0.016 |
| UA1 | 0.074 | 0.019 | 0.059 | 0.022 | 0.048 | 0.018 | 0.075 | 0.019 |
| UA2 | 0.102 | 0.017 | 0.142 | 0.025 | 0.103 | 0.022 | 0.153 | 0.023 |
| PD1 | 0.302 | 0.027 | 0.297 | 0.036 | 0.235 | 0.029 | 0.283 | 0.031 |
| PD2 | 0.300 | 0.032 | 0.316 | 0.046 | 0.341 | 0.048 | 0.379 | 0.042 |
| AD1 | 0.034 | 0.021 | 0.033 | 0.023 | 0.020 | 0.021 | 0.040 | 0.023 |
| AD2 | 0.136 | 0.025 | 0.112 | 0.028 | 0.126 | 0.028 | 0.186 | 0.027 |
| COMA | 1.000 | 0.000 | 1.000 | 0.000 | 1.000 | 0.000 | 1.000 | 0.000 |
| gamma1 | 1.066 | 0.235 | 0.518 | 0.447 | 0.464 | 0.269 | 0.531 | 0.176 |
| gamma2 | 1.489 | 0.107 | 1.472 | 0.126 | 1.535 | 0.138 | 1.604 | 0.117 |
Next, we conduct block-specific bootstrap method with different sample sizes. We estimate the SE’s of 10 preference weights and 2 gammas for various sample sizes. We re-sample with replacement for reps number of times and utilize doParallel parallel computing function.
# setup parallel processing
#detectCores()
cl <- makeCluster(detectCores() - 1, type="PSOCK")
registerDoParallel(cl)
r <- matrix(0,reps,416) # results matrix
start_time <- Sys.time() # track computational time
for (i in 1:reps) { # About 14 hours
boot_list <- list()
for (j in 1:4) {
boot <- cluster.resample(logit_it[logit_it$subjectindex == j, ],
cluster.name = "id", size = NA, newID = TRUE)
# create an newID randomizes the order of the rows of each block
boot$newID <- as.numeric(factor(boot$newID[sample(nrow(boot))],
levels = unique(boot$newID)))
boot_list[[j]] <- boot
}
boot <- do.call(rbind, boot_list) # combines the four blocks
param0 <- rbind(t(model_list[[1]]$estimate[1:13]),
t(model_list[[2]]$estimate[1:13]),
t(model_list[[3]]$estimate[1:13]),
t(model_list[[4]]$estimate[1:13])) # starting values
clusterExport(cl=cl, c('boot','llf_g','method_ml',
'iteration_ml','tol_ml','param0'))
r[i,1:416] <- t(matrix(unlist(parLapply(cl, 1:4, doBlock_SS)),ncol=1))
end_time <- Sys.time()
cat(paste(i, ifelse(i %% 10 == 0,paste(" ",round((end_time - start_time),
digits=2), "\n")," ")))
}
stopCluster(cl)boot231016c <-as.data.frame(r) # save the results
# Create a table with the SE for various columns and 13 rows
# 1 full sample
# 2 to 5 each block
# 6 to 13 eight sample sizes for block 1
# 14 to 21 eight sample sizes for block 2
# 22 to 29 eight sample sizes for block 3
# 30 to 37 eight sample sizes for block 4
sd1 <- sapply(boot231016a, sd)
sd2 <- sapply(boot231016b, sd)
sd3 <- sapply(boot231016c, sd)
sd <- cbind(t(sd1),t(sd2),t(sd3))
sd <- matrix(sd,nrow=13)
sd_var <- c("Sample Size","MO1","MO2","SC1","SC2","UA1","UA2","PD1","PD2","AD1","AD2")
sd_size <- c(20,30,40,50,75,100,125,150)
full_sd <- sd[,1:5]
b1_sd <- data.frame(round(t(rbind(matrix(sd_size,nrow=1),sd[,6:13]))[,1:11],
digits = 3))
b2_sd <- data.frame(round(t(rbind(matrix(sd_size,nrow=1),sd[,14:21]))[,1:11],
digits = 3))
b3_sd <- data.frame(round(t(rbind(matrix(sd_size,nrow=1),sd[,22:29]))[,1:11],
digits = 3))
b4_sd <- data.frame(round(t(rbind(matrix(sd_size,nrow=1),sd[,30:37]))[,1:11],
digits = 3))
colnames(b1_sd)<-sd_var
colnames(b2_sd)<-sd_var
colnames(b3_sd)<-sd_var
colnames(b4_sd)<-sd_varkable(b1_sd, align = "l", caption = "**Sample Size and Standard Errors for Block 1**",
row.names = F, escape = FALSE, centering = T)| Sample Size | MO1 | MO2 | SC1 | SC2 | UA1 | UA2 | PD1 | PD2 | AD1 | AD2 |
|---|---|---|---|---|---|---|---|---|---|---|
| 20 | 0.058 | 0.063 | 0.046 | 0.049 | 0.052 | 0.051 | 0.142 | 0.129 | 0.054 | 0.064 |
| 30 | 0.035 | 0.038 | 0.031 | 0.033 | 0.032 | 0.036 | 0.061 | 0.065 | 0.037 | 0.041 |
| 40 | 0.031 | 0.033 | 0.028 | 0.029 | 0.028 | 0.031 | 0.051 | 0.056 | 0.032 | 0.036 |
| 50 | 0.028 | 0.029 | 0.026 | 0.027 | 0.026 | 0.028 | 0.044 | 0.050 | 0.029 | 0.033 |
| 75 | 0.024 | 0.024 | 0.022 | 0.023 | 0.023 | 0.023 | 0.036 | 0.042 | 0.025 | 0.029 |
| 100 | 0.022 | 0.022 | 0.020 | 0.021 | 0.021 | 0.020 | 0.032 | 0.037 | 0.023 | 0.027 |
| 125 | 0.021 | 0.021 | 0.019 | 0.020 | 0.020 | 0.018 | 0.030 | 0.035 | 0.022 | 0.026 |
| 150 | 0.020 | 0.020 | 0.018 | 0.019 | 0.020 | 0.017 | 0.028 | 0.033 | 0.021 | 0.025 |
kable(b2_sd, align = "l", caption = "**Sample Size and Standard Errors for Block 2**",
row.names = F, escape = FALSE, centering = T)| Sample Size | MO1 | MO2 | SC1 | SC2 | UA1 | UA2 | PD1 | PD2 | AD1 | AD2 |
|---|---|---|---|---|---|---|---|---|---|---|
| 20 | 0.143 | 0.222 | 0.068 | 0.151 | 0.089 | 0.178 | 0.442 | 0.515 | 0.079 | 0.188 |
| 30 | 0.087 | 0.066 | 0.040 | 0.058 | 0.048 | 0.079 | 0.163 | 0.194 | 0.044 | 0.069 |
| 40 | 0.071 | 0.063 | 0.035 | 0.059 | 0.040 | 0.079 | 0.144 | 0.159 | 0.038 | 0.067 |
| 50 | 0.052 | 0.043 | 0.030 | 0.040 | 0.033 | 0.054 | 0.095 | 0.108 | 0.032 | 0.049 |
| 75 | 0.042 | 0.035 | 0.026 | 0.034 | 0.028 | 0.044 | 0.074 | 0.089 | 0.026 | 0.040 |
| 100 | 0.032 | 0.030 | 0.023 | 0.027 | 0.024 | 0.033 | 0.053 | 0.065 | 0.024 | 0.035 |
| 125 | 0.027 | 0.026 | 0.021 | 0.024 | 0.022 | 0.029 | 0.043 | 0.053 | 0.022 | 0.031 |
| 150 | 0.024 | 0.025 | 0.020 | 0.022 | 0.022 | 0.026 | 0.038 | 0.048 | 0.022 | 0.029 |
kable(b3_sd, align = "l", caption = "**Sample Size and Standard Errors for Block 3**",
row.names = F, escape = FALSE, centering = T)| Sample Size | MO1 | MO2 | SC1 | SC2 | UA1 | UA2 | PD1 | PD2 | AD1 | AD2 |
|---|---|---|---|---|---|---|---|---|---|---|
| 20 | 0.053 | 0.082 | 0.048 | 0.052 | 0.042 | 0.061 | 0.118 | 0.157 | 0.052 | 0.071 |
| 30 | 0.036 | 0.051 | 0.035 | 0.038 | 0.030 | 0.043 | 0.073 | 0.106 | 0.036 | 0.047 |
| 40 | 0.032 | 0.048 | 0.031 | 0.034 | 0.027 | 0.039 | 0.067 | 0.100 | 0.033 | 0.043 |
| 50 | 0.028 | 0.042 | 0.028 | 0.032 | 0.025 | 0.035 | 0.056 | 0.084 | 0.030 | 0.039 |
| 75 | 0.024 | 0.037 | 0.024 | 0.027 | 0.022 | 0.031 | 0.047 | 0.073 | 0.026 | 0.034 |
| 100 | 0.021 | 0.030 | 0.022 | 0.024 | 0.020 | 0.027 | 0.038 | 0.059 | 0.024 | 0.031 |
| 125 | 0.019 | 0.028 | 0.021 | 0.023 | 0.019 | 0.025 | 0.034 | 0.053 | 0.023 | 0.029 |
| 150 | 0.018 | 0.026 | 0.020 | 0.021 | 0.018 | 0.023 | 0.030 | 0.050 | 0.022 | 0.027 |
kable(b4_sd, align = "l", caption = "**Sample Size and Standard Errors for Block 4**",
row.names = F, escape = FALSE, centering = T)| Sample Size | MO1 | MO2 | SC1 | SC2 | UA1 | UA2 | PD1 | PD2 | AD1 | AD2 |
|---|---|---|---|---|---|---|---|---|---|---|
| 20 | 0.051 | 0.081 | 0.046 | 0.047 | 0.045 | 0.062 | 0.141 | 0.176 | 0.046 | 0.086 |
| 30 | 0.051 | 0.083 | 0.048 | 0.048 | 0.047 | 0.063 | 0.148 | 0.161 | 0.048 | 0.090 |
| 40 | 0.028 | 0.041 | 0.026 | 0.027 | 0.029 | 0.037 | 0.051 | 0.073 | 0.033 | 0.041 |
| 50 | 0.025 | 0.036 | 0.024 | 0.026 | 0.027 | 0.034 | 0.045 | 0.066 | 0.030 | 0.037 |
| 75 | 0.021 | 0.031 | 0.021 | 0.022 | 0.022 | 0.029 | 0.038 | 0.056 | 0.026 | 0.033 |
| 100 | 0.020 | 0.028 | 0.018 | 0.020 | 0.020 | 0.026 | 0.036 | 0.052 | 0.024 | 0.031 |
| 125 | 0.019 | 0.027 | 0.017 | 0.019 | 0.019 | 0.025 | 0.033 | 0.048 | 0.023 | 0.029 |
| 150 | 0.018 | 0.026 | 0.017 | 0.018 | 0.018 | 0.024 | 0.032 | 0.047 | 0.022 | 0.028 |
write.csv(b1_sd, file = "b1_sd_240218.csv", row.names = F)
write.csv(b2_sd, file = "b2_sd_240218.csv", row.names = F)
write.csv(b3_sd, file = "b3_sd_240218.csv", row.names = F)
write.csv(b4_sd, file = "b4_sd_240218.csv", row.names = F)
write.csv(full_sd, file = "full_sd_240218.csv", row.names = F)
write.csv(g2_model[["estimate"]], file = "kaizen_evidence.csv", row.names = F)
save(b1_sd,b2_sd,b3_sd,b4_sd,full_sd,g2_model,model_list,boot231016a,
boot231016b,boot231016c,file = "kaizen_saved.rda")