# public psychemedia / bowelCancer.Rnw Created 2011-11-01

Sweave script documenting the production of a funnel plot

bowelCancer.Rnw
 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 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 \documentclass[a4paper]{article}\SweaveOpts{echo=FALSE, keep.source=TRUE}\usepackage{a4wide}\usepackage{color}\usepackage{hyperref} \begin{document}\section{Example of self-documenting data journalism notes} This is an example of using Sweave to combine code and output from the R statistical programming environment and the LaTeX document processing environment to generate a self-documenting script in which the actual code used to do stats and generate statistical graphics is displayed along the charts it directly produces. \subsection{Getting Started...}The aim is to try to replicate a graphic included by Ben Goldacre in his article \emph{DIY statistical analysis: experience the thrill of touching real data}\footnote{\url{http://www.guardian.co.uk/commentisfree/2011/oct/28/bad-science-diy-data-analysis}}.  << echo = T >>=# The << echo = T >>= identifies an R code region;# echo=T means run the code, and print what happens when it's run # In the code area, lines beginning with a # are comment lines and are not executed #First, we need to load in the XML library that contains the scraper functionlibrary(XML) #Now we scrape the tablesrcURL='http://www.guardian.co.uk/commentisfree/2011/oct/28/bad-science-diy-data-analysis'cancerdata=data.frame( readHTMLTable( srcURL, which=1, header=c('Area','Rate','Population','Number') ) ) #The @ symbol on its own at the start of a line marks the end of a code block@ The format is simple: {\tt readHTMLTable(url,which=TABLENUMBER)} ({\tt TABLENUMBER} is used to extract the N'th table in the page.) The header part labels the columns (the data pulled in from the HTML table itself contains all sorts of clutter). We can inspect the data we've imported as follows: << echo = T >>=#Look at the whole table (the whole table is quite long,# so don't disply it/comment out the command for now instead.#cancerdata#If you are using RStudio, you can inspect the data using the command: View(cancerdata))#Look at the column headersnames(cancerdata)#If there's a problem with the header command not doing it's stuff, # uncomment the next two lines to fix it and check the fix...#names(cancerdata)<-c('Area','Rate','Population','Number')#names(cancerdata)#Look at the first 10 rowshead(cancerdata)#Look at the last 10 rowstail(cancerdata)#What sort of datatype is in the Number column?class(cancerdata$Number)@ The last line, {\tt class(cancerdata\$Number)}, identifies the data as type \emph{factor}. In order to do stats and plot graphs, we need the Number, Rate and Population columns to contain actual numbers. (Factors organise data according to categories; when the table is loaded in, the data is loaded in as strings of characters; rather than seeing each number as a number, it's identified as a category.) The   << echo=T >>=#Convert the numerical columns to a numeric datatypecancerdata$Rate = as.numeric(levels(cancerdata$Rate)[as.numeric(cancerdata$Rate)])cancerdata$Population = as.numeric(levels(cancerdata$Population)[as.integer(cancerdata$Population)])cancerdata$Number = as.numeric(levels(cancerdata$Number)[as.integer(cancerdata$Number)])#Just check it worked…class(cancerdata$Number)class(cancerdata$Rate)class(cancerdata$Population)head(cancerdata)@ We can now plot the data as a simple scatterplot using the {\tt plot} command (figure \ref{fig:simpleplot}) or we can add a title to the graph and tweak the axis labels (figure \ref{fig:simpleplot2}). \begin{figure}\begin{center}<>=#Plot the Number of deaths by the Populationplot(Number ~ Population, data=cancerdata)@\end{center}\caption{Vanilla scatter plot}\label{fig:simpleplot}\end{figure} \begin{figure}\begin{center}<>=#Plot the Number of deaths by the Population.#Add in a title (main) and tweak the y-axis label (ylab).plot(Number ~ Population, data=cancerdata, main='Bowel Cancer Occurrence by Population', ylab='Number of deaths')@\end{center}\caption{Vanilla scatter plot}\label{fig:simpleplot2}\end{figure} The {\tt plot} command is great for generating quick charts. If we want a bit more control over the charts we produce, the {\tt ggplot2} library is the way to go. (ggplot2 isn't part of the standard R bundle, so you'll need to install the package yourself if you haven't already installed it. In RStudio, find the Packages tab, click Install Packages, search for ggplot2 and then install it, along with its dependencies...). You can see the sort of chart ggplot creates out of the box in figure \ref{fig:simpleggplot}. \begin{figure}\begin{center}<>=require(ggplot2)#Plot the Number of deaths by the Populationp=ggplot(cancerdata)+geom_point(aes(x=Population, y=Number))print(p)@\end{center}\caption{A rather prettier plot}\label{fig:simpleggplot}\end{figure} \newpage\subsection{Generating the Funnel Plot} Doing a bit of searching for the funnel plot'' chart type used to display the data in Goldacre's article, I came across a post on Cross Validated, the Stack Overflow/Stack Exchange site dedicated to statistics related Q\&A: \emph{How to draw funnel plot using ggplot2 in R?}\footnote{\url{http://stats.stackexchange.com/questions/5195/how-to-draw-funnel-plot-using-ggplot2-in-r/5210\#5210}} The meta-analysis answer seemed to produce the similar chart type, so I had a go at cribbing the code, with confidence limits set at the 95\% and 99.9\% levels. Note that I needed to do a couple of things: \begin{enumerate}\item work out what values to use where! I did this by looking at the ggplot code to see what was plotted. p was on the y-axis and should be used to present the death rate. The data provides this as a rate per 100,000, so we need to divide by 100, 000 to make it a rate in the range 0..1. The x-axis is the population.\item change the range and width of samples used to create the curves\item change the y-axis range.\end{enumerate} You can see the result in figure \ref{fig:simpleggplot}. \begin{figure}\begin{center}<>=#TH: funnel plot code from:#stats.stackexchange.com/questions/5195/how-to-draw-funnel-plot-using-ggplot2-in-r/5210#5210#TH: Use our cancerdatanumber=cancerdata$Population#TH: The rate is given as a 'per 100,000' value, so normalise itp=cancerdata$Rate/100000 p.se <- sqrt((p*(1-p)) / (number))df <- data.frame(p, number, p.se, Area=cancerdata\$Area) ## common effect (fixed effect model)p.fem <- weighted.mean(p, 1/p.se^2) ## lower and upper limits for 95% and 99.9% CI, based on FEM estimator#TH: I'm going to alter the spacing of the samples used to generate the curvesnumber.seq <- seq(1000, max(number), 1000)number.ll95 <- p.fem - 1.96 * sqrt((p.fem*(1-p.fem)) / (number.seq))number.ul95 <- p.fem + 1.96 * sqrt((p.fem*(1-p.fem)) / (number.seq))number.ll999 <- p.fem - 3.29 * sqrt((p.fem*(1-p.fem)) / (number.seq))number.ul999 <- p.fem + 3.29 * sqrt((p.fem*(1-p.fem)) / (number.seq))dfCI <- data.frame(number.ll95, number.ul95, number.ll999, number.ul999, number.seq, p.fem) ## draw plot#TH: note that we need to tweak the limits of the y-axisfp <- ggplot(aes(x = number, y = p), data = df) +geom_point(shape = 1) +geom_line(aes(x = number.seq, y = number.ll95), data = dfCI) +geom_line(aes(x = number.seq, y = number.ul95), data = dfCI) +geom_line(aes(x = number.seq, y = number.ll999, linetype = 2), data = dfCI) +geom_line(aes(x = number.seq, y = number.ul999, linetype = 2), data = dfCI) +geom_hline(aes(yintercept = p.fem), data = dfCI) +xlab("Population") + ylab("Bowel cancer death rate") + theme_bw() #Automatically set the maximum y-axis value to be just a bit larger than the max data valuefp=fp+scale_y_continuous(limits = c(0,1.1*max(p))) #Label the outlier pointfp=fp+geom_text(aes(x = number, y = p,label=Area),size=3,data=subset(df,p>0.0003)) print(fp)@\end{center}\caption{A rather prettier plot}\label{fig:funnelplot}\end{figure} \end{document}

Please sign in to comment on this gist.

Something went wrong with that request. Please try again.