Plotting Odds Ratios (aka a forrestplot) with ggplot2 –

Hi,

if you like me work in medical research, you have to plot the results of multiple logistic regressions every once in a while. As I have not yet found a great solution to make these plots I have put together the following short skript. Do not expect too much, it’s more of a reminder to my future self than some mind-boggling new invention. The code can be found below the resulting figure looks like this:

fig_1_odds

Here comes the code. It takes the model and optionally a title as an input and generates the above plot.

 

 

plot_odds<-function(x, title = NULL){
tmp<-data.frame(cbind(exp(coef(x)), exp(confint(x))))
odds<-tmp[-1,]
names(odds)<-c(‘OR’, ‘lower’, ‘upper’)
odds$vars<-row.names(odds)
ticks<-c(seq(.1, 1, by =.1), seq(0, 10, by =1), seq(10, 100, by =10))

ggplot(odds, aes(y= OR, x = reorder(vars, OR))) +
geom_point() +
geom_errorbar(aes(ymin=lower, ymax=upper), width=.2) +
scale_y_log10(breaks=ticks, labels = ticks) +
geom_hline(yintercept = 1, linetype=2) +
coord_flip() +
labs(title = title, x = ‘Variables’, y = ‘OR’) +
theme_bw()
}

 

P.s. I know about ggplots “annotation_logticks” but they messed up my graphics, also it is not very often that ORs span more than three orders of magnitude. If they do consider playing with ggplots function or update the line beginning with “ticks <- ” in the above example

Update 29-01-2013: I replaced the nasty ” as they resulted in some nasty copy-past errors…

Reverse research – Extracting data from graphs with WebPlotDigitizer

As long as we all are still waiting for these large-scale data sharing mechanisms in place we will have to use every bit we can can our hand on. For data that is only represented as a graph – instead of tables or in the text – this means that we need some kind of graph digitizer.

The best one I found is the webplotdigitizer developed by Ankit Rohatgi. While I generally prefer desktop versions of programs, because they support reproducible research – this tools is just too easy to use that a further search would be necessary. Together with the video tutorials this is a one-stop shop for your data-extraction needs.

Congratulations to the developer for this excellent example of one task one tool!

 

Open Science – The great leaps forward

For everyone out there who believes that science (i.e. the methods, data and results) should be more open, hear the good news: You are not alone! A whole series of papers published in the worlds most highly ranked journals is on your side – or are they?

Duncal Hull compiled an excellent list of resources.

My personal favorite is the following:

Screenshot of paywall for an editorial that is supposedly about open access

 

 

 

 

fig_3_spice

When Venn diagrams are not enough – Visualizing overlapping data with Social Network Analysis in R

I recently thought about ways to visualize medications and their co-occurences in a group of children. As long as you want to visualize up to  4 different medications you can simply use Venn diagrams. There is a very nice R-package to generate these kind of graphics for you (for a  description see: Chen and Boutros, 2011). But this is of little help here.

The problem I faced involved 29 different medications and 50 children. So my data was stored in a table with 29 columns – one for each medication – and 50 rows – one for each child, so that the cells indicate whether or not the child took the medication.

M <- matrix(sample(0:1, 1450, replace=TRUE, prob=c(0.9,0.1)), nc=29)

The Solution – Social Network Analysis

There are a several R-packages to analyze and visualize social network data – I will focus on “igraph” in this post. The problem I had was that I was not – and probably I am still not –  familiar with the concepts and nomenclature of this field. The key to using the data described above in terms of network analysis was understanding that such data is called an affiliation matrix, where individuals are affiliated with certain events. As “igraph” likes adjacency matrices, where every column and row represents a different node – in our case a medication. The diagonal gives the number of times a medication was given (more information can be found on Daizaburo Shizuka site).

We transform an affilition matrix into an adjacency matrix in R simply by:

adj=M%*%t(M)

Now we can make a first bare-minimum plot:

require(igraph)
g=graph.adjacency(adj,mode=”undirected”, weighted=TRUE,diag=FALSE)
summary(g)
plot(g, main=”The bare minimum”)

 

Adding information and spicing it up a notch

In all likelihood You want to add at least three kinds of  information:

  1. Labels for the nodes
  2. Size of the nodes to represent the total number of events, aka medications
  3. Size of the links to represent the overlap between medications

name<-sample(c(LETTERS, letters, 1:99), 29, replace=TRUE)
number<-diag(adj)*5+5
width<-(E(g)$weight/2)+1
plot(g, main=”A little more information”, vertex.size=number,vertex.label=name,edge.width=width)

 

The “igraph” package lets you adopt quite a few parameters so you should consult with the manual. I only changed some of the colors, layout, fonts, etc.

plot(g, main=”Spice it up a notch”, vertex.size=number, vertex.label=name, edge.width=width, layout=layout.lgl, vertex.color=”red”, edge.color=”darkgrey”, vertex.label.family =”sans”, vertex.label.color=”black”)

 


Here is just the code:

?View Code RSPLUS
require(igraph)
setwd("~/Desktop/")
 
# Generate example data
M <- matrix(sample(0:1, 1450, replace=TRUE, prob=c(0.9,0.1)), nc=29)
 
# Transform matrices
adj=M%*%t(M)
 
# Make a simple plot
g<-graph.adjacency(adj,mode="undirected", weighted=TRUE,diag=FALSE)
summary(g)
plot(g, main="The bare minimum")
 
# Add more information
name<-sample(c(LETTERS, letters, 1:99), 29, replace=TRUE)
number<-diag(adj)*5+5
width<-(E(g)$weight/2)+1
 
plot(g, main="A little more information", vertex.size=number,vertex.label=name,edge.width=width)
 
# Adjust some plotting parameters
plot(g, main="Spice it up a notch", vertex.size=number, vertex.label=name, edge.width=width, layout=layout.lgl, vertex.color="red", edge.color="darkgrey", vertex.label.family ="sans", vertex.label.color="black")

Easy online audio player: Mp3 player

Sometimes you would like to play an audio to your participants in an online survey. But not every type of audio is working perfectly in standard online survey software and sometimes you would just like to link to your audio data on your own server. Here you could use the Flash Mp3-Player.

There are different Versions of the player available. For my tasks even the mini player was sufficient so far. How to present an audio online with this player (and so showing it in an online survey as well) is described here.

Pro:
It is really easy and the given examples are very helpful.
The code generator on the website is an additional helpful tool.

Con:
From my point of view there are no cons.

test

Visualizing GIS data with R and Open Street Map

In this post I way to share with you some code to use Openstreetmap – maps as a backdrop for a data visualization. We will use the RgoogleMaps-package for R. In the following I will show you how to make this graph.



1. Download the map

I wanted to take a closer look at an area around my former neighborhood, which is in Bochum, Germany.

lat_c<-51.47393
lon_c<-7.22667
bb<-qbbox(lat = c(lat_c[1]+0.01, lat_c[1]-0.01), lon = c(lon_c[1]+0.03, lon_c[1]-0.03))

Once this is done, you can download the corresponding Openstreetmap tile with the following line.

OSM.map<-GetMap.OSM(lonR=bb$lonR, latR=bb$latR, scale = 20000, destfile=”bochum.png”)

2. Add some points to the graphic

Now your second step will most likely be adding points to the map. I choose the following two.

lat <- c(51.47393, 51.479021)
lon <- c(7.22667, 7.222526)
val <- c(10, 100)

As the R-package was mainly build for google-maps, the coordinates need to be adjusted by hand. I made the following functions, that take the min and max value from the downloaded map.

lat_adj<-function(lat, map){(map$BBOX$ll[1]-lat)/(map$BBOX$ll[1]-map$BBOX$ur[1])}
lon_adj<-function(lon, map){(map$BBOX$ll[2]-lon)/(map$BBOX$ll[2]-map$BBOX$ur[2])

Now you can add some points to the map. If you want them to mean anything it may be handy to specify an alpha-level and change some aspects of the points, e.g. size, color, alpha corresponding to some variable of interest.

PlotOnStaticMap(OSM.map, lat = lat_adj(lat, OSM.map), lon = lon_adj(lon, OSM.map), col=rgb(200,val,0,85,maxColorValue=255),pch=16,cex=4)

Here is the full code:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
require(RgoogleMaps)
 
#define the part of the world you want to plot. Here the area around my former home.
lat_c<-51.47393
lon_c<-7.22667
bb<-qbbox(lat = c(lat_c[1]+0.01, lat_c[1]-0.01), lon = c(lon_c[1]+0.03, lon_c[1]-0.03))
 
# download the tile from OSM
OSM.map<-GetMap.OSM(lonR=bb$lonR, latR=bb$latR, scale = 20000, destfile="bochum.png")
image(OSM.map)
#Add some coordinates
lat<- c(51.47393, 51.479021)
lon<- c(7.22667, 7.222526)
val <- c(0, 255)
 
#function to adjust the coordinates
lat_adj<-function(lat, map){(map$BBOX$ll[1]-lat)/(map$BBOX$ll[1]-map$BBOX$ur[1])}
lon_adj<-function(lon, map){(map$BBOX$ll[2]-lon)/(map$BBOX$ll[2]-map$BBOX$ur[2])}
 
PlotOnStaticMap(OSM.map, lat = lat_adj(lat, OSM.map), lon = lon_adj(lon, OSM.map), col=rgb(255,0, val,90,maxColorValue=255),pch=16,cex=4)
 
dev.print(jpeg,"test.jpeg", width=1204, height=644, units="px")
yoda

The four steps to publication-grade graphics in R

For many, the main reason to use R is to generate really good-looking or at least informative graphics. However, while it is easy to find information on how to make an individual plot, it can take some time to find out how to get them out into the world. Here is my four-step program to turning your plot into a graphic-file.

In the following I will use my present favorite plot from here as an example.

 

1. Set your options

R allows you to set many general options for your plots, e.g. the margins and whether or not a box should be drawn around most of which are the documentation here.

My favorites are:

  • mfrow: To combine several plots into one (not necessary for the exaple).
  • mar: To control the margins of the plot (not necessary for the exaple).
  • las: To rotate the axis-labels (not necessary for the example)
?View Code RSPLUS
par(mar=c(2,0,2,2))

2. Make your plot

Well this part is the most heterogeneous, just take a peek at the gallery to get some inspiration, or dive into ggplot2 for a very comprehensive graphic-framework that also helps you to add legends.

?View Code RSPLUS
pie(c(1,1), labels="", col=c("black", "white"), main="Your options according to Yoda", init.angle=90)

3. Add a legend

R has a built-in function to add legends. The full documentation can be found here,  The options I use almost every time are:

  • x,y: To tell R where to put the legend. Usually I use the  name for the location (e.g. “top left”), instead of x and y-coordinates.
  • legend: To add some descriptions for the colors/line-types/shadings.
  • fill: To select the colors or alternatively “lty” for the line-type
  • bty: To get rid of the box around the legend
?View Code RSPLUS
legend("right", c("do", "do not", "try"), fill=c("black", "white", "gold"), bty="n", cex=1.4)

4. Save it to a file

R has various options to save files, as documented here. I most often save them as png, as the file-size for tiffs is extremely large at the same quality. This allows you to set the options

  • filename: well something with a  ”.png” at the end
  • width and height: To control the scale and of the image.
  • units: To
  • resolution: To Journals love images with at least 300 dpi.
  • bg: To have non-transparent background simply use “white”.
?View Code RSPLUS
dev.print(png, "yoda.png", width=8, height=6, units="in", res=300, bg="white")

Enjoy the complete script

 

?View Code RSPLUS
par(mar=c(2,0,2,2))
pie(c(1,1), labels="", col=c("black", "white"), main="Your options according to Yoda", init.angle=90)
legend("right", c("do", "do not", "try"), fill=c("black", "white", "gold"), bty="n", cex=1.4)
dev.print(png, "yoda.png", width=8, height=6, units="in", res=300, bg="white")
Para muchos el principal motivo para usar R es generar gráficos de muy buen aspecto o al menos informativos. Encontrar la información necesaria para hacer un diagrama individual puede ser fácil pero averiguar como presentar los diagramas al mundo puede ser un largo proceso. Aquí les presento mi programa, compuesto de 4 pasos, para convertir sus diagramas en archivos gráficos.A continuación usare mi grafico favorito de aquí como un ejemplo.

1. Configura tus opciones

R te permite configurar muchas opciones generales para tus diagramas, por ejemplo los márgenes y si deseas o no que una “caja” rodee la mayor parte de la documentación  aquí
Mis favoritos son:* mfrow:  Para combinar varios diagramas en un (no es necesario para el ejemplo).
* mar: Para controlar los márgenes del diagrama (no es necesario para el ejemplo)
* las: Para rotar los ejes-etiquetas (no es necesario para el ejemplo).
?View Code RSPLUS
par(mar=c(2,0,2,2))

2. Haz tu diagrama

Esta parte es la mas heterogénea. Puedes visitar nuestra galería  para inspirarte o recorrer ggplot2  para ver una estructura grafica (graphic-framework) muy completa que además te permita agregar notas.

?View Code RSPLUS
pie(c(1,1), labels="", col=c("black", "white"), main="Your options according to Yoda", init.angle=90)

3. Agregar una nota

R ha incorporado funciones para agregar notas. Puedes encontrar la documentación completa aquí. Las opciones que uso generalmente son:

  • x,y:  Para decirle a R donde colocar la nota. Generalmente uso el nombre para la localización (por ejemplo: “top left”) en lugar de coordenadas x-y.
  • legend: Para agregar algunas descripciones sobre  colors/line-types/shadings (colores/tipo de líneas/sombreado).
  • fill: Para elegir los colores o “lty” para el tipo de línea.
  • bty: Para deshacerse de la “caja” alrededor de las notas.

legend(“right”, c(“do”, “do not”, “try”), fill=c(“black”, “white”, “gold”), bty=”n”, cex=1.4)

?View Code RSPLUS
legend("right", c("do", "do not", "try"), fill=c("black", "white", "gold"), bty="n", cex=1.4)

4. Guardarlo en un archivo

R tiene varias opciones para guardar archivos, tal como esta documentado aquí. Por lo general guardo los archivos en formato png ya que los archivos tiffs son extremadamente pesados y de igual calidad. Esto te permite configurar las siguientes opciones:

  • filename: El nombre del archivos con un “.png” al final
  • width and height: Para controlar la escala de la imagen.
  • units: Para
  • resolution: Los Journals adoran las imágenes con al menos 300 dpi.
  • bg: Para tener un fondo no transparente debes usar simplemente “white”.
?View Code RSPLUS
dev.print(png, "yoda.png", width=8, height=6, units="in", res=300, bg="white")

Disfruta el Script completo!

 

?View Code RSPLUS
par(mar=c(2,0,2,2))
pie(c(1,1), labels="", col=c("black", "white"), main="Your options according to Yoda", init.angle=90)
legend("right", c("do", "do not", "try"), fill=c("black", "white", "gold"), bty="n", cex=1.4)
dev.print(png, "yoda.png", width=8, height=6, units="in", res=300, bg="white")

The latest thing I like about Zotero – Zotpress wordpress-integration

I started working with Endnote about 8 years ago as part of my first job at uni. A big chunk of my work was sorting and filing papers for my Prof. Even though it was a great way to become familiar with all the different topics, I always dreamed about all the other things I could have done during these hours. However, as long as the gold-standard for time-efficient referencing was the ability to query the social science citation index from within Endnote and only if you had a pricey subscription – these were only dreams.

With Zotero these dreams come finally true. Zotero is a plugin for Firefox that helps you organize your academic literature. It can be installed even when you do not have administrator privileges as e.g. at your work-PC. I will not talk about all the nice features of Zotero as importing existing databases, quickly adding references while browsing, or sharing libraries with your colleagues, workmates and friends.

Today I want to focus on an addon Zotpress that allows you to insert your whole library, a collection or papers tagged witha specific keay-word on your wordpress blog.The advantages are:

  • You only have to update one library, which you can also use to write your next paper.
  • As you can only show papers authored by specific authors, you can easily manage and genrate overlapping citation lists for a whole workgroup that specializes in different topics.

How to use Zotpress

To use it, you first need to have Zotero, add some references to your library, and a wordpress blog.

1. Install zotpress

Just go to the “plugin” menu in your dashboard and search for “zotpress”. Follow the white rabbit and you are done!

2. Connect to your Zotero library

This is the only step that took me some time, as I could not figure out how to generate a “public key”. This can be done right here on the zotero webpage. On this page, you can also see your userID. These two pieces of information then have to be entered in the wordpress-zotpress plugin page. Click on zotpress and then “Zotero Accounts”. And your done!

3. Add a bibliography to a post

There are many ways to add a bibliography to a post. You can either use in text citations that will also be correctly formatted according to APA, or you can use the Zotpress shortcode to enter a bibliography at a specific place.

The easiest way to test your setup is to just put a whole collection into you post. To find out the name of a specific collection go to the “citations” tab in the zotpress plugin preferences. If you enter the following line in brackets, it will print your collection – in this case my 2011 publications – sorted by date of publication. You can also select only five articles from the collection (by adding limit = “5″), or only a specific year (by adding year = “2010″).

zotpress collection=”COLLECTION_ID” sortby=”date” year = “2011″

Sorry, there's no items to display.

 

Simulation studies in R – Using all cores and other tips

After working more seriously with simulations I noticed some updates were necessary to my previous setup. Most notably are the following three:

  • It is very handy to explicitly call the different scenarios instead of using nested loops
  • Storing intermediate results in single files obliviates the need to rerun an almost finished but crashed analysis and seperates very clearly the data-generation from analysis part.
  • Using all availible cores can speed up the processing time, but may render the simulation not reproducible.

So here is my new simulation-study sceleton, that consists of five parts:

  1. Praeamble: Load all the functions that are required
  2. Simulation-function: This is the part, that will most likely be much more complicated in your case. Define the steps that will be repeated for different scenarios. The parameters of this function will be filled in by the scencarios.
  3. Scenario-Description: Explicitly show the range of values that should be passed to the Simulation-Function
  4. Run the analysis: Here you pass all the scenario-descriptions to your simulation-function. Either do this on one or all availible cores. In any case you should set a random seed to make the simulatino reproducible.
  5. Analyze the outputs: Not shown here but You propabely

Here is the complete script:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
# 1. Praeamble
setwd("c:/temp/")
require(doSNOW)
require(rlecuyer)
 
# 2. Simulation-function
sim_fun<-function(a,b,c){
results<-matrix(NA, 1000,4)
for(i in 1:1000){
#a=1; b=2;c=3;i=1
results[i,1:3]<-cbind(a, b, c)
results[i,4]<-mean(rnorm(100))#THIS MAY BE MORE COMPLEX FOR YOU HEHE!
}
write.table(results, file=paste(a,"_", b,"_", c, "_res.csv"))  
}
 
# 3. Scenario-Description
a<-seq(10, 100, 20)
b<-seq(20, 100, 30)
c<-seq(30, 200, 40)
scenarios<-expand.grid(a, b, c)
 
# 4.a Run the analysis on one core
set.seed(29012001)
for(i in 1:length(scenarios[,1])){sim_fun(scenarios[i,1], scenarios[i,2], scenarios[i,3])}
 
 
# 4.b Run the analysis on all availible Cores
cluster<-makeCluster(4, type = "SOCK")
clusterSetupRNG(cluster, seed = 29012001) 
registerDoSNOW(cluster)
 
foreach(i= 1:length(scenarios[,1])) %dopar% {sim_fun(scenarios[i,1], scenarios[i,2], scenarios[i,3])}
 
# compare the time
system.time(
for(i in 1:length(scenarios[,1])){sim_fun(scenarios[i,1], scenarios[i,2], scenarios[i,3])}
)
 
system.time(
foreach(i= 1:length(scenarios[,1])) %dopar% {sim_fun(scenarios[i,1], scenarios[i,2], scenarios[i,3])}
)

There are also other tutorials on how to run simulations in R. The one I liked most was Roger Koenkers’ “A simple protocoll for simulations in R” (accessible here) that relies more heavily on R’s built in features to solve some of the problems.

 

Keep ssh-scripts running without a connection – “screen” for Linux

I started offloading some of the heavy-lifting and the more time-consuming jobs to my server. While I am more and more comfortable of using the ssh terminal, the mayor problem was that scripts only worked as long as I was logged in. Which goes completely against the whole idea of laying off work to the server. After some goggling I found the small tool “screen” which is available for linux. A more comprehensive introduction can be found here, but after installing I only used three commands: “screen” to start; “screen -ls” to list the open terminals after reconnection, and “screen -r XXX” to reconnect to a running terminal.

Installation

With my Linux (ubuntu) server all I had to do to install it was a simple
$ sudo apt-get install screen

Use

After that you can start screen by simply typing:
$ screen

Now start your time-consuming job and just quit the terminal whenever you can’t stand waiting any longer. Whenever you want to have a look again just connect via ssh, and get a list of open screens:

$ screen -ls

this gave me:
3467.pts-0.MYSERVER (09/14/2011 08:59:46 PM) (Detached)
2238.pts-0.MYSERVER (09/14/2011 08:45:13 PM) (Attached)
2 Sockets in /var/run/screen/S-gerhi.

To reconnect to the first of these still running terminals, you just enter:
$ screen -r 3467.pts-0.MYSERVER

 

Dropbox – deine Dateien, überall

dropbox is a Web-based Service, that enables you in having access to your files at different places and in sharing them with others. It uses file synchronisation between different computers and users and represents an online data backup.

To use dropbox you have to install the dropbox client. Afterwards a new folder (the dropbox) will be generated in which you can store your files. You have acces to these data by other computers (on which you installed dropbox) as well as you can retrieve them from the Dropbox-Website.

Pro

  • You can create new folders in the dropbox, which you can release for your friends or associates (who have installed dropbox themselves). Now you can share and exchange data or work on projects together.
  • The file synchronisation runs on a server, who saves the data as well. So the involved computers don´t have to be in operation at the same time. The data will be updated, as soon as an involved Computer goes online.
  • Dropbox enables the manual adjusting of the bandwidth, so it won´t engross your hole internet connection.
  • by the defined folders “photos” and “public” you can share files with persons who are not using dropbox
  • Dropbox works with the following operation systems: Windows, Mac, Linux, iPad, iPhone, Android and BlackBerry.

Cons

  • The free version contains an initial storage capacity of 2GB. But by readvertising new members or other actions (like proofing you are a member of an university by your uni emailadress) you can increase your storage capacity.
  • The server performs a AES256-encoding on the files. The users are not able to attach their own code. But there is the opportunity e.g. to synchronize files with Dropbox, that are encoded by TrueCrypt.

Dump MySQL to CSV using R

Based on a related post on one of my favorite python-lists I remembered, that I wrote a similar snipplet some time ago.

So if you want to dump your whole MySQL database to csv-files you can recycle the following code:

?Download mysql2cvs.R
1
2
3
4
5
6
7
8
9
require(RMySQL)
m<-MySQL()
summary(m)
con<-dbConnect(m, dbname = "YOURDB", host="localhost", port=8889, user="YOURUSER", pass="YOURPASS", unix.sock="/Applications/MAMP/tmp/mysql/mysql.sock") # in case you are using MAMP 
tables<-dbListTables(con)
 
for (i in 1 : length(tables)){
temp<-(dbReadTable(con, tables[i]))
write.table(temp, tables[i], row.names=F)}

This also is a great way to use my new source-code plugin (WP-CodeBox)
Enjoy!