d.dat <- read.csv(file="initial-final-data.csv") # read in inital and final values ## Function adapted from Everitt B: An R and S-Plus Companion ## to Multivariate Analysis. London: Springer 2005. ## http://biostatistics.iop.kcl.ac.uk/publications/everitt/ bvbox<-function(a, d = 7, lin.col="black", lin.lwd="1.5") { # used for bivariate IQR (Figure 2B) # a is a data matrix # d is constant (usually 7) p <- length(a[1, ]) m1 <- mean(a[, 1]) m2 <- mean(a[, 2]) s1 <- sqrt(var(a[, 1])) s2 <- sqrt(var(a[, 2])) r <- cor(a[, 1:2])[1, 2] x <- (a[, 1] - m1)/s1 y <- (a[, 2] - m2)/s2 e <- sqrt((x * x + y * y - 2 * r * x * y)/(1 - r * r)) e2 <- e * e em <- median(e) emax <- max(e[e2 < d * em * em]) r1 <- em * sqrt((1 + r)/2) r2 <- em * sqrt((1 - r)/2) theta <- ((2 * pi)/360) * seq(0, 360, 3) xp <- m1 + (r1 * cos(theta) + r2 * sin(theta)) * s1 yp <- m2 + (r1 * cos(theta) - r2 * sin(theta)) * s2 r1 <- emax * sqrt((1 + r)/2) r2 <- emax * sqrt((1 - r)/2) theta <- ((2 * pi)/360) * seq(0, 360, 3) xpp <- m1 + (r1 * cos(theta) + r2 * sin(theta)) * s1 ypp <- m2 + (r1 * cos(theta) - r2 * sin(theta)) * s2 maxxl <- max(xpp) minxl <- min(xpp) maxyl <- max(ypp) minyl <- min(ypp) b1 <- (r * s2)/s1 a1 <- m2 - b1 * m1 y1 <- a1 + b1 * minxl y2 <- a1 + b1 * maxxl b2 <- (r * s1)/s2 a2 <- m1 - b2 * m2 x1 <- a2 + b2 * minyl x2 <- a2 + b2 * maxyl maxx <- max(c(a[, 1], xp, xpp, x1, x2)) minx <- min(c(a[, 1], xp, xpp, x1, x2)) maxy <- max(c(a[, 2], yp, ypp, y1, y2)) miny <- min(c(a[, 2], yp, ypp, y1, y2)) lines(xp, yp, lwd = lin.lwd, col = lin.col) } # makes plot similar to Figure 2A #t0 and t1 are both an n-by-2 matrix; t0 are initial values for two outcome variables # and t1 are final values for the same two variables. vec1 <- function(t0, t1, xlab="", ylab="", ...){ plot(t0, xlab, ylab, pch=16, col="grey50", xlim=c(range(c(t0[,1], t1[,1]))), ylim=c(range(c(t0[,2], t1[,2])))) arrows(t0[,1], t0[,2], t1[,1], t1[,2], length=0.1, col="grey50") points(mean(t0[,1]), mean(t0[,2]), pch=16, col="black", cex=1.5) arrows(mean(t0[,1]), mean(t0[,2]),mean(t1[,1]), mean(t1[,2]), col="black", length=0.15, lwd=3) } vec1(d.dat[,4:5], d.dat[,7:8]) # makes plot similar to Figure 2B vec2 <- function(t0, t1, xlab="", ylab="", ...){ rx <- extendrange(c(t0[,1], t1[,1])) ry <- extendrange(c(t0[,2], t1[,2])) plot(t0, xlab, ylab, pch=16, col="blue", xlim=rx, ylim=ry, xaxs="i", yaxs="i") arrows(t0[,1], t0[,2], t1[,1], t1[,2], length=0.1, col="grey50") points(x=c(mean(t0[,1]), mean(t1[,1])), y=c(ry[1],ry[1]), pch=c(21,23), cex=1.5, xpd=TRUE, bg=c("blue","red")) #proj mean x points(x=c(rx[1],rx[1]), y=c(mean(t0[,2]), mean(t1[,2])), pch=c(21,23), cex=1.5, xpd=TRUE, bg=c("blue","red"))#proj mean y bvbox(t0, lin.col="blue") bvbox(t1, lin.col="red") } vec2(d.dat[,4:5], d.dat[,7:8]) # makes plot similar to Figure 2C vec3 <- function(t0, t1, xlab="", ylab="", ...){ hull.t0 <- chull(t0) #makes convex hull of data hull.t1 <- chull(t1) rx <- extendrange(c(t0[,1], t1[,1])) ry <- extendrange(c(t0[,2], t1[,2])) par(tcl=0.3, las=1, mar=c(4.5,4.8,4,5.5)) plot(t0, xlab, ylab, xlim=rx, ylim=ry, pch=16, col="black", bty="n", xaxt="n", yaxt="n", type="n") polygon(t0[,1][hull.t0],t0[,2][hull.t0], border=rgb(0,0,1,0.4), col=rgb(0,0,1,0.4)) polygon(t1[,1][hull.t1],t1[,2][hull.t1], border=rgb(1,0,0,0.4), col=rgb(1,0,0,0.4)) points(t0, pch=16, col="black") axis(1,at=c(min(t0[,1]),max(t0[,1])), labels=c(round(min(t0[,1]),0),round(max(t0[,1]),0)), col="blue") mtext("Initial", side=1, line=2.5, at=mean(range(t0[,1])), col="blue") axis(4,at=c(min(t0[,2]),max(t0[,2])), labels=c(round(min(t0[,2]),0), round(max(t0[,2]),0)), col="blue") mtext("Initial", side=4, line=1, at=mean(range(t0[,2])), col="blue") axis(3,at=c(min(t1[,1]),max(t1[,1])), labels=c(round(min(t1[,1]),0),round(max(t1[,1]),0)), col="red", line=-1) mtext("Final", side=3, line=1.5, at=mean(range(t1[,1])), col="red") axis(2,at=c(min(t1[,2]),max(t1[,2])), labels=c(round(min(t1[,2]),0),round(max(t1[,2]),0)), col="red") mtext("Final", side=2, line=1, at=mean(range(t1[,2])), col="red") arrows(t0[,1], t0[,2], t1[,1], t1[,2], length=0.1, col="black") } vec3(d.dat[,4:5], d.dat[,7:8]) # makes plot similar to Figure 4 vec4 <- function(t0, t1, xlab="", ylab="", ...){ x.diff <- (t1-t0)[,1] y.diff <- (t1-t0)[,2] poly <- cbind(x.diff,y.diff) poly <- rbind(poly,c(0,0)) hull <- chull(poly) plot.new() plot.window(xlim=range(x.diff), ylim=range(y.diff)) polygon(poly[hull,], border="grey90", col="grey90") arrows(0,0, x.diff, y.diff, length=0.1, col="grey50") arrows(0,0,mean(x.diff), mean(y.diff), col="black", length=0.15, lwd=3) axis(1, pos=0) axis(2, pos=0) } vec4(d.dat[,4:5], d.dat[,7:8]) ##################### # 3D plots library(rgl) # makes plot similar to Figure 5A (need to rotate to obtain same orientation) vec.3d <- function(t0, tf){ # t0: an n by 3 matrix for initial values (X,Y,Z) # tf: an n by 3 matrix for final values (X,Y,Z) bg3d("white") aspect3d(1,1,3) rgl.spheres(t0[,1:3], radius=1, color="blue", alpha=0.5) #plots initial points for(i in 1:17){ # plots lines for each subject segments3d(x=c(t0[i,1],tf[i,1]), y=c(t0[i,2],tf[i,2]), z=c(t0[i,3], tf[i,3]), color="red") } grid3d(c("x","y","z+")) axis3d('x-', labels=TRUE, size=2, color="black", alpha=1) axis3d('y',labels=TRUE, pos=c(120,NA,50), size=2, color="black", alpha=1) axis3d('z+', labels=TRUE, size=2, color="black", alpha=1) mtext3d("X", edge="x", at=60, line=2, color="black", alpha=1) mtext3d("Y", edge="y", pos=c(110,NA,50), line=2.75, color="black", alpha=1) mtext3d("Z", edge="z+", line=2, color="black", alpha=1) } vec.3d(t0=d.dat[,4:6], tf=d.dat[,7:9]) ############################# # path plot # makes plot similar to Figure 6 all.dat <- read.csv(file="hd-data.csv") # read data again all.dat <- na.omit(all.dat) # removes NA values from data pred.uhdrs <- NULL # vector to store predicted UHDRS values pred.tap <- NULL # vector to store predicted taping values for(i in 1:17){ # loop through each patient, fit loess and get predicted values mod1 <- loess(avg.tap~month, data=all.dat, degree=1, span=0.9, subset=patient==i) mod2 <- loess(uhdrs~month, data=all.dat, degree=1, span=0.9, subset=patient==i) pred.tap <- append(pred.tap, predict(mod1)) pred.uhdrs <- append(pred.uhdrs, predict(mod2)) } path.plot <- function(x, y, case, subset, points=FALSE, ...){ # x: vector of numeric values for one variable # y: vector of numeric values for another variable # case: subject id (must be numeric) # subset: logical expression indicating elements used to plot the initial point at t=0 if(is.numeric(case) == TRUE) { # test to make sure case is numeric id <- length(unique(case)) } else stop("'case' must be numeric.\n") plot(y~x, subset=subset, xlim=range(x), ylim=range(y), pch=16) # plot initial points for(i in 1:id){ # draws path for each patient lines(x[case==i],y[case==i]) } if (points==TRUE) points(y~x, cex=0.2) # draws points along path for time stamp } path.plot(x=pred.tap, y=pred.uhdrs, case=as.numeric(all.dat$patient), subset=all.dat$month==0, points=TRUE) ########################## # Supplementary Figures library(lattice); library(Hmisc) # load required packages all.dat <- read.csv(file="hd-data.csv") # read data all.dat$patient <- as.factor(all.dat$patient) # make patients into a factor all.dat$patient <- reorder.factor(all.dat$patient, # order patients based on median tapping score all.dat$avg.tap, median, na.rm=TRUE) # Supplementary Figure 1 xyplot(avg.tap~(month/12)|patient,data=all.dat, ylab="Number of taps in 30 sec", xlab="Time (years)", panel=function(x,y){ panel.xyplot(x,y) panel.loess(x,y, degree=1, col="red", span=0.9,lwd=2) }) all.dat$patient <- reorder.factor(all.dat$patient, # order patients based on median UHDRS score all.dat$uhdrs, median, na.rm=TRUE) # Supplementary Figure 2 xyplot(uhdrs~(month/12)|patient,data=all.dat, ylab="UHDRS motor score", xlab="Time (years)", pty="s", panel=function(x,y){ panel.xyplot(x,y) panel.loess(x,y, degree=1, col="red", span=0.9, lwd=2) })