R/diffmean.neq.test.R

In LearningStats: Elemental Descriptive and Inferential Statistics

diffmean.neq.test <- function(x, n=NULL, y, m=NULL, Scx=NULL, Scy=NULL, alternative,
                              alpha = 0.05, plot = TRUE, lwd = 1) {

  if(!is.null(Scx)){if(is.na(Scx)) stop("'sc1' must be a valid number")}
  if(!is.null(Scy)){if(is.na(Scy)) stop("'sc2' must be a valid number")}

  if (length(x) != 1L){bar_x <- mean(x);Scx=sd(x);n=length(x)}else{bar_x <- x}
  if (length(y) != 1L){bar_y <- mean(y);Scy=sd(y);m=length(y)}else{bar_y <- y}

  DNAME <- paste0(deparse(substitute(x)), ", with sample size ", deparse(substitute(n)), ", and ",
                  deparse(substitute(y)), ", with sample size ", deparse(substitute(m)))
  NVAL <- 0

  # Statistic and pvalue
  g <- round(((Scx^2/n + Scy^2/m)^2 / ((Scx^2/n)^2/(n-1) + (Scy^2/m)^2/(m-1))) - 2, 0)
  if (g < 1) stop("not enough observations")

  STATISTIC <- (bar_x - bar_y) / sqrt(Scx^2/n + Scy^2/m)
  if (alternative == "two.sided")
    PVALUE <- 2 * pt(abs(STATISTIC), df = g, lower.tail = FALSE)
  else
    PVALUE <- pt(STATISTIC, df = g, lower.tail = (alternative == "less"))
  # Reject Region
  RR <- paste0("RR = ", switch(alternative,
                               two.sided = paste0("(-\U221E, ", round(qt(alpha/2, g), 5), "] U [",
                                                  round(qt(1 - alpha/2, g), 5), ", +\U221E)"),
                               greater = paste0("[", round(qt(1 - alpha, g), 5), ", +\U221E)"),
                               less = paste0("(-\U221E, ", round(qt(alpha, g), 5), "]")))

  # Plot
  if (plot) {
    ## Plot statistic distribution
    curve(dt(x, df = g), from = min(-3, -abs(STATISTIC) - 0.1), to = max(3, abs(STATISTIC) + 0.1),
          main = bquote(T ~ "follows" ~ "T"[.(g)]), axes = FALSE, xlab = "", ylab = "", lwd = lwd) # bquote(T ~ "\U2208" ~ "t"[.(n)])
    u <- par("usr") # x0, x1, y0, y1
    rect(u[1], 0, u[2], u[4])
    axis(2)
    legend("topright", c("p-value", "RR"), bty = "n", pch = c(22,NA), lty = c(NA,1), lwd = c(1,2),
           col = c("blue", "red"), pt.bg = adjustcolor('blue', alpha.f = 0.25), pt.cex = 2, seg.len = 1, cex = 1)
    if (alternative == "two.sided") {
      abline(h = 0, lwd = lwd + 1)
      lines(c(u[1], qt(alpha / 2, g)), c(0,0), col = "red", lwd = lwd + 1)
      lines(c(qt(1 - alpha / 2, g), u[2]), c(0,0), col = "red", lwd = lwd + 1)
      axis(1, pos = 0, col = NA, col.ticks = 1,
           at     = c(0, STATISTIC, -1*STATISTIC),
           labels = c(0, expression('T'[obs]), expression(-'T'[obs])))
      segments(x0 = c(qt(alpha / 2, g), qt(1 - alpha / 2, g)), y0 = -u[4]*0.015, y1 = u[4]*0.015, col = "red", lwd = lwd + 1)
      segments(x0 = c(qt(alpha / 2, g), qt(1 - alpha / 2, g)), y0 = rep(c(-u[4]*0.015, u[4]*0.015), each = 2),
               x1 = c(qt(alpha / 2, g) - u[2]*0.015, qt(1 - alpha / 2, g) + u[2]*0.015),
               rep(c(-u[4]*0.015, u[4]*0.015), each = 2), col = "red", lwd = lwd + 1)
      if ((abs(STATISTIC) - qt(1 - alpha / 2, g)) > 0.3) {
        axis(1, pos = 0, col = NA, col.ticks = NA,
             at     = c(qt(alpha / 2, g), qt(1 - alpha / 2, g)),
             labels = c(expression(-'t'[1-alpha/2]), expression('t'[1-alpha/2])))
        mtext("=", side = 1, line = 1.6, at = c(qt(alpha / 2, g), qt(1 - alpha / 2, g)), las = 2)
        mtext(round(qt(alpha / 2, g), 2), side = 1, line = 2.5, at = qt(alpha / 2, g))
        mtext(round(qt(1 - alpha / 2, g), 2), side = 1, line = 2.5, at = qt(1 - alpha / 2, g))
      }

      # pvalue != 0
      if (PVALUE > .Machine$double.eps) {
        # Statistic right tail
        segments(x0 = qt(1 - PVALUE / 2, g), y0 = 0,
                 x1 = qt(1 - PVALUE / 2, g), y1 = dt(qt(1 - PVALUE / 2, g), g),
                 col = 'blue', lwd = 1)
        x_vector <- seq(qt(1 - PVALUE / 2, g), 4, length = 100)
        y_vector <- dt(x_vector, g)
        polygon(c(x_vector, rev(x_vector)), c(y_vector, rep(0, length(y_vector))),
                col = adjustcolor('blue', alpha.f = 0.25), border = NA)
        mtext("=", side = 1, line = 1.6, at = qt(1 - PVALUE / 2, g), las = 2)
        mtext(round(qt(1 - PVALUE / 2, g), 2), side = 1, line = 2.5, at = qt(1 - PVALUE / 2, g))

        # Statistic left tail
        segments(x0 = qt(PVALUE / 2, g), y0 = 0,
                 x1 = qt(PVALUE / 2, g), y1 = dt(qt(PVALUE / 2, g), g),
                 col = 'blue', lwd = 1)
        x_vector <- seq(-4, qt(PVALUE / 2, g), length = 100)
        y_vector <- dt(x_vector, g)
        polygon(c(x_vector, rev(x_vector)), c(y_vector, rep(0, length(y_vector))),
                col = adjustcolor('blue', alpha.f = 0.25), border = NA)
        mtext("=", side = 1, line = 1.6, at = qt(PVALUE / 2, g), las = 2)
        mtext(round(qt(PVALUE / 2, g), 2), side = 1, line = 2.5, at = qt(PVALUE / 2, g))
      }

    } else {

      if (alternative == "less") {
        abline(h = 0, lwd = lwd + 1)
        lines(c(u[1], qt(alpha, g)), c(0,0), col = "red", lwd = lwd + 1)
        axis(1, pos = 0, col = NA, col.ticks = 1,
             at     = c(0, STATISTIC),
             labels = c(0, expression('T'[obs])))
        segments(x0 = qt(alpha, g), y0 = -u[4]*0.015, y1 = u[4]*0.015, col = "red", lwd = lwd + 1)
        segments(x0 = qt(alpha, g), y0 = c(-u[4]*0.015, u[4]*0.015),
                 x1 = qt(alpha, g) - u[2]*0.015, c(-u[4]*0.015, u[4]*0.015), col = "red", lwd = lwd + 1)
        if (abs(STATISTIC - qt(alpha, g)) > 0.3) {
          axis(1, pos = 0, col = NA, col.ticks = NA, at = qt(alpha, g), labels = expression(-'t'[1-alpha]))
          mtext("=", side = 1, line = 1.6, at = qt(alpha, g), las = 2)
          mtext(round(qt(alpha, g), 2), side = 1, line = 2.5, at = qt(alpha, g))
        }

      } else {
        abline(h = 0, lwd = lwd + 1)
        lines(c(qt(1 - alpha, g), u[2]), c(0,0), col = "red", lwd = lwd + 1)
        axis(1, pos = 0, col = NA, col.ticks = 1,
             at       = c(0, STATISTIC),
             labels   = c(0, expression('T'[obs])))
        segments(x0 = qt(1 - alpha, g), y0 = -u[4]*0.015, y1 = u[4]*0.015, col = "red", lwd = lwd + 1)
        segments(x0 = qt(1 - alpha, g), y0 = c(-u[4]*0.015, u[4]*0.015),
                 x1 = qt(1 - alpha, g) + u[2]*0.015, c(-u[4]*0.015, u[4]*0.015), col = "red", lwd = lwd + 1)
        if (abs(STATISTIC - qt(1 - alpha, g)) > 0.3) {
          axis(1, pos = 0, col = NA, col.ticks = NA, at = qt(1 - alpha, g), labels = expression('t'[1-alpha]))
          mtext("=", side = 1, line = 1.6, at = qt(1 - alpha, g), las = 2)
          mtext(round(qt(1 - alpha, g), 2), side = 1, line = 2.5, at = qt(1 - alpha, g))
        }
      }

      # pvalue != 0
      if (PVALUE > .Machine$double.eps) {
        # Statistic
        segments(x0 = STATISTIC, y0 = 0, x1 = STATISTIC, y1 = dt(STATISTIC, g),
                 col = 'blue', lwd = 1)
        if (alternative == "less") {
          x_vector <- seq(-4, STATISTIC, length = 100)
        } else {
          x_vector <- seq(STATISTIC, 4, length = 100)
        }
        y_vector <- dt(x_vector, g)
        polygon(c(x_vector, rev(x_vector)), c(y_vector, rep(0, length(y_vector))),
                col = adjustcolor('blue', alpha.f = 0.25), border = NA)
        mtext("=", side = 1, line = 1.6, at = STATISTIC, las = 2)
        mtext(round(STATISTIC, 2), side = 1, line = 2.5, at = STATISTIC)
      }
    }
  }

  ##---------------------------------------------
  METHOD <- "Test for the difference between the means of two Normal populations with different unknown variances"
  DISTNAME <- "\U2208 T_g"
  STATFORMULA <- "(bar.x\U2081 - bar.x\U2082) / sqrt((Sc\U2081\U00B2 / n\U2081) + (Sc\U2082\U00B2 / n\U2082)), \t
                  g = ((Sc\U2081\U00B2/n\U2081) + (Sc\U2082\U00B2/n\U2082))\U00B2 /
                  (((Sc\U2081\U00B2/n\U2081)\U00B2/(n\U2081-1)) + ((Sc\U2082\U00B2/n\U2082)\U00B2/(n\U2082-1))) - 2"
  ESTIMATE <- setNames(bar_x - bar_y, "\U03BC\U2081 - \U03BC\U2082")
  PARAMETER <- g
  names(NVAL) <- names(ESTIMATE)
  names(PARAMETER) <- "df"
  names(STATISTIC) <- "T"
  RVAL <- list(
    statistic = STATISTIC,
    parameter = PARAMETER,
    p.value = as.numeric(PVALUE),
    estimate = ESTIMATE,
    null.value = NVAL,
    alternative = alternative,
    method = METHOD,
    data.name = DNAME,
    alpha = alpha,
    dist.name = DISTNAME,
    statformula = STATFORMULA,
    reject.region = RR
  )
  class(RVAL) <- c("lstest", "htest")
  return(RVAL)
}