R/PlotDiff.R

In gfaity/ReachStroke: Stroke reaching's Statistics

#' Compute difference between two data set
#'
#' Compute difference between Ref and ToComp
#' 1° left panel = Density of normality
#' 2° right panel = BoxPlot - parametric (t student)
#'
#' @param Ref         Datas of reference system
#' @param ToComp      Datas of system to compare
#' @param U           "%"
#' @param RefT        Name of reference system
#' @param ToCompT     Name of system to compare
#' @param Cond        Main title
#' @param V           Current Variable
#' @param VT          Current Variable Name
#' @param Unit        Unit to use
#' @param Sym.Main
#' @param Flagqqplot  TRUE to print qqplot of each variable
#'
#' @return qqplot (normality)
#'
#' @examples
#'  V <- Variables [iVariables]
#'  VT <- VariablesNames [iVariables]
#'  U <- "%"
#'  RefT <- RefNameToRead
#'  ToCompT <- ToCompNameToRead
#'  UnitCurr <- VariablesUnit[iVariables]
#'  Ref <- subset(D, App == RefName)[,V]
#'  ToComp <- subset(D, App == ToCompName)[,V]
#'  Cond <- paste(ToCompT, "VS", RefT)
#'  Flag$qqplot <- TRUE
#'
#'  PlotDiff(DRef, DToComp, U, RefT, ToCompT, Cond, V, VT, UnitCurr, Flag$qqplot)
#'
#' @export

PlotDiff <- function(Ref,
                     ToComp,
                     U,
                     RefT,
                     ToCompT,
                     Cond,
                     V,
                     VT,
                     UnitCurr,
                     Flagqqplot,
                     Sym.Main = "") {

  #################

  par.pty = par()$pty     # save current pty (to restore it in the end)
  par(pty = "s")          # make axis square

  #################

  #Initialisation
  Alpha.Risk = 0.05

  # display some information in the consols
  cat(sprintf("\n--- %s : %s \n", Cond, V))

  #################

  #1.A Test normality
  #Shapiro-Wilk Test
  #if SW$p.value < Alpha.Risk then the distribution isn't normal
  SW_Ref <- shapiro.test(Ref)
  # SW_Ref$p.value
  SW_ToComp <- shapiro.test(ToComp)
  # SW_ToComp$p.value

  if (SW_Ref$p.value > Alpha.Risk & SW_ToComp$p.value > Alpha.Risk) {
    normality = T
    SW_Result = "These results show an idea of normality."
  } else {
    normality = F
    SW_Result = "These results show a violation of normality."
  }

  if (SW_Ref$p.value > Alpha.Risk) {
    SW_Ref_mod = ">"
  } else {
    SW_Ref_mod = "<"
  }

  if (SW_ToComp$p.value > Alpha.Risk) {
    SW_ToComp_mod = ">"
  } else {
    SW_ToComp_mod = "<"
  }

  txt = paste("Datas were (formally) tested for violation of normality with Shapiro-Wilk Test :"
              ,"\n", V, RefT, ":", round(SW_Ref$p.value,3), SW_Ref_mod, Alpha.Risk
              ,"\n", V, ToCompT, ":", round(SW_ToComp$p.value,3), SW_ToComp_mod, Alpha.Risk
              ,"\n", SW_Result
              ,"\n")
  cat (txt)

  #################

  #1.B Test density (graphical normality)
  #Compare two distributions with sm
  #Compute some informations
  Ref_char = replicate(length(Ref), 0)
  ToComp_char = replicate(length(ToComp), 1)
  # Ref_char = replicate(length(Ref), "Ref")
  # ToComp_char = replicate(length(ToComp), "ToComp")
  All_char = c(Ref_char,ToComp_char)
  All_val = c(Ref,ToComp)
  All_dat = cbind(All_char,All_val)

  #################

  #plot density to identify normality

  #if Variable name contains "Delta" then replace it by greek letter
  #Delta means difference End - Beg
  if (grepl("Delta", VT, fixed = TRUE)) {
    newVT <- gsub("Delta","", VT, fixed=TRUE)
    Xlabel = bquote(Delta ~ .(newVT) ~ .("(") * .(UnitCurr) * .(")"))
    mainT = bquote(bold(Delta ~ .(paste(newVT, ":", Cond)) ))
  } else {
    Xlabel = bquote(.(VT) ~ .("(") * .(UnitCurr) * .(")"))
    mainT = bquote(bold(.(paste(VT,":", Cond)) ))
  }

  sm.density.compare(All_dat[,2], All_dat[,1]
    , lty = c(1,1)
    , col = c("red" , "blue")
    , xlab = Xlabel
    , ylab = "Density"
  )

  title(mainT)
  txt = paste("Normality SW =",normality,"\n")
  mtext(txt, side=4, cex = 0.8, line = 1.5)#at=0.05) #inset=c(-0.5,0))

  legend("topright", c(RefT,ToCompT)
     , lty = c(1,1)
     #, title = "Condition"
     , col = c("red" , "blue")
  )

  cat(paste("Normality SW =",normality,"\n"))

  #################

  if (Flagqqplot) {
    #1.C qqPlot to verify normality
    # Solution 1 - simple but 2 plots
    qqnorm(Ref, main = paste (VT, ":", RefT,"Q-Q Plot"))

    qqnorm(ToComp, main = paste (VT, ":", ToCompT,"Q-Q Plot"))

    #Solution 2 - all in 1 plot
    # Doesn't work: it print in the big A4 pdf page, not in the reserved space
    # #function to qqplot
    # gg <- data.frame(Dat = c(Ref,ToComp), Condition = c(Ref_char,ToComp_char))
    # Ref_char = replicate(length(Ref), RefT)
    # ToComp_char = replicate(length(ToComp), ToCompT)
    #
    # ###With ggplot
    # ggploting <- ggplot(gg, aes(sample = Dat, colour = factor(Condition))) +
    #   stat_qq() +
    #   stat_qq_line()
    #
    # title("Q-Q Plot")
    # print (ggploting)
  }

  #################

  #2.A Non parametric test
  #comparaison, matched samples (whithin) = Wilcoxon
  w = wilcox.test(Ref, ToComp, paired = T)
  #print(w)

  Ref.median = round(median(Ref),2)
  ToComp.median = round(median(ToComp),2)

  #Writing some conclusion
  if (w$p.value < Alpha.Risk) {
    conclusion = 'statistically significantly different'
    conclusionPlot = "Significative difference (NP)"
  }else{
    conclusion = 'not statistically significantly different'
    conclusionPlot = "No significative difference (NP)"
  }

  A = 'The report in APA A '
  txt = paste("\n",A, w$method,'indicated that the median test ranks were '
              ,conclusion, ', ranks Z = '
              ,w$statistic, ',p =',round(w$p.value,4)
              ,', median', V, RefT, ' = ', Ref.median
              ,', median', V, ToCompT, ' = ', ToComp.median
              ,'.'
              ,"\n")
  txt_boxplot = paste(conclusionPlot
                      ,"ranks Z = ", w$statistic, ', p =',round(w$p.value,4)
                      # ,"\n"
                      # ,'median', RefT, ' = ', Ref.median
                      # ,'median', ToCompT, ' = ', ToComp.median
  )
  cat(txt)

  #################

  if (normality == T) {
    #2.B Parametric test
    #comparaison, matched samples (whithin) = t test (paired)
    tp = t.test(Ref, ToComp, paired = T)
    # print(tp)

    Ref.mean = round(mean(Ref),2)
    ToComp.mean = round(mean(ToComp),2)
    Ref.sd = round(sd(Ref),2)
    ToComp.sd = round(sd(ToComp),2)

    #Writing some conclusion
    if (tp$p.value < Alpha.Risk) {
      conclusion = 'statistically significantly different'
      conclusionPlot = "Significative difference (P)"
    }else{
      conclusion = 'not statistically significantly different'
      conclusionPlot = "No significative difference (P)"

    }

    A = 'The report in APA A '
    txt = paste("\n",A,tp$method,'indicated that the mean tests were'
                , conclusion
                ,', t =', round(tp$statistic,2)
                ,', df =', tp$parameter
                ,', p =', round(tp$p.value,4)
                ,', mean', V, RefT, ' = ', Ref.mean
                ,', SD', V, RefT, ' = ', Ref.sd
                ,', mean', V, ToCompT, ' = ', ToComp.mean
                ,', SD', V, ToCompT, ' = ', ToComp.sd
                ,',95% confident interval = ', round(tp$conf.int[1],2), 'to', round(tp$conf.int[2],2)
                ,'.'
                ,"\n")
    txt_boxplot = paste(txt_boxplot, "\n"
                        ,conclusionPlot
                        ,"t = ", round(tp$statistic,2), ', p =',round(tp$p.value,4)
    )
    cat(txt)
  }

  #################

  #3. Creation of a boxplot to visualize differences

  #if Variable name contains "Delta" then replace it by greek letter
  #Delta means difference End - Beg
  if (grepl("Delta", VT, fixed = TRUE)) {
    newVT <- gsub("Delta","", VT, fixed=TRUE)
    Ylabel = bquote(Delta ~ .(paste(newVT)) ~ .("(") * .(UnitCurr) * .(")"))
  } else {
    Ylabel = bquote(.(VT) ~ .("(") * .(UnitCurr) * .(")"))
  }

  boxplot(
    Ref,
    ToComp,
    col = "blue",
    ylab =  Ylabel,
    xlab = paste (Cond, "\n", txt_boxplot)
  )


  main = paste("Size", RefT, "=", length(Ref)
    ,"\n"
    ,"Size", ToCompT, "=", length(ToComp)
  )
  #cex.main = 0.8
  #main = c(font.main,cex.main)
  title(main,cex.main=1,line=0.3)
  # mtext("%s : %s ",txt_boxplot, side=1)

  par(pty = par.pty)  # restore original pty

  #Return ggplot (qqplot)
  #return(ggploting)
  #grid.arrange(ggploting)
}