# -------------------------------------------------------------------- # 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") }