Estimation of functional data observations using fPCA, with applications to DTI tractography data

Sahil Seth, Dept. of Biostatistics, Johns Hopkins University

 

Help

Contact:

Sahil Seth

Dept of Biostatistics

Johns Hopkins Univ

sseth4 [at] jhu.edu

R script to run some examples:


# script to create a dataset

source("http://biostat.jhsph.edu/~saseth/R/mkDat.R")

# some plotting and other functions

source("http://biostat.jhsph.edu/~saseth/R/funcs.R")


## making a sample dataset

dat <- mkDat(I=25,noVis=5,D=30)

YFinal <- dat$Y

resFinal <- LFPCA(YFinal,dat$subject,dat$time)


################# Segment ##################

## remove some segments

fracMiss <- 0.5 #fraction of functions with missing segments

minLen <- 4 # min length of the segment

maxLen <- 10 # max length of the segment, should be less than D

rmSeg <- removeDat(YFinal,type="seg",fracMiss,minLen,maxLen,sim)

resMiss <- LFPCA(rmSeg$missDat, subject1, time1)


rmFuncs <- rmSeg$missFuncs

tsubs <- subject1[rmFuncs]

tvis <- visit1[rmFuncs]

ttim <- time1[rmFuncs]


resEst <- predict(resFinal,tsubs,ttim,tvis)

resPred <- predict(resMiss,tsubs,ttim,tvis)


## plot them in the working dir

plotFuncs(model=resFinal,est=resEst,

          pred=resPred,file=”missSegs.pdf"

          ,type="seg",miss=rm,plotVis=TRUE)


## get rmpse for them

R2YEst <- getR2(YFinal[rmFuncs,],resEst$YHat)$R2

R2YPred <- getR2(YFinal[rmFuncs,],resPred$YHat)$R2

R2EstPred <- getR2(resEst$YHat,resPred$YHat)$R2


############ Tract ######################

## remove some tracts

rm <- removeDat(dat=YFinal,type="tract",fracMiss=fracMiss,

                sim=sim,

                subject=subject1,minVis=2)

## dat=YFinal;type="tract";fracMiss=fracMiss;sim=sim;

##                 subject=subject1;minVis=2

rmFuncs <- rm$missFuncs

tsubs <- subject1[rmFuncs]

tvis <- visit1[rmFuncs]

ttim <- time1[rmFuncs]


resMiss <- LFPCA( YFinal[-rmFuncs,], subject1[-rmFuncs],

                 time1[-rmFuncs])


## plot them in the working dir

resEst <- predict(resFinal,tsubs,ttim,tvis)

resPred <- predict(resMiss,tsubs,ttim)


plotFuncs(model=resFinal,est=resEst,

          pred=resPred,file=”missTracts.pdf",

          ,type="tract",miss=rm,plotVis=FALSE)


## get rmpse for them

## both should be in the same order

R2YEst <- getR2(YFinal[rmFuncs,],resEst$YHat)$R2

## this gave spline error

R2YPred <- getR2(YFinal[rmFuncs,],resPred$YHat)$R2

R2EstPred <- getR2(resEst$YHat,resPred$YHat)$R2

Step 1: Install and load these dependencies:

require(mgcv)                           #gamm function

require(Matrix)                         #used in simulations

library(orthopolynom)                   #used in simulations


Step 2: Install the predictions package as follows:

Download it and cd to the directory

[http://biostat.jhsph.edu/~saseth/R/predLFPCA_1.0.tar.gz]:

R CMD INSTALL predLFPCA_1.0.tar.gz


Step 3:

Load the package:

library(predLFPCA)