MethodOfThresholds <- function(G1, G2, #Category proportions, high to low.
epsilon = .00001, #To avoid unfittable values of 0 or 1
deltaS = 0, #Starting value of \Delta for iterations, equal means, see above
rhoS = 1, #Starting value of \rho for iterations, equal SDs, see above,
sum = T, #Whether the numbers are part of a range that sums to 1
rnd = 3,
DG = c("GDelta", "Cd")){ #Glass' Delta and SDR or Cohen's d and pooled SD. Hedge's *g* can also be added, and I might later, but it would require vectors of N's for sum = F, so I haven't attempted it yet
DG = ifelse(is.na(DG), "GDelta", DG)
colMax <- function(data) sapply(as.data.frame(data), max, na.rm = TRUE)
if (sum == T) {G1C = cumsum(G1); G2C = cumsum(G2)} else {G1C = G1; G2C = G2}
if (sum == T) {G1M = mapply('/', G1C, colMax(G1C)); G2M = mapply('/', G2C, colMax(G2C))} else {G1M = G1C; G2M = G2C}
GDF = rbind(as.numeric(unlist(G1M)), as.numeric(unlist(G2M)))
GDF[which(GDF == 1)] = GDF[which(GDF == 1)] - epsilon
GDF[which(GDF == 0)] = GDF[which(GDF == 0)] + epsilon
P_G1 = GDF[1,]
P_G2 = GDF[2,]
f0 = function(t){ #Standard normal density function
1/sqrt(2*pi)*exp(-t^2/2)}
thresholds = rep(0, length(P_G1))
for(i in 1:length(P_G1)){
F0_temp = function(t){
-P_G1[i] + integrate(f0, t, Inf)$value}
thresholds[i] = uniroot(F0_temp, lower = -10, upper = 10)$root}
F_fit = function(t, delta, rho){
f0_full = function(x, mu = delta, sigma = rho){
1/sqrt(2*pi*sigma^2)*exp(-(x-mu)^2/(2*sigma^2))}
outp = rep(0, length(t))
for(i in 1:length(outp)){
outp[i] = integrate(f0_full,t[i],Inf)$value}
return(outp)}
fitLS1 = nls(P_G2 ~ F_fit(thresholds, delta, rho), start = list(delta = deltaS, rho = rhoS))
deltaF = unname(coef(fitLS1)[1]) #Mean difference (in SD units)
rhoF = unname(coef(fitLS1)[2]) #Standard deviation ratio
out = c(deltaF, rhoF)
thresholds = rep(0, length(P_G1))
for(i in 1:length(P_G1)){
F0_temp = function(t){
-P_G2[i] + integrate(f0, t, Inf)$value}
thresholds[i] = uniroot(F0_temp, lower = -10, upper = 10)$root}
fitLS2 = nls(P_G1 ~ F_fit(thresholds, delta, rho), start = list(delta = -deltaS, rho = 1/rhoS))
deltaOF = unname(coef(fitLS2)[1])
rhoOF = unname(coef(fitLS2)[2])
rhoOF = 1/rhoOF
deltaOF = -deltaOF*rhoOF
out = c((deltaF + deltaOF)/2, (rhoF + rhoOF)/2)
#suppressWarnings(if (DG == "Cd") {out[2] = sqrt((1 + out[2]^2)/2); out[1] = out[1]/out[2]}) #Warnings suppressed out of vectorization laziness
suppressWarnings(if (DG == "Cd") {out[3] = sqrt((1 + out[2]^2)/2); out[1] = out[1]/out[3]}) #Pooled SD, then mean difference/PSD
out = round(out[1:2], rnd)
names(out) = suppressWarnings(if (DG == "GDelta") {c("Glass' Delta", "SD Ratio")} else {c("Cohen's d", "SD Ratio")})
return(out)}Here’s code to compute Cohen’s d for comparisons of White and Black LSAT test-takers in the years 2019 (X1) through 2023 (X5). Levels in order are percentages with scores of <140, 140-144, 145-149, 150-154, 155-159, 160-165, and 165+.
W1 <- c(0.023673469, 0.04952381, 0.113741497, 0.197823129, 0.278639456, 0.227755102, 0.108843537)
W2 <- c(0.025876461, 0.045353367, 0.106010017, 0.202003339, 0.283806344, 0.225375626, 0.111574847)
W3 <- c(0.018368321, 0.041507634, 0.084208015, 0.171994275, 0.263120229, 0.244274809, 0.176526718)
W4 <- c(0.019420059, 0.035115722, 0.081670657, 0.166799681, 0.271348763, 0.25698324, 0.168661878)
W5 <- c(0.025215252, 0.042435424, 0.095325953, 0.167589176, 0.253690037, 0.234932349, 0.180811808)
A1 <- c(0.059557572, 0.071469087, 0.119115145, 0.205331821, 0.236528644, 0.197390811, 0.11060692)
A2 <- c(0.052246604, 0.067398119, 0.129571578, 0.200104493, 0.245036573, 0.173458725, 0.132183908)
A3 <- c(0.042837402, 0.064025795, 0.104099493, 0.17779825, 0.238599724, 0.2192538, 0.153385537)
A4 <- c(0.03870043, 0.059722886, 0.102723364, 0.162924032, 0.231724797, 0.222169135, 0.182035356)
A5 <- c(0.052554028, 0.058447937, 0.10805501, 0.158153242, 0.214636542, 0.212180747, 0.195972495)
"White-Asian"## [1] "White-Asian"MethodOfThresholds(W1, A1, DG = "Cd")## Cohen's d SD Ratio
## 0.143 1.138MethodOfThresholds(W2, A2, DG = "Cd")## Cohen's d SD Ratio
## 0.135 1.157MethodOfThresholds(W3, A3, DG = "Cd")## Cohen's d SD Ratio
## 0.163 1.100MethodOfThresholds(W4, A4, DG = "Cd")## Cohen's d SD Ratio
## 0.094 1.172MethodOfThresholds(W5, A5, DG = "Cd")## Cohen's d SD Ratio
## 0.065 1.166H1 <- c(0.047904192, 0.107784431, 0.167664671, 0.25748503, 0.19760479, 0.167664671, 0.053892216)
H2 <- c(0.075144509, 0.092485549, 0.184971098, 0.231213873, 0.219653179, 0.167630058, 0.028901734)
H3 <- c(0.052884615, 0.096153846, 0.110576923, 0.221153846, 0.221153846, 0.211538462, 0.086538462)
H4 <- c(0.047169811, 0.051886792, 0.146226415, 0.202830189, 0.245283019, 0.193396226, 0.113207547)
H5 <- c(0.036144578, 0.102409639, 0.144578313, 0.210843373, 0.210843373, 0.168674699, 0.126506024)
"White-Hispanic"## [1] "White-Hispanic"MethodOfThresholds(W1, H1, DG = "Cd")## Cohen's d SD Ratio
## 0.429 1.035MethodOfThresholds(W2, H2, DG = "Cd")## Cohen's d SD Ratio
## 0.506 1.011MethodOfThresholds(W3, H3, DG = "Cd")## Cohen's d SD Ratio
## 0.397 1.044MethodOfThresholds(W4, H4, DG = "Cd")## Cohen's d SD Ratio
## 0.343 1.060MethodOfThresholds(W5, H5, DG = "Cd")## Cohen's d SD Ratio
## 0.361 1.074B1 <- c(0.1179941, 0.156342183, 0.203539823, 0.218289086, 0.168141593, 0.097345133, 0.038348083)
B2 <- c(0.134986226, 0.143250689, 0.228650138, 0.201101928, 0.154269972, 0.090909091, 0.046831956)
B3 <- c(0.12605042, 0.121848739, 0.18907563, 0.210084034, 0.180672269, 0.121848739, 0.050420168)
B4 <- c(0.10046729, 0.175233645, 0.170560748, 0.203271028, 0.191588785, 0.112149533, 0.046728972)
B5 <- c(0.158808933, 0.148883375, 0.17369727, 0.196029777, 0.168734491, 0.104218362, 0.049627792)
"White-Black"## [1] "White-Black"MethodOfThresholds(W1, B1, DG = "Cd")## Cohen's d SD Ratio
## 0.786 1.102MethodOfThresholds(W2, B2, DG = "Cd")## Cohen's d SD Ratio
## 0.839 1.161MethodOfThresholds(W3, B3, DG = "Cd")## Cohen's d SD Ratio
## 0.829 1.093MethodOfThresholds(W4, B4, DG = "Cd")## Cohen's d SD Ratio
## 0.889 1.092MethodOfThresholds(W5, B5, DG = "Cd")## Cohen's d SD Ratio
## 0.885 1.104Ab1 <- c(0.14084507, 0.159624413, 0.154929577, 0.178403756, 0.215962441, 0.112676056, 0.037558685)
Ab2 <- c(0.130232558, 0.106976744, 0.181395349, 0.223255814, 0.2, 0.106976744, 0.051162791)
Ab3 <- c(0.088607595, 0.126582278, 0.147679325, 0.185654008, 0.202531646, 0.135021097, 0.113924051)
Ab4 <- c(0.100917431, 0.091743119, 0.155963303, 0.178899083, 0.256880734, 0.160550459, 0.055045872)
Ab5 <- c(0.133928571, 0.111607143, 0.160714286, 0.209821429, 0.196428571, 0.107142857, 0.080357143)
"White-Aboriginal"## [1] "White-Aboriginal"MethodOfThresholds(W1, Ab1, DG = "Cd")## Cohen's d SD Ratio
## 0.691 1.162MethodOfThresholds(W2, Ab2, DG = "Cd")## Cohen's d SD Ratio
## 0.658 1.124MethodOfThresholds(W3, Ab3, DG = "Cd")## Cohen's d SD Ratio
## 0.553 1.226MethodOfThresholds(W4, Ab4, DG = "Cd")## Cohen's d SD Ratio
## 0.623 1.058MethodOfThresholds(W5, Ab5, DG = "Cd")## Cohen's d SD Ratio
## 0.705 1.118https://rpubs.com/JLLJ/LSATSEL
https://rpubs.com/JLLJ/MMTMICH22
https://rpubs.com/JLLJ/threshold
sessionInfo()## R version 4.3.1 (2023-06-16 ucrt)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 19045)
##
## Matrix products: default
##
##
## locale:
## [1] LC_COLLATE=English_United States.utf8
## [2] LC_CTYPE=English_United States.utf8
## [3] LC_MONETARY=English_United States.utf8
## [4] LC_NUMERIC=C
## [5] LC_TIME=English_United States.utf8
##
## time zone: America/Chicago
## tzcode source: internal
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## loaded via a namespace (and not attached):
## [1] digest_0.6.33 R6_2.5.1 fastmap_1.1.1 xfun_0.39
## [5] cachem_1.0.8 knitr_1.43 htmltools_0.5.5 rmarkdown_2.23
## [9] cli_3.6.1 sass_0.4.7 jquerylib_0.1.4 compiler_4.3.1
## [13] rstudioapi_0.15.0 tools_4.3.1 evaluate_0.21 bslib_0.5.0
## [17] yaml_2.3.7 rlang_1.1.1 jsonlite_1.8.7