# ------------------------------------------------------
# Program: VectorFields1.S  
#  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

graphsheet(height=6.4,width=7.5)
plot(c(0,70), c(0,100), type="n" ,
   xlab="Age",
   ylab="Percent Correct")
lines(c(0:200)/3, tX1, type="l", lty=1)
lines(c(0:200)/3, tX2, type="l", lty=1)
lines(c(0:200)/3, tX3, type="l", lty=1)


# ------------------------------------------------------
#  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

graphsheet(height=6.4,width=7.5)
plot(c(0,201)/3, c(0,100), type="n" ,
   xlab="Age",
   ylab="Percent Correct")
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)
}


# ------------------------------------------------------
#  Plot the three columns as if they were a panel study with three
#  waves of measurement spaced 15 days apart

graphsheet(height=6.4,width=7.5)
plot(c(0,201)/3, c(0,100), type="n" ,
   xlab="Age",
   ylab="Percent Correct")
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)
}



# ------------------------------------------------------
#  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=slopeLoess$fitted.values)

# ------------------------------------------------------
#  Calculate a slope field of the fitted slopes

smoothX <- seq(0, 75, by=6)
smoothY <- seq(0, 100, 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 slope field of the fitted slopes

graphsheet(height=6.4,width=7.5)
plot(c(0,220)/3, c(0,101), type="n" ,
   xlab="Age",
   ylab="Percent Correct")
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))
}


# ------------------------------------------------------
#  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)), "c:/class/psych344509/ageData.dat", ncolumns=4)


# ------------------------------------------------------
#  Quit the program here