# Illustration of the R program for implementing the methodology proposed in # "Modeling and Analysis of Functional Method Comparison Data" # by G. S. R. de Silva, L. N. Rathnayake and P. K. Choudhary. # Download the program: # Right click on the link "R Code" and save the # associated file as FMCDA.R in the R working directory. # Requires: # The program requires "mgcv", "mvtnorm" and "MASS" packages. # They can be obtained from the CRAN website of R. # Illustration using core body temperature data: # Step 1. Launch R. # Step 2. Change the working directory to the directory where # FMCDA.R file is stored. # Step 3. Run source("FMCDA.R"). This will also load # the required packages (provided they are already installed) # Step 4. Run the following commands in R source("FMCDA.R") # Illustration 1 # Note: The files Tes_a.txt and Tre_a.txt containing the actual data # are not being provided as we don't have the requisite permission. These # files have data stored as 90 (number of time points) x 12 (number of subjects) # tables. The program can be used for any data set that has the structure of # these data. tes_a <- t(read.table("Tes_a.txt")) # data file not provided tre_a <- t(read.table("Tre_a.txt")) # data file not provided str(tes_a) # num [1:12, 1:90] 37.1 37.1 37.1 36.7 37 ... # - attr(*, "dimnames")=List of 2 # ..$ : chr [1:12] "V1" "V2" "V3" "V4" ... # ..$ : NULL str(tre_a) # num [1:12, 1:90] 37.9 37.2 37.4 37 37.3 ... # - attr(*, "dimnames")=List of 2 # ..$ : chr [1:12] "V1" "V2" "V3" "V4" ... # ..$ : NULL y1 <- tes_a y2 <- tre_a Y <- cbind(y1,y2) nboot <- 500 # number of bootstrap repetitions alpha <- 0.05 # 1 - level of confidence I <- dim(y1)[1] # number of subjects D <- dim(y2)[2] # number of grid points argvals <- 1:D # estimate parameters mpace.est <- MPACE(Y, argvals, D, I) # Note: Here (nbasis, pve.biv, prob, makePD) are at their default values (10, 0.99, 0.90, FALSE) upace.est <- UPACE(y1, y2, argvals, D, I) # Note: Here (nbasis, pve.univ, pve.biv, prob, makePD) are at their default values (10, 0.99, 0.99, 0.90, FALSE) set.seed(7) # get bootstrap confidence intervals boot.result <- bootstrap(nboot, y1, y2, argvals, D, I) # Note: Here (nbasis, pve.univ, pve.biv, alpha, prob, makePD) are at their default values (10, 0.99, 0.99, 0.05, 0.90, FALSE) # The results stored in boot.result have been reported in Section 5 of the manuscript. # Illustration 2 # Note: The files y1.txt and y2.txt containing the simulated data # are being provided. These files have data stored as # 50 (number of time points) x 60 (number of subjects) tables. y1 <- as.matrix(read.table("y1.txt")) # data file provided y2 <- as.matrix(read.table("y2.txt")) # data file provided Y <- cbind(y1,y2) nboot <- 500 # number of bootstrap repetitions alpha <- 0.05 # 1 - level of confidence I <- dim(y1)[1] # number of subjects D <- dim(y2)[2] # number of grid points argvals <- 1:D # estimate parameters mpace.est <- MPACE(Y, argvals, D, I) # Note: Here (nbasis, pve.biv, prob, makePD) are at their default values (10, 0.99, 0.90, FALSE) upace.est <- UPACE(y1, y2, argvals, D, I) # Note: Here (nbasis, pve.univ, pve.biv, prob, makePD) are at their default values (10, 0.99, 0.99, 0.90, FALSE) set.seed(1234) # get bootstrap confidence intervals boot.result <- bootstrap(nboot, y1, y2, argvals, D, I) # Note: Here (nbasis, pve.univ, pve.biv, alpha, prob, makePD) are at their default values (10, 0.99, 0.99, 0.05, 0.90, FALSE) # Results: names(mpace.est) # "mean" "mean.diff" "evalues" "efunctions" "npc" # "ltau.sq" "ltau.sq.ratio" "fit" "scores" "zcor.coef" # "ZCCC" "log.tdi" mpace.est$evalues # 31.4261577 15.1709883 8.2279972 5.9498773 0.5168110 0.4697129 mpace.est$ltau.sq # -1.38732 -1.38732 head(mpace.est$log.tdi) # 0.1505816 0.1505816 0.1505816 0.1505816 0.1505816 0.1505816 head(mpace.est$ZCCC) # 0.8628706 0.7970728 0.7687018 0.7728092 0.7898173 0.7984797 names(upace.est) # "mean" "mean.diff" "evalues" "efunctions" "npc" # "ltau.sq" "ltau.sq.ratio" "fit" "scores" "zcor.coef" # "ZCCC" "log.tdi" upace.est$evalues # 30.634223 14.433924 7.514765 5.221867 upace.est$ltau.sq # -1.313705 -1.313705 head(upace.est$log.tdi) # 0.1873835 0.1873835 0.1873835 0.1873835 0.1873835 0.1873835 head(upace.est$ZCCC) # 0.7947647 0.7415159 0.7181287 0.7197505 0.7310968 0.7355361 names(boot.result) # "mu1.univ.pci" "mu2.univ.pci" "mu1.univ.ci" "mu2.univ.ci" # "mean.diff.univ.pci" "mean.diff.univ.ci" "ltau.sq1.univ.ci" "ltau.sq2.univ.ci" # "ltau.sq.ratio.univ.ci" "zcor.coef.univ.pci" "zccc.univ.pci" "ltdi.univ.pci" # "zcor.coef.univ.ci" "zccc.univ.ci" "ltdi.univ.ci" "mu1.biv.pci" # "mu2.biv.pci" "mu1.biv.ci" "mu2.biv.ci" "mean.diff.biv.pci" # "mean.diff.biv.ci" "ltau.sq1.biv.ci" "ltau.sq2.biv.ci" "ltau.sq.ratio.biv.ci" # "zcor.coef.biv.pci" "zccc.biv.pci" "ltdi.biv.pci" "zcor.coef.biv.ci" # "zccc.biv.ci" "ltdi.biv.ci" head(boot.result$zccc.biv.pci) # lcl 0.3040576 0.2803214 0.258041 0.2376141 0.2194386 0.203925 head(boot.result$zccc.univ.pci) # lcl 0.231412 0.2159941 0.20101 0.1869621 0.174391 0.1638637