Tag Archives: R

Evolution of a logistic regression

In my last post I showed how one can easily summarize the outcome of a logistic regression. Here I want to show how this really depends on the data-points that are used to estimate the model. Taking a cue from the evolution of a correlation I have plotted the estimated Odds Ratios (ORs) depending on the number of included participants. The result is bad news for those working with small (< 750 participants) data-sets.

“eval_reg” Function to estimate model parameters for subsets of data

eval_reg<-function(model){
mod<-model
dat<-mod$data[sample(nrow(mod$data)),]
vars<-names(coef(mod))
est<-data.frame(matrix(nrow=nrow(dat), ncol=length(vars)))
pb <- txtProgressBar(min = 50, max = nrow(dat), style = 3)

for(i in 50:nrow(dat)){
try(boot_mod<-update(mod, data=dat[1:i,]))
try(est[i,]<-exp(coef(boot_mod)))
setTxtProgressBar(pb, i)
}
est$mod_nr<-1:length(dat[,1])
names(est)<-c(vars, ‘mod_nr’)
return(est)
}

As I randomized the order of data you can run it again and again to arrive at an even deeper mistrust as some of the resulting permutations will look like they stabilize earlier. On the balance you need to set the random-number seed to make it reproducible.

Run and plot the development

set.seed(29012001)

mod_eval<-eval_reg(gp_mod)

tmp<-melt(mod_eval,id=’mod_nr’)
tmp2<-tmp[tmp$variable!='(Intercept)',]

ticks<-c(seq(.1, 1, by =.1), seq(0, 10, by =1), seq(10, 100, by =10))

ggplot(tmp2, aes(y=value, x = mod_nr, color = variable)) +
geom_line() +
geom_hline(y=1, linetype=2) +
labs(title = ‘Evolution of logistic regression’, y = ‘OR’, x = ‘number of participants’) +
scale_y_log10(breaks=ticks, labels = as.character(ticks)) +
theme_bw()

Update 29-01-2013:

I added my definition of the ticks on the log-scale. The packages needed are ggplot2 and reshape.

Plotting Odds Ratios (aka a forrestplot) with ggplot2 –

8 Replies

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:

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…

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

5 Replies

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:

Labels for the nodes
Size of the nodes to represent the total number of events, aka medications
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")

Visualizing GIS data with R and Open Street Map

5 Replies

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:

^?Download visualizing_data_on_osm.R

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")

The four steps to publication-grade graphics in R

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")

Simulation studies in R – Using all cores and other tips

1 Reply

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:

Praeamble: Load all the functions that are required
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.
Scenario-Description: Explicitly show the range of values that should be passed to the Simulation-Function
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.
Analyze the outputs: Not shown here but You propabely

Here is the complete script:

^?Download simulation_in_r_update.R

# 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.

Dump MySQL to CSV using R

GIS on a shoestring – Getting traveltimes from google

2 Replies

The analysis of geospatial information is currently a big trend in medicine and public health. Even though some may want to convince you that this can only be achieved with the latest and most expensive software, I am not convinced. First, analysis of spatial data dates back to at least 1856 when John Snow investigated Cholera-outbreaks in London. Second, as I try to demonstrate today some very interesting analysis and data can be retreived essentially for free.

While I have already made a post on how to plot freely availible geospatial data in R in a previous post , this post will show you how to use Python to access the google maps database and gather e.g. travel times and distances to/from various locations with known zip-codes.

Please note that this is my first Python skript. So it will certainly not meet the high standards you might have developed based on previous posts. On the up-side, you will get the baby step instructions.

Update 2011/07/03: A much more user-friendly version of the script that adds guis to select a proper csv-file, containing start and end-adressess and to store the results can be found here. If you are afraid of Python, you can use the stand-alone Mac app “batchtimer” that basically contains all files necessary from here.

Continue reading →

Visualizing small-scale paired data – combining boxplots, stripcharts, and confidence-intervals in R

2 Replies

Sometimes when working with small paired data-sets it is nice to see/show all the data in a structured form. For example when looking at pre-post comparisons, connected dots are a natural way to visualize which data-points belong together. In R this can be easily be combined with boxplots expressing the overall distribution of the data. This also has the advantage of beeing more true to non-normal data that is not correctly represented by means +/- 95%CI. I have not come up against a good tutorial of how to do such a plot (although the right hand plot borrows heavily from this post on the excellent R-mailing list), so in the post you will find the code to generate such a graph in R.

Continue reading →

RStudio the missing link between your brain and statistics

Simulation studies in R – Reproducing MacCullums et al. 2002 “Effects of variable dichotomization”

2 Replies

[en]
I recently came across an excellent paper “On the Practice of Dichotomization of Quantitiative Variables” by MacCallum and colleagues (2002) . As I use ANOVAs a lot in my research, it really got me thinking about the whole issue. Even though I have no great idea for an innovative simulation study, you might have one. If you read through this post, you will notice that it’s really simple – at least the technical part.

I will only explain the two-variable scenario. But the setup is basically the same for more complex, e.g. two-variable, setups. Let’s start with their small numerical example before turning to the simulation study.

Continue reading →

Visualizing geographic data with R

1 Reply

[de]

Viele Daten, die man erhebt, sind örtlich gebunden. In R kann man solche Daten sehr schön mit dem Package sp darstellen. Alles, was man braucht sind:

Kartenmaterial: Also Beschreibungen der Kanten. Eine sehr gute Datenbank hierfür ist die GADM database of Global Administrative Areas. Für Deutschland gibt es dort eine Karte aller Kreise, die als R-Datenobjekt gespeichert sind – nett oder?
Daten, die einen lokalen Bezug haben. Viele gibt es in Genesis Datenbank des Statistikportals dem gemeinsamen Datenangebot der Statistischen Bundesämter. Ich hab’ mir vorgenommen einfach mal die Anzahl der Krankenhausbetten zu plotten. Hier ist das Ergebnis:

Continue reading →

Sustainable Research

Improving research by promoting Free and Open Source Software

Tag Archives: R

Evolution of a logistic regression

Plotting Odds Ratios (aka a forrestplot) with ggplot2 –

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

The Solution – Social Network Analysis

Adding information and spicing it up a notch

Visualizing GIS data with R and Open Street Map

1. Download the map

2. Add some points to the graphic

The four steps to publication-grade graphics in R

1. Set your options

2. Make your plot

3. Add a legend

4. Save it to a file

Enjoy the complete script

1. Configura tus opciones

2. Haz tu diagrama

3. Agregar una nota

4. Guardarlo en un archivo

Disfruta el Script completo!

Simulation studies in R – Using all cores and other tips

Dump MySQL to CSV using R

GIS on a shoestring – Getting traveltimes from google

Visualizing small-scale paired data – combining boxplots, stripcharts, and confidence-intervals in R

RStudio the missing link between your brain and statistics

Simulation studies in R – Reproducing MacCullums et al. 2002 “Effects of variable dichotomization”

Continue reading →

Visualizing geographic data with R