############# Frieman test ############# DrugResponse <- matrix(c(5.40, 5.50, 5.55, 5.85, 5.70, 5.75, 5.20, 5.60, 5.50, 5.55, 5.50, 5.40, 5.90, 5.85, 5.70, 5.45, 5.55, 5.60, 5.40, 5.40, 5.35, 5.45, 5.50, 5.35, 5.25, 5.15, 5.00, 5.85, 5.80, 5.70, 5.25, 5.20, 5.10, 5.65, 5.55, 5.45, 5.60, 5.35, 5.45, 5.05, 5.00, 4.95, 5.50, 5.50, 5.40, 5.45, 5.55, 5.50, 5.55, 5.55, 5.35, 5.45, 5.50, 5.55, 5.50, 5.45, 5.25, 5.65, 5.60, 5.40, 5.70, 5.65, 5.55, 6.30, 6.30, 6.25), nrow = 22, byrow = TRUE, dimnames = list(1 : 22, c("DrugA", "DrugB", "DrugC"))) DrugResponse install.packages("reshape") library(reshape) DR2 <- melt(DrugResponse) out <- aov(value~X2, data = DR2) shapiro.test(resid(out)) boxplot(value~X2, data= DR2) friedman.test(DrugResponse) DRA <- DR2[DR2$X2 != "DrugA", ] t.test(value~X2, data= DRA, paired = TRUE) DRB <- DR2[DR2$X2 != "DrugB", ] t.test(value~X2, data= DRB, paired = TRUE) DRC <- DR2[DR2$X2 != "DrugC", ] t.test(value~X2, data= DRC, paired = TRUE) 0.05/3 friedman.test.with.post.hoc <- function(formu, data, to.print.friedman = T, to.post.hoc.if.signif = T, to.plot.parallel = T, to.plot.boxplot = T, signif.P = .05, color.blocks.in.cor.plot = T, jitter.Y.in.cor.plot =F) { # formu is a formula of the shape: Y ~ X | block # data is a long data.frame with three columns: [[ Y (numeric), X (factor), block (factor) ]] # Note: This function doesn't handle NA's! In case of NA in Y in one of the blocks, then that entire block should be removed. # Loading needed packages if(!require(coin)) { print("You are missing the package 'coin', we will now try to install it...") install.packages("coin") library(coin) } if(!require(multcomp)) { print("You are missing the package 'multcomp', we will now try to install it...") install.packages("multcomp") library(multcomp) } if(!require(colorspace)) { print("You are missing the package 'colorspace', we will now try to install it...") install.packages("colorspace") library(colorspace) } # get the names out of the formula formu.names <- all.vars(formu) Y.name <- formu.names[1] X.name <- formu.names[2] block.name <- formu.names[3] if(dim(data)[2] >3) data <- data[,c(Y.name,X.name,block.name)] # In case we have a "data" data frame with more then the three columns we need. This code will clean it from them... # Note: the function doesn't handle NA's. In case of NA in one of the block T outcomes, that entire block should be removed. # stopping in case there is NA in the Y vector if(sum(is.na(data[,Y.name])) > 0) stop("Function stopped: This function doesn't handle NA's. In case of NA in Y in one of the blocks, then that entire block should be removed.") # make sure that the number of factors goes with the actual values present in the data: data[,X.name ] <- factor(data[,X.name ]) data[,block.name ] <- factor(data[,block.name ]) number.of.X.levels <- length(levels(data[,X.name ])) if(number.of.X.levels == 2) { warning(paste("'",X.name,"'", "has only two levels. Consider using paired wilcox.test instead of friedman test"))} # making the object that will hold the friedman test and the other. the.sym.test <- symmetry_test(formu, data = data, ### all pairwise comparisons teststat = "max", xtrafo = function(Y.data) { trafo( Y.data, factor_trafo = function(x) { model.matrix(~ x - 1) %*% t(contrMat(table(x), "Tukey")) } ) }, ytrafo = function(Y.data){ trafo(Y.data, numeric_trafo = rank, block = data[,block.name] ) } ) # if(to.print.friedman) { print(the.sym.test) } if(to.post.hoc.if.signif) { if(pvalue(the.sym.test) < signif.P) { # the post hoc test The.post.hoc.P.values <- pvalue(the.sym.test, method = "single-step") # this is the post hoc of the friedman test # plotting if(to.plot.parallel & to.plot.boxplot) par(mfrow = c(1,2)) # if we are plotting two plots, let's make sure we'll be able to see both if(to.plot.parallel) { X.names <- levels(data[, X.name]) X.for.plot <- seq_along(X.names) plot.xlim <- c(.7 , length(X.for.plot)+.3) # adding some spacing from both sides of the plot if(color.blocks.in.cor.plot) { blocks.col <- rainbow_hcl(length(levels(data[,block.name]))) } else { blocks.col <- 1 # black } data2 <- data if(jitter.Y.in.cor.plot) { data2[,Y.name] <- jitter(data2[,Y.name]) par.cor.plot.text <- "Parallel coordinates plot (with Jitter)" } else { par.cor.plot.text <- "Parallel coordinates plot" } # adding a Parallel coordinates plot matplot(as.matrix(reshape(data2, idvar=X.name, timevar=block.name, direction="wide")[,-1]) , type = "l", lty = 1, axes = FALSE, ylab = Y.name, xlim = plot.xlim, col = blocks.col, main = par.cor.plot.text) axis(1, at = X.for.plot , labels = X.names) # plot X axis axis(2) # plot Y axis points(tapply(data[,Y.name], data[,X.name], median) ~ X.for.plot, col = "red",pch = 4, cex = 2, lwd = 5) } if(to.plot.boxplot) { # first we create a function to create a new Y, by substracting different combinations of X levels from each other. subtract.a.from.b <- function(a.b , the.data) { the.data[,a.b[2]] - the.data[,a.b[1]] } temp.wide <- reshape(data, idvar=X.name, timevar=block.name, direction="wide") #[,-1] wide.data <- as.matrix(t(temp.wide[,-1])) colnames(wide.data) <- temp.wide[,1] Y.b.minus.a.combos <- apply(with(data,combn(levels(data[,X.name]), 2)), 2, subtract.a.from.b, the.data =wide.data) names.b.minus.a.combos <- apply(with(data,combn(levels(data[,X.name]), 2)), 2, function(a.b) {paste(a.b[2],a.b[1],sep=" - ")}) the.ylim <- range(Y.b.minus.a.combos) the.ylim[2] <- the.ylim[2] + max(sd(Y.b.minus.a.combos)) # adding some space for the labels is.signif.color <- ifelse(The.post.hoc.P.values < .05 , "green", "grey") boxplot(Y.b.minus.a.combos, names = names.b.minus.a.combos , col = is.signif.color, main = "Boxplots (of the differences)", ylim = the.ylim ) legend("topright", legend = paste(names.b.minus.a.combos, rep(" ; PostHoc P.value:", number.of.X.levels),round(The.post.hoc.P.values,5)) , fill = is.signif.color ) abline(h = 0, col = "blue") } list.to.return <- list(Friedman.Test = the.sym.test, PostHoc.Test = The.post.hoc.P.values) if(to.print.friedman) {print(list.to.return)} return(list.to.return) } else { print("The results where not significant, There is no need for a post hoc test") return(the.sym.test) } } # Original credit (for linking online, to the package that performs the post hoc test) goes to "David Winsemius", see: # http://tolstoy.newcastle.edu.au/R/e8/help/09/10/1416.html } friedman.test.with.post.hoc(value~X2|X1, DR2) ################ Two way repeated measures ANOVA ######### acne <- read.csv("10_rmanova.csv",header = TRUE) acL <- reshape(acne, direction = "long", varying=3:6, sep ="") acL$group <- factor(acL$group) acL$id <- factor(acL$id) acL$time <- factor(acL$time) interaction.plot(acL$time, acL$group, acL$month) acD <- aov(month~group*time +Error(id), data=acL) summary(acD) acL0 <- acL[acL$time=="0",] acL1 <- acL[acL$time=="1",] acL3 <- acL[acL$time=="3",] acL6 <- acL[acL$time=="6",] t.test(month~group, data=acL0) t.test(month~group, data=acL1) t.test(month~group, data=acL3) t.test(month~group, data=acL6) ########## Corelation test ############## install.packages("UsingR") library(UsingR) data("galton") str(galton) View(galton) plot(child~parent, data =galton) cor.test(galton$child, galton$parent) out = lm(child~parent, data= galton) summary(out) # y = 23.94153 + 0.64629 x X abline(out, col ="red") plot(jitter(child,5)~jitter(parent,5), galton) sunflowerplot(galton) install.packages("SwissAir") library(SwissAir) data("AirQual") View(AirQual) str(AirQual) Ox <- AirQual [ ,c("ad.O3","lu.O3","sz.O3")] + AirQual[, c("ad.NOx","lu.NOx","sz.NOx")]- AirQual[, c("ad.NO", "lu.NO","sz.NO")] names(Ox) <- c("ad","lu","sz") plot(lu~sz, data= Ox) install.packages("hexbin") library(hexbin) bin = hexbin(Ox$lu, Ox$sz, xbins = 50) plot(bin) smoothScatter(Ox$lu, Ox$sz) install.packages("IDPmisc") library(IDPmisc) iplot(Ox$lu, Ox$sz) par(mfrow=c(1,1))