### Cohen's d test
library(effsize)
try_with_time_limit <- function(expr, cpu = Inf, elapsed = Inf){
y <- try({setTimeLimit(cpu, elapsed); expr}, silent = TRUE)
if(inherits(y, "try-error")) stop("Operation timed out") else y
}
generate_data <- function(n,m,stdev){
x <- rnorm(n,m,stdev)
sd.adj = stdev/sd(x)
x <- x * sd.adj
m.adj = m - mean(x)
x <- x + m.adj
return(x)
}
context("Cohen d")
test_that("Two samples with large difference", {
set.seed(52)
x = rnorm(100,mean=10)
y = rnorm(100,mean=12)
d = (c(x,y))
f = rep(c("A","B"),each=100)
eff.d = cohen.d(d,f)
expect_lt(eff.d$conf.int[1], -2)
expect_gt(eff.d$conf.int[2], -2)
})
test_that("Two samples with large difference with Hedges G", {
set.seed(52)
x = rnorm(100,mean=10)
y = rnorm(100,mean=12)
d = (c(x,y))
f = rep(c("A","B"),each=100)
eff.g = cohen.d(d,f,hedges.correction=TRUE)
expect_lt(eff.g$conf.int[1], -2)
expect_gt(eff.g$conf.int[2], -2)
})
test_that("Two independent samples with Hedges G, example from Borenstein et al, Introduction to Meta-Analysis", {
set.seed(321)
x1 = generate_data(50,103,5.5)
x2 = generate_data(50,100,4.5)
eff.g = cohen.d(x1,x2,hedges.correction=TRUE)
expect_equal(eff.g$estimate, 0.5924, tolerance=.0001)
})
test_that("Two samples from same population", {
set.seed(54)
d <- rnorm(200)
f <- rep(c(1,2),100)
expect_warning( eff.d <- cohen.d(d ~ f) )
expect_lt(eff.d$conf.int[1], 0)
expect_gt(eff.d$conf.int[2], 0)
})
# test_that("Paired measures", {
# delta = c(1.73, 1.06, 2.03, 1.40, 0.95, 1.13, 1.41, 1.73, 1.63, 1.56) - 1
# set.seed(50)
# a = delta
# set.seed(50)
# b = delta + runif(10)
#
# eff.d = cohen.d(a, b, paired=TRUE)
# expect_equal(abs(eff.d$estimate),1.42, tolerance=.01)
# })
test_that("Paired measures w/NA", {
set.seed(60)
x = rnorm(10,mean=10)
y = rnorm(10,mean=12)
x[4] <- NA
eff.d = cohen.d(x,y,paired=TRUE,na.rm=TRUE)
expect_equal(abs(eff.d$estimate),2.8, tolerance=.01)
})
test_that("Paired measures with multiple NA", {
eff.d = cohen.d(c(0.2,0.4,0.8,0.9,NA),
c(0.12,0.14,0.18,0.119,NA),
paired = T, na.rm = T)
expect_equal(abs(eff.d$estimate),1.55, tolerance=.01)
})
test_that("Non centrality parameter",{
# From Confidence Intervals on Effect Size, David C. Howell
# Data from: Adams, Wright, and Lohr (1996)
set.seed(22)
Homophobic = generate_data(35,24,sqrt(148.87))
Nonhomophobic = generate_data(29,16.5,sqrt(139.16))
eff.d = cohen.d(Homophobic, Nonhomophobic, noncentral = TRUE)
expect_equal(as.numeric(eff.d$estimate[1]),0.62,tolerance=0.01)
expect_equal(as.numeric(eff.d$conf.int[1]),0.117,tolerance=0.001)
expect_equal(as.numeric(eff.d$conf.int[2]),1.125,tolerance=0.001)
})
test_that("Two samples with large negative difference and noncentral", {
set.seed(52)
x = rnorm(100,mean=10)
y = rnorm(100,mean=12)
d = (c(x,y))
f = rep(c("A","B"),each=100)
eff.d = cohen.d(d,f,noncentral = TRUE)
expect_lt(eff.d$conf.int[1], -2)
expect_gt(eff.d$conf.int[2], -2)
})
# issue #23 noncentral error (infinite loop)
test_that("Two samples with large negative difference and noncentral", {
a = c(9.81605624621576, 8.93891560898168, 9.05436620537713, 6.01771305071382,
8.98575772172126, 7.2595761413452, 13.7537304479165, 9.89861975363729,
7.37246838418731, 8.33888742066376)
b = c(8.62798051757027, 5.25390514295654, 12.0869496333763, 7.20553319908951,
6.99017188584104, 8.22363490828756, 8.05485287987422, 6.50534569228974,
8.68871033652186, 4.62897740315003)
res = try_with_time_limit(cohen.d(a,b,hedges.correction = FALSE, noncentral =TRUE),1)
expect_equal(as.numeric(res$estimate),0.63236,tolerance = .0001)
})
#issue #27 Cohen.d gives wrong value when data is not arranged by f
test_that("Order of factor values does not affect result",{
set.seed(25)
d.data <- data.frame(group = factor(sample(c(1,2), 20, replace = T, prob = c(.33,.67))),
value = rnorm(60,100,15))
group1 <- d.data$value[d.data$group==1]
group2 <- d.data$value[d.data$group==2]
d.data.arranged <- d.data[order(d.data$group),]
r1 <- cohen.d(d.data$value ~ d.data$group)
r2 <- cohen.d(d.data.arranged$value ~ d.data.arranged$group)
expect_equal(as.numeric(r1$estimate),as.numeric(r2$estimate))
})
#spin-off of issue #27 Cohen.d gives wrong value when data is not arranged by f
test_that("When inverting control and treatment effsize just change sign",{
set.seed(7)
group1 <- rnorm(20,100,15)
group2 <- rnorm(40,100,15)
r1 <- cohen.d(group1,group2)
r2 <- cohen.d(group2,group1)
expect_equal(r1$estimate,-r2$estimate)
})
# issue #28 confidence interval with non-central distribution for paired data
test_that("confidence interval with non-central distribution for paired data",{
# data from https://www.uvm.edu/%7Edhowell/methods7/Supplements/Confidence%20Intervals%20on%20Effect%20Size.pdf
moon.data = c(1.73, 1.06, 2.03, 1.40, 0.95, 1.13, 1.41, 1.73, 1.63, 1.56)
g1 = rep(1,length(moon.data))
# x2 = seq(1,1.9,by=0.1)
# x1 = moon.data+x2
res = cohen.d(moon.data,g1,
paired = TRUE,
noncentral = TRUE,
conf.level = 0.95
)
expect_equal(as.numeric(res$conf.int[1]),0.4907785,tolerance = .000001)
expect_equal(as.numeric(res$conf.int[2]),2.3358769,tolerance = .000001)
})
# issue #28 confidence interval with non-central distribution for two samples
test_that("confidence interval with non-central distribution for two samples",{
# data from https://www.uvm.edu/%7Edhowell/methods7/Supplements/Confidence%20Intervals%20on%20Effect%20Size.pdf
set.seed(537)
hf = generate_data(35,24,sqrt(148.87))
nhf = generate_data(29,16.5,sqrt(139.16))
res= cohen.d(hf,nhf,noncentral=TRUE)
expect_equal(as.numeric(res$conf.int[1]),0.117,tolerance = .001)
expect_equal(as.numeric(res$conf.int[2]),1.126,tolerance = .001)
})
# issue #?? unpaired noncentral error (infinite loop)
test_that("Two samples paired formula paired and noncentral", {
data(sleep)
res = try_with_time_limit(cohen.d(extra ~ group,
data = sleep,
noncentral =TRUE),1)
expect_equal(as.numeric(res$estimate),-0.8321,tolerance = .0001)
})
# issue #32 paired Hedges correction
# Using the example from
# Michael Borenstein, L. V. Hedges, J. P. T. Higgins and H. R. Rothstein
# Introduction to Meta-Analysis.
test_that("Two samples paired pooled with Hedges correction", {
set.seed(7828)
x1 = generate_data(50,103,5.5)
d = generate_data(50,3,5.5)
x2 = x1 - d
res.d = cohen.d(x1, x2, paired=TRUE)
res.g = cohen.d(x1, x2, paired=TRUE, hedges.correction = TRUE)
expect_equal(as.numeric(res.d$estimate),0.4225,tolerance = .0001)
expect_equal(as.numeric(res.d$var),0.0131,tolerance = .0001)
expect_equal(as.numeric(res.g$estimate),0.4160,tolerance = .0001)
expect_equal(as.numeric(res.g$var),0.0127,tolerance = .0001)
})
# from issue #34
test_that("Two samples paired normal cohen d", {
x1 = c(131,124,130,105,102,120,101,120)
x2 = c(133,134,129,112,106,125,106,129)
res.d = cohen.d(x2, x1, paired=TRUE)
expect_equal(as.numeric(res.d$estimate),0.422,tolerance = .001)
})
test_that("Single sample with null mu",{
moon.data = c(1.73, 1.06, 2.03, 1.40, 0.95, 1.13, 1.41, 1.73, 1.63, 1.56)
res.d = cohen.d( ~ moon.data,data=data.frame(moon.data=moon.data))
expect_equal(as.numeric(res.d$estimate),4.294,tolerance = .001)
})
test_that("Single sample with non-null mu",{
moon.data <- c(1.73, 1.06, 2.03, 1.40, 0.95, 1.13, 1.41, 1.73, 1.63, 1.56)
res.d = cohen.d( ~ moon.data, mu=1,data=data.frame(moon.data=moon.data))
expect_equal(as.numeric(res.d$estimate),1.359,tolerance = .001)
})
test_that("Cohen pooled false",{
Before = c(45L, 52L, 63L, 68L, 57L, 55L, 60L, 59L)
After = c(49L, 50L, 70L, 71L, 53L, 61L, 62L, 67L)
res.d = cohen.d(After, Before, pooled=FALSE)
expect_equal(as.numeric(res.d$estimate),0.429,tolerance = .001)
})
test_that("Cohen numeric first arg with own class",{
d = c (2,3,4,5,6,7)
f = factor(rep(c("a","b"),each=3))
class(d) <- "Cheese"
res.d = cohen.d(d, f)
expect_equal(as.numeric(res.d$estimate),-3,tolerance = .1)
})
test_that("Same results independent on data's sorting order",{
## From Issue #48
#Setup some dummy data
set.seed(1234)
numSubjects <- 20
subject <- 1:numSubjects
before <- runif(numSubjects)
after <- before - 0.35*runif(length(before))
df <- data.frame(
subject=rep(subject,2),
group=c(
rep("before",numSubjects),
rep("after",numSubjects)
),
value=c(before, after)
)
##################################
valueBefore <- df[df$group == "before",]$value
valueAfter <- df[df$group == "after",]$value
# d1 == d2
d1 <- effsize::cohen.d(valueAfter, valueBefore, paired=T)
d2 <- expect_warning(effsize::cohen.d(value~group, data=df, paired=T))
expect_equal(as.numeric(d1$estimate),as.numeric(d2$estimate),tolerance = .001)
##################################
# Sorted by subjectId
df <- df[order(df$subject),]
valueBefore <- df[df$group == "before",]$value
valueAfter <- df[df$group == "after",]$value
d3 <- effsize::cohen.d(valueAfter, valueBefore, paired=T)
d4 <- expect_warning(effsize::cohen.d(value~group, data=df, paired=T))
expect_equal(as.numeric(d3$estimate),as.numeric(d4$estimate),tolerance = .001)
##################################
# Sorted by group
df <- df[order(df$group),]
valueBefore <- df[df$group == "before",]$value
valueAfter <- df[df$group == "after",]$value
d5 <- effsize::cohen.d(valueAfter, valueBefore, paired=T)
d6 <- expect_warning(effsize::cohen.d(value~group, data=df, paired=T))
expect_equal(as.numeric(d5$estimate),as.numeric(d6$estimate),tolerance = .001)
##################################
d7 <- expect_silent(effsize::cohen.d(value~group | Subject(subject), data=df, paired=T))
expect_equal(as.numeric(d6$estimate),as.numeric(d7$estimate),tolerance = .001)
expect_error(effsize::cohen.d(value~group | Subject(subject), data=df[c(-1),], paired=T))
})
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.