diff --git a/DESCRIPTION b/DESCRIPTION index 99df558..8218208 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: rigr Title: Regression, Inference, and General Data Analysis Tools in R -Version: 1.0.7 +Version: 1.0.8 Authors@R: c(person("Amy D", "Willis", email = "adwillis@uw.edu", role = c("aut", "cre"), comment = c(ORCID = "0000-0002-2802-4317")), person("Taylor", "Okonek", role = "aut"), person("Charles J", "Wolock", role = "aut"), diff --git a/R/proptest.R b/R/proptest.R index c8b2e00..ebe335a 100644 --- a/R/proptest.R +++ b/R/proptest.R @@ -187,10 +187,19 @@ proptest <- function (var1, var2 = NULL, by = NULL, exact = FALSE, binom <- stats::binom.test(x = sum(var1), n = length(var1), p = null.hypoth, alternative = alternative, conf.level = conf.level) se <- sqrt(est1*(1-est1)/length(var1)) - cil <- as.numeric(format(min(binom$conf.int), - digits = digits)) - cih <- as.numeric(format(max(binom$conf.int), - digits = digits)) + if (alternative == "two.sided") { + cil <- as.numeric(format(min(binom$conf.int), + digits = digits)) + cih <- as.numeric(format(max(binom$conf.int), + digits = digits)) + } else { + ci_test <- stats::binom.test(x = sum(var1), n = length(var1), p = null.hypoth, + alternative = "two.sided", conf.level = conf.level) + cil <- as.numeric(format(min(ci_test$conf.int), + digits = digits)) + cih <- as.numeric(format(max(ci_test$conf.int), + digits = digits)) + } se <- as.numeric(format(se, digits = digits)) pval <- as.numeric(format(binom$p.value, digits = digits)) diff --git a/R/proptesti.R b/R/proptesti.R index cffe917..01a6f9a 100644 --- a/R/proptesti.R +++ b/R/proptesti.R @@ -114,10 +114,18 @@ proptesti <- function(x1, n1, x2 = NULL, n2 = NULL, exact = FALSE, zstat <- NULL pval <- as.numeric(format(test$p.value, digits = digits)) - cil <- as.numeric(format(min(test$conf.int), - digits = digits)) - cih <- as.numeric(format(max(test$conf.int), - digits = digits)) + if (alternative == "two.sided") { + cil <- as.numeric(format(min(test$conf.int), + digits = digits)) + cih <- as.numeric(format(max(test$conf.int), + digits = digits)) + } else { + ci_test <- stats::binom.test(x1, n1, p = null.hypoth, alternative = "two.sided", conf.level = conf.level) + cil <- as.numeric(format(min(ci_test$conf.int), + digits = digits)) + cih <- as.numeric(format(max(ci_test$conf.int), + digits = digits)) + } } else { test <- stats::prop.test(x1, n1, p = null.hypoth, alternative = alternative, conf.level = conf.level, correct = correct) @@ -129,13 +137,14 @@ proptesti <- function(x1, n1, x2 = NULL, n2 = NULL, exact = FALSE, digits = digits)) cih <- as.numeric(format(est1 + stats::qnorm(cl)*se1, digits = digits)) - if (alternative == "less") { - cil <- as.numeric(format(min(test$conf.int), - digits = digits)) - } else if (alternative == "greater") { - cih <- as.numeric(format(max(test$conf.int), - digits = digits)) - } + # as of issue 166, setting all CI's to be two sided even if test is one sided + # if (alternative == "less") { + # cil <- as.numeric(format(min(test$conf.int), + # digits = digits)) + # } else if (alternative == "greater") { + # cih <- as.numeric(format(max(test$conf.int), + # digits = digits)) + # } } est1 <- as.numeric(format(est1, digits = digits)) se1 <- as.numeric(format(se1, digits = digits)) @@ -144,6 +153,13 @@ proptesti <- function(x1, n1, x2 = NULL, n2 = NULL, exact = FALSE, "]", sep = "")), ncol = 5) colnames(printMat) <- c("Variable", "Obs", "Mean", "Std. Error", paste0(conf.level*100, "% CI")) rownames(printMat) <- "" + # bound ci in [0, 1] + if (cil < 0) { + cil <- 0 + } + if (cih > 1) { + cih <- 1 + } } else { twosamp <- TRUE est <- c(x1/n1, x2/n2, x1/n1- x2/n2) @@ -176,6 +192,8 @@ proptesti <- function(x1, n1, x2 = NULL, n2 = NULL, exact = FALSE, names(printMat) <- c("Group", "Obs", "Mean", "Std. Err.", paste(conf.level * 100, "% CI", sep = "")) row.names(printMat) <- c("", " ", " ") + cil[1:2] <- c(ifelse(cil[1] < 0, 0, cil[1]), ifelse(cil[2] < 0, 0, cil[2])) + cih[1:2] <- c(ifelse(cih[1] > 1, 1, cih[1]), ifelse(cih[2] > 1, 1, cih[2])) } par <- c(null.hypoth = null.hypoth, alternative = alternative, conf.level = conf.level, exact = exact, twosamp = twosamp, diff --git a/tests/testthat/test_proptest.R b/tests/testthat/test_proptest.R index e3d70f8..917b11f 100644 --- a/tests/testthat/test_proptest.R +++ b/tests/testthat/test_proptest.R @@ -265,9 +265,13 @@ p2 <- binom.test(sum(a), length(a), alternative = "less") test_that("proptest() returns correct numbers for one-sample exact test, left-sided", { expect_s3_class(p1, "proptest") expect_equal(p1$pval, p2$p.value, tolerance = 1e-2) # p-value - expect_equal(as.numeric(strsplit(substr(p1$tab[[6]], start = 2, stop = nchar(p1$tab[[6]])-1), ", ")[[1]]), - p2$conf.int[1:2], - tolerance = 1e-2) # conf int + # as of issue #166, make all CIs two sided + # expect_equal(as.numeric(strsplit(substr(p1$tab[[6]], start = 2, stop = nchar(p1$tab[[6]])-1), ", ")[[1]]), + # p2$conf.int[1:2], + # tolerance = 1e-2) # conf int + expect_false(all.equal(as.numeric(strsplit(substr(p1$tab[[6]], start = 2, stop = nchar(p1$tab[[6]])-1), ", ")[[1]]), + p2$conf.int[1:2], + tolerance = 1e-2) == TRUE) expect_equal(p1$tab[[1]], "a") # var name expect_equal(as.numeric(p1$tab[[2]]), length(a)) # n obs expect_equal(as.numeric(p1$tab[[3]]), sum(is.na(a))) # NAs @@ -290,9 +294,14 @@ p2 <- binom.test(sum(a), length(a), alternative = "greater") test_that("proptest() returns correct numbers for one-sample exact test, right-sided", { expect_s3_class(p1, "proptest") expect_equal(p1$pval, p2$p.value, tolerance = 1e-2) # p-value - expect_equal(as.numeric(strsplit(substr(p1$tab[[6]], start = 2, stop = nchar(p1$tab[[6]])-1), ", ")[[1]]), - p2$conf.int[1:2], - tolerance = 1e-2) # conf int + expect_equal(p1$pval, p2$p.value, tolerance = 1e-2) # p-value + # as of issue #166, make all CIs two sided + # expect_equal(as.numeric(strsplit(substr(p1$tab[[6]], start = 2, stop = nchar(p1$tab[[6]])-1), ", ")[[1]]), + # p2$conf.int[1:2], + # tolerance = 1e-2) # conf int + expect_false(all.equal(as.numeric(strsplit(substr(p1$tab[[6]], start = 2, stop = nchar(p1$tab[[6]])-1), ", ")[[1]]), + p2$conf.int[1:2], + tolerance = 1e-2) == TRUE) expect_equal(p1$tab[[1]], "a") # var name expect_equal(as.numeric(p1$tab[[2]]), length(a)) # n obs expect_equal(as.numeric(p1$tab[[3]]), sum(is.na(a))) # NAs diff --git a/tests/testthat/test_proptesti.R b/tests/testthat/test_proptesti.R index f55b7d4..fc250bc 100644 --- a/tests/testthat/test_proptesti.R +++ b/tests/testthat/test_proptesti.R @@ -86,8 +86,8 @@ test_that("proptesti() returns correct numbers for one-sample test, left-sided", expect_equal(abs(p1$zstat), sqrt(p2$statistic[[1]]), tolerance = 1e-2) # test statistic expect_equal(p1$pval, p2$p.value, tolerance = 1e-2) # p-value expect_equal(as.numeric(strsplit(substr(p1$tab[[5]], start = 2, stop = nchar(p1$tab[[5]])-1), ", ")[[1]]), - c(0, - #p2$estimate[[1]] - 1.96*sqrt(p2$estimate[[1]]*(1-p2$estimate[[1]])/length(a)), + # as of issue 166 make all CIs two sided + c(p2$estimate[[1]] - 1.96*sqrt(p2$estimate[[1]]*(1-p2$estimate[[1]])/length(a)), p2$estimate[[1]] + 1.96*sqrt(p2$estimate[[1]]*(1-p2$estimate[[1]])/length(a))), tolerance = 1e-2) # conf int expect_equal(p1$tab[[1]], "var1") # var name @@ -112,9 +112,9 @@ test_that("proptesti() returns correct numbers for one-sample test, right-sided" expect_equal(abs(p1$zstat), sqrt(p2$statistic[[1]]), tolerance = 1e-2) # test statistic expect_equal(p1$pval, p2$p.value, tolerance = 1e-2) # p-value expect_equal(as.numeric(strsplit(substr(p1$tab[[5]], start = 2, stop = nchar(p1$tab[[5]])-1), ", ")[[1]]), + # as of issue 166 make all CIs two sided c(p2$estimate[[1]] - 1.96*sqrt(p2$estimate[[1]]*(1-p2$estimate[[1]])/length(a)), - 1), - #p2$estimate[[1]] + 1.96*sqrt(p2$estimate[[1]]*(1-p2$estimate[[1]])/length(a))), + p2$estimate[[1]] + 1.96*sqrt(p2$estimate[[1]]*(1-p2$estimate[[1]])/length(a))), tolerance = 1e-2) # conf int expect_equal(p1$tab[[1]], "var1") # var name expect_equal(as.numeric(p1$tab[[2]]), length(a)) # n obs @@ -183,9 +183,13 @@ p2 <- binom.test(sum(a), length(a), alternative = "less") test_that("proptesti() returns correct numbers for one-sample test, exact, left-sided", { expect_s3_class(p1, "proptesti") expect_equal(p1$pval, p2$p.value, tolerance = 1e-2) # p-value - expect_equal(as.numeric(strsplit(substr(p1$tab[[5]], start = 2, stop = nchar(p1$tab[[5]])-1), ", ")[[1]]), - p2$conf.int[1:2], - tolerance = 1e-2) # conf int + # as of issue #166, make all CIs two sided + # expect_equal(as.numeric(strsplit(substr(p1$tab[[5]], start = 2, stop = nchar(p1$tab[[5]])-1), ", ")[[1]]), + # p2$conf.int[1:2], + # tolerance = 1e-2) # conf int + expect_false(all.equal(as.numeric(strsplit(substr(p1$tab[[5]], start = 2, stop = nchar(p1$tab[[5]])-1), ", ")[[1]]), + p2$conf.int[1:2], + tolerance = 1e-2) == TRUE) expect_equal(p1$tab[[1]], "var1") # var name expect_equal(as.numeric(p1$tab[[2]]), length(a)) # n obs expect_equal(as.numeric(p1$tab[[3]]), p2$estimate[[1]], tolerance = 3) # estimate of mean @@ -205,9 +209,13 @@ p2 <- binom.test(sum(a), length(a), alternative = "greater") test_that("proptesti() returns correct numbers for one-sample test, exact, right-sided", { expect_s3_class(p1, "proptesti") expect_equal(p1$pval, p2$p.value, tolerance = 1e-2) # p-value - expect_equal(as.numeric(strsplit(substr(p1$tab[[5]], start = 2, stop = nchar(p1$tab[[5]])-1), ", ")[[1]]), - p2$conf.int[1:2], - tolerance = 1e-2) # conf int + # as of issue #166, make all CIs two sided + # expect_equal(as.numeric(strsplit(substr(p1$tab[[5]], start = 2, stop = nchar(p1$tab[[5]])-1), ", ")[[1]]), + # p2$conf.int[1:2], + # tolerance = 1e-2) # conf int + expect_false(all.equal(as.numeric(strsplit(substr(p1$tab[[5]], start = 2, stop = nchar(p1$tab[[5]])-1), ", ")[[1]]), + p2$conf.int[1:2], + tolerance = 1e-2) == TRUE) expect_equal(p1$tab[[1]], "var1") # var name expect_equal(as.numeric(p1$tab[[2]]), length(a)) # n obs expect_equal(as.numeric(p1$tab[[3]]), p2$estimate[[1]], tolerance = 3) # estimate of mean diff --git a/tests/testthat/test_regress.R b/tests/testthat/test_regress.R index 86c12dc..7c37097 100644 --- a/tests/testthat/test_regress.R +++ b/tests/testthat/test_regress.R @@ -367,7 +367,7 @@ test_that("regress() returns same output as coxph() for fnctl = hazard", { }) -mri2 <- read.table("https://rct-design.com/TeachingMaterials/Datasets/mri.txt", header = T) +mri2 <- mri mri2$obstime_yrs <- mri2$obstime/365.25 mri2$ldlcat <- cut(mri2$ldl, breaks=c(0, 70, 100, 130, 160, 190, 250), right=FALSE) mri2$surv <- Surv(mri2$obstime_yrs, mri2$death)