# Fancy Rugs in Regression Plots

Rugplots along the axes show the distribution of the underlying data in regression model plots. This is particulary useful in connection with additive (nonparametric) models where the plotted smooth function is the exclusive representation of the model in order to assess how much data contributed to the model fit at the different values of the exlanatory variable.

The custom plot.gam() function includes the possibility of such rugs and pointwise conficence intervalls by default.

Adding quartiles to the rugs requires some customization, though. I included the complete code to produce the above plot underneath.

The example for this GAM is borrowed from the excellent book of Alain Zuur et. al. Mixed Effects Models and Extensions in Ecology with R p.55ff. The ISIT data to run the code above is included in the R package AED which can be downloaded from the books website.

Note: the package AED is needed for the example dataset only. It is NOT necessary to use the example code on ones own dataset.

The main points concerning the rugs and quantile lables on the x-axis are:

1. Plot the coordinate system without lables
2. plot( ... , axes = FALSE )

3. Plot the x-axis with x-lables
4. axis(side = 1 , line = 0.3 , at = 0:5*1000 , tick = TRUE)

5. Plot rugs – the jitter() is necessary since a lot of datapoints sit on the same values of SampleDepth – so shake them a bit.
6. axis(side = 1 , line = -0.9 , at = jitter(ISIT$SampleDepth) , labels = F , tick = T , tcl = 0.8 , lwd.ticks = 0.1 , lwd = 0) 7. Print the lables “1Q”, “Median” and “3Q” or whatever you like to call them on the right position. The line and padj parameter set the position of the text and cex.axis the textsize. 8. axis(side = 1 , line = -0.8 , at = fivenum(ISIT$SampleDepth)[2:4], lwd = 0 , tick = F , labels = c("1Q","median","3Q"), cex.axis = 0.7, col.axis = "black" , padj = -2.8)

9. Plot thick tickmarks crossing through the rug cloud at 1Q, median and 3Q
10. axis(side = 1 , line = -0.8 , at = fivenum(ISIT$SampleDepth)[2:4], lwd = 0 , tick = T, tcl = 1.1 , lwd.ticks = 1 , col.ticks = "black", labels = FALSE) 11. and finally short thick tickmarks under the text touching the x-axis 12. axis(side = 1 , line = 0.3, at = fivenum(ISIT$SampleDepth)[2:4], lwd = 0 , tick = T, tcl = 0.2 , lwd.ticks = 1 , col.ticks = "black", labels = FALSE)

Here goes the complete code:
library(mgcv)
library(AED)
data(ISIT)
#
# Fit a univariate GAM model
model <- gam(Sources ~ s(SampleDepth) , data = ISIT)
fit <- predict(model, se = T)$fit se <- predict(model, se = T)$se.fit
lcl <- fit - 1.96 * se
ucl <- fit + 1.96 * se
#
# open a jpeg
jpeg("FancyRugs.jpg" , width=400, height=400)
#
# set plotting options: 1 plot per page, horizontal labels and textsize
par(mfrow = c(1,1) , las = 1 , cex = 1)
#
# plot coordinatesystem and labels
plot(0 , bty = "n" , type = "n" , xlim = c(0,5000) , ylim = c(-10,50) , xlab = "Depth (m)" , ylab = expression(paste("Number of sources (" , m^-3 , ")")) , axes = FALSE)
#
title(main="Association between number of sources of\nbioluminescent organisms and ocean depth" , cex.main = 0.8)
#
## _____ X-AXIS ______
# x-axis values
axis(side = 1 , line = 0.3 , at = 0:5*1000 , tick = TRUE)
#
# rugs at datapoints
axis(side = 1 , line = -0.9 , at = jitter(ISIT$SampleDepth) , labels = F , tick = T , tcl = 0.8 , lwd.ticks = 0.1 , lwd = 0) # # labels at 1Q, median and 3Q axis(side = 1 , line = -0.8 , at = fivenum(ISIT$SampleDepth)[2:4], lwd = 0 , tick = F , labels = c("1Q","median","3Q"), cex.axis = 0.7, col.axis = "black" , padj = -2.8)
#
# tick marks at 1Q, median and 3Q
axis(side = 1 , line = 0.3, at = fivenum(ISIT$SampleDepth)[2:4], lwd = 0 , tick = T, tcl = 0.2 , lwd.ticks = 1 , col.ticks = "black", labels = FALSE) # axis(side = 1 , line = -0.8 , at = fivenum(ISIT$SampleDepth)[2:4], lwd = 0 , tick = T, tcl = 1.1 , lwd.ticks = 1 , col.ticks = "black", labels = FALSE)
#
## _____ Y-AXIS ______
# y-axis values
axis(side = 2 , at = 0:5*10)
#
# rugs at datapoints
axis(side = 2 , line = -0.9 , at = jitter(ISIT$Sources) , labels = F , tick = T , tcl = 0.8 , lwd.ticks = 0.1 , lwd = 0) # # labels at 1Q, median and 3Q axis(side = 2 , line = -0.7 , at = fivenum(ISIT$Sources)[2:4], lwd = 0 , tick = F , labels = c("1Q","median","3Q"), cex.axis = 0.7, col.axis = "black")
#
# thicker tick marks at 1Q, median and 3Q
axis(side = 2 , line = 0.3, at = fivenum(ISIT$Sources)[2:4], lwd = 0 , tick = T, tcl = 0.3 , lwd.ticks = 1 , col.ticks = "black", labels = FALSE , padj = -2) axis(side = 2 , line = -0.7 , at = fivenum(ISIT$Sources)[2:4], lwd = 0 , tick = T, tcl = 1.1 , lwd.ticks = 1 , col.ticks = "black" , labels = FALSE)
#
# horizontal line marking the intercept = mean(Sources) (for univariate model only)
abline(h=mean(ISIT$Sources), lty=3) # # Scatterplot lines(ISIT$SampleDepth , ISIT$Source , type = "p" , cex = 0.4 , lwd = 0.2 , col = "grey") # # plot main figure lines(ISIT$SampleDepth[order(ISIT$SampleDepth)] , fit[order(ISIT$SampleDepth)] , col = "black" , lwd = 2)
#
# plot lower confidence limit (lcl)
lines(ISIT$SampleDepth[order(ISIT$SampleDepth)] , lcl[order(ISIT$SampleDepth)] , col = "grey" , lwd = 1) # # plot upper confidence limit (ucl) lines(ISIT$SampleDepth[order(ISIT$SampleDepth)] , ucl[order(ISIT$SampleDepth)] , col = "grey" , lwd = 1)
#
# closing the jpg file
dev.off()

# Install R Package XML with Ubuntu amd64

install.packages("XML")

failed while complaining

Cannot find xml2-config

Seth Falcon pointed out library libxml2-dev – Development files for the GNOME XML library – where missing.

sudo aptitude install libxml2-dev

fixed the problem.

# Installing RODBC with Ubuntu/Debian amd64

Trying to install RODBC in Ubuntu with
install.packages("RODBC")
failed throwing an error message

A glance at the r-help showed that it had to do with something called unixODBC – an ODBC driver manager.

The package was installed, but not the development package and thus not the headers which R complained about.

Again a not-so-obvious-for-the-newbee-Linux-Unix-Shell-goblish thing. The fix is
sudo aptitude install unixodbc-dev

# Open Access .mdb Files with RODBC

Getting data into R can be done by reading colon separated files (.csv) via the read.table() function. It is also possible to access databases directly and send SQL queries directly from R to the database. This has some advantages: Using Sweave the queries get documented in the analysis report, variable formats are retained.

To install the RODBC package:
install.packages("RODBC")

Open a database connection to an Microsoft Access database file, e.g. “MyDataBase.mdb” sitting in the Folder “C:\ MyPath\MyDataBase.mdb”:
channel <- odbcConnectAccess("C:/MyPath/MyDataBase")
note that the Windows backslashes “\” become slashes “/” in R and the extension “.mdb” is omitted.

Getting the database table “MyTable” into the R dataframe “R.Table”
R.Table <- sqlQuery( channel , paste ("select * from MyTable"))
MyTable can also be a sql query in the Access database.

# Create LaTex tables from R output

Creating reasonable layouted LaTeX tables from R output was easier then expected. I should have googled it long ago…

install.packages("xtable")

Lets say you created a tabular output in R called “tab1”, e.g by doing:
data(CO2)
tab1 <- with(CO2, table(Treatment , Type))
tab1

your text output in R would look like

Type
Treatment Quebec Mississippi
nonchilled 21 21
chilled 21 21

Now you would like this or whatever table or data.frame as a nice LaTeX-table, the only thing to do is:
library(xtable)
xtable(tab1)

and the output will be:

% latex table generated in R 2.9.0 by xtable 1.5-5 package
% Wed Jul 08 16:20:54 2009
\begin{table}[ht]
\begin{center}
\begin{tabular}{rrr}
\hline
& Quebec & Mississippi \\
\hline
nonchilled & 21 & 21 \\
chilled & 21 & 21 \\
\hline
\end{tabular}
\end{center}
\end{table}

If you are using Sweave the usage becomes
<< echo = FALSE , results = tex >>
library(xtable)
data(CO2)
with(CO2, xtable(table(Treatment , Type)))
@

and a the result of the R run is a LaTeX document. Another post will give a hint how to paste R graphics into the same Sweave or LaTeX document… later…

The full usage of xtable() is
xtable(x, caption=NULL, label=NULL, align=NULL, digits=NULL, display=NULL)
xtable(table , caption = "My table caption")
and labels
xtable(table , label = " MyLaTeXlable")
to the LaTeX tables.

# Install R packages without internet connection

You might have no internet connection or as in my case you have one, but the firewall prevents anything beside Windows IE to access the net. Grrr.

Ok. It takes more time, but is possible. Lets have an example:
Installing the package ‘gam’ for R-2.9 on Windows being in Norway (with cran.ii.uib.no as the nearest and fastest “Comprehensive R Archive Network”.

1. Find another place with internet connection or as in my case, the browser works, so use the browser.
2. Find the package on the repository corresponding to your operating system and your version of R in our above example here
3. klick on the link gam_1.0.zip, download it and save it (e.g. using an USB disk). Remember the path to the file (in our example in the folder H:/DATA/).
4. Open R and write install.packages("H:/DATA/gam_1.0.zip", repos=NULL) at the command line.

Thats it.

# Generate LaTeX tables of descriptive statistics

Using a LOT of time on fiddling together tables of descriptive statistics manually in R did not inspire me to look in the CRAN repositories for a R function doing exactly this… yes a bit stubborn… want to find out myself…

So what I had to do was:

1. handcode the descriptives for each variable, something like prop.table(xtabs( ~ Variable + Group1 + Group2), c(2,3)),
2. calculate some statistics for each variable (chi^2, one-way-ANOVA, Fisher-exact..) and attach it to the descriptives,
3. run the code,
4. copy the plain text output from R back into the editor and finally
5. add LaTeX or CSV syntax to get the desired results.

Anyway this is over now. A solution exists, of course (thx Kjetil… again):

The package ‘reporttools’ of Kaspar Rufibach.
install.packages("reporttools", dependecies=TRUE)

The most important functions are well:

tableContinuous( vars = c(bmi, ejectionfraction, systolicBP, diastolicBP) , group = sex , subset = significantstenosis , print.pval = "anova")

and

tableNominal() for nominal variables.

The output is LaTeX and it is possible to specify table captions and lables to the tables in the function call. I will give it a try inside Sweave… but first I have to get it to my Linux machine … which is blocked by the corporate firewall … after a recent Mircrosoft powergrab … after a Virusattack (Conficker) … after running XP unpached.

But thats another story

# Difficulties installing rgl

rgl is a R package for three-dimensional visualisation using OpenGL. The package provides functions implementing a new graphics device suitable for visualisation of R objects in three dimensions using the OpenGL libraries.

It can be installed from the Ubuntu repostitories with

sudo aptitude install r-cran-rgl

or on all platforms inside R with

install.packages("rgl")

BUT:

It depends on some GL libraries installed, which do not get installed by default. So you might expect unsuccessful installation, with an error message like

The solution is to install the missing library manually with

sudo aptitude install libglu1-mesa-dev

… one of those inconveniences preventing mainstream users switching to OpenSource software – it seems to me.

# How2plot nicer GAM curves

Generalized additive models visualize potential non-linear associations between a predictor and a response variable.

The default plotting method produces clean plots for all covariates in the model (or a selection) but: They do not have presentation quality by any means in terms of colors, understandable axes-labels, scaling, etc.

This example shows a customization variant for a additive regression model with one covariate. The goal was to display the absolute value of the response variable on the y-axis and not the “difference from intercept” which is default.

This is only meaningful for a single covariate in the model.

## 1 The default plot.gam() method

MyGAM1<- with(MyData[MyData$Strata==1,], gam(Y ~ s(Covariate))) MyGAM0<- with(MyData[MyData$Strata==0,], gam(Y ~ s(Covariate)))

par( mfcol=c(1,2))
plot(MyGAM0)
plot(MyGAM1)

This is the resulting plot:

## 2 The fancy way

1. Extract the values of the model response from the GAM object:

response1 <- predict(MyGAM1, type="response", se.fit=T)
response0 <- predict(MyGAM0, type="response", se.fit=T)

2. Print the response values against the covariate (note: this works just with one covariate)

par(mfcol=c(1,1))

plot(0, type="n", bty="n", main="Fancy GAM plot", xlab="MyCovariate", ylab="MyResponse", lwd=3,ylim=c(0,60), xlim=c(0,200))
legend("bottomright", bty="n", lwd=5, col=c("green","red"), legend=c("Strata = 0", "Strata = 1"))

lines(sm.spline(MyGAM1$model$Covariate , response1$fit) , lwd = 3 , col = "red") lines(sm.spline(MyGAM1$model$Covariate , response1$fit+1.96*response1$se) , lty = 3 , lwd = 2 , col = "red") lines(sm.spline(MyGAM1$model$Covariate , response1$fit-1.96*response1$se) , lty = 3 , lwd = 2 , col = "red") lines(sm.spline(MyGAM0$model$Covariate , response0$fit) , lwd = 3 , col = "green")
lines(sm.spline(MyGAM0$model$Covariate, response0$fit + 1.96 * response0$se) , lty = 3 , lwd = 2, col = "green")
lines(sm.spline(MyGAM0$model$Covariate, response0$fit - 1.96 * response0$se) , lty = 3 , lwd = 2 , col = "green")

abline(h=gam.dm1$coefficients[1], lty=2, lwd=1, col="red") abline(h=gam.dm0$coefficients[1], lty=2, lwd=1, col="green")

Update:I have written a much more detailed static page about the additive COX model: http://rforge.org/plothr/
The page has a download link to the function plotHR() which does all the fuzz. It is extensively commented. It should be easy to understand the syntax and modify it for individual purposes.

Therneau et al. refer to the proportional hazards model or COX-regression model as “the workhorse of regression analysis for censored data”. They show how to implement the additive form of this model in SAS and S-pluss; already mentioned by Hastie and Tibshirany in 1986 when introducing Generalized Additive Models (GAM).

I found modelling the functional form of the covariates in a regression model for rightcensored survival times with smoothing splines extremely useful. And the implementation is absolutely straightforward in R.

The only thing needed is the installation of the R-libraries “survival” and “pspline”:

install.packages("pspline")
and
install.packages("survival")

In the following code I will refer to a dataset “MyData” with a binary status variable “death” and a time-to-event variable “days2death”.
The status variable “death” should be (not necessarily) 1 if the event of interesst occured to the subject and “days2death” gives then the time to this event.

Viualizing the functional form of a covariate takes the following steps:

1. create the survival object of interesst
2. fit a proportional hazards model with smoothing splines,
3. predict the functional form of the covariate of interesst and
4. plot it!

Note that there is the termplot() function in R which gives you the GAM plots after the modelfit, so step 3 would not be necessary – BUT: it has a bug and fails plotting a single covariate; and it does not allow all to much customizing.

This is the R code to achieve the analysis:

## 1 Create survival object:

surv.death <- Surv(MyData$days2death, MyData$death)

## 2 Fit proportional hazards model with smoothing splines for continuous covariates:

library(survival)
library(pspline)
pham.fit <- coxph( surv.death ~ pspline(EF, df=4) + pspline(Age, df=4) + strata (Sex, df=4) , data = MyData)

The model above includes the continuous covariates “EF” (ejection fraction) and “Age” and stratifies for “Sex”.

## 3 Produce the fitted smoothing spline for the first covariate in the above model formula with standard errors

predicted <- predict(pham.fit , type = "terms" , se.fit = TRUE , terms = 1)
“terms=1” refers to “pspline(EF,df=4)”

## 4 Plot it

First plotting axes and labels
plot(0 , xlab="Ejection Fraction" , ylab = "Hazard Ratio" , main = "All-cause Death" , type = "n" , xlim=c(0,100) , ylim=c(0,3))
the range of values on the x-axis (“xlim=c(0,100)”) is chosen manually for this specific covariate; of course it is possible to use something like ylim = c( 0 , max(MyData$EF) ). Now plot the fitted smoothing spline using the lines() function: lines( sm.spline(MyData$EF , exp(predicted$fit)) , col = "red" , lwd = 0.8) Note that the term prediction gives log-hazard-ratios; therefore exp(predicted$fit) is plotted against the values of the covariate. The sm.spline() function is necessary since the points of the plot appear in random order and density, according to the underlying dataset; a plain lines() function would produce just a chaotic pattern. Alternative:
plot(MyData$EF , exp(predicted$fit) , col = "red" , cex = 0.2)
produces a scattered plot that reflects the distribution of the underlying data – I do prefer adding a rug-plot on the bottom of the graph to illustrate this (see under).

… upper and lower confidence limits with dashed thinner lines

lines(sm.spline(MyData$EF , exp(predicted$fit + 1.96 * predicted$se)) , col = "orange" , lty = 2 , lwd = 0.4) and lines(sm.spline(MyData$EF , exp(predicted$fit - 1.96 * predicted$se)) , col = "orange" , lty = 2 , lwd = 0.4)

… a tiny horizontal line at hazard level 1, do see where the confidence limits cross:
abline( h = 1 , col = "lightgrey" , lty = 2 , lwd = 0.4)

… tiny tickmarks on the x-axes to reflect the distribution of the underlying data:
axis( side = 1 , at = MyData$EF, labels = F , tick = T , tcl = 0.4 , lwd.ticks = 0.1) … and some fancy red tickmarks to mark minimum, lower hinge, median, upper hinge and maximum of the covariate in the dataset: axis( side = 1 , at = fivenum(MyData$EF), labels = F , tick = T , tcl = -0.2 , lwd.ticks = 1 , col.ticks = "red")

Thats it!

## 4b) The easy way (works ONLY with MORE then 1 continous covariate) – predicting the terms can be omitted:

termplot(pham.fit, se=T, rug=T)

Resulting in …