Skip to content

Instantly share code, notes, and snippets.

@johncolby
Created June 1, 2011 20:16
Show Gist options
  • Save johncolby/1003205 to your computer and use it in GitHub Desktop.
Save johncolby/1003205 to your computer and use it in GitHub Desktop.
An example showing how to plot longitudinal data in R using base graphics and ggplot2
library(plyr)
library(ggplot2)
# Simulate data
numSubs = 40
x = runif(numSubs, 0, 30)
x = c(x, x + abs(rnorm(numSubs, 2)))
y = 2 + x - 0.02*x^2 + rnorm(2*numSubs, 0, 1)
data = data.frame(ID = 1:numSubs, Age = x, FA = y, Time = factor(rep(c(1,2), each=numSubs)), Gender = sample(c('Male', 'Female'), numSubs, replace=T))
# Fit quadratic model
quad.fit = lm(FA ~ Age + I(Age^2), data=data)
pred.ages = with(data, data.frame(Age = seq(min(Age), max(Age), length.out=100)))
quad.pred = pred.ages
quad.pred$FA = predict(quad.fit, quad.pred)
# Fit exponential model
exp.fit = nls(FA ~ SSasymp(Age, Asym, R0, lrc), data=data)
exp.pred = pred.ages
exp.pred$FA = predict(exp.fit, exp.pred)
# Plot with base graphics
dev.new(width=5, height=5)
colorpicks = c('pink', 'skyblue')
with(data, plot(Age, FA, pch=c(1,19)[Time], col=colorpicks[as.numeric(Gender)], bty='l', main='Longitudinal example', xlab='Age', ylab='FA'))
d_ply(data, 'ID', function(x) lines(x$Age, x$FA, lty=3, col='#666666'))
lines(quad.pred)
lines(exp.pred, lty=2)
legend('bottomright', bty='n', col=c('black', 'black', colorpicks), pch=c(1,19,19,19), c('Time 1', 'Time 2', 'Female', 'Male'))
# Plot with ggplot2
dev.new(width=5, height=4)
p = ggplot(data, aes(x=Age, y=FA)) +
geom_point(aes(color=Gender, shape=Time)) +
geom_line(aes(group=ID), alpha=0.3, linetype=3) +
geom_smooth(aes(group=1), color='black', method='lm', formula=y~x+I(x^2)) +
geom_smooth(aes(group=1), color='black', linetype=2, method='nls', formula=y~SSasymp(x, Asym, R0, lrc), se=F) +
scale_color_manual(values=colorpicks) +
scale_shape_manual(values=c(1,19)) +
opts(title='Longitudinal example')
p
# Facets with ggplot2
dev.new(width=8, height=4)
(p = p + facet_grid(facets=.~Gender))
# For more complex summary geoms, and other arbitrary annotations, the plyr package can be used to put together custom data frames with the info needed for plotting
ExpFit <- function(df) {
fit = nls(FA ~ SSasymp(Age, Asym, R0, lrc), data=df)
pred = with(df, data.frame(Age = seq(min(Age), max(Age), length.out=100)))
pred$FA = predict(fit, pred)
cbind(pred, t(coefficients(fit)), xpos=range(df$Age)[2], ypos=range(df$FA)[1])
}
exp.fits = ddply(data, 'Gender', ExpFit)
dev.new(width=8, height=4)
p + geom_text(aes(group=1, x=max(xpos), y=min(ypos), label=paste('list(FA[infinity]==', format(Asym, dig=4), ',FA[0]==', format(R0, dig=4), ',tau==', format(exp(lrc), dig=4), ')')), data=exp.fits[!duplicated(exp.fits$Asym),], parse=T, hjust=1, vjust=0, size=3.5)
# Or as an equation...
dev.new(width=8, height=4)
p + geom_text(aes(group=1, x=max(xpos), y=min(ypos), label=paste('FA==', format(Asym, dig=4), format(R0-Asym, dig=4), '*italic(e)^{-', format(exp(lrc), dig=4), '*Age}')), data=exp.fits[!duplicated(exp.fits$Asym),], parse=T, hjust=1, vjust=0, size=3.5)
@johncolby
Copy link
Author

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment