R/tt.1group.R

In lessR: Less Code with More Comprehensive Results

.OneGroup  <-
function(Y, Ynm, mu=NULL, n=NULL, m=NULL, s=NULL, brief, bw1,
         from.data, conf_level, alternative, digits_d, mmd, msmd,
         Edesired, paired, graph, xlab, line_chart, show_title, quiet,
         full_curve=TRUE, pdf_file, width, height, fill=NULL, ...) {

  # get lab_x_size  lab_y_size
  lab_size <- getOption("lab_size")
  lab_x_size <- getOption("lab_x_size")
  lab_x_size <- ifelse(is.null(lab_x_size), lab_size, lab_x_size)
  adj <- .RSadj(lab_size=lab_x_size); lab_x_size <- adj$lab_size

  # get variable labels if exist plus axes labels
  gl <- .getlabels(xlab=NULL, ylab=xlab, main=NULL, lab_x_size=lab_x_size,
                   graph.win=FALSE)  # # graphics window not yet set-up
  x.name <- gl$yn; x.lbl <- gl$yl; x.lab <- gl$yb
  main.lab <- gl$mb
  sub.lab <- gl$sb

  # get variable label if exists (redundant with above
  gl <- .getlabels(graph.win=FALSE)  # graphics window not yet set-up
  y.lbl <- gl$yl
  if ( (!is.null(y.lbl))  && !quiet) {
    cat("Response Variable:  ", Ynm, ", ", as.character(y.lbl), sep="", "\n")
    cat("\n")
  }

  if (!brief && !quiet) cat("\n------ Describe ------\n\n")

  if (from.data) {
    n <- sum(!is.na(Y))
    n.miss <- sum(is.na(Y))
    Y <- na.omit(Y)
    m <- mean(Y)
    s <- sd(Y)
    v <- var(Y)
  }
  else {
    v <- s^2
  }

  clpct <- paste(toString(round((conf_level)*100, 2)), "%", sep="")

  dig.smr.d <- ifelse (from.data, digits_d, digits_d - 1)

  if (!quiet) {
    if (Ynm != "Y") cat(Ynm,  ": ", sep="")
    if (from.data) cat(" n.miss = ", n.miss, ",  ", sep="") 
    cat("n = ", n, ",   mean = ", .fmt(m, dig.smr.d), 
        ",  sd = ", .fmt(s, dig.smr.d), sep="", "\n")

    if
      (brief) cat("\n")
    else  {
      if (from.data) {

        cat("\n\n------ Normality Assumption ------\n\n")
        # Normality
        if (n > 30) {
          cat("Sample mean assumed normal because n > 30, so no test needed.",
              sep="", "\n")
        }
        else {
          cat("Null hypothesis is a normal distribution", sep="")
          if (Ynm != "Y") cat(" of ", Ynm, ".", sep="") else cat(".")
          cat("\n")
          if (n > 2 && n < 5000) {
            nrm1 <- shapiro.test(Y)
            W.1 <- round(nrm1$statistic,min(4,digits_d+1))
            p.val1 <- round(nrm1$p.value,min(4,digits_d))
            cat(nrm1$method, ":  W = ", W.1, ",  p-value = ", p.val1,
                sep="", "\n")
          }
          else
            cat("Sample size out of range for Shapiro-Wilk normality test.",
                "\n")
        }  
      } 
    }

    if (!brief) cat("\n\n------ Infer ------\n\n")
  }  # end !quiet

  # t-test
  if (alternative == "two_sided")
    alt <- "two.sided"
  else
    alt <- alternative


  if (!is.null(mu)) m.dist <- m - mu
  df <- n - 1
  sterr <- s * sqrt(1/n)
  if (alternative == "two_sided")
    tcut <- qt((1-conf_level)/2, df=df, lower.tail=FALSE)
  else if (alternative == "less")
    tcut <- qt(1-conf_level, df=df, lower.tail=FALSE)
  else if (alternative == "greater")
    tcut <- qt(1-conf_level, df=df, lower.tail=FALSE)

  if (from.data) {
    if (!is.null(mu)) mu.null <- mu else mu.null <- 0
    ttest <- t.test(Y, conf.level=conf_level, alternative=alt, mu=mu.null)
    df <- ttest$parameter
    lb <- ttest$conf.int[1]
    ub <- ttest$conf.int[2]
    E <- (ub-lb)/2
    if (!is.null(mu)) {
      tvalue <- ttest$statistic
      pvalue <- ttest$p.value
    }
  }  # end from.data
  else {
    E <- tcut * sterr
    lb <- m - E
    ub <- m + E
    if (!is.null(mu)) {
      tvalue <- m.dist/sterr
      if (alternative == "two_sided")
        pvalue <- 2 * pt(abs(tvalue), df=df, lower.tail=FALSE)
      else if (alternative == "less")
        pvalue <- pt(tvalue, df=df, lower.tail=TRUE)
      else if (alternative == "greater")
        pvalue <- pt(tvalue, df=df, lower.tail=FALSE)
    }
  }     
  # difference from mu and standardized mean difference
  if (!is.null(mu)) {
    mdiff <- m - mu
    smd <- abs(mdiff/s)
  }

  if (!quiet) {
    cat("t-cutoff for ", clpct, " range of variation: tcut = ", .fmt(tcut,3), "\n") 
    cat("Standard Error of Mean: SE = ", .fmt(sterr), "\n\n")

    if (!is.null(mu)) {
    if (alt != "two.sided") 
      cat("\nAlternative hypothesis: Population mean difference is", alt, 
          "than", mu, "\n")
      cat("Hypothesized Value H0: mu =", mu, "\n")
      txt <- "Hypothesis Test of Mean:  t-value = "
      cat(txt, .fmt(tvalue,3), ",  df = ", df, ",  p-value = ", .fmt(pvalue,3),
          sep="", "\n\n")
    }
    cat("Margin of Error for ", clpct, " Confidence Level:  ", .fmt(E),
        sep="", "\n")
    txt <- " Confidence Interval for Mean:  "
    cat(clpct, txt, .fmt(lb), " to ", .fmt(ub), sep="", "\n")

    # difference from mu and standardized mean difference
    if (!is.null(mu)) {
      if (!brief) cat("\n\n------ Effect Size ------\n\n") else cat("\n")
      cat("Distance of sample mean from hypothesized:  " , .fmt(mdiff), "\n",
          "Standardized Distance, Cohen's d:  ", .fmt(smd),
          sep="", "\n")
    }

    # needed sample size from Edesired
    if (!is.null(Edesired)) {
      zcut <- qnorm((1-conf_level)/2)
      ns <- ((zcut*s)/Edesired)^2 
      n.needed <- ceiling(1.132*ns + 7.368) 
      cat("\n\n------ Needed Sample Size ------\n\n")
      if (Edesired > E) {
        cat("Note: Desired margin of error,", .fmt(Edesired), 
           "is worse than what was obtained,", .fmt(E), "\n\n") 
      }
      cat("Desired Margin of Error: ", .fmt(Edesired), "\n")
      cat("\n")
      cat("For the following sample size there is a 0.9 probability of obtaining\n")
      cat("the desired margin of error for the resulting 95% confidence interval.\n")
      cat("-------\n")
      cat("Needed sample size: ", n.needed, "\n")
      cat("\n")
      cat("Additional data values needed: ", n.needed-n, "\n")
    }
  }  # end !quiet


  # graphs
  if (graph) {

    # keep track of the number of plots in this routine, see if manage graphics
    plt.i <- 0
    plt.title  <- character(length=0)
    manage.gr <- .graphman()
 
    if (is.null(pdf_file)) {
      if (manage.gr) {
        n.win <- 1
        if (paired) n.win <- n.win + 1
        if (line_chart) n.win <- n.win + 1
        if (n.win > 0) {
          .graphwin(n.win, width, height)
          i.win <- 2   # start first graphics window on 3
          orig.params <- par(no.readonly=TRUE)
          on.exit(par(orig.params))
        }
      }
    }

    if (line_chart) {
      if (!is.null(pdf_file))
        pdf(file=paste("LineChart_",Ynm,".pdf",sep=""), width=width, height=height)

      if (manage.gr) {
        i.win <- i.win + 1 
        dev.set(which=i.win)
      }

      plt.i <- plt.i + 1
      plt.title[plt.i] <- paste("Sequentially Plotted Data for", Ynm)

      Y.df = data.frame(Y)
      x.df <- data.frame(1:nrow(Y.df))
      .plt.main(x.df, Y.df, segments=TRUE, size=.9,
                xlab="Index", ylab=NULL, main=plt.title[plt.i],
                use_plotly=FALSE)

      if (!is.null(pdf_file)) {
        dev.off()
        .showfile(paste("LineChart_", Ynm, ".pdf", sep=""),
                  paste("line chart of", Ynm))
      }

    }


    if (!is.null(mu)) {

      if (manage.gr) {
        i.win  <- i.win + 1 
        dev.set(which=i.win)
      }

      plt.i <- plt.i + 1
      plt.title[plt.i] <- "One-Group Plot"

      if (paired) x.lab <- "Differences of Matched Pairs"

      if (!is.null(pdf_file)) {
        if (!grepl(".pdf", pdf_file))
          pdf_file <- paste(pdf_file, ".pdf", sep="")
        .opendev(pdf_file, width, height)
      }

      .OneGraph(Y, bw1, Ynm, digits_d, brief,
           n, m, mu, mdiff, s, smd, mmd, msmd,
           clpct, tvalue, pvalue, ub, lb, x.lab, alt,
           show_title, full_curve=full_curve, quiet,
           fill=fill)

      if (!is.null(pdf_file)) {
        dev.off()
        .showfile(pdf_file, paste("density plot of", Ynm))
      }
    }

    else {  # mu=NULL: CI only — show density

      if (manage.gr) {
        i.win <- i.win + 1
        dev.set(which=i.win)
      }

      plt.i <- plt.i + 1
      plt.title[plt.i] <- "Density Plot"

      if (!is.null(pdf_file)) {
        if (!grepl(".pdf", pdf_file))
          pdf_file <- paste(pdf_file, ".pdf", sep="")
        .opendev(pdf_file, width, height)
      }

      if (!grepl("gray", getOption("theme"))) {
        clr.v  <- .color_range(.get_fill(), 1)[1]
        c1     <- as.vector(col2rgb(clr.v) / 255)
        col.1  <- rgb(c1[1], c1[2], c1[3])
        col.1d <- rgb(c1[1], c1[2], c1[3], alpha=.2)
      } else {
        col.1  <- rgb(.40, .40, .40)
        col.1d <- rgb(.40, .40, .40, .20)
      }
      col.tx <- getOption("lab_color")
      if (!is.null(fill)) {
        c1 <- as.vector(col2rgb(fill[1]) / 255)
        col.1  <- rgb(c1[1], c1[2], c1[3])
        col.1d <- rgb(c1[1], c1[2], c1[3], alpha=.2)
      }

      dYA   <- suppressWarnings(density(Y, bw1))
      min.x <- min(dYA$x)
      max.x <- max(dYA$x)
      max.y <- max(dYA$y) * 1.1

      xl <- if (length(x.lab) == 1) x.lab else x.lab[4]
      if (nchar(xl) > 52) xl <- paste(substr(xl, 1, 50), "...")

      orig.params <- par(no.readonly=TRUE)
      on.exit(par(orig.params), add=TRUE)
      par(bg=getOption("window_fill"))
      par(mar=c(4.1, 1.5, 4, .4), mgp=c(3, .6, 0), cex=.8)

      plot.new()
      plot.window(xlim=c(min.x, max.x), ylim=c(0, max.y))
      ax <- .axes_dim()
      axT1 <- pretty(c(min.x, max.x))
      .axes(NULL, NULL, axT1, NULL)
      box(col=getOption("panel_color"))
      adj <- .RSadj(axis_size=ax$axis_x_size)
      title(xlab=xl, col.lab=col.tx, cex.lab=adj$axis_size * 1.2)

      polygon(c(min(dYA$x), dYA$x, max(dYA$x)),
              c(0, dYA$y, 0), col=col.1d, border=NA)
      segments(min(dYA$x), 0, max(dYA$x), 0, col=col.1)
      lwd.b <- if (.Platform$OS == "windows") 2 else 1.75
      lines(dYA, col=col.1, lty="solid", lwd=lwd.b)
      y.m <- approx(dYA$x, dYA$y, xout=m)$y
      lines(c(m, m), c(0, y.m), lty="solid", lwd=.85,
            col="gray50")

      if (show_title) {
        mtext("One-Group Density Plot", side=3, line=2.4,
              font=2, col=col.tx)
        mtext(bquote(paste("  n=", .(n),
          "  m=", .(.fmt(m, digits_d)),
          "  s=", .(.fmt(s, digits_d)))),
          side=3, line=1.2, cex=.9, adj=0, col=col.tx)
        mtext(bquote(paste(
          "  ", .(clpct), " CI:  ",
          .(.fmt(lb, 3)), " to ", .(.fmt(ub, 3)))),
          side=3, line=0.1, cex=.9, adj=0, col=col.tx)
      }

      if (.allow.interactive()) {
        plt <- suppressWarnings(dn.plotly(
          x          = Y,
          x_name     = Ynm,
          fill       = col.1,
          x_lab      = xl,
          y_lab      = "Density",
          ax         = NULL,
          bw         = bw1,
          full_curve = full_curve,
          style_opts = .plotly_style()
        ))
        print(plt)
      }

      if (!is.null(pdf_file)) {
        dev.off()
        .showfile(pdf_file, paste("density plot of", Ynm))
      }
    }

    return(list(i=plt.i, ttl=plt.title))
  }  # end if graph

}  # End One Group