R/t.diff.R

In alemiani/explora: Explore and Visualize Data

t.diff <- function(dataset, group, bonferroni.correction = FALSE){
  nums <- sapply(dataset, is.numeric)
  H <- dataset[ , nums]
  drops <- as.character(substitute(group)[[3]])
  H[,c(drops)] <- list(NULL)

  if(is.factor(group)) {
    gg <- levels(as.factor(group))
    lev <- (as.numeric(as.factor(group)))
    name.g.1 <- paste(substitute(gg))[[1]]
    name.g.2 <- paste(substitute(gg))[[2]]
  } else {
    low.factor  <- paste('low.', substitute(group), sep = '')[[3]]
    high.factor <- paste('High.', substitute(group), sep = '')[[3]]
    name <- NA
    name[group < median(group, na.rm = T)] <- low.factor
    name[group > median(group, na.rm = T)] <- high.factor
    group <- as.factor(name)
    name.g.1 <- high.factor
    name.g.2 <- low.factor
  }


  df <- data.frame(
    N.g1           = unlist(lapply(H, function(x) td.test(x, group, verbose = FALSE)$N.g1)),
    mean.1         = round(unlist(lapply(H, function(x) td.test(x, group, verbose = FALSE)$mean.1)),2),
    sd.1           = round(unlist(lapply(H, function(x) td.test(x, group, verbose = FALSE)$sd.1)),2),
    N.g2           = unlist(lapply(H, function(x) td.test(x, group, verbose = FALSE)$N.g2)),
    mean.2         = round(unlist(lapply(H, function(x) td.test(x, group, verbose = FALSE)$mean.2)),2),
    sd.2           = round(unlist(lapply(H, function(x) td.test(x, group, verbose = FALSE)$sd.2)),2),
    var_statistic  = round(unlist(lapply(H, function(x) td.test(x, group, verbose = FALSE)$var.statistic)),2),
    var_p.value    = unlist(lapply(H, function(x) td.test(x, group, verbose = FALSE)$var.p.value)),
    t_statistic    = round(unlist(lapply(H, function(x) td.test(x, group, verbose = FALSE)$t.statistic)),2),
    t_p.value      = unlist(lapply(H, function(x) td.test(x, group, verbose = FALSE)$t.p.value)),
    effectSize_d   = unlist(lapply(H, function(x) td.test(x, group, verbose = FALSE)$effectSize_d))
  )
  df$var_p.value <- symnum(df$var_p.value,legend = F, corr = FALSE, na = FALSE,
                           cutpoints = c(0, 0.001, 0.01, 0.05, 0.1, 1), symbols = c("***", "**", "*", ".", " "))
  if(bonferroni.correction == TRUE) {
    N <- nrow(df)
    df$t_p.value <- df$t_p.value * N
    df$p.values <- df$t_p.value
    df$t_p.value <- symnum(df$t_p.value,legend = F, corr = FALSE, na = FALSE,
                           cutpoints = c(0, 0.001, 0.01, 0.05, 0.1, max(df$t_p.value, na.rm = TRUE)),
                           symbols = c("***", "**", "*", ".", " "))

  } else {
    df$p.values <- df$t_p.value
    df$t_p.value <- symnum(df$t_p.value,legend = F, corr = FALSE, na = FALSE,
                           cutpoints = c(0, 0.001, 0.01, 0.05, 0.1, max(df$t_p.value, na.rm = TRUE)),
                           symbols = c("***", "**", "*", ".", " "))

  }



  names(df)[names(df) == 'N.g1'] <- paste('N_', name.g.1, sep = '')
  names(df)[names(df) == 'N.g2'] <- paste('N_', name.g.2, sep = '')
  names(df)[names(df) == 'mean.1'] <- paste('M_', name.g.1, sep = '')
  names(df)[names(df) == 'mean.2'] <- paste('M_', name.g.2, sep = '')
  names(df)[names(df) == 'sd.1'] <- paste('sd_', name.g.1, sep = '')
  names(df)[names(df) == 'sd.2'] <- paste('sd_', name.g.2, sep = '')

  df <- df[order(-abs(df$effectSize_d)),]


  par(mar=c(10, 4, 2, 2) + 0.1)
  Title <- paste(name.g.1, 'vs' ,name.g.2, sep = ' ')
  plot(abs(df$effectSize_d),las=3, xaxt = "n", xlab = '', ylab = '')
  mtext(side = 2, line = -2, 'Effect size D')
  title(Title, line = -2)
  labels <- rownames(df[order(-abs(df$effectSize_d)),])
  lines(abs(df$effectSize_d))
  par(new = T)
  p.values <- df$p.values
  plot(p.values, axes=F, xlab=NA, ylab=NA, col = 'purple', ylim = c(0,0.2))
  abline(h = 0.05, col = 'purple', lty = 3, pch = 20)
  axis(side = 4)
  mtext(side = 4, line = -2, 'p.values')
  axis(1, at = 1:length(labels),labels = labels,las=3)
  df$p.values <- NULL
  par(mfrow=c(1,1))

  return(df)
}

# swiss$Factor.2 <- as.factor(c(rep('Factor1', 20), rep('Factor2', 27)))
# t.diff(dataset = swiss,
#        group = swiss$Factor.2,
#        bonferroni.correction = F)