diff --git a/.Rbuildignore b/.Rbuildignore index 2ac6917..54148bd 100644 --- a/.Rbuildignore +++ b/.Rbuildignore @@ -14,7 +14,7 @@ ^build/js$ ^temp$ ^.*\.jmo$ -^Validaion$ +^Validation$ ^NEWS\.Rmd$ ^jmvpower$ ^ANCOVA$ diff --git a/.github/workflows/R-CMD-check.yaml b/.github/workflows/R-CMD-check.yaml index b3483f0..7e35d1a 100644 --- a/.github/workflows/R-CMD-check.yaml +++ b/.github/workflows/R-CMD-check.yaml @@ -43,4 +43,5 @@ jobs: with: extra-packages: rcmdcheck - - uses: r-lib/actions/check-r-package@v1 \ No newline at end of file + - uses: r-lib/actions/check-r-package@v1 + diff --git a/.gitignore b/.gitignore index 7bbef1e..c751c58 100644 --- a/.gitignore +++ b/.gitignore @@ -11,3 +11,5 @@ temp inst/doc .DS_Store +NUL +tests/testthat/Rplots.pdf diff --git a/DESCRIPTION b/DESCRIPTION index 54207bc..bc54190 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: Superpower Title: Simulation-Based Power Analysis for Factorial Designs -Version: 0.1.2 +Version: 0.2.0 Authors@R: c(person(given = "Aaron", family = "Caldwell", role = c("aut", "cre"), @@ -26,7 +26,7 @@ URL: https://aaroncaldwell.us/SuperpowerBook/ BugReports: https://github.com/arcaldwell49/Superpower/issues License: MIT + file LICENSE Encoding: UTF-8 -RoxygenNote: 7.1.1 +RoxygenNote: 7.1.2 Imports: mvtnorm, MASS, @@ -39,7 +39,6 @@ Imports: dplyr, magrittr, tidyselect, - Hmisc, tidyr Suggests: knitr, @@ -47,6 +46,8 @@ Suggests: pwr, testthat, covr, - jmvcore + jmvcore, + spelling VignetteBuilder: knitr +Language: en-US diff --git a/NAMESPACE b/NAMESPACE index 5837cfc..67d2fcb 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -4,12 +4,16 @@ S3method(confint,sim_result) S3method(emmeans_power,data.frame) S3method(emmeans_power,emmGrid) S3method(emmeans_power,summary_em) +S3method(plot,ancova_power) S3method(plot,design_aov) S3method(plot,opt_alpha) S3method(plot,sim_result) +S3method(print,ancova_power) S3method(print,design_aov) S3method(print,opt_alpha) S3method(print,sim_result) +export(ANCOVA_analytic) +export(ANCOVA_contrast) export(ANOVA_compromise) export(ANOVA_design) export(ANOVA_exact) @@ -25,12 +29,12 @@ export(optimal_alpha) export(p_standardized) export(plot_power) export(power.ftest) +export(power_oneway_ancova) export(power_oneway_between) export(power_oneway_within) export(power_standardized_alpha) export(power_threeway_between) export(power_twoway_between) -import(Hmisc) import(emmeans) import(ggplot2) importFrom(MASS,mvrnorm) @@ -50,7 +54,14 @@ importFrom(magrittr,"%>%") importFrom(magrittr,'%>%') importFrom(reshape2,melt) importFrom(stats,as.formula) +importFrom(stats,contr.sum) +importFrom(stats,contrasts) +importFrom(stats,dbeta) +importFrom(stats,df) +importFrom(stats,dt) importFrom(stats,median) +importFrom(stats,model.matrix) +importFrom(stats,optim) importFrom(stats,optimize) importFrom(stats,p.adjust) importFrom(stats,pf) @@ -61,7 +72,9 @@ importFrom(stats,pt) importFrom(stats,qf) importFrom(stats,qnorm) importFrom(stats,qt) +importFrom(stats,qtukey) importFrom(stats,sd) +importFrom(stats,terms) importFrom(stats,uniroot) importFrom(tidyr,expand_grid) importFrom(tidyselect,matches) diff --git a/R/ANCOVA_analytic.R b/R/ANCOVA_analytic.R new file mode 100644 index 0000000..0eb2606 --- /dev/null +++ b/R/ANCOVA_analytic.R @@ -0,0 +1,561 @@ +#' Power Calculations for Factorial ANCOVAs +#' +#' Complete power analyses for ANCOVA omnibus tests and contrasts. This function does not support within subjects factors. +#' +#' @param design Output from the ANOVA_design function +#' @param mu Vector specifying mean for each condition +#' @param n Sample size in each condition +#' @param sd Standard deviation for all conditions (or a vector specifying the sd for each condition) +#' @param r2 Coefficient of Determination of the model with only the covariates +#' @param n_cov Number of covariates +#' @param alpha_level Alpha level used to determine statistical significance +#' @param beta_level Type II error probability (power/100-1) +#' @param cmats List of matrices for specific contrasts of interest +#' @param label_list An optional list to specify the factor names and condition (recommended, if not used factors and levels are indicated by letters and numbers). +#' @param design_result Output from the ANOVA_design function +#' @param round_up Logical indicator (default = TRUE) for whether to round up sample size calculations to nearest whole number +#' +#' @return +#' One, or two, data frames containing the power analysis results from the power analysis for the omnibus ANCOVA (main_results) or contrast tests (contrast_results). +#' In addition, every F-test (aov_list and con_list) is included in a list of power.htest results. +#' Lastly, a (design_param) list containing the design parameters is also included in the results. +#' @examples +#' # Simple 2x3 ANCOVA +#' +#' ANCOVA_analytic( +#' design = "2b*3b", +#' mu = c(400, 450, 500, +#' 400, 500, 600), +#' n_cov = 3, +#' sd = 100, +#' r2 = .25, +#' alpha_level = .05, +#' beta_level = .2, +#' round_up = TRUE +#' ) +#' @section References: +#' Shieh, G. (2020). Power analysis and sample size planning in ANCOVA designs. Psychometrika, 85(1), 101-120. +#' @importFrom stats uniroot pf df qf contr.sum dt dbeta qtukey model.matrix terms optim contrasts +#' @export +#' + + +ANCOVA_analytic <- function(design, + mu, + n = NULL, + sd, + r2 = NULL, + n_cov, + alpha_level = Superpower_options("alpha_level"), + beta_level = NULL, + cmats = list(), + label_list = NULL, + design_result = NULL, + round_up = TRUE) { + + METHOD <- "Power Calculation for ANCOVA" + METHOD2 <- "Power Calculation for ANCOVA contrast" + TYPE = "Exact" + if(!is.null(design_result)){ + mu = design_result$mu + sd = design_result$sd + n = design_result$n + + + factornames = design_result$factornames + frml3 = as.formula(paste(gsub("\\+", "*",design_result$frml2),collapse=" ")) + labelnameslist = design_result$labelnameslist + design = design_result$design + } else{ + + #Check String for an acceptable digits and factor (w or b) + if (grepl("^(\\d{1,3}(w|b)\\*){0,2}\\d{1,3}(w|b)$", design, ignore.case = FALSE, perl = TRUE) == FALSE) { + stop("Problem in the design argument: must input number of levels as integer (2-999) and factor-type (between or within) as lower case b (between) or w (within)") + } + + #Ensure sd is greater than 0 + if (any(sd <= 0) || length(sd) != 1) { + stop("Standard deviation (sd) is less than or equal to zero; input a single value greater than zero") + } + + #Ensure, if single correlation is input, that it is between 0 and 1 + if(!is.null(r2)){ + if ((r2 <= 0) | r2 >= 1 ) { + stop("Coefficient of Determination must be greater than -1 and less than 1") + } + } + + + factor_levels <- as.numeric(strsplit(design, "\\D+")[[1]]) + if(prod(factor_levels) != length(mu)){ + stop("Length of means does not match design.") + } + labelnames = NULL + if (is.null(label_list)) { + label_list = list() + for (i1 in 1:length(factor_levels)){ + label_list1 = NULL + labelnames <- append(labelnames,paste(paste(letters[i1]), sep = "")) + + for (i2 in 1:factor_levels[i1]){ + labelnames <- append(labelnames,paste(paste(letters[i1]), paste(i2), sep = "")) + label_list1 = c(label_list1,paste(paste(letters[i1]), paste(i2), sep = "")) + } + + label_list[[i1]] = as.vector(label_list1) + names(label_list)[i1] = paste(paste(letters[i1]), sep = "") + + } + + } else{ + for(i in 1:length(label_list)){ + + labelnames = append(labelnames, names(label_list)[i]) + labelnames = append(labelnames, label_list[[i]]) + } + + } + + labelnameslist = label_list + factornames = names(label_list) + + } + + if( missing(mu) || missing(sd) || missing(n_cov)){ + stop("mu, sd, and n_cov are missing and must be provided") + } + + n_grp = length(mu) + + if(!is.null(n)){ + if(length(n)==1){ + nvec = rep(n,length(mu)) + } else { + nvec = n + } + + if(length(nvec) != length(mu)){ + stop("Length of n and mu do not match") + } + } else { + nvec = NULL + } + factorlist = list() + cons = cmats + if(length(factornames) == 1){ + con_df = data.frame(a = 1:length(labelnameslist[[1]])) + con_df$a = as.factor(con_df$a) + contrasts(con_df$a) = "contr.sum" + con_mat = model.matrix(~a,con_df) + colnames(con_mat) = attr(con_mat, "assign") + cmat_a = t(as.matrix(con_mat[,which(colnames(con_mat)==1)])) + cmats = list(a = cmat_a) + factorlist = factornames + } else if(length(factornames) == 2){ + con_df = as.data.frame(expand.grid(b = 1:length(labelnameslist[[2]]), + a = 1:length(labelnameslist[[1]]))) + con_df$a = as.factor(con_df$a) + con_df$b = as.factor(con_df$b) + contrasts(con_df$a) = "contr.sum" + contrasts(con_df$b) = "contr.sum" + #colnames(con_df) = design_result$factornames + con_mat = model.matrix(~a*b,con_df) + colnames(con_mat) = attr(con_mat, "assign") + cmat_a = t(as.matrix(con_mat[,which(colnames(con_mat)==1)])) + cmat_b = t(as.matrix(con_mat[,which(colnames(con_mat)==2)])) + cmat_ab = t(as.matrix(con_mat[,which(colnames(con_mat)==3)])) + cmats = list(a = cmat_a, b = cmat_b, ab=cmat_ab) + factorlist = c(factornames[1], + factornames[2], + paste0(factornames[1], ":", + factornames[2])) + } else { + con_df = as.data.frame(expand.grid(c = 1:length(labelnameslist[[3]]), + b = 1:length(labelnameslist[[2]]), + a = 1:length(labelnameslist[[1]]))) + con_df$a = as.factor(con_df$a) + con_df$b = as.factor(con_df$b) + con_df$c = as.factor(con_df$c) + contrasts(con_df$a) = "contr.sum" + contrasts(con_df$b) = "contr.sum" + contrasts(con_df$c) = "contr.sum" + #colnames(con_df) = design_result$factornames + con_mat = model.matrix(~a*b*c,con_df) + colnames(con_mat) = attr(con_mat, "assign") + cmat_a = t(as.matrix(con_mat[,which(colnames(con_mat)==1)])) + cmat_b = t(as.matrix(con_mat[,which(colnames(con_mat)==2)])) + cmat_c = t(as.matrix(con_mat[,which(colnames(con_mat)==3)])) + cmat_ab = t(as.matrix(con_mat[,which(colnames(con_mat)==4)])) + cmat_ac = t(as.matrix(con_mat[,which(colnames(con_mat)==5)])) + cmat_bc = t(as.matrix(con_mat[,which(colnames(con_mat)==6)])) + cmat_abc = t(as.matrix(con_mat[,which(colnames(con_mat)==7)])) + cmats = list(a = cmat_a, b = cmat_b, c = cmat_c, + ab = cmat_ab, ac = cmat_ac, bc = cmat_bc, + abc = cmat_abc) + factorlist = c(factornames[1], + factornames[2], + factornames[3], + paste0(factornames[1], ":", + factornames[2]), + paste0(factornames[1], ":", + factornames[3]), + paste0(factornames[2], ":", + factornames[3]), + paste0(factornames[1], ":", + factornames[2], ":", + factornames[3])) + } + #cmat = t(contr.sum(length(nvec))) # Contrast matrix + # Create matrix for Wald statistic (W star; eq 10 from Shieh) + pow_res = list() + pow_res2 = list() + con_res = list() + con_res2 = list() + if(grepl("w", design)){ + stop("Design contains a within subject factor. \n This is not supported by this function at this time") + } + + if (!is.null(beta_level)){ + pow = 1 - beta_level + } + + p.body = quote({ + pow_anc_meth( + cmat, + mu, + nvec, + n_cov, + r2, + sd, + alpha_level + )}) + + if (is.null(beta_level)){ + for(i1 in 1:length(cmats)){ + cmat = cmats[[i1]] + res <- eval(p.body) + beta_level <- 1-res$pow + pow_res[[i1]] = list( + cmat = cmat, + mu = mu, + nvec = nvec, + n_cov = n_cov, + r2 = res$r2, + sd = sd, + alpha_level = res$alpha_level, + pow = res$pow, + beta_level = beta_level, + num_df = res$num_df, + den_df = res$den_df, + N_tot = res$N_tot + ) + } + names(pow_res) = names(cmats) + + if(length(cons) != 0){ + for(i1 in 1:length(cons)){ + cmat = cons[[i1]] + res <- eval(p.body) + beta_level <- 1-res$pow + con_res[[i1]] = list( + cmat = cmat, + mu = mu, + nvec = nvec, + n_cov = n_cov, + r2 = res$r2, + sd = sd, + alpha_level = res$alpha_level, + pow = res$pow, + beta_level = beta_level, + num_df = res$num_df, + den_df = res$den_df, + N_tot = res$N_tot + ) + } + names(con_res) = names(cons) + } + + } else if(is.null(r2)) { + + for(i1 in 1:length(cmats)){ + cmat = cmats[[i1]] + + r2 = optim(par=.5, + fn=function(r2){ abs(eval(p.body)$pow - pow)}, + #c(.001, .999), + upper = .999, + lower = .001, + method= "Brent")$par + # Uniroot produces error + #r2 <- uniroot(function(r2) eval(p.body)$pow - + # pow, c(1e-10, 1 - 1e-10))$root + res <- eval(p.body) + beta_level <- 1-res$pow + pow_res[[i1]] = list( + cmat = cmat, + mu = mu, + nvec = res$nvec, + n_cov = n_cov, + r2 = res$r2, + sd = sd, + alpha_level = res$alpha_level, + pow = res$pow, + beta_level = res$beta_level, + num_df = res$num_df, + den_df = res$den_df, + N_tot = res$N_tot + ) + } + names(pow_res) = names(cmats) + + if(length(cons) != 0){ + for(i1 in 1:length(cons)){ + cmat = cons[[i1]] + r2 <- uniroot(function(r2) eval(p.body)$pow - + pow, c(1e-10, 1 - 1e-10))$root + res <- eval(p.body) + beta_level <- 1-res$pow + con_res[[i1]] = list( + cmat = cmat, + mu = mu, + nvec = nvec, + n_cov = n_cov, + r2 = res$r2, + sd = sd, + alpha_level = res$alpha_level, + pow = res$pow, + beta_level = beta_level, + num_df = res$num_df, + den_df = res$den_df, + N_tot = res$N_tot + ) + } + names(con_res) = names(cons) + } + + } else if (is.null(alpha_level)) { + + for(i1 in 1:length(cmats)){ + cmat = cmats[[i1]] + alpha_level <- uniroot(function(alpha_level) eval(p.body)$pow - + pow, c(1e-10, 1 - 1e-10))$root + res <- eval(p.body) + beta_level <- 1-res$pow + pow_res[[i1]] = list( + cmat = cmat, + mu = mu, + nvec = nvec, + n_cov = n_cov, + r2 = res$r2, + sd = sd, + alpha_level = res$alpha_level, + pow = res$pow, + beta_level = beta_level, + num_df = res$num_df, + den_df = res$den_df, + N_tot = res$N_tot + ) + } + names(pow_res) = names(cmats) + + if(length(cons) != 0){ + for(i1 in 1:length(cons)){ + cmat = cons[[i1]] + alpha_level <- uniroot(function(alpha_level) eval(p.body)$pow - + pow, c(1e-10, 1 - 1e-10))$root + res <- eval(p.body) + beta_level <- 1-res$pow + con_res[[i1]] = list( + cmat = cmat, + mu = mu, + nvec = nvec, + n_cov = n_cov, + r2 = res$r2, + sd = sd, + alpha_level = res$alpha_level, + pow = res$pow, + beta_level = beta_level, + num_df = res$num_df, + den_df = res$den_df, + N_tot = res$N_tot + ) + } + names(con_res) = names(cons) + } + + } else if (is.null(n)){ + + p.body2 <- quote({ + pow_anc_meth(cmat, + mu, + nvec = rep(N_tot / length(mu), length(mu)), + n_cov, + r2, + sd, + alpha_level) + }) + + for(i1 in 1:length(cmats)){ + cmat = cmats[[i1]] + + N_tot <- optim(par=(2*length(mu)+n_cov), + fn=function(N_tot){ abs(eval(p.body2)$pow - pow)}, + c(2*length(mu)+n_cov, 1000000000), + control=list(warn.1d.NelderMead = FALSE))$par + res <- eval(p.body2) + beta_level <- 1-res$pow + if(round_up == TRUE && (!(N_tot%%1==0) || !any(res$nvec%%1==0))) { + N_tot = length(mu)*ceiling(N_tot / length(mu)) + + res <- eval(p.body2) + + } + pow_res[[i1]] = list( + cmat = cmat, + mu = mu, + nvec = res$nvec, + n_cov = res$n_cov, + r2 = res$r2, + sd = sd, + alpha_level = res$alpha_level, + pow = res$pow, + beta_level = res$beta_level, + num_df = res$num_df, + den_df = res$den_df, + N_tot = res$N_tot + ) + } + names(pow_res) = names(cmats) + + if(length(cons) != 0){ + for(i1 in 1:length(cons)){ + cmat = cons[[i1]] + N_tot <- optim(par=(2*length(mu)+n_cov), + fn=function(N_tot){ abs(eval(p.body2)$pow - pow)}, + c(2*length(mu)+n_cov, 1000000000), + control=list(warn.1d.NelderMead = FALSE))$par + res <- eval(p.body2) + beta_level <- 1-res$pow + if(round_up == TRUE && (!(N_tot%%1==0) || !any(res$nvec%%1==0))) { + N_tot = length(mu)*ceiling(N_tot / length(mu)) + + res <- eval(p.body2) + + } + con_res[[i1]] = list( + cmat = cmat, + mu = mu, + nvec = nvec, + n_cov = n_cov, + r2 = res$r2, + sd = sd, + alpha_level = res$alpha_level, + pow = res$pow, + beta_level = beta_level, + num_df = res$num_df, + den_df = res$den_df, + N_tot = res$N_tot + ) + } + names(con_res) = names(cons) + } + + + }else { + stop("internal error: exactly one of n, r2, beta_level, and alpha_level must be NULL") + } + + df_pow = data.frame(factor = character(), + N_tot = integer(), + n_cov = integer(), + r2 = double(), + alpha_level = double(), + beta_level = double(), + power = double()) + + for(i1 in 1:length(pow_res)){ + df_pow[i1,]$factor = factorlist[i1] + df_pow[i1,]$N_tot = pow_res[[i1]]$N_tot + df_pow[i1,]$n_cov = pow_res[[i1]]$n_cov + df_pow[i1,]$r2 = pow_res[[i1]]$r2 + df_pow[i1,]$alpha_level = pow_res[[i1]]$alpha_level + df_pow[i1,]$beta_level = pow_res[[i1]]$beta_level + df_pow[i1,]$power = pow_res[[i1]]$pow*100 + + pow_res2[[i1]] = structure(list(dfs = c(pow_res[[i1]]$num_df, + pow_res[[i1]]$den_df), + N = pow_res[[i1]]$N_tot, + n = pow_res[[i1]]$nvec, + n_cov = pow_res[[i1]]$n_cov, + mu = pow_res[[i1]]$mu, + sd = pow_res[[i1]]$sd, + r2 = pow_res[[i1]]$r2, + alpha_level = pow_res[[i1]]$alpha_level, + beta_level = pow_res[[i1]]$beta_level, + power = pow_res[[i1]]$pow*100, + type = TYPE, + method = METHOD), class = "power.htest") + } + names(pow_res2) = names(pow_res) + rownames(df_pow) = NULL + + if(length(cons) == 0){ + con_res = NULL + con_res2 = NULL + df_con = NULL + } else{ + df_con = data.frame(contrast = character(), + N_tot = integer(), + n_cov = integer(), + r2 = double(), + alpha_level = double(), + beta_level = double(), + power = double()) + + for(i1 in 1:length(con_res)){ + df_con[i1,]$contrast = names(con_res)[i1] + df_con[i1,]$N_tot = con_res[[i1]]$N_tot + df_con[i1,]$n_cov = con_res[[i1]]$n_cov + df_con[i1,]$r2 = con_res[[i1]]$r2 + df_con[i1,]$alpha_level = con_res[[i1]]$alpha_level + df_con[i1,]$beta_level = con_res[[i1]]$beta_level + df_con[i1,]$power = con_res[[i1]]$pow*100 + + con_res2[[i1]] = structure(list(dfs = c(con_res[[i1]]$num_df, + con_res[[i1]]$den_df), + N = con_res[[i1]]$N_tot, + n = con_res[[i1]]$nvec, + n_cov = con_res[[i1]]$n_cov, + contrast = con_res[[i1]]$cmat, + mu = con_res[[i1]]$mu, + sd = con_res[[i1]]$sd, + r2 = con_res[[i1]]$r2, + alpha_level = con_res[[i1]]$alpha_level, + beta_level = con_res[[i1]]$beta_level, + power = con_res[[i1]]$pow*100, + type = TYPE, + method = METHOD2), class = "power.htest") + } + names(con_res2) = names(con_res) + rownames(df_con) = NULL + } + + + structure( + list( + main_results = df_pow, + aov_list = pow_res2, + contrast_results = df_con, + con_list = con_res2, + design_params = list( + design = design, + mu = mu, + sd = sd, + label_list = label_list + ) + ), + class = "ancova_power" + ) + + +} diff --git a/R/ANCOVA_contrast.R b/R/ANCOVA_contrast.R new file mode 100644 index 0000000..824dd5a --- /dev/null +++ b/R/ANCOVA_contrast.R @@ -0,0 +1,228 @@ +#' Power Calculations for ANCOVA Contrasts +#' +#' Complete power analyses for specific ANCOVA contrasts. This function does not support within subjects factors. +#' +#' @param cmat Matrix of the specific contrasts of interest +#' @param mu Vector specifying mean for each condition +#' @param n Sample size in each condition +#' @param sd Standard deviation for all conditions (or a vector specifying the sd for each condition) +#' @param r2 Coefficient of Determination of the model with only the covariates +#' @param n_cov Number of covariates +#' @param alpha_level Alpha level used to determine statistical significance +#' @param beta_level Type II error probability (power/100-1) +#' @param round_up Logical indicator (default = TRUE) for whether to round up sample size calculations to nearest whole number +#' +#' @return +#' Object of class "power.htest", a list of the arguments (including the computed one) augmented with method and note elements. +#' @examples +#' ANCOVA_contrast(cmat = c(-1,1), +#' n = 15, +#' mu = c(0,1), +#' sd = 1, +#' r2 = .2, +#' n_cov = 1) +#' @section References: +#' Shieh, G. (2020). Power analysis and sample size planning in ANCOVA designs. Psychometrika, 85(1), 101-120. +#' @importFrom stats uniroot pf df qf contr.sum dt dbeta qtukey model.matrix terms optim +#' @export + +ANCOVA_contrast <- function(cmat, + mu, + n = NULL, + sd, + r2 = NULL, + n_cov, + alpha_level = Superpower_options("alpha_level"), + beta_level = NULL, + round_up = TRUE) { + + METHOD <- "Power Calculation for ANCOVA contrast" + TYPE = "Exact" + + if(is.vector(cmat)){ + cmat = matrix(cmat, nrow=1) + } + + if(ncol(cmat) != length(mu)){ + stop("cmat must match the length of mu") + } + + #Ensure, if single correlation is input, that it is between 0 and 1 + if(!is.null(r2)){ + if ((r2 <= 0) | r2 >= 1 ) { + stop("Coefficient of Determination must be greater than -1 and less than 1") + } + } + + #Ensure sd is greater than 0 + if (any(sd <= 0) || length(sd) != 1) { + stop("Standard deviation (sd) is less than or equal to zero; input a single value greater than zero") + } + + if( missing(mu) || missing(sd) || missing(n_cov)){ + stop("mu, sd, and n_cov are missing and must be provided") + } + + n_grp = length(mu) + + if(!is.null(n)){ + if(length(n==1)){ + nvec = rep(n,length(mu)) + } else {nvec = n} + + if(length(nvec) != length(mu)){ + stop("Length of N and mu do not match!") + } + } else { + nvec = NULL + } + + if (!is.null(beta_level)){ + pow = 1 - beta_level + } + + p.body = quote({ + pow_anc_meth(cmat, + mu, + nvec, + n_cov, + r2, + sd, + alpha_level) + }) + + if (is.null(beta_level)){ + res <- eval(p.body) + beta_level <- 1-res$pow + con_res = list( + cmat = cmat, + mu = mu, + nvec = nvec, + n_cov = n_cov, + r2 = res$r2, + sd = sd, + alpha_level = res$alpha_level, + pow = res$pow, + beta_level = beta_level, + num_df = res$num_df, + den_df = res$den_df, + N_tot = res$N_tot + ) + + + + } else if(is.null(r2)) { + + r2 = optim( + par = .5, + fn = function(r2) { + abs(eval(p.body)$pow - pow) + }, + #c(.001, .999), + upper = .999, + lower = .001, + method = "Brent" + )$par + + #r2 <- uniroot(function(r2) eval(p.body)$pow - + # pow, c(1e-10, 1 - 1e-10))$root + res <- eval(p.body) + beta_level <- 1-res$pow + con_res = list( + cmat = cmat, + mu = mu, + nvec = nvec, + n_cov = n_cov, + r2 = res$r2, + sd = sd, + alpha_level = res$alpha_level, + pow = res$pow, + beta_level = beta_level, + num_df = res$num_df, + den_df = res$den_df, + N_tot = res$N_tot + ) + + } else if (is.null(alpha_level)) { + + alpha_level <- uniroot(function(alpha_level) eval(p.body)$pow - + pow, c(1e-10, 1 - 1e-10))$root + res <- eval(p.body) + beta_level <- 1-res$pow + con_res = list( + cmat = cmat, + mu = mu, + nvec = nvec, + n_cov = n_cov, + r2 = res$r2, + sd = sd, + alpha_level = res$alpha_level, + pow = res$pow, + beta_level = beta_level, + num_df = res$num_df, + den_df = res$den_df, + N_tot = res$N_tot + ) + + + } else if (is.null(n)){ + + p.body2 <- quote({ + pow_anc_meth(cmat, + mu, + nvec = rep(N_tot / length(mu), length(mu)), + n_cov, + r2, + sd, + alpha_level) + }) + + N_tot <- optim(par=(2*length(mu)+n_cov), + fn=function(N_tot){ abs(eval(p.body2)$pow - pow)}, + c(2*length(mu)+n_cov, 1000000000), + control=list(warn.1d.NelderMead = FALSE))$par + res <- eval(p.body2) + beta_level <- 1-res$pow + if(round_up == TRUE && (!(N_tot%%1==0) || !any(res$nvec%%1==0))) { + N_tot = length(mu)*ceiling(N_tot / length(mu)) + + res <- eval(p.body2) + + } + con_res = list( + cmat = cmat, + mu = mu, + nvec = nvec, + n_cov = n_cov, + r2 = res$r2, + sd = sd, + alpha_level = res$alpha_level, + pow = res$pow, + beta_level = beta_level, + num_df = res$num_df, + den_df = res$den_df, + N_tot = res$N_tot + ) + + + + }else { + stop("internal error: exactly one of n, r2, beta_level, and alpha_level must be NULL") + } + + structure(list(dfs = c(con_res$num_df, + con_res$den_df), + N = con_res$N_tot, + n = con_res$nvec, + n_cov = con_res$n_cov, + contrast = con_res$cmat, + mu = con_res$mu, + sd = con_res$sd, + r2 = con_res$r2, + alpha_level = con_res$alpha_level, + beta_level = con_res$beta_level, + power = con_res$pow*100, + type = TYPE, + method = METHOD), + class = "power.htest") +} diff --git a/R/ANCOVA_factorial_power.R b/R/ANCOVA_factorial_power.R new file mode 100644 index 0000000..76a9f0e --- /dev/null +++ b/R/ANCOVA_factorial_power.R @@ -0,0 +1,78 @@ +pow_anc_meth = function(cmat, + mu, # Vector of Group Means + nvec, # Group Sample Sizes + n_cov, # Number of Covariates + r2, # Coefficient of Determination + sd, # SD UNADJUSTED + alpha_level # Alpha level + ){ + + n_grp = length(mu) + mu = mu + var_e = (sd^2*(1-r2)) + + numint <- 2000 + dd <- 1e-5 + coevec <- c(1, + rep(c(4, 2), numint / 2 - 1), + 4, 1) + bl <- dd + bu <- 1 - dd + intl <- (bu - bl) / numint + bvec <- bl + intl * (0:numint) + cmu <- cmat %*% matrix(mu, n_grp, 1) + # Derive other details + num_df <- nrow(cmat) + N_tot <- sum(nvec) + qmat <- diag(N_tot / nvec) + # Matrix multiplication to get lower case gamma squared (eq. 23 from Shieh) + l_gamma2 <- + t(cmu) %*% solve(cmat %*% qmat %*% t(cmat)) %*% cmu / var_e + + N_tot <- sum(nvec) + den_df <- N_tot - n_grp - n_cov + dfx <- den_df + 1 + b <- n_cov / dfx + # Get critical F-statistic + fcrit <- qf(1 - alpha_level, num_df, den_df) + # Solution differs by number of covariates + if (n_cov == 1) { + tl <- qt(dd, dfx) + tu <- qt(1 - dd, dfx) + intl <- (tu - tl) / numint + tvec <- tl + intl * (0:numint) + wtpdf <- (intl / 3) * coevec * dt(tvec, dfx) + + pow <- sum(wtpdf * pf( + fcrit, + num_df, + den_df, + c(N_tot * l_gamma2) / (1 + b * tvec ^ 2), + lower.tail = FALSE + )) + } + else { + wbpdf <- (intl / 3) * coevec * dbeta(bvec, dfx / 2, n_cov / 2) + pow <- sum(wbpdf * pf(fcrit, + num_df, + den_df, + c(N_tot * l_gamma2) * bvec, + lower.tail = FALSE) + ) + } + + return(list(pow = pow, + num_df = num_df, + den_df = den_df, + N_tot = N_tot, + beta_level = 1-pow, + alpha_level = alpha_level, + cmat = cmat, + mu = mu, + nvec = nvec, + n_cov = n_cov, + r2 = r2, + sd = sd)) +} + + diff --git a/R/ANOVA_design.R b/R/ANOVA_design.R index 7fc90a7..0af831d 100644 --- a/R/ANOVA_design.R +++ b/R/ANOVA_design.R @@ -4,7 +4,8 @@ #' @param mu Vector specifying mean for each condition #' @param sd standard deviation for all conditions (or a vector specifying the sd for each condition) #' @param r Correlation between dependent variables (single value or matrix) -#' @param labelnames Optional vector to specifying factor and condition names (recommended, if not used factors and levels are indicated by letters and numbers) +#' @param label_list An optional list to specify the factor names and condition (recommended, if not used factors and levels are indicated by letters and numbers). +#' @param labelnames Optional vector to specifying factor and condition names. This parameter is deprecated and will be overridden by input from label_list. #' @param plot Should means plot be printed (defaults to TRUE) #' @return Returns single list with simulated data, design, design list, factor names, formulas for ANOVA, means, sd, correlation, sample size per condition, correlation matrix, covariance matrix, design string, labelnames, labelnameslist, factor names, meansplot #' @@ -35,7 +36,8 @@ #' ## with a mean pattern of 1, 0, 1, 0, conditions labeled 'condition' and #' ## 'voice', with names for levels of "cheerful", "sad", and "human", "robot" #' ANOVA_design(design = "2w*2w", n = 40, mu = c(1, 0, 1, 0), sd = 2, r = 0.8, -#' labelnames = c("condition", "cheerful", "sad", "voice", "human", "robot")) +#' label_list= list(condition = c("cheerful", "sad"), +#' voice = c("human", "robot"))) #' @section Warnings: #' Varying the sd or r (e.g., entering multiple values) violates assumptions of homoscedascity and sphericity respectively #' @importFrom stats pnorm pt qnorm qt as.formula median @@ -48,9 +50,10 @@ #' @export #' -ANOVA_design <- function(design, n, mu, sd, r = 0, - labelnames = NULL, - plot = Superpower_options("plot")){ +ANOVA_design <- function(design, n, mu, sd, r = 0, + label_list = NULL, + labelnames = NULL, + plot = Superpower_options("plot")){ #Check String for an acceptable digits and factor (w or b) if (grepl("^(\\d{1,3}(w|b)\\*){0,2}\\d{1,3}(w|b)$", design, ignore.case = FALSE, perl = TRUE) == FALSE) { @@ -79,6 +82,15 @@ ANOVA_design <- function(design, n, mu, sd, r = 0, #Store factor levels (used many times in the script, calculate once) factor_levels <- as.numeric(strsplit(design, "\\D+")[[1]]) + if(!is.null(label_list)){ + + labelnames = vector() + for(i in 1:length(label_list)){ + labelnames = append(labelnames, names(label_list)[i]) + labelnames = append(labelnames, label_list[[i]]) + } + } + if (is.null(labelnames)) { for (i1 in 1:length(factor_levels)){ labelnames <- append(labelnames,paste(paste(letters[i1]), sep = "")) diff --git a/R/ANOVA_exact.R b/R/ANOVA_exact.R index e9ca3af..f0ec14c 100644 --- a/R/ANOVA_exact.R +++ b/R/ANOVA_exact.R @@ -43,7 +43,6 @@ #' @importFrom afex aov_car #' @importFrom graphics pairs #' @importFrom dplyr mutate -#' @import Hmisc #' @import emmeans #' @import ggplot2 @@ -362,7 +361,7 @@ ANOVA_exact <- function(design_result, meansplot2 = meansplot + geom_jitter(position = position_jitter(0.2)) + stat_summary( - fun.data = "mean_sdl", + fun.data = "smean.sdl", fun.args = list(mult = 1), geom = "crossbar", color = "red" @@ -766,7 +765,7 @@ ANOVA_exact2 <- function(design_result, meansplot2 = meansplot + stat_summary( - fun.data = "mean_sdl", + fun.data = "smean.sdl", fun.args = list(mult = 1), geom = "crossbar", color = "red" diff --git a/R/alpha_standardized.R b/R/alpha_standardized.R index 39f67f2..413b5b3 100644 --- a/R/alpha_standardized.R +++ b/R/alpha_standardized.R @@ -1,7 +1,7 @@ #' Compute standardized alpha level based on unstandardized alpha level and the number of observations N. #' @param alpha The unstandardized alpha level (e.g., 0.05), independent of the sample size. #' @param N The number of observations (e.g., the sample size) in the dataset -#' @param standardize_N The nuber of observations (e.g., the sample size) you want to use to standardize the alpha level for. Defaults to 100 (base on Good, 1982). +#' @param standardize_N The number of observations (e.g., the sample size) you want to use to standardize the alpha level for. Defaults to 100 (base on Good, 1982). #' @examples #' ## Check it yields .05 for N = 100: #' alpha_standardized(alpha = 0.05, N = 100) diff --git a/R/ancova_helper_functions.R b/R/ancova_helper_functions.R new file mode 100644 index 0000000..2c258a3 --- /dev/null +++ b/R/ancova_helper_functions.R @@ -0,0 +1,134 @@ + +# Shieh function for one way ANOVA +pwr_ancova_shieh <- function (mu, # Vector of Group Means + n, # Group Sample Sizes + n_cov, # Number of Covariates + r2, # Coefficient of Determination + sd, # SD UNADJUSTED + alpha_level # Alpha level +) { + n_grp = length(mu) + mu = mu + nvec = n + var_e = (sd^2*(1-r2)) + cmat = t(contr.sum(length(nvec))) # Contrast matrix + # Create matrix for Wald statistic (W star; eq 10 from Shieh) + numint <- 2000 + dd <- 1e-5 + coevec <- c(1, + rep(c(4, 2), numint / 2 - 1), + 4, 1) + bl <- dd + bu <- 1 - dd + intl <- (bu - bl) / numint + bvec <- bl + intl * (0:numint) + cmu <- cmat %*% matrix(mu, n_grp, 1) + # Derive other details + num_df <- nrow(cmat) + N_tot <- sum(nvec) + qmat <- diag(N_tot / nvec) + # Matrix multiplication to get lower case gamma squared (eq. 23 from Shieh) + l_gamma2 <- t(cmu) %*% solve(cmat %*% qmat %*% t(cmat)) %*% cmu / var_e + + N_tot <- sum(nvec) + den_df <- N_tot - n_grp - n_cov + dfx <- den_df + 1 + b <- n_cov / dfx + # Get critical F-statistic + fcrit <- qf(1 - alpha_level, num_df, den_df) + # Solution differs by number of covariates + if (n_cov == 1) { + tl <- qt(dd, dfx) + tu <- qt(1 - dd, dfx) + intl <- (tu - tl) / numint + tvec <- tl + intl * (0:numint) + wtpdf <- (intl / 3) * coevec * dt(tvec, dfx) + + pow <- sum(wtpdf * pf( + fcrit, + num_df, + den_df, + c(N_tot * l_gamma2) / (1 + b * tvec ^ 2), + lower.tail = FALSE + )) + } + else { + wbpdf <- (intl / 3) * coevec * dbeta(bvec, dfx / 2, n_cov / 2) + pow <- sum(wbpdf * pf(fcrit, + num_df, + den_df, + c(N_tot * l_gamma2) * bvec, + lower.tail = FALSE)) + } + + + return(pow) + +} + +# Keppel (PASS) function for one way ANOVA +pwr_ancova_keppel = function(mu, # Vector of Group Means + n, # Group Sample Sizes + n_cov, # Number of Covariates + r2, # Coefficient of Determination + sd, # SD UNADJUSTED + alpha_level # Alpha level +){ + num_df = length(mu) - 1 + N_tot = sum(n) + + pow = pf(qf(alpha_level, + num_df, + (N_tot - length(mu) - n_cov ), + lower.tail = FALSE), + num_df, + (N_tot - length(mu) - n_cov ), + (anc_var_m(n,mu)/(sd^2*(1-r2))) * (N_tot / length(mu)) * length(mu), + lower.tail = FALSE) + + return(pow) +} + +# helper function for adding +anc_var_m = function(n, mu) { + N = sum(n) + + G = length(mu) + + mu_bar = sum((n/N) * mu) + + N_bar = N/G + + var_m = sum((n/N)*(mu-mu_bar)^2) + + return(var_m) + +} + +# selection of type of method for ANCOVA power + +pwr_method <- function(mu, + n, + n_cov, + r2, + sd, + alpha_level, type) { + switch(type, + exact = pwr_ancova_shieh(mu, + n, + n_cov, + r2, + sd, + alpha_level), + approx = pwr_ancova_keppel(mu, + n, + n_cov, + r2, + sd, + alpha_level)) +} + + + + + diff --git a/R/data_generation.R b/R/data_generation.R deleted file mode 100644 index 4801c15..0000000 --- a/R/data_generation.R +++ /dev/null @@ -1,52 +0,0 @@ -#' @importFrom MASS mvrnorm -#' @importFrom reshape2 melt dcast -#' @importFrom stats lm confint -#' @importFrom magrittr '%>%' -#' @importFrom dplyr select mutate everything - -ancova_gen <- function(post_diff, sd, r, alpha_level, conf_level, grp_n) { - cor_mat = matrix(c(1,r, - r,1), - nrow = 2) - - #Need this to avoid "undefined" global error from occuring - id <- variable <- NULL - - cov_mat = cor_mat*(sd)^2 - - df_1 <- as.data.frame(mvrnorm(n = grp_n, - mu = c(0, post_diff), - Sigma = cov_mat)) - colnames(df_1) = c("y_pre", "y_post") - - df_1$id = rownames(df_1) - df_1$id = paste0("g1_",df_1$id) - - df_2 <- as.data.frame(mvrnorm(n = grp_n, - mu = c(0, 0), - Sigma = cov_mat)) - colnames(df_2) = c("y_pre", "y_post") - - df_2$id = rownames(df_2) - df_2$id = paste0("g2_",df_2$id) - - df_sim = rbind(df_1,df_2) - - df_sim <- melt(df_sim, id.vars = c("id")) %>% - mutate(group = ifelse(grepl("g1",id),"group1","group2"), - time = ifelse(grepl("pre",variable),"pre","post")) - - tmp = dcast(df_sim, id+group ~ time, value.var = "value") - tmp$diff = tmp$post-tmp$pre - - - fit_lm <- lm(post ~ pre + group, tmp) - width_ac <- diff(confint(fit_lm, level = (conf_level))[3,])[[1]]/2 - p_ac = as.numeric(summary(fit_lm)[["coefficients"]][, "Pr(>|t|)"][3]) - fit_aov <- lm(diff ~ group, tmp) - width_aov <- diff(confint(fit_aov, level = (conf_level))[2,])[[1]]/2 - p_aov = as.numeric(summary(fit_aov)[["coefficients"]][, "Pr(>|t|)"][2]) - - ancova_df = data.frame(width_ac = width_ac, p_ac = p_ac, - width_aov = width_aov, p_aov = p_aov) -} \ No newline at end of file diff --git a/R/globals.R b/R/globals.R index 5cc8e61..8d0b815 100644 --- a/R/globals.R +++ b/R/globals.R @@ -1,3 +1,3 @@ # Globals -utils::globalVariables(c("monte_gen","exact_gen","lambda","cohen_f")) \ No newline at end of file +utils::globalVariables(c("monte_gen","exact_gen","lambda","cohen_f","contrasts<-")) \ No newline at end of file diff --git a/R/methods.ancova_power.R b/R/methods.ancova_power.R new file mode 100644 index 0000000..e8323b6 --- /dev/null +++ b/R/methods.ancova_power.R @@ -0,0 +1,74 @@ +#' Methods for ancova_power objects +#' +#' Methods defined for objects returned from the ANCOVA_analytic function. +#' +#' @param x object of class \code{ancova_power} as returned from one of the simulation functions in Superpower. +#' @param ... further arguments passed through, see description of return value +#' @return +#' \describe{ +#' \item{\code{print}}{Prints short summary of the simulation result} +#' \item{\code{plot}}{Returns a meansplot of from the defined design} +#' } +#' @name ancova_power-methods + + +### methods for ancova_power + +#' @rdname ancova_power-methods +#' @method print ancova_power +#' @export + +print.ancova_power <- function(x,...){ + + if (!is.null(x$main_results)) { + pow_tab1 = data.frame(`Total N`= x$main_results$N_tot, + `Covariates`= x$main_results$n_cov, + `r2`= x$main_results$r2, + `Alpha Level`= x$main_results$alpha_level, + `Beta Level`= x$main_results$beta_level, + `Power` = x$main_results$power, + check.names = FALSE) + row.names(pow_tab1) = x$main_results$factor + + cat("Power Analysis Results for ANCOVA") + cat("\n") + print(pow_tab1, digits = 4) + cat("\n") + } + if (!is.null(x$contrast_results)) { + pow_tab2 = data.frame(`Total N`= x$contrast_results$N_tot, + `Covariates`= x$contrast_results$n_cov, + `r2`= x$contrast_results$r2, + `Alpha Level`= x$contrast_results$alpha_level, + `Beta Level`= x$contrast_results$beta_level, + `Power` = x$contrast_results$power, + check.names = FALSE) + row.names(pow_tab2) = x$contrast_results$factor + cat("Power Analysis Results for ANCOVA contrasts") + cat("\n") + print(pow_tab2, digits = 4) + } + + +} + +#' @rdname ancova_power-methods +#' @method plot ancova_power +#' @import ggplot2 +#' @export + +plot.ancova_power <- function(x,...){ + if(is.null(x$design_params$design)){ + stop("No plotting method when design is not provided.") + } + + plot(ANOVA_design(design = x$design_params$design, + n = 40, + mu = x$design_params$mu, + sd = x$design_params$sd, + r = 0, + label_list = x$design_params$label_list, + plot = FALSE)) + +} + diff --git a/R/misc_functions.R b/R/misc_functions.R index 9fbbc55..355b2b8 100644 --- a/R/misc_functions.R +++ b/R/misc_functions.R @@ -13,4 +13,18 @@ ci_binom = function(centx,n,conf_level=.95){ ucl <- (p1 + p2)/p3 res = c(lcl,ucl) return(res) +} + + + +smean.sdl = function (x, mult = 1, na.rm = TRUE) { + if (na.rm) + x <- x[!is.na(x)] + n <- length(x) + if (n == 0) + return(c(Mean = NA, Lower = NA, Upper = NA)) + xbar <- sum(x)/n + sd <- sqrt(sum((x - xbar)^2)/(n - 1)) + c(Mean = xbar, Lower = xbar - mult * sd, Upper = xbar + + mult * sd) } \ No newline at end of file diff --git a/R/oneway_ancova.R b/R/oneway_ancova.R new file mode 100644 index 0000000..f7c542f --- /dev/null +++ b/R/oneway_ancova.R @@ -0,0 +1,189 @@ +#' Power Calculations for a one-way ANCOVA +#' +#' Compute power of ANCOVA omnibus test (power_oneway_ancova) or contrast (power_oneway_ancova) for one-way (single factor), between subjects designs. +#' +#' @param n Sample size in each condition. +#' @param mu Vector specifying mean for each condition. +#' @param n_cov Number of covariates. +#' @param r2 Coefficient of determination (r^2) of the combined covariates. +#' @param sd Standard deviation for all conditions (residual SD without covariate adjustment). +#' @param alpha_level Alpha level used to determine statistical significance. +#' @param beta_level Type II error probability (power/100-1) +#' @param round_up Logical indicator for whether to round up the sample size(s) to a whole number. Default is TRUE. +#' @param type Sets the method for estimating power. "exact" will use the Shieh (2020) approach while "approx" will use the Keppel (1991) approach. +#' +#' @return +#' dfs = degrees of freedom, +#' N = Total sample size, +#' n = Sample size per group/condition, +#' n_cov = Number of covariates, +#' mu = Mean for each condition, +#' sd = Standard deviation, +#' r2 = Coefficient of determination of combined covariates. +#' alpha_level = Type 1 error probability, +#' beta_level = Type 2 error probability, +#' power = Power of test (1-beta_level\*100%), +#' type = Method (Shieh or Keppel) for estimating power +#' +#' @examples +#' # Example from Table 1 Shieh 2020 +#' power_oneway_ancova(mu = c(400, 450, 500), n = c(21,21,21), +#' r2 = .1^2, sd = 100) +#' @section References: +#' Keppel, G. (1991). Design and Analysis A Researcher's Handbook. 3rd Edition. Prentice Hall. Englewood Cliffs, New Jersey. See pages 323 - 324. +#' Shieh, G. (2017). Power and sample size calculations for contrast analysis in ANCOVA. Multivariate behavioral research, 52(1), 1-11. +#' Shieh, G. (2020). Power analysis and sample size planning in ANCOVA designs. Psychometrika, 85(1), 101-120. +#' @importFrom stats uniroot pf df qf contr.sum dt dbeta qtukey optim +#' @export +#' + + +power_oneway_ancova <- function(n = NULL, + mu = NULL, + n_cov = 1, + r2 = NULL, + sd = 1, + alpha_level = Superpower_options('alpha_level'), + beta_level = NULL, + round_up = TRUE, + type = "exact"){ + + if (!is.null(alpha_level) && !is.numeric(alpha_level) || any(0 > + alpha_level | alpha_level > 1)) { + stop(sQuote("alpha_level"), " must be numeric in [0, 1]") + } + + if (!is.null(beta_level) && !is.numeric(beta_level) || any(0 > beta_level | + beta_level > 1)) { + stop(sQuote("beta_level"), " must be numeric in [0, 1].") + } + + if(is.null(mu) || length(mu) <= 1){ + stop("mu cannot be NULL or have length of 1. Please provide group means.") + } + + if(!is.null(n) && length(n) != length(mu)){ + stop("length of n and mu must match.") + } + + if(length(sd) > 1){ + stop("sd can only have length of 1.") + } + + if (!is.null(beta_level)){ + pow = 1 - beta_level + } + #var_e = sd^1*(1-r2) + + if(type != "exact" && type != "approx"){ + stop("type of power calculation can only be exact or approx") + } + + + if(is.null(n_cov) || !is.numeric(n_cov) || length(n_cov) > 1) { + stop("Please provide number of covariates (n_cov).") + } + + if(!is.null(r2) && (r2 >= 1 || r2 <= 0)){ + stop(sQuote("r2"), " must be numeric in [0, 1]") + } + + if(!is.null(n)) { + N_tot = sum(n) + } else { + N_tot = NULL + } + + #internal error: exactly one of n, r2, beta_level, and alpha_level must be NULL + if(sum(c(is.null(n), + is.null(r2), + is.null(beta_level), + is.null(alpha_level))) != 1){ + stop("internal error: exactly one of n, r2, beta_level, and alpha_level must be NULL") + } + + num_df = length(mu) - 1 + + + p.body <- quote({ + pwr_method(mu, + n, + n_cov, + r2, + sd, + alpha_level, + type) + }) + + if (!is.null(beta_level)){ + pow = 1 - beta_level + } + + if (is.null(beta_level)){ + pow <- eval(p.body) + } else if(is.null(r2)) { + r2 <- uniroot(function(r2) eval(p.body) - + pow, c(1e-10, 1 - 1e-10))$root + } else if (is.null(alpha_level)) { + alpha_level <- uniroot(function(alpha_level) eval(p.body) - + pow, c(1e-10, 1 - 1e-10))$root + } else if (is.null(n)){ + + p.body2 <- quote({ + pwr_method(mu, + n = rep(N_tot/length(mu), length(mu)), + n_cov, + r2, + sd, + alpha_level, + type) + }) + + + N_tot <- optim(par=(2*length(mu)+n_cov), + fn=function(N_tot){ abs(eval(p.body2) - pow)}, + c(2*length(mu)+n_cov, 1000000000), + control=list(warn.1d.NelderMead = FALSE))$par + n_1 = N_tot / length(mu) + n = rep(n_1, length(mu)) + }else { + stop("internal error: exactly one of n, r2, beta_level, and alpha_level must be NULL") + } + + if(round_up == TRUE && (!(N_tot%%1==0) || !any(n%%1==0))) { + N_tot = length(mu)*ceiling(N_tot / length(mu)) + n_1 = N_tot / length(mu) + n = rep(n_1, length(mu)) + + pow = pwr_method(mu = mu, + n = n, + n_cov = n_cov, + r2 = r2, + sd = sd, + alpha_level = alpha_level, + type = type) + beta_level = 1-pow + + } + + + den_df = N_tot - length(mu) - n_cov + #sd_m = sqrt(anc_var_m(n,mu)) + #sd_e = sqrt((sd^2*(1-r2))) + power_final = pow * 100 + beta_level = 1 - pow + METHOD <- "Power Calculation for 1-way ANCOVA" + TYPE = type + structure(list(dfs = c(num_df, den_df), + N = N_tot, + n = n, + n_cov = n_cov, + mu = mu, + sd = sd, + r2 = r2, + alpha_level = alpha_level, + beta_level = beta_level, + power = power_final, + type = TYPE, + method = METHOD), class = "power.htest") +} diff --git a/R/p_standardized.R b/R/p_standardized.R index 348fda4..b8d3d87 100644 --- a/R/p_standardized.R +++ b/R/p_standardized.R @@ -1,7 +1,7 @@ #' Compute standardized alpha level based on unstandardized alpha level and the number of observations N. #' @param p The observed p-value. #' @param N The number of observations (e.g., the sample size) in the dataset -#' @param standardize_N The nuber of observations (e.g., the sample size) you want to use to standardize the alpha level for. Defaults to 100 (base on Good, 1982). +#' @param standardize_N The number of observations (e.g., the sample size) you want to use to standardize the alpha level for. Defaults to 100 (base on Good, 1982). #' @examples #' ## Check it yields .05 for N = 100: #' p_standardized(p = 0.05, N = 100) diff --git a/README.Rmd b/README.Rmd index 81e3a46..cdc446b 100644 --- a/README.Rmd +++ b/README.Rmd @@ -53,14 +53,14 @@ devtools::install_github("arcaldwell49/Superpower") ``` ## ANOVA_design function -Currently the ANOVA_design function can create designs up three factors, for both within, between, and mixed designs. It requires the following input: design, n, mu, sd, r, and optionally allows you to set labelnames. +Currently the ANOVA_design function can create designs up three factors, for both within, between, and mixed designs. It requires the following input: design, n, mu, sd, r, and optionally allows you to set label_list. 1. design: string that specifies the design (see below). 2. n: the sample size for each between subject condition. 3. mu: a vector with the means for each condition. 4. sd: the population standard deviation. Assumes homogeneity of variances (only one standard deviation can be provided). 5. r: the correlation for within designs (or 0 for between designs). -6. labelnames: This is an optional vector of words that indicates factor names and level names (see below). +6. label_list: This is an optional named list of words that indicates factor names and level names (see below). 7. A final optional setting is to specify if you want to output a plot or not (plot = TRUE or FALSE) ### Specifying the design using 'design' @@ -79,9 +79,9 @@ Note that for each cell in the design, a mean must be provided. Thus, for a "2b* Means need to be entered in the correct order. ANOVA_design outputs a plot so you can check if you entered all means as you intended. Always carefully check if the plot that is generated matches your expectations. -The general principle is that the code generates factors, indicated by the factor names you entered in the labelnames variable, (i.e., *condition* and *time*). Levels are indicated by factor names and levels (e.g., control_time1, control_time2, control_time3, etc). +The general principle is that the code generates factors, indicated by the factor names you entered in the label_list variable, (i.e., *condition* and *time*). Levels are indicated by factor names and levels (e.g., control_time1, control_time2, control_time3, etc). -If your design has just one factor, just enter the means in the same order as the labelnames (see below). For more factors, note the general pattern in the example below. Means are entered in the following order for a 3 factors design (each with 2 levels): +If your design has just one factor, just enter the means in the same order as the label_list (see below). For more factors, note the general pattern in the example below. Means are entered in the following order for a 3 factors design (each with 2 levels): 1. a1 b1 c1 2. a1 b1 c2 @@ -107,7 +107,7 @@ design_result <- ANOVA_design(design = "2b*2w*2b", ### Specifying label names -To make sure the plots and tables with simulation results are easy to interpret, it really helps to name all factors and levels. You can enter the labels in the 'labelnames' variable. You can also choose not to specify names. Then all factors are indicated by letters (a, b, c) and all levels by numbers (a1, a2, a3). +To make sure the plots and tables with simulation results are easy to interpret, it really helps to name all factors and levels. You can enter the labels in the 'label_list' variable. You can also choose not to specify names. Then all factors are indicated by letters (a, b, c) and all levels by numbers (a1, a2, a3). For the 2x3 design we have been using as an example, where there are 2 factors (condition and time of measurement), the first with 2 levels (placebo vs. medicine) and the second with three levels (time1, time2, and time3) we would enter the labels as follows: @@ -255,7 +255,9 @@ design_result <- ANOVA_design(design = "2b*2w", mu = c(1.03, 1.41, 0.98, 1.01), sd = 1.03, r = 0.8, - labelnames = c("voice", "human", "robot", "emotion", "cheerful", "sad")) + label_list = list("voice" = c( "human", "robot"), + "emotion" = c( "cheerful", "sad")) +) power_result_vig_1 <- readRDS(file = "vignettes/sim_data/power_result_vig_1.rds") @@ -282,13 +284,13 @@ design <- "2b" n <- 100 mu <- c(24, 26.2) sd <- 6.4 -labelnames <- c("condition", "control", "pet") # +label_list1 <- list("condition" =c ("control", "pet") ) design_result <- ANOVA_design(design = design, n = n, mu = mu, sd = sd, - labelnames = labelnames) + label_list = label_list1) power_result_vig_2 <- readRDS(file = "vignettes/sim_data/power_result_vig_2.rds") @@ -331,7 +333,7 @@ design_result <- ANOVA_design(design = "2b", n = 100, mu = c(24, 26.2), sd = 6.4, - labelnames = c("condition", "control", "pet")) + label_list = list("condition"= c("control", "pet"))) ANOVA_exact(design_result)$main_results$power # power of 67.7 is a bit low. Let's increase it a bit to n = 150 to see if we are closer to our goal of 90% power. @@ -340,7 +342,7 @@ design_result <- ANOVA_design(design = "2b", n = 150, mu = c(24, 26.2), sd = 6.4, - labelnames = c("condition", "control", "pet"), + label_list = list("condition"= c("control", "pet")), plot = FALSE) ANOVA_exact(design_result)$main_results$power @@ -350,7 +352,7 @@ design_result <- ANOVA_design(design = "2b", n = 175, mu = c(24, 26.2), sd = 6.4, - labelnames = c("condition", "control", "pet"), + label_list = list("condition"= c("control", "pet")), plot = FALSE) ANOVA_exact(design_result)$main_results$power @@ -360,7 +362,7 @@ design_result <- ANOVA_design(design = "2b", n = 180, mu = c(24, 26.2), sd = 6.4, - labelnames = c("condition", "control", "pet"), + label_list = list("condition"= c("control", "pet")), plot = FALSE) ANOVA_exact(design_result)$main_results$power @@ -392,7 +394,7 @@ design_result <- ANOVA_design(design = "2b", n = 180, mu = c(24, 26), sd = 6.4, - labelnames = c("condition", "control", "pet"), + label_list = list("condition"= c("control", "pet")), plot = FALSE) plot_power(design_result, max_n = 250, @@ -406,7 +408,7 @@ design_result <- ANOVA_design(design = "2b", n = 180, mu = c(24, 26), sd = 6.8, - labelnames = c("condition", "control", "pet"), + label_list = list("condition"= c("control", "pet")), plot = FALSE) plot_power(design_result, min_n = 10, max_n = 250, diff --git a/README.md b/README.md index 48c39a0..b3743e2 100644 --- a/README.md +++ b/README.md @@ -38,7 +38,7 @@ devtools::install_github("arcaldwell49/Superpower") Currently the ANOVA\_design function can create designs up three factors, for both within, between, and mixed designs. It requires the following input: design, n, mu, sd, r, and optionally allows you to set -labelnames. +label\_list. 1. design: string that specifies the design (see below). 2. n: the sample size for each between subject condition. @@ -46,7 +46,7 @@ labelnames. 4. sd: the population standard deviation. Assumes homogeneity of variances (only one standard deviation can be provided). 5. r: the correlation for within designs (or 0 for between designs). -6. labelnames: This is an optional vector of words that indicates +6. label\_list: This is an optional named list of words that indicates factor names and level names (see below). 7. A final optional setting is to specify if you want to output a plot or not (plot = TRUE or FALSE) @@ -78,12 +78,12 @@ plot so you can check if you entered all means as you intended. Always carefully check if the plot that is generated matches your expectations. The general principle is that the code generates factors, indicated by -the factor names you entered in the labelnames variable, (i.e., +the factor names you entered in the label\_list variable, (i.e., *condition* and *time*). Levels are indicated by factor names and levels (e.g., control\_time1, control\_time2, control\_time3, etc). If your design has just one factor, just enter the means in the same -order as the labelnames (see below). For more factors, note the general +order as the label\_list (see below). For more factors, note the general pattern in the example below. Means are entered in the following order for a 3 factors design (each with 2 levels): @@ -112,7 +112,7 @@ works. To make sure the plots and tables with simulation results are easy to interpret, it really helps to name all factors and levels. You can enter -the labels in the ‘labelnames’ variable. You can also choose not to +the labels in the ‘label\_list’ variable. You can also choose not to specify names. Then all factors are indicated by letters (a, b, c) and all levels by numbers (a1, a2, a3). @@ -413,7 +413,9 @@ design_result <- ANOVA_design(design = "2b*2w", mu = c(1.03, 1.41, 0.98, 1.01), sd = 1.03, r = 0.8, - labelnames = c("voice", "human", "robot", "emotion", "cheerful", "sad")) + label_list = list("voice" = c( "human", "robot"), + "emotion" = c( "cheerful", "sad")) +) ``` ![](README_files/figure-gfm/unnamed-chunk-8-1.png) @@ -485,13 +487,13 @@ design <- "2b" n <- 100 mu <- c(24, 26.2) sd <- 6.4 -labelnames <- c("condition", "control", "pet") # +label_list1 <- list("condition" =c ("control", "pet") ) design_result <- ANOVA_design(design = design, n = n, mu = mu, sd = sd, - labelnames = labelnames) + label_list = label_list1) ``` ![](README_files/figure-gfm/unnamed-chunk-9-1.png) @@ -575,7 +577,7 @@ design_result <- ANOVA_design(design = "2b", n = 100, mu = c(24, 26.2), sd = 6.4, - labelnames = c("condition", "control", "pet")) + label_list = list("condition"= c("control", "pet"))) ``` ![](README_files/figure-gfm/unnamed-chunk-12-1.png) @@ -601,7 +603,7 @@ design_result <- ANOVA_design(design = "2b", n = 150, mu = c(24, 26.2), sd = 6.4, - labelnames = c("condition", "control", "pet"), + label_list = list("condition"= c("control", "pet")), plot = FALSE) ANOVA_exact(design_result)$main_results$power @@ -624,7 +626,7 @@ design_result <- ANOVA_design(design = "2b", n = 175, mu = c(24, 26.2), sd = 6.4, - labelnames = c("condition", "control", "pet"), + label_list = list("condition"= c("control", "pet")), plot = FALSE) ANOVA_exact(design_result)$main_results$power @@ -647,7 +649,7 @@ design_result <- ANOVA_design(design = "2b", n = 180, mu = c(24, 26.2), sd = 6.4, - labelnames = c("condition", "control", "pet"), + label_list = list("condition"= c("control", "pet")), plot = FALSE) ANOVA_exact(design_result)$main_results$power @@ -712,7 +714,7 @@ design_result <- ANOVA_design(design = "2b", n = 180, mu = c(24, 26), sd = 6.4, - labelnames = c("condition", "control", "pet"), + label_list = list("condition"= c("control", "pet")), plot = FALSE) plot_power(design_result, max_n = 250, @@ -736,7 +738,7 @@ design_result <- ANOVA_design(design = "2b", n = 180, mu = c(24, 26), sd = 6.8, - labelnames = c("condition", "control", "pet"), + label_list = list("condition"= c("control", "pet")), plot = FALSE) plot_power(design_result, min_n = 10, max_n = 250, diff --git a/Validation/ANCOVA_oneway_sims.Rmd b/Validation/ANCOVA_oneway_sims.Rmd new file mode 100644 index 0000000..d28ef42 --- /dev/null +++ b/Validation/ANCOVA_oneway_sims.Rmd @@ -0,0 +1,208 @@ +--- +title: "Factorial ANCOVA Validation: One Way Design" +output: github_document +--- + +```{r setup, include=TRUE} +knitr::opts_chunk$set(echo = TRUE) +library(simstudy) +library(data.table) +library(Superpower) +library(car) +library(testthat) +set.seed(61036) +``` + +# Simulation 1 + +```{r onewayrho1} + +nsim = 5000 +res_group = vector() +res_aov = vector() + +for(i in 1:nsim){ + # Generate outcome and covariate matrix + dt <- genCorData( + 21 * 3, + mu = c(0, 400), + sigma = 100, + rho = .1, + corstr = "cs", + cnames = c("cov", "y1") + ) + +# generate random treatment assignment +study1 = trtAssign(dt, + n = 3, + balanced = TRUE, + grpName = "group") + +# add treatment effect +study1[, y := ifelse(group == 3, y1+100, + ifelse(group == 2, y1+50, + y1))] +# factorize group +study1[, grp := as.factor(group)] + +a = as.data.frame(Anova(lm(y ~ cov + grp, data = study1), + type ='III')) + +res_group = append(a["grp",]$`Pr(>F)`, + res_group) + +a2 = as.data.frame(Anova(lm(y ~ grp, data = study1), + type ='III')) + +res_aov = append(a2["grp",]$`Pr(>F)`, + res_aov) + +} + +print(mean(res_group < .05)) + +``` + +```{r} + +test_that("Simulation 1", { + # Both should be about .81 + # Check ANCOVA result + expect_equal(mean(res_group < .05), + ANCOVA_analytic(design = "3b", + mu = c(400,450,500), + n = 21, + sd = 100, + r2 = .1^2, + n_cov = 1, + alpha_level = .05, + beta_level = NULL)$main_results$power/100, + tolerance = .01) + # Check ANOVA result + expect_equal(mean(res_aov < .05), .8148, + tolerance = .01) +}) + +``` + +# Simulation 2 + +```{r} +res_group2 = vector() + +for(i in 1:nsim){ + # Generate outcome and covariate matrix + dt <- genCorData( + 48, + mu = c(0, 400), + sigma = 100, + rho = .5, + corstr = "cs", + cnames = c("cov", "y1") + ) + + # generate random treatment assignment + study1 = trtAssign(dt, + n = 3, + balanced = TRUE, + grpName = "group") + + # add treatment effect + study1[, y := ifelse(group == 3, y1+100, + ifelse(group == 2, y1+50, + y1))] + # factorize group + study1[, grp := as.factor(group)] + + a = as.data.frame(Anova(lm(y ~ cov + grp, data = study1), + type ='III')) + + res_group2 = append(a["grp",]$`Pr(>F)`, + res_group2) + +} + +``` + +```{r} + +test_that("Simulation 2",{ + # Again should be ~80% + # Check ANCOVA result + expect_equal(mean(res_group2 < .05), + ANCOVA_analytic(design = "3b", + mu = c(400,450,500), + n = 48/3, + sd = 100, + r2 = .5^2, + n_cov = 1, + alpha_level = .05, + beta_level = NULL)$main_results$power/100, + tolerance = .01) + +}) + +print(mean(res_group2 < .05)) +``` + +# Simulation 3 + +```{r} +res_group3 = vector() + +for(i in 1:nsim){ + # Generate outcome and covariate matrix + dt <- genCorData( + 18, + mu = c(0, 400), + sigma = 100, + rho = .9, + corstr = "cs", + cnames = c("cov", "y1") + ) + + # generate random treatment assignment + study1 = trtAssign(dt, + n = 3, + balanced = TRUE, + grpName = "group") + + # add treatment effect + study1[, y := ifelse(group == 3, y1+100, + ifelse(group == 2, y1+50, + y1))] + # factorize group + study1[, grp := as.factor(group)] + + a = as.data.frame(Anova(lm(y ~ cov + grp, data = study1), + type ='III')) + + res_group3 = append(a["grp",]$`Pr(>F)`, + res_group3) + +} + +print(mean(res_group3 < .05)) +``` + +```{r} + +test_that("Simulation 3",{ + # Again should be ~80% + # Check ANCOVA result + expect_equal(mean(res_group3 < .05), + ANCOVA_analytic(design = "3b", + mu = c(400,450,500), + n = 6, + sd = 100, + r2 = .9^2, + n_cov = 1, + alpha_level = .05, + beta_level = NULL)$main_results$power/100, + tolerance = .01) + +}) + +``` + + diff --git a/Validation/ANCOVA_oneway_sims.md b/Validation/ANCOVA_oneway_sims.md new file mode 100644 index 0000000..55d788e --- /dev/null +++ b/Validation/ANCOVA_oneway_sims.md @@ -0,0 +1,211 @@ +Factorial ANCOVA Validation +================ + +``` r +knitr::opts_chunk$set(echo = TRUE) +library(simstudy) +library(data.table) +library(Superpower) +``` + + ## Registered S3 methods overwritten by 'lme4': + ## method from + ## cooks.distance.influence.merMod car + ## influence.merMod car + ## dfbeta.influence.merMod car + ## dfbetas.influence.merMod car + +``` r +library(car) +``` + + ## Loading required package: carData + +``` r +library(testthat) +set.seed(61036) +``` + +# Simulation 1 + +``` r +nsim = 5000 +res_group = vector() +res_aov = vector() + +for(i in 1:nsim){ + # Generate outcome and covariate matrix + dt <- genCorData( + 21 * 3, + mu = c(0, 400), + sigma = 100, + rho = .1, + corstr = "cs", + cnames = c("cov", "y1") + ) + +# generate random treatment assignment +study1 = trtAssign(dt, + n = 3, + balanced = TRUE, + grpName = "group") + +# add treatment effect +study1[, y := ifelse(group == 3, y1+100, + ifelse(group == 2, y1+50, + y1))] +# factorize group +study1[, grp := as.factor(group)] + +a = as.data.frame(Anova(lm(y ~ cov + grp, data = study1), + type ='III')) + +res_group = append(a["grp",]$`Pr(>F)`, + res_group) + +a2 = as.data.frame(Anova(lm(y ~ grp, data = study1), + type ='III')) + +res_aov = append(a2["grp",]$`Pr(>F)`, + res_aov) + +} +``` + +``` r +test_that("Simulation 1", { + # Both should be about .81 + # Check ANCOVA result + expect_equal(mean(res_group < .05), + ANCOVA_analytic(design = "3b", + mu = c(400,450,500), + n = 21, + sd = 100, + r2 = .1^2, + n_cov = 1, + alpha_level = .05, + beta_level = NULL)$main_result$power/100, + tolerance = .01) + # Check ANOVA result + expect_equal(mean(res_aov < .05), .8148, + tolerance = .01) +}) +``` + + ## Test passed 🎉 + +# Simulation 2 + +``` r +res_group2 = vector() + +for(i in 1:nsim){ + # Generate outcome and covariate matrix + dt <- genCorData( + 48, + mu = c(0, 400), + sigma = 100, + rho = .5, + corstr = "cs", + cnames = c("cov", "y1") + ) + + # generate random treatment assignment + study1 = trtAssign(dt, + n = 3, + balanced = TRUE, + grpName = "group") + + # add treatment effect + study1[, y := ifelse(group == 3, y1+100, + ifelse(group == 2, y1+50, + y1))] + # factorize group + study1[, grp := as.factor(group)] + + a = as.data.frame(Anova(lm(y ~ cov + grp, data = study1), + type ='III')) + + res_group2 = append(a["grp",]$`Pr(>F)`, + res_group2) + +} +``` + +``` r +test_that("Simulation 2",{ + # Again should be ~80% + # Check ANCOVA result + expect_equal(mean(res_group2 < .05), + ANCOVA_analytic(design = "3b", + mu = c(400,450,500), + n = 48/3, + sd = 100, + r2 = .5^2, + n_cov = 1, + alpha_level = .05, + beta_level = NULL)$main_result$power/100, + tolerance = .01) + +}) +``` + + ## Test passed 😀 + +# Simulation 3 + +``` r +res_group3 = vector() + +for(i in 1:nsim){ + # Generate outcome and covariate matrix + dt <- genCorData( + 18, + mu = c(0, 400), + sigma = 100, + rho = .9, + corstr = "cs", + cnames = c("cov", "y1") + ) + + # generate random treatment assignment + study1 = trtAssign(dt, + n = 3, + balanced = TRUE, + grpName = "group") + + # add treatment effect + study1[, y := ifelse(group == 3, y1+100, + ifelse(group == 2, y1+50, + y1))] + # factorize group + study1[, grp := as.factor(group)] + + a = as.data.frame(Anova(lm(y ~ cov + grp, data = study1), + type ='III')) + + res_group3 = append(a["grp",]$`Pr(>F)`, + res_group3) + +} +``` + +``` r +test_that("Simulation 3",{ + # Again should be ~80% + # Check ANCOVA result + expect_equal(mean(res_group3 < .05), + ANCOVA_analytic(design = "3b", + mu = c(400,450,500), + n = 6, + sd = 100, + r2 = .9^2, + n_cov = 1, + alpha_level = .05, + beta_level = NULL)$main_result$power/100, + tolerance = .01) + +}) +``` + + ## Test passed 🌈 diff --git a/Validation/ANCOVA_threeway_sims.Rmd b/Validation/ANCOVA_threeway_sims.Rmd new file mode 100644 index 0000000..3fe2f86 --- /dev/null +++ b/Validation/ANCOVA_threeway_sims.Rmd @@ -0,0 +1,137 @@ +--- +title: "Factorial ANCOVA Validation: Two Way Design" +output: github_document +--- + +```{r setup, include=TRUE} +knitr::opts_chunk$set(echo = TRUE) +library(simstudy) +library(data.table) +library(Superpower) +library(car) +library(testthat) +set.seed(19256) +``` + +# Simulation 1 + +```{r } + +nsim = 5000 +res_2x2x2a = vector() +res_2x2x2b = vector() +res_2x2x2ab = vector() +res_2x2x2abc = vector() +res_aov = vector() + +for(i in 1:nsim){ + # Generate outcome and covariate matrix + dt <- genCorData( + 10 * 8, + mu = c(0, 0), + sigma = 2.5, + rho = .33, + corstr = "cs", + cnames = c("cov", "y1") + ) + +# generate random treatment assignment +study1 = trtAssign(dt, + n = 8, + balanced = TRUE, + grpName = "group") + +study1[, a := ifelse(group %in% c(1:4),1,2)] + +study1[, b := ifelse(group %in% c(2,4,6,8),2,1)] +study1[, c := ifelse(group %in% c(3,4,7,8),2,1)] +# add treatment effect +study1[, y := ifelse(a == 1, y1, y1+1)] +# factorize group +study1[, a := as.factor(a)] +study1[, b := as.factor(b)] +study1[, c := as.factor(c)] +study1[, group := as.factor(group)] + +study1[,.('means' = mean(y),'sd' = sd(y)), by = .(group)] +a = suppressMessages({suppressWarnings({ + as.data.frame(afex::aov_car(y ~ cov + a * b * c + Error(1 / id), + factorize = FALSE, + include_aov = FALSE, + data = study1)$anova_table) + }) }) + +a2 = suppressMessages({suppressWarnings({ + as.data.frame(afex::aov_car( + y ~ a * b *c + Error(1 / id), + factorize = FALSE, + include_aov = FALSE, + data = study1 + )$anova_table) + }) }) +res_aov = append(a2["a",]$`Pr(>F)`, + res_aov) + +res_2x2x2a = append(a["a",]$`Pr(>F)`, + res_2x2x2a) +res_2x2x2b = append(a["b",]$`Pr(>F)`, + res_2x2x2b) +res_2x2x2ab = append(a["a:b",]$`Pr(>F)`, + res_2x2x2ab) +res_2x2x2abc = append(a["a:b:c",]$`Pr(>F)`, + res_2x2x2ab) + +} + +``` + +```{r} +des1 = ANOVA_design(design = "2b*2b*2b", + n = 10, + mu = c(1,1,1,1,0,0,0,0), sd = 2.5) +ex1 = ANOVA_exact2(des1, verbose=FALSE) + + +test_that("Simulation 1", { + # check ANOVA result + expect_equal(mean(res_aov < .05), ex1$main_results$power[2]/100, + tolerance = .02) + # Both should be about .81 + # Check ANCOVA result + expect_equal(mean(res_2x2x2a < .05), + ANCOVA_analytic(design = "2b*2b*2b", + mu = c(1,1,1,1,0,0,0,0), + n = 10, + sd = 2.5, + r2 = .33^2, + n_cov = 1, + alpha_level = .05, + beta_level = NULL)$main_results$power[1]/100, + tolerance = .01) + + expect_equal(mean(res_2x2x2b < .05), + ANCOVA_analytic(design = "2b*2b*2b", + mu = c(1,1,1,1,0,0,0,0), + n = 10, + sd = 2.5, + r2 = .33^2, + n_cov = 1, + alpha_level = .05, + beta_level = NULL)$main_results$power[2]/100, + tolerance = .01) + + expect_equal(mean(res_2x2x2ab < .05), + ANCOVA_analytic(design = "2b*2b*2b", + mu = c(1,1,1,1,0,0,0,0), + n = 10, + sd = 2.5, + r2 = .33^2, + n_cov = 1, + alpha_level = .05, + beta_level = NULL)$main_results$power[7]/100, + tolerance = .01) + +}) + +``` + diff --git a/Validation/ANCOVA_twoway_sims.Rmd b/Validation/ANCOVA_twoway_sims.Rmd new file mode 100644 index 0000000..862b30c --- /dev/null +++ b/Validation/ANCOVA_twoway_sims.Rmd @@ -0,0 +1,373 @@ +--- +title: "Factorial ANCOVA Validation: Two Way Design" +output: github_document +--- + +```{r setup, include=TRUE} +knitr::opts_chunk$set(echo = TRUE) +library(simstudy) +library(data.table) +library(Superpower) +library(car) +library(testthat) +set.seed(19256) +``` + +# Simulation 1 + +```{r } + +nsim = 5000 +res_2x2a = vector() +res_2x2b = vector() +res_2x2ab = vector() +res_aov = vector() + +for(i in 1:nsim){ + # Generate outcome and covariate matrix + dt <- genCorData( + 20 * 4, + mu = c(0, 0), + sigma = 2.5, + rho = .1, + corstr = "cs", + cnames = c("cov", "y1") + ) + +# generate random treatment assignment +study1 = trtAssign(dt, + n = 4, + balanced = TRUE, + grpName = "group") + +study1[, a := ifelse(group == 1 | group == 2,1,2)] + +study1[, b := ifelse(group == 1 | group == 3,1,2)] +# add treatment effect +study1[, y := ifelse((group == 1 | group == 3), y1+1, y1)] +# factorize group +study1[, a := as.factor(a)] +study1[, b := as.factor(b)] +study1[, group := as.factor(group)] + +study1[,.('means' = mean(y),'sd' = sd(y)), by = .(group)] +a = as.data.frame(anova(lm(y ~ cov + a*b, data = study1))) +a2 = as.data.frame(anova(lm(y ~ a*b, data = study1))) +res_aov = append(a2["b",]$`Pr(>F)`, + res_aov) + +res_2x2a = append(a["a",]$`Pr(>F)`, + res_2x2a) +res_2x2b = append(a["b",]$`Pr(>F)`, + res_2x2b) +res_2x2ab = append(a["a:b",]$`Pr(>F)`, + res_2x2ab) + +} + +mean(res_2x2a < .05) +mean(res_2x2b < .05) +mean(res_2x2ab < .05) + +``` + +```{r} +des1 = ANOVA_design(design = "2b*2b", n = 20, mu = c(1, 0, 1, 0), sd = 2.5) +ex1 = ANOVA_exact2(des1, verbose=FALSE) + + +test_that("Simulation 1", { + # check ANOVA result + expect_equal(mean(res_aov < .05), ex1$main_results$power[2]/100, + tolerance = .02) + # Both should be about .81 + # Check ANCOVA result + expect_equal(mean(res_2x2a < .05), + ANCOVA_analytic(design = "2b*2b", + mu = c(1,0,1,0), + n = 20, + sd = 2.5, + r2 = .1^2, + n_cov = 1, + alpha_level = .05, + beta_level = NULL)$main_results$power[1]/100, + tolerance = .01) + + expect_equal(mean(res_2x2b < .05), + ANCOVA_analytic(design = "2b*2b", + mu = c(1,0,1,0), + n = 20, + sd = 2.5, + r2 = .1^2, + n_cov = 1, + alpha_level = .05, + beta_level = NULL)$main_results$power[2]/100, + tolerance = .01) + +}) + +``` + +# Simulation 2 + +```{r } + +res_2x3a = vector() +res_2x3b = vector() +res_2x3ab = vector() +res_aov = vector() + +for(i in 1:nsim){ + # Generate outcome and covariate matrix + dt <- genCorData( + 33 * 6, + mu = c(0, 0), + sigma = 1.5, + rho = .4, + corstr = "cs", + cnames = c("cov", "y1") + ) + +# generate random treatment assignment +study1 = trtAssign(dt, + n = 6, + balanced = TRUE, + grpName = "group") + +study1[, a := ifelse(group == 1 | group == 2 | group == 3, 1, 2)] + +study1[, b := ifelse(group == 1 | group == 4,1, + ifelse(group == 2 | group == 5, + 2, 3))] +# add treatment effect +study1[, y := ifelse(group == 2, y1+1, + ifelse(group == 3, y1+2, + ifelse(group == 5, y1+1.5, + ifelse(group == 6, y1+3,y1))))] +# factorize group +study1[, a := as.factor(a)] +study1[, b := as.factor(b)] +study1[, group := as.factor(group)] + +#library(ggplot2) +#ggplot(study1, aes(x=b,y=y,group=group)) + geom_boxplot() + facet_wrap(~a) + +#a = as.data.frame(anova(lm(y ~ cov + a*b, data = study1))) + +a = suppressMessages({suppressWarnings({ + as.data.frame(afex::aov_car(y ~ cov + a * b + Error(1 / id), + factorize = FALSE, + include_aov = FALSE, + data = study1)$anova_table) + }) }) + +a2 = suppressMessages({suppressWarnings({ + as.data.frame(afex::aov_car( + y ~ a * b + Error(1 / id), + factorize = FALSE, + include_aov = FALSE, + data = study1 + )$anova_table) + }) }) + +res_aov = append(a2["a",]$`Pr(>F)`, + res_aov) + +res_2x3a = append(a["a",]$`Pr(>F)`, + res_2x3a) +res_2x3b = append(a["b",]$`Pr(>F)`, + res_2x3b) +res_2x3ab = append(a["a:b",]$`Pr(>F)`, + res_2x3ab) + +} + +mean(res_2x3a < .05) +mean(res_2x3b < .05) +mean(res_2x3ab < .05) + +``` + +```{r} + +des1 = ANOVA_design(design = "2b*3b", n = 33, + mu = c(0,1,2,0,1.5,3), + sd = 1.5) +ex1 = ANOVA_exact2(des1, verbose=FALSE) + + +test_that("Simulation 2", { + + expect_equal(mean(res_aov < .05), + ex1$main_results$power[1]/100, + tolerance = .015) + + # Both should be about .81 + # Check ANCOVA result + expect_equal(mean(res_2x3a < .05), + ANCOVA_analytic(design = "2b*3b", + mu = c(0,1,2,0,1.5,3), + n = 33, + sd = 1.5, + r2 = .4^2, + n_cov = 1, + alpha_level = .05, + beta_level = NULL)$main_results$power[1]/100, + tolerance = .01) + + expect_equal(mean(res_2x3b < .05), + ANCOVA_analytic(design = "2b*3b", + mu = c(0,1,2,0,1.5,3), + n = 33, + sd = 1.5, + r2 = .4^2, + n_cov = 1, + alpha_level = .05, + beta_level = NULL)$main_results$power[2]/100, + tolerance = .01) + + expect_equal(mean(res_2x3ab < .05), + ANCOVA_analytic(design = "2b*3b", + mu = c(0,1,2,0,1.5,3), + n = 33, + sd = 1.5, + r2 = .4^2, + n_cov = 1, + alpha_level = .05, + beta_level = NULL)$main_results$power[3]/100, + tolerance = .01) + +}) +``` + +# Simulation 3 + +```{r } + +res_3x3a = vector() +res_3x3b = vector() +res_3x3ab = vector() +res_aov = vector() + +for(i in 1:nsim){ + # Generate outcome and covariate matrix + dt <- genCorData( + 17 * 9, + mu = c(0, 0), + sigma = 1.78, + rho = .2, + corstr = "cs", + cnames = c("cov", "y1") + ) + +# generate random treatment assignment +study1 = trtAssign(dt, + n = 9, + balanced = TRUE, + grpName = "group") + +study1[, a := ifelse(group == 1 | + group == 2 | + group == 3, 1, + ifelse(group == 4 | + group == 5 | + group == 6, 2,3))] + +study1[, b := ifelse(group == 1 | + group == 4 | + group == 7, 1, + ifelse(group == 2 | + group == 5 | + group == 8, 2,3))] +# add treatment effect +study1[, y := ifelse(a == 1, y1+1, + ifelse(a == 2, y1+1.5, y1+2))] +# factorize group +study1[, a := as.factor(a)] +study1[, b := as.factor(b)] +study1[, group := as.factor(group)] + +a = suppressMessages({suppressWarnings({ + as.data.frame(afex::aov_car(y ~ cov + a * b + Error(1 / id), + factorize = FALSE, + include_aov = FALSE, + data = study1)$anova_table) + }) }) + +a2 = suppressMessages({suppressWarnings({ + as.data.frame(afex::aov_car( + y ~ a * b + Error(1 / id), + factorize = FALSE, + include_aov = FALSE, + data = study1 + )$anova_table) + }) }) + +res_aov = append(a2["a",]$`Pr(>F)`, + res_aov) + +res_3x3a = append(a["a",]$`Pr(>F)`, + res_3x3a) +res_3x3b = append(a["b",]$`Pr(>F)`, + res_3x3b) +res_3x3ab = append(a["a:b",]$`Pr(>F)`, + res_3x3ab) + +} + +mean(res_3x3a < .05) +mean(res_3x3b < .05) +mean(res_3x3ab < .05) + +``` + +```{r} + +test_that("Simulation 3", { + + expect_equal(mean(res_aov < .05), + .72, + tolerance = .015) + + # Both should be about .81 + # Check ANCOVA result + expect_equal(mean(res_3x3a < .05), + ANCOVA_analytic( + design = "3b*3b", + mu = c(1, 1, 1, 1.5, 1.5, 1.5, 2, 2, 2), + n = 17, + sd = 1.78, + r2 = .2 ^ 2, + n_cov = 1, + alpha_level = .05, + beta_level = NULL)$main_results$power[1]/100, + tolerance = .015) + + expect_equal(mean(res_3x3b < .05), + ANCOVA_analytic( + design = "3b*3b", + mu = c(1, 1, 1, 1.5, 1.5, 1.5, 2, 2, 2), + n = 17, + sd = 1.78, + r2 = .2 ^ 2, + n_cov = 1, + alpha_level = .05, + beta_level = NULL + )$main_results$power[2]/100, + tolerance = .01) + + expect_equal(mean(res_3x3ab < .05), + ANCOVA_analytic( + design = "3b*3b", + mu = c(1, 1, 1, 1.5, 1.5, 1.5, 2, 2, 2), + n = 17, + sd = 1.78, + r2 = .2 ^ 2, + n_cov = 1, + alpha_level = .05, + beta_level = NULL)$main_results$power[3]/100, + tolerance = .01) + +}) +``` + diff --git a/Validation/ANCOVA_twoway_sims.md b/Validation/ANCOVA_twoway_sims.md new file mode 100644 index 0000000..366a514 --- /dev/null +++ b/Validation/ANCOVA_twoway_sims.md @@ -0,0 +1,380 @@ +Factorial ANCOVA Validation: Two Way Design +================ + +``` r +knitr::opts_chunk$set(echo = TRUE) +library(simstudy) +library(data.table) +library(Superpower) +``` + + ## Registered S3 methods overwritten by 'lme4': + ## method from + ## cooks.distance.influence.merMod car + ## influence.merMod car + ## dfbeta.influence.merMod car + ## dfbetas.influence.merMod car + +``` r +library(car) +``` + + ## Loading required package: carData + +``` r +library(testthat) +set.seed(19256) +``` + +# Simulation 1 + +``` r +nsim = 5000 +res_2x2a = vector() +res_2x2b = vector() +res_2x2ab = vector() +res_aov = vector() + +for(i in 1:nsim){ + # Generate outcome and covariate matrix + dt <- genCorData( + 20 * 4, + mu = c(0, 0), + sigma = 2.5, + rho = .1, + corstr = "cs", + cnames = c("cov", "y1") + ) + +# generate random treatment assignment +study1 = trtAssign(dt, + n = 4, + balanced = TRUE, + grpName = "group") + +study1[, a := ifelse(group == 1 | group == 2,1,2)] + +study1[, b := ifelse(group == 1 | group == 3,1,2)] +# add treatment effect +study1[, y := ifelse((group == 1 | group == 3), y1+1, y1)] +# factorize group +study1[, a := as.factor(a)] +study1[, b := as.factor(b)] +study1[, group := as.factor(group)] + +study1[,.('means' = mean(y),'sd' = sd(y)), by = .(group)] +a = as.data.frame(anova(lm(y ~ cov + a*b, data = study1))) +a2 = as.data.frame(anova(lm(y ~ a*b, data = study1))) +res_aov = append(a2["b",]$`Pr(>F)`, + res_aov) + +res_2x2a = append(a["a",]$`Pr(>F)`, + res_2x2a) +res_2x2b = append(a["b",]$`Pr(>F)`, + res_2x2b) +res_2x2ab = append(a["a:b",]$`Pr(>F)`, + res_2x2ab) + +} +``` + +``` r +des1 = ANOVA_design(design = "2b*2b", n = 20, mu = c(1, 0, 1, 0), sd = 2.5) +``` + +![](ANCOVA_twoway_sims_files/figure-gfm/unnamed-chunk-2-1.png) + +``` r +ex1 = ANOVA_exact2(des1, verbose=FALSE) + + +test_that("Simulation 1", { + # check ANOVA result + expect_equal(mean(res_aov < .05), ex1$main_results$power[2]/100, + tolerance = .02) + # Both should be about .81 + # Check ANCOVA result + expect_equal(mean(res_2x2a < .05), + ANCOVA_analytic(design = "2b*2b", + mu = c(1,0,1,0), + n = 20, + sd = 2.5, + r2 = .1^2, + n_cov = 1, + alpha_level = .05, + beta_level = NULL)$main_results$power[1]/100, + tolerance = .01) + + expect_equal(mean(res_2x2b < .05), + ANCOVA_analytic(design = "2b*2b", + mu = c(1,0,1,0), + n = 20, + sd = 2.5, + r2 = .1^2, + n_cov = 1, + alpha_level = .05, + beta_level = NULL)$main_results$power[2]/100, + tolerance = .01) + +}) +``` + + ## Test passed 🥳 + +# Simulation 2 + +``` r +res_2x3a = vector() +res_2x3b = vector() +res_2x3ab = vector() +res_aov = vector() + +for(i in 1:nsim){ + # Generate outcome and covariate matrix + dt <- genCorData( + 33 * 6, + mu = c(0, 0), + sigma = 1.5, + rho = .4, + corstr = "cs", + cnames = c("cov", "y1") + ) + +# generate random treatment assignment +study1 = trtAssign(dt, + n = 6, + balanced = TRUE, + grpName = "group") + +study1[, a := ifelse(group == 1 | group == 2 | group == 3, 1, 2)] + +study1[, b := ifelse(group == 1 | group == 4,1, + ifelse(group == 2 | group == 5, + 2, 3))] +# add treatment effect +study1[, y := ifelse(group == 2, y1+1, + ifelse(group == 3, y1+2, + ifelse(group == 5, y1+1.5, + ifelse(group == 6, y1+3,y1))))] +# factorize group +study1[, a := as.factor(a)] +study1[, b := as.factor(b)] +study1[, group := as.factor(group)] + +#library(ggplot2) +#ggplot(study1, aes(x=b,y=y,group=group)) + geom_boxplot() + facet_wrap(~a) + +#a = as.data.frame(anova(lm(y ~ cov + a*b, data = study1))) + +a = suppressMessages({suppressWarnings({ + as.data.frame(afex::aov_car(y ~ cov + a * b + Error(1 / id), + factorize = FALSE, + include_aov = FALSE, + data = study1)$anova_table) + }) }) + +a2 = suppressMessages({suppressWarnings({ + as.data.frame(afex::aov_car( + y ~ a * b + Error(1 / id), + factorize = FALSE, + include_aov = FALSE, + data = study1 + )$anova_table) + }) }) + +res_aov = append(a2["a",]$`Pr(>F)`, + res_aov) + +res_2x3a = append(a["a",]$`Pr(>F)`, + res_2x3a) +res_2x3b = append(a["b",]$`Pr(>F)`, + res_2x3b) +res_2x3ab = append(a["a:b",]$`Pr(>F)`, + res_2x3ab) + +} +``` + +``` r +des1 = ANOVA_design(design = "2b*3b", n = 33, + mu = c(0,1,2,0,1.5,3), + sd = 1.5) +``` + +![](ANCOVA_twoway_sims_files/figure-gfm/unnamed-chunk-4-1.png) + +``` r +ex1 = ANOVA_exact2(des1, verbose=FALSE) + + +test_that("Simulation 2", { + + expect_equal(mean(res_aov < .05), + ex1$main_results$power[1]/100, + tolerance = .015) + + # Both should be about .81 + # Check ANCOVA result + expect_equal(mean(res_2x3a < .05), + ANCOVA_analytic(design = "2b*3b", + mu = c(0,1,2,0,1.5,3), + n = 33, + sd = 1.5, + r2 = .4^2, + n_cov = 1, + alpha_level = .05, + beta_level = NULL)$main_results$power[1]/100, + tolerance = .01) + + expect_equal(mean(res_2x3b < .05), + ANCOVA_analytic(design = "2b*3b", + mu = c(0,1,2,0,1.5,3), + n = 33, + sd = 1.5, + r2 = .4^2, + n_cov = 1, + alpha_level = .05, + beta_level = NULL)$main_results$power[2]/100, + tolerance = .01) + + expect_equal(mean(res_2x3ab < .05), + ANCOVA_analytic(design = "2b*3b", + mu = c(0,1,2,0,1.5,3), + n = 33, + sd = 1.5, + r2 = .4^2, + n_cov = 1, + alpha_level = .05, + beta_level = NULL)$main_results$power[3]/100, + tolerance = .01) + +}) +``` + + ## Test passed 🥇 + +# Simulation 3 + +``` r +res_3x3a = vector() +res_3x3b = vector() +res_3x3ab = vector() +res_aov = vector() + +for(i in 1:nsim){ + # Generate outcome and covariate matrix + dt <- genCorData( + 17 * 9, + mu = c(0, 0), + sigma = 1.78, + rho = .2, + corstr = "cs", + cnames = c("cov", "y1") + ) + +# generate random treatment assignment +study1 = trtAssign(dt, + n = 9, + balanced = TRUE, + grpName = "group") + +study1[, a := ifelse(group == 1 | + group == 2 | + group == 3, 1, + ifelse(group == 4 | + group == 5 | + group == 6, 2,3))] + +study1[, b := ifelse(group == 1 | + group == 4 | + group == 7, 1, + ifelse(group == 2 | + group == 5 | + group == 8, 2,3))] +# add treatment effect +study1[, y := ifelse(a == 1, y1+1, + ifelse(a == 2, y1+1.5, y1+2))] +# factorize group +study1[, a := as.factor(a)] +study1[, b := as.factor(b)] +study1[, group := as.factor(group)] + +a = suppressMessages({suppressWarnings({ + as.data.frame(afex::aov_car(y ~ cov + a * b + Error(1 / id), + factorize = FALSE, + include_aov = FALSE, + data = study1)$anova_table) + }) }) + +a2 = suppressMessages({suppressWarnings({ + as.data.frame(afex::aov_car( + y ~ a * b + Error(1 / id), + factorize = FALSE, + include_aov = FALSE, + data = study1 + )$anova_table) + }) }) + +res_aov = append(a2["a",]$`Pr(>F)`, + res_aov) + +res_3x3a = append(a["a",]$`Pr(>F)`, + res_3x3a) +res_3x3b = append(a["b",]$`Pr(>F)`, + res_3x3b) +res_3x3ab = append(a["a:b",]$`Pr(>F)`, + res_3x3ab) + +} +``` + +``` r +test_that("Simulation 3", { + + expect_equal(mean(res_aov < .05), + .72, + tolerance = .015) + + # Both should be about .81 + # Check ANCOVA result + expect_equal(mean(res_3x3a < .05), + ANCOVA_analytic( + design = "3b*3b", + mu = c(1, 1, 1, 1.5, 1.5, 1.5, 2, 2, 2), + n = 17, + sd = 1.78, + r2 = .2 ^ 2, + n_cov = 1, + alpha_level = .05, + beta_level = NULL)$main_results$power[1]/100, + tolerance = .015) + + expect_equal(mean(res_3x3b < .05), + ANCOVA_analytic( + design = "3b*3b", + mu = c(1, 1, 1, 1.5, 1.5, 1.5, 2, 2, 2), + n = 17, + sd = 1.78, + r2 = .2 ^ 2, + n_cov = 1, + alpha_level = .05, + beta_level = NULL + )$main_results$power[2]/100, + tolerance = .01) + + expect_equal(mean(res_3x3ab < .05), + ANCOVA_analytic( + design = "3b*3b", + mu = c(1, 1, 1, 1.5, 1.5, 1.5, 2, 2, 2), + n = 17, + sd = 1.78, + r2 = .2 ^ 2, + n_cov = 1, + alpha_level = .05, + beta_level = NULL)$main_results$power[3]/100, + tolerance = .01) + +}) +``` + + ## Test passed 😀 diff --git a/Validation/ANCOVA_twoway_sims_files/figure-gfm/unnamed-chunk-2-1.png b/Validation/ANCOVA_twoway_sims_files/figure-gfm/unnamed-chunk-2-1.png new file mode 100644 index 0000000..5e55f21 Binary files /dev/null and b/Validation/ANCOVA_twoway_sims_files/figure-gfm/unnamed-chunk-2-1.png differ diff --git a/Validation/ANCOVA_twoway_sims_files/figure-gfm/unnamed-chunk-4-1.png b/Validation/ANCOVA_twoway_sims_files/figure-gfm/unnamed-chunk-4-1.png new file mode 100644 index 0000000..457eac4 Binary files /dev/null and b/Validation/ANCOVA_twoway_sims_files/figure-gfm/unnamed-chunk-4-1.png differ diff --git a/docs/404.html b/docs/404.html index 1e65db5..9c17e80 100644 --- a/docs/404.html +++ b/docs/404.html @@ -71,7 +71,7 @@ Superpower - 0.1.2 + 0.2.0 @@ -93,6 +93,9 @@