Nothing
#' Creates a regresion table in APA style with bootstrap confidence intervals
#' @param ... Regression (i.e., lm) result objects. Typically, one for each block in the regression.
#' @param filename (optional) Output filename document filename (must end in .rtf or .doc only)
#' @param table.number Integer to use in table number output line
#' @param number.samples Number of samples to create for bootstrap CIs
#' @return APA table object
#' @references
#' Algina, J. Keselman, H.J. & Penfield, R.J. (2008). Note on a confidence interval for the squared semipartial correlation coefficient. Educational and Psychological Measurement, 68, 734-741.
#' @examples
#' \dontrun{
#' #Note: number.samples = 50 below.
#' # However, please use a value of 1000 or higher
#'
#' # View top few rows of goggles data set
#' # from Discovering Statistics Using R
#' set.seed(1)
#' head(album)
#'
#' # Single block example
#' blk1 <- lm(sales ~ adverts + airplay, data=album)
#' apa.reg.boot.table(blk1)
#' apa.reg.boot.table(blk1,filename="exRegTable.doc")
#'
#' # Two block example, more than two blocks can be used
#' blk1 <- lm(sales ~ adverts, data=album)
#' blk2 <- lm(sales ~ adverts + airplay + attract, data=album)
#' apa.reg.boot.table(blk1,blk2,filename="exRegBlocksTable.doc")
#'
#' # Interaction product-term test with blocks
#' blk1 <- lm(sales ~ adverts + airplay, data=album)
#' blk2 <- lm(sales ~ adverts + airplay + I(adverts * airplay), data=album)
#' apa.reg.boot.table(blk1,blk2,filename="exInteraction1.doc")
#' }
#' @export
apa.reg.boot.table<-function(...,filename=NA, table.number=NA, number.samples = 1000) {
cat("\n\napa.reg.boot.table is a beta version.\n")
prop_var_conf_level <- .95
K <- number.samples
conf.level <- .95
conf_level <- conf.level
prop.var.conf.level = .95
if (prop.var.conf.level == .90) {
prop_var_conf_level <- .90
} else {
prop_var_conf_level <- .95
}
regression_results_list <- list(...)
table_number <- table.number
is_random_predictors <- FALSE
if (is.na(filename)) {
make_file_flag=FALSE
} else {
make_file_flag=TRUE
}
L=length(regression_results_list)
is_same_criterion <- c()
first_result <- regression_results_list[[1]]
first_criterion <- colnames(first_result$model)[1]
for (i in 1:L) {
cur_result <- regression_results_list[[i]]
cur_criterion_name <- colnames(cur_result$model)[1]
is_same_criterion[i] <- first_criterion == cur_criterion_name
}
if (any(is_same_criterion == FALSE)) {
cat("apa.reg.table error:\nAll regression objects (i.e., blocks) must use the same criterion.\n")
cat("The regression objects used had different criterion variables.\n\n")
return(FALSE)
}
is_same_predictors <- c()
first_result <- regression_results_list[[1]]
first_model <- first_result$model
last_model_number_predictors <- dim(first_model)[2]
last_predictors <- colnames(first_result$model)[2:last_model_number_predictors]
n <- dim(first_result$model)[1]
for (i in 1:L) {
cur_result <- regression_results_list[[i]]
cur_model <- cur_result$model
cur_model_number_predictors <- dim(cur_model)[2]
cur_predictors <- colnames(cur_model)[2:cur_model_number_predictors]
is_same_predictors[i] <- all(intersect(last_predictors,cur_predictors) == last_predictors)
last_predictors <- cur_predictors
}
if (any(is_same_predictors==FALSE)) {
cat("apa.reg.table error:\nEach regression objects (i.e., block) must contain all of the predictors from the preceeding regression object (i.e., block).\n\n")
cat("For example:\n")
cat("block1 <- lm(y ~ a + b)\n")
cat("block2 <- lm(y ~ a + b + c)\n\n")
cat("The second block contains all of the predictors from the first block plus additional predictors.\n\n")
cat("Therefore the command below will work: \n\n")
cat("apa.reg.table(block1, block2)\n\n")
return(FALSE)
}
#get analyses for each block
block_results <- list()
L <- length(regression_results_list)
for (i in 1:L) {
cat(sprintf("Block %g: Generating %g bootstrap samples\n", i, K))
cur_result <- apa_single_boot_block(regression_results_list[[i]],is_random_predictors, K = K, prop_var_conf_level = prop_var_conf_level)
block_results[[i]] <- cur_result
}
cat(sprintf("Bootstrap for Delta RSQ in progress\n"))
R2_boot_results<- boot::boot(data = regression_results_list[[L]]$model, statistic = boot_blocks, R=K, regression_results_list=regression_results_list)
is_multiple_blocks <- FALSE
if (L>1) {
is_multiple_blocks <- TRUE
}
if (is_multiple_blocks == TRUE) {
#contain (for each boot dataset/row) R2 for each block
# must do some subtraction b2-b1 (b3-b2) etc... to get Delta's
delta_R2_data <- calculate_boot_delta_R2(R2_boot_results$t)
delta_R2_result <- R2_boot_results$t0
delta_R2_CI <- get_delta_R2_CI(delta_R2_data, prop_var_conf_level)
}
#Combine blocks
block_out_txt <- block_results[[1]]$model_details_txt
block_out_rtf <- block_results[[1]]$model_details_rtf
first_block_calculate_cor <- block_results[[1]]$calculate_cor # use later for table note
first_block_calculate_beta <- block_results[[1]]$calculate_beta # use later for table note
last_block_summary <- block_results[[1]]$model_summary_extended
last_block_lm <- regression_results_list[[1]]
if (is_multiple_blocks == TRUE) {
for (i in 2:L) {
cur_block_lm <- regression_results_list[[i]]
cur_block_summary <- block_results[[i]]$model_summary_extended
cur_block_out_txt <- block_results[[i]]$model_details_txt
cur_block_out_rtf <- block_results[[i]]$model_details_rtf
delta_R2_details <- get_delta_R2_blocks(blk2=cur_block_lm,blk1=last_block_lm,summary2=cur_block_summary,summary1=last_block_summary,n, prop_var_conf_level = prop_var_conf_level)
delta_R2_details <- update_with_delta_R2_boot_CI(delta_R2_details, delta_R2_CI[i,], prop_var_conf_level = prop_var_conf_level)
num_lines <- dim(cur_block_out_txt)[1]
cur_block_out_txt$difference[num_lines-2] <- delta_R2_details$deltaR2_txt
cur_block_out_txt$difference[num_lines-1] <- delta_R2_details$deltaR2_CI_txt
cur_block_out_rtf$difference[num_lines-2] <- delta_R2_details$deltaR2_rtf
cur_block_out_rtf$difference[num_lines-1] <- delta_R2_details$deltaR2_CI_rtf
last_block_summary <- cur_block_summary
last_block_lm <- cur_block_lm
block_out_txt <- rbind(block_out_txt,cur_block_out_txt)
block_out_rtf <- rbind(block_out_rtf,cur_block_out_rtf)
}
} else {
block_out_txt <- dplyr::select(block_out_txt, -difference)
block_out_rtf <- dplyr::select(block_out_rtf, -difference)
}
#console table
table_title <- sprintf("Regression results using %s as the criterion\n",first_criterion)
txt_column_names <- get_txt_column_names(block_out_txt)
if (prop_var_conf_level == .95) {
txt_column_names <- sub("XX","95",txt_column_names)
} else {
txt_column_names <- sub("XX","90",txt_column_names)
}
names(block_out_txt) <- txt_column_names
table_body <- block_out_txt
table_note <- get_reg_table_note_txt(first_block_calculate_cor, first_block_calculate_beta)
table_block_results <- block_results
tbl_console <- list(table_number = table_number,
table_title = table_title,
table_body = as.data.frame(table_body),
table_note = table_note,
table_block_results = table_block_results)
table_block_results <- block_results
class(tbl_console) <- "apa_table"
if (make_file_flag==TRUE) {
table_title <- sprintf("Regression results using %s as the criterion\n",first_criterion)
table_note <- get_reg_table_note_rtf(first_block_calculate_cor, first_block_calculate_beta)
#set columns widths and names
colwidths <- get_rtf_column_widths(block_out_rtf)
regression_table <- as.matrix(block_out_rtf)
rtf_column_names <- get_rtf_column_names(block_out_rtf)
if (prop_var_conf_level == .95) {
rtf_column_names <- sub("XX","95",rtf_column_names)
} else {
rtf_column_names <- sub("XX","90",rtf_column_names)
}
new_col_names <- rtf_column_names
colnames(regression_table) <- new_col_names
#Create RTF code
rtfTable <- RtfTable$new(isHeaderRow=TRUE, defaultDecimalTableProportionInternal=.15)
rtfTable$setTableContent(regression_table)
rtfTable$setCellWidthsInches(colwidths)
rtfTable$setRowSecondColumnDecimalTab(.4)
txt_body <- rtfTable$getTableAsRTF(FALSE,FALSE)
if (is_multiple_blocks==TRUE) {
write.rtf.table(filename = filename,txt.body = txt_body,table.title = table_title, table.note = table_note,landscape=TRUE,table.number=table_number)
} else {
write.rtf.table(filename = filename,txt.body = txt_body,table.title = table_title, table.note = table_note,table.number=table_number,landscape=TRUE)
}
}
return(tbl_console)
}
apa_single_boot_block<-function(cur_blk,is_random_predictors, K, prop_var_conf_level) {
#Regression data
reg_table_data <- cur_blk$model
n <- dim(reg_table_data)[1]
col_names <- names(reg_table_data)
#Determine if beta weights should be calculated
is_weighted <- "(weights)" %in% col_names
is_var_factor <- is_variable_factor(reg_table_data)
calculate_beta <- TRUE
calculate_cor <- TRUE
if (is_weighted == TRUE) {
calculate_beta <- FALSE
}
if (is_var_factor == TRUE) {
calculate_beta <- FALSE
calculate_cor <- FALSE
}
#Summary statistics
model_summary_extended <- broom::glance(cur_blk) #glance include lm constant in df
model_summary_extended$df <- summary(cur_blk)$df[1] #use summary information to get df without constant
#Regression table statistics
reg_table <- broom::tidy(cur_blk)
names(reg_table) <- c("predictor","b","SE","t","p")
#adjust df
is_intercept <- FALSE
if (reg_table$predictor[1] == "(Intercept)") {
model_summary_extended$df[1] <- model_summary_extended$df[1]-1
is_intercept <- TRUE
}
R2 <- model_summary_extended$r.squared
R2_pvalue <- model_summary_extended$p.value
df1 <- model_summary_extended$df
df2 <- model_summary_extended$df.residual
# Get bootstrap samples
boot_results <- boot::boot(data = reg_table_data, statistic = boot_equation, R = K, lm_out=cur_blk)
boot_data <- get_boot_data(boot_results = boot_results)
prop_var_conf_level_str <- sprintf("%g", round(prop_var_conf_level*100))
R2CI_boot <- get_boot_R2_CI(boot_data, conf_level = prop_var_conf_level)
R2LL <- R2CI_boot[1]
R2UL <- R2CI_boot[2]
R2_txt <- strip.leading.zero(add.sig.stars(sprintf("%1.3f",R2),R2_pvalue))
if (R2_pvalue>= .05) {
R2LL <- 0 # Algina, Keselman, & Penfield (2008) Correction
cat(" R2 CI lower bound adjusted to zero as per Algina, Keselman, & Penfield (2008)\n\n")
}
R2LL_txt <- strip.leading.zero(sprintf("%1.2f",R2LL))
R2UL_txt <- strip.leading.zero(sprintf("%1.2f",R2UL))
model_summary_txt <- sprintf("R2 = %s",R2_txt)
model_summary_CI_txt <- sprintf("%s%% CI[%s,%s]",prop_var_conf_level_str, R2LL_txt,R2UL_txt)
model_summary_rtf <- sprintf("{\\i R\\super 2 \\nosupersub} = %s",R2_txt)
model_summary_CI_rtf <- sprintf("%s%% CI[%s,%s]",prop_var_conf_level_str, R2LL_txt,R2UL_txt)
#Add b-weight CI's
b_CI_boot <- get_boot_b_CI(boot_data = boot_data, cur_blk = cur_blk, conf_level = .95)
if (!is.null(dim(b_CI_boot))){
LLb <- b_CI_boot[,c("LL")]
ULb <- b_CI_boot[,c("UL")]
} else {
LLb <- b_CI_boot[1]
ULb <- b_CI_boot[2]
}
# LLb <- b_CI_boot[,c("LL")]
# ULb <- b_CI_boot[,c("UL")]
# b_CI <- confint(cur_blk)
# LLb <- b_CI[,c("2.5 %")]
# ULb <- b_CI[,c("97.5 %")]
reg_table <- dplyr::mutate(reg_table, LLb=LLb, ULb=ULb)
reg_table <- dplyr::select(reg_table,predictor, b,LLb,ULb,SE,t,p)
#table parts
reg_table_components <- table_without_intercept_row(reg_table)
reg_table_lower <- reg_table_components$lower_table
reg_table_first <- reg_table_components$first_row #intercept row
#correlation
if (calculate_cor==TRUE) {
r_with_criterion <- correlations_with_criterion(reg_table_data)
reg_table_lower <- dplyr::inner_join(reg_table_lower,r_with_criterion,by="predictor")
}
#semi-partial correlation squared
reg_table_lower <- dplyr::mutate(reg_table_lower, sr2=(t*t)*(1-R2)/df2, LLsr2=-999,ULsr2=-999)
number_predictors <- dim(reg_table_lower[,1])[1]
# if (number_predictors>1) {
# #use delta R2 process for CI
# for (i in 1:number_predictors) {
# sr2 <- reg_table_lower$sr2[i]
# ci <- get_sr2_ci(sr2=sr2,R2=R2,n=n)
# LL <- ci$LL
# UL <- ci$UL
# reg_table_lower$LLsr2[i] <- LL
# reg_table_lower$ULsr2[i] <- UL
# }
# } else {
# #single predictor so use R2 CI
# reg_table_lower$LLsr2[1] <- R2LL
# reg_table_lower$ULsr2[1] <- R2UL
# }
sr2_CI_boot <- get_boot_sr2_CI(boot_data, cur_blk, prop_var_conf_level)
sr2_CI_boot <- trim_intercept_row(sr2_CI_boot)
if (!is.null(dim(sr2_CI_boot))){
reg_table_lower$LLsr2 <- sr2_CI_boot[,c("LL")]
reg_table_lower$ULsr2 <- sr2_CI_boot[,c("UL")]
} else {
reg_table_lower$LLsr2 <- sr2_CI_boot[1]
reg_table_lower$ULsr2 <- sr2_CI_boot[2]
}
id_is_ns <- (reg_table_lower$p >= .05)
if (sum(id_is_ns) >0) {
reg_table_lower$LLsr2[id_is_ns] <- 0
cat(" sr2 CI lower bound(s) adjusted to zero as per Algina, Keselman, & Penfield (2008)\n\n")
}
#beta
if (calculate_beta==TRUE) {
reg_table_lower <- dplyr::mutate(reg_table_lower, beta = - 999, LLbeta = -999, ULbeta = -999)
for (i in 1:number_predictors) {
p_name <- reg_table_lower$predictor[i]
b <- reg_table_lower$b[i]
#LLb <- reg_table_lower$LLb[i]
#ULb <- reg_table_lower$ULb[i]
sd_crit <- stats::sd(reg_table_data[,1],na.rm = TRUE)
sd_pred <- stats::sd(reg_table_data[,p_name],na.rm = TRUE)
beta <- convert_b_to_beta(b = b , sd_pred = sd_pred, sd_crit = sd_crit)
# LLbeta <- convert_b_to_beta(b = LLb, sd_pred = sd_pred, sd_crit = sd_crit)
# ULbeta <- convert_b_to_beta(b = ULb, sd_pred = sd_pred, sd_crit = sd_crit)
# beta <- b*(sd_pred / sd_crit)
# LLbeta <- LLb*(sd_pred / sd_crit)
# ULbeta <- ULb*(sd_pred / sd_crit)
reg_table_lower$beta[i] <- beta
# reg_table_lower$LLbeta[i] <- LLbeta
# reg_table_lower$ULbeta[i] <- ULbeta
}
beta_CI_boot <- get_boot_beta_CI(boot_data, cur_blk, .95)
beta_CI_boot <- trim_intercept_row(beta_CI_boot)
if (!is.null(dim(beta_CI_boot))){
reg_table_lower$LLbeta <- beta_CI_boot[,c("LL")]
reg_table_lower$ULbeta <- beta_CI_boot[,c("UL")]
} else {
reg_table_lower$LLbeta <- beta_CI_boot[1]
reg_table_lower$ULbeta <- beta_CI_boot[2]
}
}
#Combine with intercept
if (reg_table_components$is_intercept == TRUE) {
intercept_row <- get_empty_row(reg_table_lower)
intercept_row$predictor[1] <- reg_table_first$predictor[1]
intercept_row$b[1] <- reg_table_first$b[1]
intercept_row$LLb[1] <- reg_table_first$LLb[1]
intercept_row$ULb[1] <- reg_table_first$ULb[1]
intercept_row$t[1] <- reg_table_first$t[1]
intercept_row$p[1] <- reg_table_first$p[1]
intercept_row$SE[1] <- reg_table_first$SE[1]
model_details_extended <- rbind(intercept_row, reg_table_lower)
} else {
model_details_extended <- reg_table_lower
}
L <- dim(model_details_extended)[1]
CIb_str <- txt.ci.brackets(model_details_extended$LLb, model_details_extended$ULb, strip_zero = FALSE)
CIsr2_str <- txt.ci.brackets(model_details_extended$LLsr2, model_details_extended$ULsr2, strip_zero = TRUE)
if (calculate_beta == TRUE) {
CIbeta_str <- txt.ci.brackets(model_details_extended$LLbeta, model_details_extended$ULbeta, strip_zero = FALSE)
}
model_details <- model_details_extended[,1,drop=FALSE]
model_details$b <- add.sig.stars(sprintf("%1.2f",model_details_extended$b), model_details_extended$p)
model_details$b_CI <- CIb_str
if (calculate_beta==TRUE) {
model_details$beta <- sprintf("%1.2f",model_details_extended$beta)
model_details$beta_CI <- CIbeta_str
if (is_intercept == TRUE) {
model_details$beta[1] <- ""
model_details$beta_CI[1] <- ""
}
}
model_details$sr2 <- strip.leading.zero(sprintf("%1.2f", model_details_extended$sr2))
model_details$sr2_CI <- CIsr2_str
if (is_intercept == TRUE) {
model_details$sr2[1] <- ""
model_details$sr2_CI[1] <- ""
}
if (calculate_cor==TRUE) {
model_details$r <- strip.leading.zero(add.sig.stars(sprintf("%1.2f",model_details_extended$r), model_details_extended$r_pvalue)) #intercept row issue
product_rows <- is_product_row(model_details$predictor)
if (is_intercept == TRUE) {
model_details$r[1] <- ""
model_details$r[product_rows] <- ""
}
}
model_details_txt <- model_details
model_details_txt$summary <- ""
model_details_txt$difference <- "" #blank column to be populated later
model_details_txt <- add_row_to_model_summary(model_details_txt)
model_details_txt <- add_row_to_model_summary(model_details_txt)
model_details_txt <- add_row_to_model_summary(model_details_txt)
num_row <- dim(model_details_txt)[1]
model_details_rtf <- model_details_txt
model_details_txt$summary[num_row-2] <- model_summary_txt
model_details_txt$summary[num_row-1] <- model_summary_CI_txt
model_details_rtf$summary[num_row-2] <- model_summary_rtf
model_details_rtf$summary[num_row-1] <- model_summary_CI_rtf
output <- list()
output$model_summary_extended <- model_summary_extended
output$model_details_extended <- model_details_extended
output$model_details_txt <- model_details_txt
output$model_details_rtf <- model_details_rtf
output$calculate_beta <- calculate_beta
output$calculate_cor <- calculate_cor
return(output)
}
boot_equation <- function(data, indices, lm_out){
data2 <- data[indices, ]
#data2 <- attitude
cur_lm <- lm(formula = lm_out$call$formula, data=data2)
lm_data <- tibble::as_tibble(cur_lm$model)
cur_results <- summary(cur_lm)
cur_coef <- cur_results$coefficients
t_values <- cur_coef[,colnames(cur_coef)=="t value"]
b_values <- cur_coef[,colnames(cur_coef)=="Estimate"]
names(b_values) <- rownames(cur_coef)
names(t_values) <- rownames(cur_coef)
df2 <- cur_results$df[2]
R2 <- cur_results$r.squared
other_stats <- c(df2, R2)
names(other_stats) <- c("df2","R2")
#Get beta
sd_model <- apply(lm_data, FUN=sd, 2)
sd_criterion <- sd_model[1]
sd_predictors <- sd_model[2:length(sd_model)]
b_to_beta_multiplier <- sd_predictors/sd_criterion
colnames_predictors <- names(b_values)
is_intercept <- FALSE # Fix in original XYZZY
if (!is.null(colnames_predictors)) {
if (colnames_predictors[1]=="(Intercept)") {
is_intercept <- TRUE
}
}
if (is_intercept == TRUE) {
b_values_to_convert <- b_values[2:length(b_values)]
initial_names <- names(b_values_to_convert)
beta_values <- c(-999, b_values_to_convert*b_to_beta_multiplier)
names(beta_values) <- c("(Intercept)", initial_names)
} else {
b_values_to_convert <- b_values
beta_values <- b_values_to_convert*b_to_beta_multiplier
}
names(t_values) <- paste("t_",names(t_values), sep="")
names(b_values) <- paste("b_",names(b_values), sep="")
names(beta_values) <- paste("beta_",names(beta_values), sep="")
output <- c(t_values, b_values, beta_values, other_stats)
return(output)
}
t_df_to_sr2_df <- function(t_data, other_data) {
df2 <- other_data$df2
R2 <- other_data$R2
n_col <- dim(t_data)[2]
t_data_names <- names(t_data)
sr2_data_names <- gsub("t_","sr2_", t_data_names)
t_to_sr2_multiplier <- (1-R2)/df2
t_to_sr2_multiplier_matrix <- t_to_sr2_multiplier %o% rep(1,n_col)
t_matrix <- as.matrix(t_data)
sr2_matrix <- (t_matrix^2) * t_to_sr2_multiplier_matrix
sr2_data <- as.data.frame(sr2_matrix)
names(sr2_data) <- sr2_data_names
return(sr2_data)
}
get_boot_data <- function(boot_results) {
boot_data <- boot_results$t
boot_column_names <- names(boot_results$t0)
boot_data <- tibble::as_tibble(boot_data)
names(boot_data) <- boot_column_names
return(boot_data)
}
get_percentile_conf <- function(x, conf_level) {
K <- length(x)
lower_index <- round( ((1-conf_level)/2) * K) + 1
upper_index <- round( ((1-(1-conf_level)/2)) * K)
x_sorted <- sort(x)
LL <- x_sorted[lower_index]
UL <- x_sorted[upper_index]
return(c(LL,UL))
}
get_boot_R2_CI <- function(boot_data, conf_level) {
x <- boot_data$R2
ci_values <- get_percentile_conf(x, conf_level=conf_level)
return(ci_values)
}
get_boot_b_CI <- function(boot_data, cur_blk, conf_level) {
predictor_names <- rownames(summary(cur_blk)$coefficients)
b_names <- paste("b_", predictor_names, sep="")
b_data <- boot_data[, b_names]
matrix_out <- ci_matrix_from_data(cur_data = b_data, predictor_names = predictor_names, conf_level = conf_level)
return(matrix_out)
}
ci_matrix_from_data <- function(cur_data, predictor_names, conf_level) {
num_predictors <- length(predictor_names)
# if (!is.null(dim(cur_data))) {
# num_weights <- dim(cur_data)[2]
# } else {
# num_weights <- 1
# }
matrix_out <- matrix(NA, nrow=num_predictors, ncol=2)
rownames(matrix_out) <- predictor_names
colnames(matrix_out) <- c("LL","UL")
for (i in 1:num_predictors) {
cur_col <- unlist(cur_data[,i])
cur_ci <- get_percentile_conf(x = cur_col, conf_level=conf_level)
matrix_out[i, 1] <- cur_ci[1]
matrix_out[i, 2] <- cur_ci[2]
}
return(matrix_out)
}
get_boot_beta_CI <- function(boot_data, cur_blk, conf_level) {
lm_summary <- summary(cur_blk)
predictor_names <- rownames(lm_summary$coefficients)
beta_names <- paste("beta_", predictor_names, sep="")
beta_data <- boot_data[, beta_names]
matrix_out <- ci_matrix_from_data(cur_data = beta_data, predictor_names = predictor_names, conf_level = conf_level)
return(matrix_out)
}
get_boot_sr2_CI <- function(boot_data, cur_blk, conf_level) {
lm_summary <- summary(cur_blk)
predictor_names <- rownames(lm_summary$coefficients)
t_names <- paste("t_", predictor_names, sep="")
t_data <- boot_data[, t_names]
sr2_data <- t_df_to_sr2_df(t_data, boot_data)
matrix_out <- ci_matrix_from_data(cur_data = sr2_data, predictor_names = predictor_names, conf_level = conf_level)
return(matrix_out)
}
trim_intercept_row <- function(matrix_in) {
cur_rownames <- rownames(matrix_in)
id <- cur_rownames != "(Intercept)"
matrix_out <- matrix_in[id,]
return(matrix_out)
}
boot_blocks <- function(data, indices, regression_results_list){
data2 <- data[indices, ] # select obs. in bootstrap sample
L <- length(regression_results_list)
R2_vec <- rep(NA,L)
for (i in 1:L) {
orig_reg_result <- regression_results_list[[i]]
cur_reg_result <- lm(formula = orig_reg_result$call$formula, data=data2)
R2_vec[i] <- summary(cur_reg_result)$r.squared
}
return(R2_vec)
}
calculate_boot_delta_R2 <- function(matrix_in){
matrix_dims <- dim(matrix_in)
number_blocks <- matrix_dims[2]
K <- matrix_dims[1]
delta_R2_data <- matrix(NA, nrow=K, ncol=number_blocks)
for (i in 2:number_blocks) {
cur_R2 <- matrix_in[,i]
prev_R2 <- matrix_in[,(i-1)]
delta_value <- cur_R2 - prev_R2
delta_R2_data[,i] <- delta_value
}
return(delta_R2_data)
}
get_delta_R2_CI <- function(delta_R2_data, conf_level) {
number_blocks <- dim(delta_R2_data)[2]
block_names <- paste("block",1:number_blocks,sep="")
colnames(delta_R2_data) <- block_names
matrix_out <- ci_matrix_from_data(cur_data = delta_R2_data, predictor_names = block_names, conf_level = conf_level)
return(matrix_out)
}
update_with_delta_R2_boot_CI <- function(delta_R2_details, delta_R2_CI, prop_var_conf_level) {
pvalue <- delta_R2_details$deltaR2_pvalue
deltaR2 <- delta_R2_details$deltaR2
is_sig <- FALSE
if (pvalue<.05) {
is_sig <- TRUE
}
LL <- delta_R2_CI[1]
UL <- delta_R2_CI[2]
#correction from literature
if (is_sig == FALSE) {
LL <- 0
cat(" Delta R2 CI lower bound adjusted to zero as per Algina, Keselman, & Penfield (2008)\n\n")
}
prop_var_conf_level_str <- sprintf("%g", round(prop_var_conf_level*100))
delta_R2_details$deltaR2_CI_txt<- sprintf("%s%% CI[%s, %s]", prop_var_conf_level_str, strip.leading.zero(sprintf("%1.2f",LL)), strip.leading.zero(sprintf("%1.2f",UL)))
delta_R2_details$deltaR2_CI_rtf<- sprintf("{%s%% CI}[%s, %s]",prop_var_conf_level_str, strip.leading.zero(sprintf("%1.2f",LL)), strip.leading.zero(sprintf("%1.2f",UL)))
return(delta_R2_details)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.