JRR Tolkien wrote a story (sans hobbits) called ‘Leaf by Niggle’ that has always resonated with me. I offer you an imperfect, incomplete tree (but my roof is intact).
Suggestions for improvements are encouraged.
And here is the link tree for the document (for your easy reviewing of the offered content) :
In this post I showcase a nice bar-plot and a balloon-plot listing recommended Nutritional supplements , according to how much evidence exists for thier benefits, scroll down to see it(and click here for the data behind it) * * * * The gorgeous blog “Information Is Beautiful” recently publish an eye candy post showing a “balloon race” image (see a static version of the image here) illustrating how much evidence exists for the benefits of various Nutritional supplements (such as: green tea, vitamins, herbs, pills and so on) . The higher the bubble in the Y axis score (e.g: the bubble size) for the supplement the greater the evidence there is for its effectiveness (But only for the conditions listed along side the supplement).
There are two reasons this should be of interest to us:
This shows a fun plot, that R currently doesn’t know how to do (at least I wasn’t able to find an implementation for it). So if anyone thinks of an easy way for making one – please let me know.
The data for the graph is openly (and freely) provided to all of us on this Google Doc.
The advantage of having the data on a google doc means that we can see when the data will be updated. But more then that, it means we can easily extract the data into R and have our way with it (Thanks to David Smith’s post on the subject)
For example, I was wondering what are ALL of the top recommended Nutritional supplements, an answer that is not trivial to get from the plot that was in the original post.
In this post I will supply two plots that present the data: A barplot (that in retrospect didn’t prove to be good enough) and a balloon-plot for a table (that seems to me to be much better).
Barplot (You can click the image to enlarge it)
The R code to produce the barplot of Nutritional supplements efficacy score (by evidence for its effectiveness on the listed condition).
# loading the data
supplements.data.0 <- read.csv("http://spreadsheets.google.com/pub?key=0Aqe2P9sYhZ2ndFRKaU1FaWVvOEJiV2NwZ0JHck12X1E&output=csv")
supplements.data <- supplements.data.0[supplements.data.0[,2] >2,] # let's only look at "good" supplements
supplements.data <- supplements.data[!is.na(supplements.data[,2]),] # and we don't want any missing data
supplement.score <- supplements.data[, 2]
ss <- order(supplement.score, decreasing = F) # sort our data
supplement.score <- supplement.score[ss]
supplement.name <- supplements.data[ss, 1]
supplement.benefits <- supplements.data[ss, 4]
supplement.score.col <- factor(as.character(supplement.score))
levels(supplement.score.col) <- c("red", "orange", "blue", "dark green")
supplement.score.col <- as.character(supplement.score.col)
# mar: c(bottom, left, top, right) The default is c(5, 4, 4, 2) + 0.1.
par(mar = c(5,9,4,13)) # taking care of the plot margins
bar.y <- barplot(supplement.score, names.arg= supplement.name, las = 1, horiz = T, col = supplement.score.col, xlim = c(0,6.2),
main = c("Nutritional supplements efficacy score","(by evidence for its effectiveness on the listed condition)", "(2010)"))
axis(4, labels = supplement.benefits, at = bar.y, las = 1) # Add right axis
abline(h = bar.y, col = supplement.score.col , lty = 2) # add some lines so to easily follow each bar
Also, the nice things is that if the guys at Information Is Beautiful will update there data, we could easily run the code and see the updated list of recommended supplements.
Balloon plot So after some web surfing I came around an implementation of a balloon plot in R (Thanks to R graph gallery) There where two problems with using the command out of the box. The first one was that the colors where non informative (easily fixed), the second one was that the X labels where overlapping one another. Since there is no "las" parameter in the function, I just opened the function up, found where this was plotted and changed it manually (a bit messy, but that's what you have to do sometimes...)
Here are the result (you can click the image for a larger image):
And here is The R code to produce the Balloon plot of Nutritional supplements efficacy score (by evidence for its effectiveness on the listed condition). (it's just the copy of the function with a tiny bit of editing in line 146, and then using it)
Daniel Malter just shared on the R mailing list (link to the thread) his code for performing the Siegel-Tukey (Nonparametric) test for equality in variability. Excited about the find, I contacted Daniel asking if I could republish his code here, and he kindly replied “yes”. From here on I copy his note at full.
My goal in this post is to give an overview of Friedman’s Test and then offer R code to perform post hoc analysis on Friedman’s Test results. (The R function can be downloaded from here)
Preface: What is Friedman’s Test
Friedman test is a non-parametric randomized block analysis of variance. Which is to say it is a non-parametric version of a one way ANOVA with repeated measures. That means that while a simple ANOVA test requires the assumptions of a normal distribution and equal variances (of the residuals), the Friedman test is free from those restriction. The price of this parametric freedom is the loss of power (of Friedman’s test compared to the parametric ANOVa versions).
The hypotheses for the comparison across repeated measures are:
H0: The distributions (whatever they are) are the same across repeated measures
H1: The distributions across repeated measures are different
The test statistic for the Friedman’s test is a Chi-square with [(number of repeated measures)-1] degrees of freedom. A detailed explanation of the method for computing the Friedman test is available on Wikipedia.
Performing Friedman’s Test in R is very simple, and is by using the “friedman.test” command.
Post hoc analysis for the Friedman’s Test
Assuming you performed Friedman’s Test and found a significant P value, that means that some of the groups in your data have different distribution from one another, but you don’t (yet) know which. Therefor, our next step will be to try and find out which pairs of our groups are significantly different then each other. But when we have N groups, checking all of their pairs will be to perform [n over 2] comparisons, thus the need to correct for multiple comparisons arise. The tasks: Our first task will be to perform a post hoc analysis of our results (using R) – in the hope of finding out which of our groups are responsible that we found that the null hypothesis was rejected. While in the simple case of ANOVA, an R command is readily available (“TukeyHSD”), in the case of friedman’s test (until now) the code to perform the post hoc test was not as easily accessible. Our second task will be to visualize our results. While in the case of simple ANOVA, a boxplot of each group is sufficient, in the case of a repeated measures – a boxplot approach will be misleading to the viewer. Instead, we will offer two plots: one of parallel coordinates, and the other will be boxplots of the differences between all pairs of groups (in this respect, the post hoc analysis can be thought of as performing paired wilcox.test with correction for multiplicity).
R code for Post hoc analysis for the Friedman’s Test
The analysis will be performed using the function (I wrote) called “friedman.test.with.post.hoc”, based on the packages “coin” and “multcomp”. Just a few words about it’s arguments:
formu – is a formula object of the shape: Y ~ X | block (where Y is the ordered (numeric) responce, X is a group indicator (factor), and block is the block (or subject) indicator (factor)
data – is a data frame with columns of Y, X and block (the names could be different, of course, as long as the formula given in “formu” represent that)
All the other parameters are to allow or suppress plotting of the results.
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
}
Example
(The code for the example is given at the end of the post)
Let's make up a little story: let's say we have three types of wine (A, B and C), and we would like to know which one is the best one (in a scale of 1 to 7). We asked 22 friends to taste each of the three wines (in a blind fold fashion), and then to give a grade of 1 till 7 (for example sake, let's say we asked them to rate the wines 5 times each, and then averaged their results to give a number for a persons preference for each wine. This number which is now an average of several numbers, will not necessarily be an integer).
After getting the results, we started by performing a simple boxplot of the ratings each wine got. Here it is:
The plot shows us two things: 1) that the assumption of equal variances here might not hold. 2) That if we are to ignore the "within subjects" data that we have, we have no chance of finding any difference between the wines.
So we move to using the function "friedman.test.with.post.hoc" on our data, and we get the following output:
$Friedman.Test
Asymptotic General Independence Test
data: Taste by
Wine (Wine A, Wine B, Wine C)
stratified by Taster
maxT = 3.2404, p-value = 0.003421
$PostHoc.Test
Wine B - Wine A 0.623935139
Wine C - Wine A 0.003325929
Wine C - Wine B 0.053772757
The conclusion is that once we take into account the within subject variable, we discover that there is a significant difference between our three wines (significant P value of about 0.0034). And the posthoc analysis shows us that the difference is due to the difference in tastes between Wine C and Wine A (P value 0.003). and maybe also with the difference between Wine C and Wine B (the P value is 0.053, which is just borderline significant).
Plotting our analysis will also show us the direction of the results, and the connected answers of each of our friends answers:
Here is the code for the example:
source("https://www.r-statistics.com/wp-content/uploads/2010/02/Friedman-Test-with-Post-Hoc.r.txt") # loading the friedman.test.with.post.hoc function from the internet
### Comparison of three Wine ("Wine A", "Wine B", and
### "Wine C") for rounding first base.
WineTasting <- data.frame(
Taste = 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),
Wine = factor(rep(c("Wine A", "Wine B", "Wine C"), 22)),
Taster = factor(rep(1:22, rep(3, 22))))
with(WineTasting , boxplot( Taste ~ Wine )) # boxploting
friedman.test.with.post.hoc(Taste ~ Wine | Taster ,WineTasting) # the same with our function. With post hoc, and cool plots
If you find this code useful, please let me know (in the comments) so I will know there is a point in publishing more such code snippets...
Inviting other R lovers (as myself) to find fun things to do with this (and similar) dataset.
The Data
The data set is comprised of 51 rows, one for each US states with the two variables (columns):
CTR – The CTR means “Click Through Rate” and is from chitika data base and collected from over two random days in January (a total of 31,667,158 total impressions), and is from the full range of Internet users (they don’t have traditional demographic data – every impression is completely anonymous).
Percent of the population who graduated college.
Super basic analysis and plot
This data presents a stunning -0.63 correlation between the two measurements. Hinting that “The Educated are Harder to Advertise To” (as the original post suggested). The data can be easily visualized using a scatter plot:
Created using just a few lines of R code:
aa <- read.table("https://www.r-statistics.com/wp-content/uploads/2010/02/State_CTR_Date.txt", sep = "t", header = T)
aa[,2:3] <- aa[,2:3] * 100
plot(aa[,2] ~ aa[,3], sub = paste("Correlation: ", round(cor(aa[,2], aa[,3]), 2)),
main = "Scatter plot of %CTR VS %College_Grad per State",
xlab = "%College_Grad per State",
ylab = "%CTR per State"
)
abline(lm(aa[,2] ~ aa[,3]), col = "blue")
My conclusion from the analysis
I was asked in the comments (by Eyal) to add my own conclusions to the analysis. Does higher intelligence imply lower chances of clicking ads, my answer (under the present data) is simple "I don't know". The only real conclusion I can make of the data is that there might be a point in checking this effect in a more rigorous way (which I am sure is already being done).
What should we have done in order to know? When doing scientific research, we often ask ourselves how sure are we of our results. The rule of thumb for this type of question is called "the pyramid of evidence". It is a way to organize various ways of getting "information" about the world, in an hierarchy of reliability. Here is a picture of this pyramid:
We can see that the most reliable source is a systematic review of randomized controlled trials. In our case, that would mean having controlled experiments where you take groups of people with different levels of "intelligence" (how would you measure that?), and check their CTR (click through rates) on banner ads. This should be done in various ways, correcting for various confounders , and later the results and conclusions (from several such experiments) should be systematically reviewed by experts on the subject.
All of this should be done in order to make a real assessment of the underlying question - how does smarts effects banner clicking. And the reason we need all of this work is because of what is said in the title of the next section:
"Correlation does not imply causation" is a phrase used in science and statistics to emphasize that correlation between two variables does not automatically imply that one causes the other (though it does not remove the fact that correlation can still be a hint, whether powerful or otherwise). The opposite belief, correlation proves causation, is a logical fallacy by which two events that occur together are claimed to have a cause-and-effect relationship.
But a much clearer explenation of it was given by the following XKCD comic strip:
Next step: other resources to play with
The motivation for my post is based on this digg post trying to hint how Religiousness is connected to "negative" things such as crimes, poverty and so on. That post was based on the following links:
If someone is motivated, he/she can extract that data and combine it with the current provided data.
In conclusion: this simplistic dataset, combined with other data resources, provides opportunity for various fun demonstrations of pairs correlation plots and of nice spatial plots (of states colored by their matching variable). It is a good opportunity to emphasize (to students, friends and the like) that "Correlation does not imply causation!". And finally - If you are an R lover/blogger and feel like playing with this - please let me know 🙂 .
I just noticed a google buzz from Jeroen ooms, with a Youtube video titled “RApache Hello World + POST arguments + catching errors.”
In this ~7 min video tutorial, Jeroen shares with us:
How to write “Hello World” in a website using RApache.
How to extract arguments from a form submited by the website visitor (and then inserting it into an “rnorm” function so to control the output). And finally,
How to catch an error in case of an invalid argument on an R Web Application.
Thank you Jeroen for a very simple, step by step, tutorial:
In case you have a self hosted WordPress blog, and you wish to show your R code in it, how would you do it?
The simplest solution would be to just paste the code as plain text, which will look like this:
x <- rnorm(100, mean = 2, sd = 3) plot(x, xlab = “index”, main = “Example code”)
But if you would like to help our readers orient themselves inside your code by giving different colors to different commands in the code (a.k.a: syntax highlighting). So it would like something like this:
x <- rnorm(100, mean = 2, sd = 3) # Creating a vector
plot(x, xlab = "index", main = "Example code") # Plotting it
WP-Syntax provides clean syntax highlighting using GeSHi -- supporting a wide range of popular languages (including R). It supports highlighting with or without line numbers and maintains formatting while copying snippets of code from the browser.
But there is a problem. The current WP-Syntax version is using an old version of GeSHi, and only the newer version (currently GeSHi version 1.0.8.6) includes support for R syntax. In order to solve this I patched the plugin and I encourage you to download (the fixed version of) WP-Syntax from here, which will allow you to highlight your R code.
Usage
After installing (and activating) the plugin, in order to add R code to your post you will need to: 1) Only work in HTML mode (not the Visual mode). Or else, the code you will paste will be messed up. 2) Put your code between the <pre> tag, like this:
(Note: make sure that you rewrite the " - so it will work.)
<pre lang="rsplus" line="1">
...Your R code here...
</pre>
Final note: R Syntax highlight in other ways
If you wish to have R syntax higlight inside an HTML file, I encourage you can have a look at the highlight package, by Romain Francois.
If you want to higlight your R syntax inside wordpress.com, here is a blog post by Erik Iverson showing how to do that using Emacs.
p.s: If you have a blog in which you write about R, please let me know about it in the comments (Or just join R-bloggers.com) - I'd love to follow you 🙂
(The R code for Barnard’s exact test is at the end of the article, and you could also just download it from here, or from github)
Barnards exact test - p-value based on the nuisance parameter
About Barnard’s exact test
About half a year ago, I was studying various statistical methods to employ on contingency tables. I came across a promising method for 2×2 contingency tables called “Barnard’s exact test“. Barnard’s test is a non-parametric alternative to Fisher’s exact test which can be more powerful (for 2×2 tables) but is also more time-consuming to compute (References can be found in the Wikipedia article on the subject).
When comparing Fisher’s and Barnard’s exact tests, the loss of power due to the greater discreteness of the Fisher statistic is somewhat offset by the requirement that Barnard’s exact test must maximize over all possible p-values, by choice of the nuisance parameter, π. For 2 × 2 tables the loss of power due to the discreteness dominates over the loss of power due to the maximization, resulting in greater power for Barnard’s exact test. But as the number of rows and columns of the observed table increase, the maximizing factor will tend to dominate, and Fisher’s exact test will achieve greater power than Barnard’s.
About the R implementation of Barnard’s exact test
After finding about Barnard’s test I was sad to discover that (at the time) there had been no R implementation of it. But last week, I received a surprising e-mail with good news. The sender, Peter Calhoun, currently a graduate student at the University of Florida, had implemented the algorithm in R. Peter had found my posting on the R mailing list (from almost half a year ago) and was so kind as to share with me (and the rest of the R community) his R code for computing Barnard’s exact test. Here is some of what Peter wrote to me about his code:
On a side note, I believe there are more efficient codes than this one. For example, I’ve seen codes in Matlab that run faster and display nicer-looking graphs. However, this code will still provide accurate results and a plot that gives the p-value based on the nuisance parameter. I did not come up with the idea of this code, I simply translated Matlab code into R, occasionally using different methods to get the same result. The code was translated from:
Trujillo-Ortiz, A., R. Hernandez-Walls, A. Castro-Perez, L. Rodriguez-Cardozo. Probability Test. A MATLAB file. URL
My goal was to make this test accessible to everyone. Although there are many ways to run this test through Matlab, I hadn’t seen any code to implement this test in R. I hope it is useful for you, and if you have any questions or ways to improve this code, please contact me at [email protected]
Duncan Murdoch just posted a youtube video presenting an animation clip of a 3d rgl object.
Duncan even went further and wrote an explanation on how he made the video:
here are the steps I used: 1. Design a shape to be displayed, and then play with the animation functions to make it change over time. Use play3d to do it live in R, movie3d to write the individual frames of the movie to .png files. 2. Use the ffmpeg package (not an R package, a separate project at http://ffmpeg.org) to convert the .png files to an .mp4 file. The individual frames totalled about 1 GB; the compressed movie is about 45 MB. 3. Upload to Youtube. I’m not a musician, so I had to use one of their licensed background tracks, I couldn’t write my own. I spent a lot of time picking one and then adjusting the timing of the video to compensate. Each render/upload cycle at full resolution took about an hour and a half. It’s a lot faster to render in a smaller window with fewer frames per second, but it’s still tedious. It’s easier to synchronize if you actually have a copy of the music locally, but Youtube doesn’t let you download their music. So the timing isn’t perfect, but it’s good enough for me!
