#' 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)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.