# --------------------------------------------------------------------
# Example differential equations analysis of Widowhood Project data.
# 
# This example can be run in the free version of the R software
# which can be downloaded from http://cran.r-project.org

# --------------------------------------------------------------------
# Load the mixed effects library and set the output width

library(nlme)
options(width=80)

# --------------------------------------------------------------------
# Read the data from disk as a list.

theDataList <- scan("test2.dat", list(id=0, day=0, mhi=0), skip=0)

# --------------------------------------------------------------------
# Create a matrix of these data

rawMatrix <- cbind(theDataList$id, theDataList$day, theDataList$mhi)

# --------------------------------------------------------------------
# Create a dataframe of these data and print summary statistics

rawFrame <- data.frame(theDataList)
summary(rawFrame)

# --------------------------------------------------------------------
# Calculate a vector of ID numbers and the total number of subjects

uniqueSubjects <- unique(rawFrame$id)
maxSubjects <- length(uniqueSubjects)

# --------------------------------------------------------------------
# Set the number of occasions to skip when creating the embedding

theTau <- 11

# --------------------------------------------------------------------
# Set the time scale to 1 day (interval of time between occasions)

theInterval <- 1

# --------------------------------------------------------------------
# Create matrices to store results

idMHI <- matrix(NA, maxSubjects, 100)
xMHI <- matrix(NA, maxSubjects, 100)
dxMHI <- matrix(NA, maxSubjects, 100)
d2xMHI <- matrix(NA, maxSubjects, 100)


# --------------------------------------------------------------------
# Loop over all subjects and calculate local linear approximation 
#   derivatives for the Mental Health Inventory aggregate variable.

for (i in c(1:maxSubjects)) {
    theSubject <- uniqueSubjects[i]
    cat ("\n Subject = ", theSubject, sep="")
    tSelect <- rawFrame$id == theSubject
    theData <- rawFrame$mhi[tSelect]
    theTime <- rawFrame$day[tSelect]
    maxSamples <- length(theData)
    theFit <- lsfit(theTime,theData)
    t0 <- theFit$residuals		
    t1 <- t0[1:(maxSamples-(theTau*2))]
    t2 <- t0[(1+theTau):(maxSamples-theTau)]
    t3 <- t0[(1+(2*theTau)):maxSamples]
    tSelect <- !is.na(t1) & !is.na(t2) & !is.na(t3)
    x <- t2[tSelect]
    dx <- (t3[tSelect]-t1[tSelect]) / (2*theTau*theInterval)
    d2x <- (t3[tSelect] - (2 * t2[tSelect]) + t1[tSelect]) / 
            (theTau*theTau*theInterval*theInterval)
    idMHI[i,1:length(x)] <- theSubject
    xMHI[i,1:length(x)] <- x
    dxMHI[i,1:length(x)] <- dx
    d2xMHI[i,1:length(x)] <- d2x
}

# --------------------------------------------------------------------
# Turn output matrices into vectors

tid <- c(idMHI)
tx <- c(xMHI)
tdx <- c(dxMHI)
td2x <- c(d2xMHI)

# --------------------------------------------------------------------
# Remove all data triplets with missing values for MHI

tSelect <- !is.na(tx)

# --------------------------------------------------------------------
# 

rawFrameTau11 <- data.frame(id=tid[tSelect],
                  mhi=tx[tSelect], 
                  dmhi=tdx[tSelect], 
                  d2mhi=td2x[tSelect] )
summary(rawFrameTau11)

# --------------------------------------------------------------------
# 

tregTau3a <- lme(d2mhi ~ dmhi + mhi - 1, 
                 random= ~ mhi + dmhi - 1 | id, 
                 data=rawFrameTau11)
summary(tregTau3a)

# --------------------------------------------------------------------
# 

tcoefTau3a <- coef(tregTau3a)

# --------------------------------------------------------------------
# 

1 - (var(tregTau3a$residuals[,2]) / var(rawFrameTau11$d2mhi))

# --------------------------------------------------------------------
# 

for (tSubject in uniqueSubjects) {
        tSelect <- rawFrameTau11$id == tSubject
        cat(tSubject, "  r^2=", 
        1 - (var(tregTau3a$residuals[tSelect,2]) / 
             var(rawFrameTau11$d2mhi[tSelect])), 
        "N=", length(tregTau3a$residuals[tSelect,2]), "\n")
}