# Supplementary analyses of the Mansoura University Department of Obstetrics and Gynaecology papers.
# cf. https://steamtraen.blogspot.com/2021/10/a-catastrophic-failure-of-peer-review.html
# By Nick Brown, October 2021.
# License: CC-0

library(rpsychi)

# Function to display the possible ranges of the F or t statistic from a one- or two-way ANOVA.
f_range <- function (m, s, n, title=FALSE, show.t=FALSE, dp.p=-1, labels=c()) {
  m.ok <- m
  if (class(m.ok) == "matrix") {
    func <- ind.twoway.second
    useF <- c(3, 2, 4)
    default_labels <- c("col F", "row F", "inter F")
  }
  else {
    m.ok <- matrix(m)
    func <- ind.oneway.second
    useF <- 1
    default_labels <- c("F")
    if (show.t) {
      default_labels <- c("t")
    }
  }
  
  # Determine how many DPs to use from input numbers, if not specified
  dp <- dp.p
  if (dp.p == -1) {
    dp <- 0
    numbers <- c(m, s)
    for (i in numbers) {
      if (i != round(i, 0)) {
        dp <- max(dp, 1)
        j <- i * 10
        if (j != round(j, 0)) {
          dp <- max(dp, 2)
        }
      }
    }
  }
  
  if (length(labels) == 0) {
    labels <- default_labels
  }
  
  # Calculate the nominal test statistic(s) (i.e., assuming no rounding error)
  f.nom <- func(m=m.ok, sd=s, n=n)$anova.table$F

  # We correct for rounding in reported numbers by allowing for the maximum possible rounding error.
  # For the maximum F estimate, we subtract .005 from all SDs; for minimum F estimate, we add .005.
  # We then add or subtract .005 to every mean, in all possible permutations.
  # (".005" is an example, based on 2 decimal places of precision.)
  delta <- (0.1 ^ dp) / 2    #typically 0.005
  tiny.s <- delta / 100      #minimum SD, to prevent divide-by-zero errors

  s.hi <- pmax(s - delta, tiny.s)
  s.lo <- pmax(s + delta, tiny.s)

  # Initialise maximum and minimum F statistics to unlikely values.
  f.hi <- rep(-1, length(useF))
  f.lo <- rep(999999, length(useF))
  f.hi <- f.nom
  f.lo <- f.nom
  
  # Generate every possible combination of +/- maximum rounding error to add to each mean.
  l <- length(m.ok)
  rawcomb <- combn(rep(c(-delta, delta), l), l)
  comb <- rawcomb[,!duplicated(t(rawcomb))]
  
  # Generate every possible set of test statistics within the bounds of rounding error,
  #  and retain the largest and smallest.
  for (i in 1:ncol(comb)) {
    m.adj <- m.ok + comb[,i]

    if (all((abs(m.adj - m.adj[1])) < 1e-14)) {
      adj.f.hi <- 0
      adj.f.lo <- 0
    }
    else {
      adj.f.hi <- func(m=m.adj, sd=s.hi, n=n)$anova.table$F
      adj.f.lo <- func(m=m.adj, sd=s.lo, n=n)$anova.table$F
    }

    f.hi <- pmax(f.hi, adj.f.hi)
    f.lo <- pmin(f.lo, adj.f.lo)
  }
  
  if (show.t) {
    f.nom <- sqrt(f.nom)
    f.hi <-  sqrt(f.hi)
    f.lo <-  sqrt(f.lo)
  }
  
  if (title != FALSE) {
    cat(title)
  }
  
  sp <- " "
  fdp <- 3     # number of DPs for F or t statistic
  dpf <- paste("%.", fdp, "f", sep="")
  for (i in 1:length(useF)) {
    j <- useF[i]
    cat(sp, labels[i], ": ", sprintf(dpf, f.nom[j]),
        " (min=", sprintf(dpf, f.lo[j]),
        ", max=", sprintf(dpf, f.hi[j]), ")",
        sep="")
    sp <- "  "
  }
  
  if ((dp.p == -1)  && (dp < 2)) {
    cat(" <<< dp set to", dp, "automatically")
  }
  
  cat("\n", sep="")
}

# Examples for Badawy papers:

# BAD202 (10.1016/j.fertnstert.2008.01.016) height
f_range(m=c(155.30, 159.10), s=c(5.21, 5.34), n=c(110, 102), show.t=TRUE)
# t: 5.243 (min=5.224, max=5.262)

# BAD205 (10.1016/j.fertnstert.2006.02.097) height
f_range(m=c(168.0, 165.3), s=c(4.9, 5.1), n=c(404, 400), show.t=TRUE)
# t: 7.655 (min=7.298, max=8.019) <<< dp set to 1 automatically

# BAD207 (10.1016/j.fertnstert.2007.08.034) age
f_range(m=c(23.8, 25.3), s=c(3.1, 2.9), n=c(115, 101), show.t=TRUE)
# t: 3.657 (min=3.357, max=3.966) <<< dp set to 1 automatically

# ANOVA analyses for spreadsheet:

# FAW110 (10.1007/s00404-007-0527-x)
n_faw110 <- c(57, 53, 50)

f_range(m=c(28.7, 29.1, 29.4), s=c(5.4, 6.1, 7.1), n=n_faw110)                #T1 Age, p=NS
# F: 0.172 (min=0.123, max=0.233) <<< dp set to 1 automatically
# *** p not NS, as was claimed by authors

f_range(m=c(3.2, 3.9, 3.7), s=c(1.3, 1.5, 2.4), n=n_faw110)                   #T1 Miscarriages, p=NS
# F: 2.284 (min=1.553, max=3.218) <<< dp set to 1 automatically
# *** p not NS, as was claimed by authors

f_range(m=c(38.3, 37, 37.1), s=c(1.8, 2.2, 2.4), n=n_faw110)                  #T2 EGA_at_birth, p=NS
# F: 6.334 (min=5.126, max=7.738) <<< dp set to 1 automatically
# *** p not NS, as was claimed by authors

f_range(m=c(3212, 3115, 3090), s=c(547, 659, 715), n=n_faw110)                #T2 Birth_weight, p=NS
# F: 0.552 (min=0.542, max=0.563) <<< dp set to 0 automatically

f_range(m=c(9.2, 8.93, 8.32), s=c(0.89, 0.72, 0.98), n=n_faw110)              #T2 Apgar_1min, p=NS
# F: 14.147 (min=13.651, max=14.658)
# *** p not NS, as was claimed by authors

f_range(m=c(9.86, 9.76, 9.1), s=c(0.41, 0.36, 0.78), n=n_faw110)              #T2 Apgar_5min, p=NS
# F: 30.215 (min=28.881, max=31.603)
# *** p not NS, as was claimed by authors

f_range(m=c(38.1, 36.8, 36.9), s=c(2.1, 2.3, 2.1), n=n_faw110)                #T2 EGA_at_SVD, p=NS
# F: 6.139 (min=4.970, max=7.496) <<< dp set to 1 automatically
# *** p not NS, as was claimed by authors

f_range(m=c(38.6, 37.9, 37.5), s=c(1.9, 1.5, 1.4), n=n_faw110)                #T2 EGA_at_CS, p=NS
# F: 6.325 (min=4.854, max=8.122) <<< dp set to 1 automatically
# *** p not NS, as was claimed by authors

f_range(m=c(11.8, 14.2, 13.5), s=c(0.8, 0.5, 0.1), n=c(11, 8, 26))            #T4 EGA_at_loss, p=NS
# F: 79.741 (min=61.230, max=103.682) <<< dp set to 1 automatically
# *** p not NS, as was claimed by authors

# BAD209 (10.1016/s1472-6483(10)61046-2)
n_bad209 <- c(58, 61, 60)

f_range(m=c(28.80, 31.80, 29.70), s=c(0.30, 0.60, 0.50), n=n_bad209)                   #T1 Age, p=NS
# F: 602.097 (min=462.556, max=796.486) <<< dp set to 1 automatically
# *** p not NS, as was claimed by authors

f_range(m=c(0.50, 0.60, 0.60), s=c(0.10, 0.20, 0.30), n=n_bad209)                      #T1 Parity, p=NS
# F: 4.167 (min=0.000, max=16.668) <<< dp set to 1 automatically 

f_range(m=c(2.80, 3.10, 4.10), s=c(0.30, 0.30, 0.40), n=n_bad209)                      #T1 Infertility, p=NS
# F: 242.723 (min=155.160, max=391.126) <<< dp set to 1 automatically
# *** p not NS, as was claimed by authors

f_range(m=c(24.10, 23.10, 25.10), s=c(3.20, 3.30, 3.10), n=n_bad209)                   #T1 BMI, p=NS
# F: 5.902 (min=5.167, max=6.722) <<< dp set to 1 automatically
# *** p not NS, as was claimed by authors

f_range(m=c(4.60, 4.90, 3.90), s=c(0.70, 0.50, 0.50), n=n_bad209)                      #T1 FSH, p=NS
# F: 48.443 (min=32.764, max=71.164) <<< dp set to 1 automatically
# *** p not NS, as was claimed by authors

f_range(m=c(4.80, 5.20, 4.70), s=c(0.50, 0.40, 0.40), n=n_bad209)                      #T1 LH, p=NS
# F: 22.386 (min=11.165, max=42.102) <<< dp set to 1 automatically
# *** p not NS, as was claimed by authors

f_range(m=c(4.20, 5.80, 6.90), s=c(0.40, 0.30, 0.30), n=n_bad209)                      #T2 Total_follicles, p=<0.05
# F: 963.285 (min=675.957, max=1429.742) <<< dp set to 1 automatically
# *** p not NS, as was claimed by authors

f_range(m=c(1.60, 2.10, 4.10), s=c(0.10, 0.10, 0.20), n=n_bad209)                      #T2 Follicles>14mm, p=<0.05
# Error in delta.upper + dfb : non-numeric argument to binary operator
# Numbers are too extreme (F > many thousands), but at least p is indeed <0.05

f_range(m=c(1.00, 1.40, 3.40), s=c(0.00, 0.10, 0.30), n=n_bad209)                      #T2 Follicles≥18mm, p=<0.05
# F: 2932.325 (min=2677.531, max=3198.998) <<< dp set to 1 automatically

f_range(m=c(12.40, 11.60, 9.60), s=c(0.30, 0.50, 0.40), n=n_bad209)                    #T2 Stimulation, p=<0.05
# F: 734.646 (min=541.553, max=1021.638) <<< dp set to 1 automatically

f_range(m=c(4.40, 4.10, 4.30), s=c(0.50, 0.50, 0.40), n=n_bad209)                      #T2 Pretreatment_endometrial_thickness, p=NS
# F: 6.348 (min=2.211, max=14.791) <<< dp set to 1 automatically
# *** p not NS, as was claimed by authors

f_range(m=c(8.50, 8.20, 10.20), s=c(0.30, 0.30, 0.40), n=n_bad209)                     #T2 Endometrial_thickness, p=<0.05
# F: 615.833 (min=419.782, max=937.737) <<< dp set to 1 automatically

f_range(m=c(257.10, 284.00, 302.00), s=c(68.20, 91.30, 71.30), n=n_bad209)             #T2 Oestradiol, p=NS
# F: 4.967 (min=4.937, max=4.996) <<< dp set to 1 automatically
# *** p not NS, as was claimed by authors

f_range(m=c(7.10, 9.10, 11.30), s=c(0.80, 1.20, 1.30), n=n_bad209)                     #T2 Progesterone, p=<0.05
# F: 205.963 (min=180.036, max=236.293) <<< dp set to 1 automatically

# BAD210 (10.1080/00016340802638199)
n_bad210 <- c(269, 107, 420, 200)

f_range(m=c(25.80, 24.30, 23.50, 25.00), s=c(3.20, 2.90, 2.60, 1.90), n=n_bad210)      #T1 Age; reported F=1.09, p=0.67
# F: 42.736 (min=37.367, max=48.721) <<< dp set to 1 automatically

f_range(m=c(0.40, 0.30, 0.60, 0.30), s=c(0.12, 0.18, 0.20, 0.11), n=n_bad210)          #T1 Parity; reported F=1.1, p=0.62
# F: 212.656 (min=186.356, max=242.489)

f_range(m=c(3.20, 4.10, 3.70, 3.30), s=c(1.22, 1.61, 1.21, 0.11), n=n_bad210)          #T1 Infertility; reported F=1.2, p=0.08
# F: 22.250 (min=21.400, max=23.126)

f_range(m=c(158.30, 155.10, 157.20, 160.10), s=c(5.12, 4.20, 5.90, 3.51), n=n_bad210)  #T1 Height; reported F=2.6, p=0.55
# F: 26.186 (min=25.991, max=26.382)

f_range(m=c(80.30, 79.10, 75.10, 74.10), s=c(5.42, 4.22, 4.13, 5.11), n=n_bad210)      #T1 Weight; reported F=6.4, p=0.07
# F: 97.667 (min=97.100, max=98.237)

f_range(m=c(30.10, 29.10, 27.20, 26.20), s=c(2.91, 3.12, 2.33, 3.44), n=n_bad210)      #T1 BMI; reported F=7.3, p=0.08
# F: 93.076 (min=92.181, max=93.978)

f_range(m=c(37.10, 37.30, 38.00, 38.10), s=c(2.30, 2.70, 1.80, 1.80), n=n_bad210)      #T3 Gestational_age; reported F=n/a, p=NS
# F: 14.492 (min=10.914, max=18.801) <<< dp set to 1 automatically


f_range(m=c(2916, 3011, 2975, 3005), s=c(354, 401, 410, 398), n=n_bad210)              #T3 Birthweight; reported F=n/a, p=NS
# F: 2.641 (min=2.573, max=2.710) <<< dp set to 0 automatically
1-pf(2.573,3 , (sum(n_bad210) - 4))
# [1] 0.05280154

# Fisher's exact test analyses for spreadsheet:

fisher.test(matrix(c(56, 114, 62, 108), ncol=2))	#10.1080/01443610802042688
fisher.test(matrix(c(23, 147, 28, 142), ncol=2))	#10.1080/01443610802042688
fisher.test(matrix(c(30, 140, 28, 142), ncol=2))	#10.1080/01443610802042688
fisher.test(matrix(c(10, 160, 12, 158), ncol=2))	#10.1080/01443610802042688
fisher.test(matrix(c(17, 153, 11, 159), ncol=2))	#10.1080/01443610802042688
fisher.test(matrix(c(9, 161, 6, 164), ncol=2))	#10.1080/01443610802042688
fisher.test(matrix(c(7, 163, 15, 155), ncol=2))	#10.1080/01443610802042688
fisher.test(matrix(c(2, 168, 4, 166), ncol=2))	#10.1080/01443610802042688
fisher.test(matrix(c(2, 168, 3, 167), ncol=2))	#10.1080/01443610802042688
fisher.test(matrix(c(6, 164, 9, 161), ncol=2))	#10.1080/01443610802042688
fisher.test(matrix(c(51, 119, 48, 122), ncol=2))	#10.1080/01443610802042688
fisher.test(matrix(c(10, 160, 9, 161), ncol=2))	#10.1080/01443610802042688
fisher.test(matrix(c(30, 140, 28, 142), ncol=2))	#10.1080/01443610802042688
fisher.test(matrix(c(7, 163, 3, 167), ncol=2))	#10.1080/01443610802042688
fisher.test(matrix(c(20, 150, 24, 146), ncol=2))	#10.1080/01443610802042688
fisher.test(matrix(c(16, 154, 13, 157), ncol=2))	#10.1080/01443610802042688
fisher.test(matrix(c(2, 168, 0, 170), ncol=2))	#10.1080/01443610802042688
fisher.test(matrix(c(3, 167, 1, 169), ncol=2))	#10.1080/01443610802042688

fisher.test(matrix(c(38, 16, 38, 13), ncol=2))	#10.1111/j.1447-0756.2012.02016.x
fisher.test(matrix(c(12, 42, 10, 41), ncol=2))	#10.1111/j.1447-0756.2012.02016.x
fisher.test(matrix(c(16, 38, 6, 45), ncol=2))	#10.1111/j.1447-0756.2012.02016.x
fisher.test(matrix(c(14, 40, 5, 46), ncol=2))	#10.1111/j.1447-0756.2012.02016.x
fisher.test(matrix(c(1, 53, 0, 51), ncol=2))	#10.1111/j.1447-0756.2012.02016.x
fisher.test(matrix(c(2, 52, 1, 50), ncol=2))	#10.1111/j.1447-0756.2012.02016.x

fisher.test(matrix(c(55, 61, 57, 51), ncol=2))	#10.1007/s00404-015-3767-1
fisher.test(matrix(c(61, 55, 51, 57), ncol=2))	#10.1007/s00404-015-3767-1
fisher.test(matrix(c(33, 83, 30, 78), ncol=2))	#10.1007/s00404-015-3767-1
fisher.test(matrix(c(59, 57, 50, 58), ncol=2))	#10.1007/s00404-015-3767-1
fisher.test(matrix(c(75, 41, 74, 34), ncol=2))	#10.1007/s00404-015-3767-1
fisher.test(matrix(c(28, 88, 28, 80), ncol=2))	#10.1007/s00404-015-3767-1
fisher.test(matrix(c(13, 103, 6, 102), ncol=2))	#10.1007/s00404-015-3767-1
fisher.test(matrix(c(90, 26, 80, 28), ncol=2))	#10.1007/s00404-015-3767-1
fisher.test(matrix(c(20, 96, 21, 87), ncol=2))	#10.1007/s00404-015-3767-1
fisher.test(matrix(c(6, 110, 7, 101), ncol=2))	#10.1007/s00404-015-3767-1
fisher.test(matrix(c(41, 75, 57, 51), ncol=2))	#10.1007/s00404-015-3767-1
fisher.test(matrix(c(14, 102, 32, 76), ncol=2))	#10.1007/s00404-015-3767-1
fisher.test(matrix(c(2, 114, 5, 103), ncol=2))	#10.1007/s00404-015-3767-1
fisher.test(matrix(c(18, 98, 16, 92), ncol=2))	#10.1007/s00404-015-3767-1
fisher.test(matrix(c(98, 18, 92, 16), ncol=2))	#10.1007/s00404-015-3767-1
fisher.test(matrix(c(50, 66, 39, 69), ncol=2))	#10.1007/s00404-015-3767-1
fisher.test(matrix(c(86, 30, 73, 35), ncol=2))	#10.1007/s00404-015-3767-1
fisher.test(matrix(c(32, 84, 38, 70), ncol=2))	#10.1007/s00404-015-3767-1
fisher.test(matrix(c(3, 113, 0, 108), ncol=2))	#10.1007/s00404-015-3767-1
fisher.test(matrix(c(25, 91, 21, 87), ncol=2))	#10.1007/s00404-015-3767-1
fisher.test(matrix(c(1, 115, 0, 108), ncol=2))	#10.1007/s00404-015-3767-1
fisher.test(matrix(c(6, 110, 8, 100), ncol=2))	#10.1007/s00404-015-3767-1
fisher.test(matrix(c(25, 204, 50, 160), ncol=2))	#10.1007/s00404-015-3767-1
fisher.test(matrix(c(172, 57, 168, 42), ncol=2))	#10.1007/s00404-015-3767-1
fisher.test(matrix(c(82, 147, 111, 99), ncol=2))	#10.1007/s00404-015-3767-1
fisher.test(matrix(c(166, 63, 168, 42), ncol=2))	#10.1007/s00404-015-3767-1
fisher.test(matrix(c(33, 196, 64, 146), ncol=2))	#10.1007/s00404-015-3767-1
fisher.test(matrix(c(30, 199, 56, 154), ncol=2))	#10.1007/s00404-015-3767-1

fisher.test(matrix(c(32, 7, 31, 8), ncol=2))	#10.1016/j.fertnstert.2007.12.044
fisher.test(matrix(c(25, 14, 30, 9), ncol=2))	#10.1016/j.fertnstert.2007.12.044
fisher.test(matrix(c(35, 4, 13, 26), ncol=2))	#10.1016/j.fertnstert.2007.12.044
fisher.test(matrix(c(27, 12, 10, 29), ncol=2))	#10.1016/j.fertnstert.2007.12.044
fisher.test(matrix(c(4, 35, 21, 18), ncol=2))	#10.1016/j.fertnstert.2007.12.044

fisher.test(matrix(c(22, 8, 2, 28), ncol=2))	#10.1016/j.ijgo.2015.03.043
fisher.test(matrix(c(4, 26, 2, 28), ncol=2))	#10.1016/j.ijgo.2015.03.043
fisher.test(matrix(c(4, 26, 26, 4), ncol=2))	#10.1016/j.ijgo.2015.03.043

fisher.test(matrix(c(15, 22, 9, 25), ncol=2))	#10.1016/j.jmig.2012.11.006
fisher.test(matrix(c(7, 30, 8, 26), ncol=2))	#10.1016/j.jmig.2012.11.006
fisher.test(matrix(c(13, 24, 13, 21), ncol=2))	#10.1016/j.jmig.2012.11.006
fisher.test(matrix(c(2, 35, 4, 30), ncol=2))	#10.1016/j.jmig.2012.11.006
fisher.test(matrix(c(34, 3, 25, 9), ncol=2))	#10.1016/j.jmig.2012.11.006
fisher.test(matrix(c(2, 35, 2, 32), ncol=2))	#10.1016/j.jmig.2012.11.006
fisher.test(matrix(c(1, 36, 7, 27), ncol=2))	#10.1016/j.jmig.2012.11.006

fisher.test(matrix(c(2, 52, 10, 44), ncol=2))	#10.1016/j.ejogrb.2012.09.014