R/dn.main.R

In lessR: Less Code, More Results

.dn.main <- 
function(x, 
         bw="nrd0", type="both",
         histo=TRUE, bin_start=NULL, bin_width=NULL,
         fill_hist, col.nrm, col.gen, fill_nrm, fill_gen,
         rotate_x=0, rotate_y=0, offset=0.5, 
         x.pt=NULL, xlab=NULL, main=NULL, sub=NULL,
         y_axis=FALSE, x.min=NULL, x.max=NULL,
         band=FALSE, color_rug="gray40", size_rug=0.5, quiet,
         fncl=NULL, ...)  {

  if (!is.null(x.pt)) {
    y_axis <- TRUE
    type <- "general"
  }

  col.bg <- getOption("panel_fill")
  col.box <- getOption("panel_color")

  lab_cex <- getOption("lab_cex")
  axis_cex <- getOption("axis_cex")
  col.axis <- getOption("axis_text_color")

  # get lab_x_cex  lab_y_cex
  lab_cex <- getOption("lab_cex")
  lab_x_cex <- getOption("lab_x_cex")
  lab_x_cex <- ifelse(is.null(lab_x_cex), lab_cex, lab_x_cex)
  adj <- .RSadj(lab_cex=lab_x_cex); lab_x_cex <- adj$lab_cex

  # get variable labels if exist plus axes labels
  gl <- .getlabels(xlab, ylab=NULL, main, lab_x_cex=lab_x_cex) 
  x.name <- gl$xn; x.lbl <- gl$xl; x.lab <- gl$xb
  main.lab <- gl$mb
  sub.lab <- gl$sb

  # get breaks from user supplied bin width and/or supplied start value
  # otherwise, breaks="Sturges by default
  # copied from histogram
  if (!is.null(bin_width)  || !is.null(bin_start)) {
    if (is.null(bin_start))
       bin_start <- pretty(min(x, na.rm=TRUE):max(x, na.rm=TRUE))[1]
    if (is.null(bin_width)) {
      h <- hist(x, plot=FALSE, breaks="Sturges")
      bin_width <- h$breaks[2]-h$breaks[1]
    }
    max.x <- max(x, na.rm = TRUE)
    seq.end <- max.x
    breaks <- seq(bin_start, seq.end, bin_width)
    while (max(breaks) < max.x) {
      seq.end <- seq.end + bin_width
      breaks <- seq(bin_start, seq.end, bin_width)
    }

  }
  else
    breaks="Sturges"

    cat("\n")

  if (is.null(main)) {
    orig.params <- par(no.readonly=TRUE)
    on.exit(par(orig.params))
    par(mar=c(4,4,2,2)+0.1)
  }

  # histogram calculations, no plot
  h <- hist(x, plot=FALSE, breaks)

  # no missing data
  n <- sum(!is.na(x))
  n.miss <- sum(is.na(x))
  if (n.miss > 0) x <- na.omit(x)
 
  # general density curve, no plot
  # suppress warnings about possible graphic parameters
  d.gen <- suppressWarnings(density(x, bw, ...))
  
  mx <- mean(x)

  # min and max x coordinates for graph, make symmetric
  min_dev.x <- min(d.gen$x) - mx
  max.dev.x <- max(d.gen$x) - mx
  if (abs(min_dev.x) > abs(max.dev.x)) {
    if (is.null(x.min)) x.min <- min(d.gen$x)
    if (is.null(x.max)) x.max <- mx + abs(min_dev.x)
  }
  if (abs(max.dev.x) >= abs(min_dev.x)) {
    if (is.null(x.min)) x.min <- mx - abs(max.dev.x)
    if (is.null(x.max)) x.max <- max(d.gen$x)
  }
  
  # normal density curve, no plot
  xx <- seq(x.min, x.max, length=200)
  lw  <- 2
  d.nrm <- dnorm(xx,mean(x),sd(x))

  # max y coordinate for graph
  max.y <- max(max(d.nrm), max(d.gen$y), max(h$density))

  # set margins
  margs <- .plt.marg(0, y.lab=NULL, x.lab, main, sub, lab_x_cex=lab_x_cex)
  lm <- margs$lm
  tm <- margs$tm
  rm <- margs$rm
  bm <- margs$bm

  orig.params <- par(no.readonly=TRUE)
  on.exit(par(orig.params))  
 
  if (y_axis) lm <- lm + 0.6  # allow extra room for densities
  par(bg=getOption("window_fill"))
  par(mai=c(bm, lm, tm, rm))

  
  # set up plot area
  plot(h, border="transparent", freq=FALSE,
     xlim=c(x.min,x.max), ylim=c(0,max.y),
     axes=FALSE, ann=FALSE, xlab=NULL, ylab=NULL, main=NULL, ...)

  # axis, axis ticks
  if (!y_axis)
    .axes(x.lvl=NULL, y.lvl=NULL, axTicks(1), NULL,
          rotate_x=rotate_x, rotate_y=rotate_y, offset=offset, ...)
  else
    .axes(x.lvl=NULL, y.lvl=NULL, axTicks(1), axTicks(2),
          rotate_x=rotate_x, rotate_y=rotate_y, offset=offset, ...)

  # axis value labels
  if (!y_axis) y.lab="" else y.lab="Density"
  .axlabs(x.lab, y.lab=NULL, main.lab, sub.lab,
          xy_ticks=TRUE, offset=offset, ...) 

  # colored background for plotting area
  usr <- par("usr")
  rect(usr[1], usr[3], usr[2], usr[4], col=col.bg, border=col.box,
    lwd=getOption("panel_lwd"), lty=getOption("panel_lty"))
  
  # plot the histogram
  if (histo)
    plot(h, add=TRUE, freq=FALSE, col=fill_hist, border="transparent")

  # plot the normal curve
  if (type == "normal" || type == "both") {
    if (fill_nrm == "transparent") lw <- 1.35 else lw <- 1
    polygon(c(x.min,xx,x.max), c(0,d.nrm,0), col=fill_nrm, 
            border=col.nrm, lwd=lw)
  }

  # plot the general curve
  if (type == "general" || type == "both") {
    lw <- ifelse (fill_gen == "transparent", 1.35, 1)
    polygon(d.gen, col=fill_gen, border=col.gen, lwd=lw)
  }

  if (band) rug(x, col=color_rug, lwd=size_rug)
 
  # plot the optional bar about a chosen point for general curve only
  if (!is.null(x.pt)  &&  type == "general") {
    d <- d.gen
    y.pt <- d$y[which.min(abs(x.pt-d$x))]
    xbeg <- x.pt - 0.5
    xend <- x.pt + 0.5
    xsub <- d$x[d$x > xbeg & d$x < xend]
    ysub <- d$y[d$x > xbeg & d$x < xend]
    txt <- paste("Density =", .fmt(y.pt,3), "at x =", .fmt(x.pt, 2))
    title(main=txt)
    polygon(c(xbeg, xsub, xend), c(0, ysub, 0), col="lightsteelblue")
    if (min(x) > 0) left <- -2*min(x) else left <- 2*min(x) 
    lines(c(left,x.pt), c(y.pt,y.pt), col="darkblue", lwd = 0.5)
    lines(c(x.pt,x.pt), c(0,y.pt), col="darkblue", lwd = 0.5)
  }

  # text output
  tx=""
# if (!quiet) {

    tx <- character(length = 0)
    if (getOption("suggest")) {
  
      # function call for suggestions
      fncl <- .fun_call.deparse(fncl) 
      fncl <- gsub(")$", "", fncl)  # get function call less closing ) 
      fncl <- gsub(" = ", "=", fncl)

      tx[length(tx)+1] <- ">>> Suggestions"
        tx[length(tx)+1] <- "bin_width: set the width of each bin"
      txt <- "  # histogram only"
      tx[length(tx)+1] <- paste("Histogram(", getOption("xname"),
         ")", txt, sep="")      
      txt <- "  # Violin/Box/Scatterplot (VBS) plot"
      tx[length(tx)+1] <- paste("Plot(", getOption("xname"), ")", txt,
         sep="")      
    }
    txsug <- tx
    if (length(txsug) == 0) txsug <- ""

    tx <- character(length = 0)

    tx[length(tx)+1] <- paste( "--- Bandwidth ---",
           "     for general curve:", .fmt(d.gen$bw,4), "\n")

    W <- NA; p.val <- NA
    if (type == "normal" || type == "both") {
      digits_d <- 3
      if (n > 2 && n < 5000) {
        nrm <- shapiro.test(x)
        W <- .fmt(nrm$statistic,min(4,digits_d+1))
        p.val <- .fmt(nrm$p.value,min(4,digits_d+1))
        tx[length(tx)+1] <- " "
        tx[length(tx)+1] <- paste("Null hypothesis is a normal population")
        tx[length(tx)+1] <- paste(nrm$method, ":  W = ", W, ",  p-value = ",
          p.val, sep="")
      }
      else 
        tx[length(tx)+1] <- paste("Sample size out of range for Shapiro-Wilk",
                            " normality test.")
    }

  return(list(tx=tx, txsug=txsug, 
              bw=d.gen$bw, n=n, n.miss=n.miss, W=W, pvalue=p.val))

# }  # not quiet

}