`WRS2` feature requests: Part 2

[IndrajeetPatil]

`WRS2` feature requests: Part 2

[IndrajeetPatil]

bootstrap-t method for one-sample test

trimcibt <-
  function(x,
           tr = .2,
           alpha = .05,
           nboot = 599,
           side = TRUE,
           plotit = FALSE,
           op = 1,
           nullval = 0,
           SEED = TRUE,
           prCRIT = FALSE,
           pr = TRUE,
           xlab = "") {
    #
    #  Compute a 1-alpha confidence interval for the trimmed mean
    #  using a bootstrap percentile t method.
    #
    #  The default amount of trimming is tr=.2
    #  side=T, for true,  indicates the symmetric two-sided method
    #
    #  Side=F yields an equal-tailed confidence interval
    #
    #  NOTE: p.value is reported when side=T only.
    #
    #  prCRIT=TRUE, if side=FALSE, lower and upper critical values are returned.
    # plotit=TRUE, plots an estimate of the null distribution of the Tukey--McLaughlin test statistic.
    # op=1, use adaptive kernel density estimator when plotit=TRUE
    # op=2, uses the expected frequency curve, which will have a lower execution if nboot is large.
    #
    x = WRS2:::elimna(x)
    side <- as.logical(side)
    p.value <- NA
    if (SEED)
      set.seed(2) # set seed of random number generator so that
    #   results can be duplicated.
    test <- (mean(x, tr) - nullval) / WRS2::trimse(x, tr)
    data <-
      matrix(sample(x, size = length(x) * nboot, replace = TRUE), nrow = nboot)
    data <- data - mean(x, tr)
    top <- apply(data, 1, mean, tr)
    bot <- apply(data, 1, WRS2::trimse, tr)
    tval <- top / bot
    if (plotit) {
      if (op == 1)
        akerd(tval, xlab = xlab)
      if (op == 2)
        rdplot(tval, xlab = xlab)
    }
    if (side)
      tval <- abs(tval)
    tval <- sort(tval)
    icrit <- round((1 - alpha) * nboot)
    ibot <- round(alpha * nboot / 2)
    itop <-
      nboot - ibot #altered code very slightly to correspond to recent versions of my books.
    if (!side) {
      if (prCRIT)
        print(paste("Lower crit=", tval[ibot], "Upper crit=", tval[itop]))
      if (prCRIT)
        print(paste(".025 Lower Type I=", mean(tval <= 0 - 1.96)))
      if (prCRIT)
        print(paste(".05 Lower Type I=", mean(tval <= 0 - 1.645)))
      if (prCRIT)
        print(paste(".025 Upper Type I=", mean(tval >= 1.96)))
      if (prCRIT)
        print(paste(".05 Upper Type I=", mean(tval >= 1.645)))
      trimcibt <- mean(x, tr) - tval[itop] * WRS2::trimse(x, tr)
      trimcibt[2] <- mean(x, tr) - tval[ibot] * WRS2::trimse(x, tr)
      if (pr)
        print("NOTE: p.value is computed only when side=T")
    }
    if (side) {
      if (prCRIT)
        print(paste("Symmetric Crit.val=", tval[icrit]))
      trimcibt <- mean(x, tr) - tval[icrit] * WRS2::trimse(x, tr)
      trimcibt[2] <- mean(x, tr) + tval[icrit] * WRS2::trimse(x, tr)
      p.value <- (sum(abs(test) <= abs(tval))) / nboot
    }
    list(
      estimate = mean(x, tr),
      ci = trimcibt,
      test.stat = test,
      p.value = p.value,
      n = length(x)
    )
  }


trimcibt(mtcars$wt, nullval = 3)
#> $estimate
#> [1] 3.197
#> 
#> $ci
#> [1] 2.841377 3.552623
#> 
#> $test.stat
#> [1] 1.179181
#> 
#> $p.value
#> [1] 0.2537563
#> 
#> $n
#> [1] 32

^{Created on 2020-04-04 by the reprex package (v0.3.0.9001)}

[IndrajeetPatil]

Simpler version of the same function:

trimcibt <- function(x,
                     tr = 0.2,
                     alpha = 0.05,
                     nboot = 200,
                     nullval = 0,
                     ...) {
  test <- (mean(x, tr) - nullval) / WRS2::trimse(x, tr)
  data <- matrix(sample(x, size = length(x) * nboot, replace = TRUE), nrow = nboot) - mean(x, tr)
  top <- apply(data, 1, mean, tr)
  bot <- apply(data, 1, WRS2::trimse, tr)
  tval <- sort(abs(top / bot))
  icrit <- round((1 - alpha) * nboot)
  ibot <- round(alpha * nboot / 2)
  itop <- nboot - ibot
  trimcibt <- mean(x, tr) - tval[icrit] * WRS2::trimse(x, tr)
  trimcibt[2] <- mean(x, tr) + tval[icrit] * WRS2::trimse(x, tr)
  p.value <- (sum(abs(test) <= abs(tval))) / nboot

  list(
    estimate = mean(x, tr),
    ci = trimcibt,
    test.stat = test,
    p.value = p.value,
    n = length(x)
  )
}

trimcibt(mtcars$wt, nullval = 3)
#> $estimate
#> [1] 3.197
#> 
#> $ci
#> [1] 2.857201 3.536799
#> 
#> $test.stat
#> [1] 1.179181
#> 
#> $p.value
#> [1] 0.275
#> 
#> $n
#> [1] 32

^{Created on 2021-01-09 by the reprex package (v0.3.0.9001)}

[IndrajeetPatil]

nboot argument needs to be added.
(This is painfully slow function)

library(WRS2)

wmcpAKP <- function(x, tr = .2) {
  #  Compute Algina et al measure of effect size for pairs of J dependent groups
  if (is.matrix(x) || is.data.frame(x)) x <- WRS2:::listm(x)
  J <- length(x)
  C <- (J^2 - J) / 2
  A <- matrix(NA, nrow = C, ncol = 5)
  dimnames(A) <- list(NULL, c("Group", "Group", "Effect.Size", "ci.low", "ci.up"))
  ic <- 0
  for (j in 1:J) {
    for (k in 1:J) {
      if (j < k) {
        ic <- ic + 1
        A[ic, 1] <- j
        A[ic, 2] <- k
        A[ic, 3] <- WRS2::dep.effect(x[[j]], x[[k]], tr = tr)[5]
        A[ic, 4] <- WRS2::dep.effect(x[[j]], x[[k]], tr = tr)[21]
        A[ic, 5] <- WRS2::dep.effect(x[[j]], x[[k]], tr = tr)[25]
      }
    }
  }
  A
}


wAKP.avg <- function(x, tr = .2) {
  #
  # Have J dependent groups. For each pair of groups, compute
  # AKP type measure of effect size and average the results
  #
  #  For tr=0, get Cohen d type measure
  #
  #
  a <- wmcpAKP(x, tr = tr)

  data.frame(
    "Effect.Size" = mean(a[, 3]),
    "ci.low" = mean(a[, 4]),
    "ci.up" = mean(a[, 5])
  )
}


before <- c(190, 210, 300, 240, 280, 170, 280, 250, 240, 220)
now <-    c(170, 280, 250, 240, 190, 260, 180, 200, 100, 200)
after <-  c(210, 210, 340, 190, 260, 180, 200, 220, 230, 200)

df <- data.frame(before, now, after)

wAKP.avg(df)
#>   Effect.Size     ci.low    ci.up
#> 1   0.2867495 -0.7193793 1.424093

^{Created on 2021-01-14 by the reprex package (v0.3.0.9001)}