# ------------------------------------------------------ # Program: VectorFields1.R # Author: Steven M. Boker # Date: Tue Nov 30 10:58:06 EST 2004 # # Vector Field plotting # # ------------------------------------------------------ # ------------------------------------------------------ # Create and plot two autoregressive time series tX1 <- rep(NA, 201) tX1[1] <- 90 for (i in 2:201) { tX1[i] <- .98 * tX1[i-1] } tX1 <- 100 - tX1 + rnorm(201, mean=0, sd=10) #tX1[tX1 > 100] <- 100 tX2 <- rep(NA, 201) tX2[1] <- 90 for (i in 2:201) { tX2[i] <- .99 * tX2[i-1] } tX2 <- 100 - tX2 + rnorm(201, mean=0, sd=10) #tX2[tX2 > 100] <- 100 tX3 <- rep(NA, 201) tX3[1] <- 90 for (i in 2:201) { tX3[i] <- .995 * tX3[i-1] } tX3 <- 100 - tX3 + rnorm(201, mean=0, sd=10) #tX3[tX3 > 100] <- 100 pdf("VF1ARLongPlot1.pdf", height=6.4,width=7.5) plot(c(0,70), c(0,120), type="n" , xlab="Age", ylab="Total Score", cex.lab=1.75, cex.axis=1.5, cex.main=1.75, cex=1.75, mgp=c(2.5,.75,0)) lines(c(0:200)/3, tX1, type="l", lty=1, col=1) lines(c(0:200)/3, tX2, type="l", lty=1, col=2) lines(c(0:200)/3, tX3, type="l", lty=1, col=4) dev.off() # ------------------------------------------------------ # Plot a state space plots of the Autoregressive Data # with two waves of measurement spaced 15 days apart for (theTau in seq(1, 56, by=5)) { tfile <- paste("VF1AR1StateSpaceTau", theTau, ".pdf", sep="") pdf(tfile, height=6.4,width=6.4) plot(c(0,120), c(0,120), type="n" , xlab="Test Score (t)", ylab=paste("Test Score (t+", theTau, ")", sep=""), cex.lab=1.75, cex.axis=1.5, cex.main=1.75, cex=1.75, mgp=c(2.5,.75,0)) for (i in seq(1, (length(tX1)-(theTau+1)),by=1)) { lines(tX1[i],tX1[i+theTau], type="p", col=1) lines(tX2[i],tX2[i+theTau], type="p", col=2) lines(tX3[i],tX3[i+theTau], type="p", col=4) } dev.off() } # ------------------------------------------------------ # Create a timelagged dataframe (3-D Embedding) tau <- 15 tX <- c(tX1, rep(NA,2*tau+1), tX2, rep(NA,2*tau+1), tX3) tAge <- c(1:length(tX1), rep(NA,2*tau+1), 1:length(tX2), rep(NA,2*tau+1), 1:length(tX3))/3 tLen <- length(tX) ageData1 <- data.frame(x1=tX[1:(tLen-(2*tau))], x2=tX[(1+tau):(tLen-tau)], x3=tX[(1+(2*tau)):tLen], age=tAge[1:(tLen-(2*tau))]) # ------------------------------------------------------ # Plot the first two columns as if they were a panel study with two # waves of measurement spaced 15 days apart pdf("VF1ARLagged1.pdf", height=6.4,width=7.5) plot(c(0,201)/3, c(0,120), type="n" , xlab="Age", ylab="Total Score", cex.lab=1.75, cex.axis=1.5, cex.main=1.75, cex=1.75, mgp=c(2.5,.75,0)) for (i in seq(1,length(ageData1$x1),by=2)) { if (is.na(ageData1$age[i]) | ageData1$age[i] > (68-tau)) { next } lines(c(ageData1$age[i],ageData1$age[i]+tau), c(ageData1$x1[i],ageData1$x2[i]), type="l", lty=1) } dev.off() # ------------------------------------------------------ # Plot the three columns as if they were a panel study with three # waves of measurement spaced 15 days apart pdf("VF1ARLong2.pdf", height=6.4,width=7.5) plot(c(0,201)/3, c(0,120), type="n" , xlab="Age", ylab="Total Score", cex.lab=1.75, cex.axis=1.5, cex.main=1.75, cex=1.75, mgp=c(2.5,.75,0)) for (i in seq(1,length(ageData1$x1),by=2)) { if (is.na(ageData1$age[i]) | ageData1$age[i] > (68-(2*tau))) { next } lines(c(ageData1$age[i], ageData1$age[i]+tau, ageData1$age[i]+(2*tau)), c(ageData1$x1[i], ageData1$x2[i], ageData1$x3[i]), type="l", lty=1) } dev.off() # ------------------------------------------------------ # Calculate the slope and midpoint for each pair from the # first two columns of ageData1 ageMid <- ageData1$age + (tau/2) scoreMid <- (ageData1$x1 + ageData1$x2) / 2 slopeMid <- (ageData1$x2 - ageData1$x1) / tau tSelect <- !is.na(ageMid) & !is.na(scoreMid) & !is.na(slopeMid) ageSlopeFrame <- data.frame( age = ageMid[tSelect], score = scoreMid[tSelect], slope = slopeMid[tSelect]) # ------------------------------------------------------ # Calculate a loess smoothed slope surface slopeLoess <- loess(slope ~ age * score, data=ageSlopeFrame) ageSlopeFrame2 <- data.frame(age=ageSlopeFrame$age, score=ageSlopeFrame$score, slope=ageSlopeFrame$slope, fitted=predict(slopeLoess)) # ------------------------------------------------------ # Calculate a slope field of the fitted slopes smoothX <- seq(0, 75, by=6) smoothY <- seq(0, 120, by=8) smoothP <- matrix(NA, length(smoothX), length(smoothY)) smoothZ <- matrix(NA, length(smoothX), length(smoothY)) h <- 10 for (i in 1:length(smoothX)) { x <- smoothX[i] t1x <- abs(x - ageSlopeFrame2$age)/h for (j in 1:length(smoothY)) { y <- smoothY[j] t1y <- abs(y - ageSlopeFrame2$score)/h t2 <- rep(0, length(t1x)) t2[t1x < 1 & t1y < 1] <- 1/2 smoothP[i,j] <- sum(t2)/(length(ageSlopeFrame2$age)*h) smoothZ[i,j] <- mean(ageSlopeFrame2$fitted[t1x < 1 & t1y < 1]) } } smoothP <- 100 * (smoothP/sum(smoothP)) smoothSlopeFrame <- data.frame( age=rep(smoothX,length(smoothY)), score=rep(smoothY,each=length(smoothX)), density=c(smoothP), slope=c(smoothZ)) # ----------------------------------------------------- # Plot a 3-d Perspective plot of the aggregated slopes for (i in seq(0, 360, by=30)) { tfile <- paste("VF1SmoothSlopePersp", i, ".pdf", sep="") pdf(tfile, height=6.4,width=7.5) persp(unique(smoothSlopeFrame$age), unique(smoothSlopeFrame$score), matrix(smoothSlopeFrame$slope,nrow=13,ncol=16), theta = i, phi = 25, r = 3, d = 1, shade = .25, col = 4, xlab = "Age", ylab = "Score", zlab = "Slope", cex.lab=1.25, cex.axis=1.25, cex.main=1.25, cex=1.25, mgp=c(2.5,.75,0)) dev.off() } # ----------------------------------------------------- # Plot a 3-d Perspective plot of the aggregated densities pdf("VF1SmoothDensePersp1.pdf", height=6.4,width=7.5) persp(unique(smoothSlopeFrame$age), unique(smoothSlopeFrame$score), matrix(smoothSlopeFrame$density,nrow=13,ncol=16), theta = -45, phi = 25, r = 3, d = 1, shade = .25, col = 4, xlab = "Age", ylab = "Score", zlab = "Density", cex.lab=1.25, cex.axis=1.25, cex.main=1.25, cex=1.25, mgp=c(2.5,.75,0)) dev.off() # ------------------------------------------------------ # Plot a slope field of the fitted slopes pdf("VF1ARSlope1.pdf", height=6.4,width=7.5) plot(c(0,220)/3, c(0,121), type="n" , xlab="Age", ylab="Total Score", cex.lab=1.75, cex.axis=1.5, cex.main=1.75, cex=1.75, mgp=c(2.5,.75,0)) scale <- 2.5 for (i in 1:length(smoothSlopeFrame$age)) { tx <- smoothSlopeFrame$age[i] ty <- smoothSlopeFrame$score[i] ta <- smoothSlopeFrame$slope[i] td <- smoothSlopeFrame$density[i] txinc <- (scale * td) / sqrt( 1 + (ta^2)) tyinc <- ta * txinc lines(c(tx+txinc, tx-txinc), c(ty+tyinc, ty-tyinc)) } dev.off() # ------------------------------------------------------ # Write a file suitable for use with svf tSelect <- !is.na(ageData1$x1) & !is.na(ageData1$x2) & !is.na(ageData1$age) tMatrix <- cbind(ageData1$age[tSelect], ageData1$x1[tSelect], ageData1$age[tSelect] + tau, ageData1$x2[tSelect]) write(t(round(tMatrix,4)), "ageData.dat", ncolumns=4) # ------------------------------------------------------ # Quit the program here