---
title: "fPASS: An R package for Power and Sample Size analysis (PASS) for Projection-based Two-Sample test for functional data."
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{fPASS: An R package for Power and Sample Size analysis (PASS) for Projection-based Two-Sample test for functional data.}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  eval=TRUE
)
#options(rmarkdown.html_vignette.check_title = FALSE)
```

This document will present a elaborate understanding of the `fPASS` R package and its capabilities. First we need to load the current development version of the package using the following code. The package is designed to conduct for Power and Sample Size analysis (PASS) for Projection-based Two-Sample test for functional data. Please see [Wang (2021)](https://doi.org/10.1214/21-EJS1802) and [Koner and Luo (2023)](https://arxiv.org/abs/2302.05612) for the details of the testing procedure and its advantage over the current state-of-the-art testing procedures used for longitudinal data such as mixed model and GEE.

## Installation

```{r, eval = FALSE}
# Install development version from GitHub
devtools::install_github("SalilKoner/fPASS", build_vignettes = FALSE)
```

## Usage

```{r setup, eval=TRUE, message=FALSE}
library(fPASS)
library(nlme)
library(face)
```

```{r, include=F, eval=TRUE}
st <- proc.time()
```

The key function to use for finding the power and sample size (PASS) for the test is `fPASS::PASS_Proj_Test_ufDA()`. In the following, we will see how the function can be used to conduct PASS for the Projection-based test for longitudinal and functional data under different covariance structure in order to design a randomized clinical trials (RCT). We will present the examples of using the function with different case studies. 

### Case study: Calculating PASS analysis for designing RCT with Longitudinal data with a specific visit times. 

Suppose that we want to conduct a clinical trial to test the efficacy of a treatment effect against placebo, where the visit times for each subjects are at baseline, 3 months, 6 months, and every 6 months thereafter for 24 months. Assume that for each visit (except for the baseline), there is an window of 15 days around the stipulated visit time when the patients can appear to the clinic. This means that, a patient can come to the clinic for third visit around $5.5 - 6.5$ months from baseline. Assuming that the higher the response, the better is the treatment, it is hypothesized that the effect size (mean to variance ratio) between the treatment and placebo will be $1$ unit each year, and the treatment effect follows a linear trend. Further, it is assumed that the longitudinal response trajectory for each subject has a marginal standard deviation about 5 units at each visit time and assumes a dying correlation which can be characterized a conditional auto regressive process with correlation parameter $0.5$. Suppose that after collection of data, we want to use the [Two-sample projection-based test by Wang (2021)](https://doi.org/10.1214/21-EJS1802) to assess the treatment efficacy. Assuming that we can have equal number of subjects in each of the treatment and the placebo group, we want to find the minimum sample size required to test the treatment efficacy with a power of $80\%$ at a significance level of $5\%$, using the routine `fPASS::PASS_Proj_Test_ufDA()`. 

#### Argument specifiction for `fPASS::PASS_Proj_Test_ufDA()`

The argument of the function can be obtained by running `?fPASS::PASS_Proj_Test_ufDA()` in R. The key arguments of are `mean_diff_fnm` denoting the group difference function, `obs.design` specifying the observation design and `cov.par` specifying the covariance structure of the latent trajectory.  We will demonstrate how to supply each of the parameters in the function from the statement of the problem. 

* `obs.design` and `nobs_per_subj` : The number of observation points for each subject are pre-determined, and they are not randomly varying among each other. This constitutes that it falls under `longitudinal` design, and the visit times are $0, 3, 6, 12, 18$ and $24$ months with a window of $0.5$ months for each visit other than baseline. The `obs.design` argument should be set as

```{r, eval=TRUE}
obs.design <- list("design" = "longitudinal",
                  "visit.schedule" = c(3,6,12,18,24), # pre-determined visit times, other than baseline
                  "visit.window" = 0.5) # visit window
nobs_per_subj <- length(obs.design$visit.schedule) + 1 # number of observation per subject is 6.
```

* `cov.par` and `cov.type` : The conditional auto regressive process with correlation parameter $0.5$ can be characterized by the `nlme::corCAR1()` structure available in `nlme::corClasses`. Therefore, we can characterize the correlation as

```{r, eval=TRUE}
# The function will internally create a data set
# and a factor variable Subject denoting subject id.
# You can change time and Subject by any other name as well.
cor.str <- nlme::corCAR1(0.5, form = ~ time | Subject);
# The marginal sd at each time-point is assumed to be 5. 
cov.par  <- list("var" = 25, "cor" = cor.str) 
# We must set cov.type = "ST" for any covariance 
# structure belonging to nlme::corClasses. 
cov.type <- "ST"
```

* `mean_diff_fnm` :  Assuming that the group difference function has a linear trend, and that the group difference is 1 unit every year, which means the group difference is 4 unit after two years (24 months). As mentioned in the details in the argument section, internally, the function scales the specified visit times into $[0,1]$, so that baseline is 0 and 24 months is 1, therefore, the mean function can be written as

```{r, eval=TRUE}
mean_diff_fn  <- function(t){2*t}
mean_diff_fnm <- "mean_diff_fn" # the name of mean difference function needs to
                                # be specified with the argument. 
```

* `missing_type` and `missing_percent` : We assume no missing observation at each visit, so we will set the missing value related parameters as follows:

```{r, eval=TRUE}
missing_type <- "nomiss"
missing_percent <- 0 # Has be to a number between 0 and 0.8
```

* `alloc.ratio` : Assuming that we have an equal allocation ratio of the samples in each group, 

```{r, eval=TRUE}
alloc.ratio <- c(1,1)
```

* `target.power` and `sig.level` : We want to find the sample size required to achieve a power of $80\%$ of test at significance level $5\%$, 
```{r, eval=TRUE}
target.power <- 0.8
sig.level  <- 0.05
```

* `eval_SS` : The number of independent and identically distributed (i.i.d) replication of subjects based on which the eigencomponents of the covariance process will be estimated empirically. This should be set as a large number in order to correctly estimate the true number of eigenfunctions and to get enough precision on the estimated eigencomponents. Setting a small value of `eval_SS` might lead to incorrect detection of large number of eigenfunctions, and then the power of the projection-based test will incorrectly inflated because the dimension and the magnitude of projection-vector gets arbitrarily large, as a result of projecting the mean difference function on additionally detected eigenfunctions (possibly with very small eigenvalues). Users must remember that the **computation time of the function is directly proportional to the size of `eval_SS`, so they must expect that the run time will be higher if `eval_SS` is larger, especially when `fpca.method` is set as 'face'**, which is described next item. 

```{r, eval=TRUE}
eval_SS <- 5000
```

* `fpca_method` : The module for functional principal component analysis that will be used to estimate the eigencomponents, can be any one of
'fpca.sc', in which case `fPASS::fpca.sc()` to be used or 'face' for `face::face.sparse` function to be used as the primary eigenfunction estimation routine. Note that, `fPASS::fpca.sc()` is just a copy of `refund::fpca.sc()` routine, except here the `shrinkage` scores (Yao et. al. 2005, JASA) are correctly estimated, and the issue of `NA` value for the scores when the measurement error variance are estimated to be zero, are corrected. 

```{r, eval=TRUE}
fpca_method <- "fpca.sc"
```

* `sigma2.e` : Set a small but non-zero measurement error to ensure nugget effect and stability in inversion of covariance

```{r, eval=TRUE}
sigma2.e <- 0.001
```

*  `fpca_optns` : The argument `fpca_optns` must be a named list with elements that could specified as arguments of either of `fPASS::fpca.sc()` function or  `face::face.sparse()` depending on the choice of the argument `fpca_optns`. For example, the default value of `pve` argument is $0.99$, but the user can overwrite by setting it as follows. See the details of argument in the help page for what are compatible names for the elements in the list. A higher value of percentage of variation explained might lead to detection of extra eigenfunctions with close to zero eigenvalues, which might lead to inflated power of the test for fixed sample size. 

```{r, eval=TRUE}
fpca_optns  <- list("pve" = 0.95) 
```

*  `mean_diff_add_args` : The user can specify any additional parameter in form of a named list that might need to be specified in the mean different function specified in the argument `mean_diff_fnm`. In this case, there is nothing, so it will be specified as an empty list as follows:

```{r, eval=TRUE}
mean_diff_add_args  <- list() 
```

* `nWgrid ` : The length of grid points used to estimate the eigenfunction and to approximately compute the projection $\int \mu_1(t) - \mu_2(t)\phi_k(t) dt$. We keep the default value of `nWgrid`, $201$. 

```{r, eval=TRUE}
nWgrid <- 201
```

* `nsim` : The number of samples to be generated from the alternate distribution of the Hotelling $T^2$ statistic to accurate compute the power. The default value is set at 10000, but the user can change it in the function. In our example, we have used `nsim`$=1000$ for faster computation of the vignette. 

#####  Distribution of minimum sample size using `fPASS::fpca.sc()` function:

Now that we have explained how to specify each and every arguments of `fPASS::PASS_Proj_Test_ufDA()` function, we will compute the minimum sample size required. However, as the testing procedure internally computes the eigenfunctions based on large number of samples and then uses the estimated eigenfunctions to compute the PASS, the computed power and the minimum sample size will be different due to the fact that the estimated eigenfunctions are different at each run because of sampling variation. Therefore, it is important to run the function `fPASS::PASS_Proj_Test_ufDA()` for few iterations to see how the sampling variation affects the computed power or sample size. In this example, we run the same function in parallel for about 10 times, and present the percentiles of the estimated sample size. It is recommended to run the the function at least about 50 to 100 times and use the median of these computed sample sizes obtained after a reasonable number of iteration, as a point estimate of the minimum sample size and the interquartile range (IQR) as a valid measure of range of the sample size that could be considered. From a practical perspective, depending on the budget of the trial, one can choose the minimum number of subjects that will be enrolled for the trial based on the value of IQR.

```{r, message=FALSE, results='hide'}
library(foreach)
required_sample_size <- foreach(i=1:10, .combine='c', .packages = "fPASS") %do% {
        mean_diff_fn <- mean_diff_fn
        PASS_Proj_Test_ufDA(sample_size = NULL, target.power = target.power,
                            sig.level = sig.level, nobs_per_subj = nobs_per_subj,
                            obs.design = obs.design, mean_diff_fnm = "mean_diff_fn",
                            cov.type = cov.type, cov.par = cov.par,
                            sigma2.e = sigma2.e, missing_type = missing_type,
                            missing_percent = missing_percent, eval_SS = eval_SS,
                            alloc.ratio = alloc.ratio, fpca_method = fpca_method,
                            mean_diff_add_args = mean_diff_add_args,
                            fpca_optns = fpca_optns, npc_to_use = NULL,
                            nsim = 1e3)$required_SS
}
```

```{r}
quantile(required_sample_size, probs = c(0.25,0.5,0.75), names = TRUE)
```

We have intentionally left out the discussion of `npc_to_use` argument for the `fPASS::PASS_Proj_Test_ufDA()` function. The number of eigenfunction to use to compute the PASS is not before hand. Therefore, for a specified sampling design and covariance structure of the response trajectory, it is not possible to to provide an appropriate value unless the user has an substantial information about the eigencomponents before. This is where the function `fPASS::Extract_Eigencomp_fDA()` can assist. If we look at the help page of that function, you can notice that except for one argument, all of arguments of that function matches that of `fPASS::PASS_Proj_Test_ufDA()`. Therefore, for a specified sampling design and covariance parameter we can take a look the estimated eigenfunctions and eigenvalues to have a more informed decision how many eigencomponents to use to in the main function to conduct the PASS analysis. Here is an example of concering to in our case study.

#####  Estimate eigenfunctions using `fPASS::fpca.sc()` function

```{r, message = F, results='hide'}
fpca_method <- "fpca.sc"

library(foreach)
eigenlist <- foreach(i=1:3, .packages = "fPASS") %do% {
  mean_diff_fn <- mean_diff_fn
eigencomp <- Extract_Eigencomp_fDA(nobs_per_subj = nobs_per_subj,
                      obs.design = obs.design, mean_diff_fnm = mean_diff_fnm,
                      cov.type = cov.type, cov.par = cov.par,
                      sigma2.e = sigma2.e, missing_type = missing_type,
                      missing_percent = missing_percent, eval_SS = eval_SS,
                      alloc.ratio = alloc.ratio, fpca_method = fpca_method,
                      mean_diff_add_args = mean_diff_add_args,
                      fpca_optns = fpca_optns,
                      data.driven.scores = FALSE) 
 eigencomp[c("working.grid", "est_eigenfun", "est_eigenval")]
}
```

```{r}
print(lapply(eigenlist, function(x) x$est_eigenval))
```

For a specified PVE of $95\%$, `fPASS::fpca.sc()` function detects two leading eigenfunction with non-zero eigenvalues for the covariance structure in most of the cases. Next, we see what happens by using `face::face.sparse()` function for the same PVE The `face::face.sparse()` function is much slower than the `fPASS::fpca.sc()`, so we have kept `eval_SS` argument at $500$ to not significantly increase the run time of the vignette. The user are advised to use a large value of `eval_SS` in order to get high enough precision for the estimated eigenfunction. 

#####  Estimate eigenfunctions using `face::face.sparse()` function

```{r, message = F, results='hide'}

fpca_method <- "face"
# Setting eval_SS to 500, specifically for 
# the face::face.sparse() function
# to ensure that the vignette builds within
# real time. The user must set a much larger value
# of eval_SS e.g. ranging between 2000-5000, to ensure
# enough precision accuracy for the eigenfunctions. 
eval_SS  <- 500 # increase it to 2000 for practical case

library(foreach)
# library(doParallel)
# cl <- makeCluster(detectCores()-1)
# registerDoParallel(cl)
eigenlist <- foreach(i=1:3, .packages = "fPASS") %do% {
  mean_diff_fn <- mean_diff_fn
eigencomp <- Extract_Eigencomp_fDA(nobs_per_subj = nobs_per_subj,
                      obs.design = obs.design, mean_diff_fnm = mean_diff_fnm,
                      cov.type = cov.type, cov.par = cov.par,
                      sigma2.e = sigma2.e, missing_type = missing_type,
                      missing_percent = missing_percent, eval_SS = eval_SS,
                      alloc.ratio = alloc.ratio, fpca_method = fpca_method,
                      mean_diff_add_args = mean_diff_add_args,
                      fpca_optns = fpca_optns,
                      data.driven.scores = FALSE) 
 eigencomp[c("working.grid", "est_eigenfun", "est_eigenval")]
}
# stopCluster(cl)
# matplot(eigencomp$working.grid, eigencomp$est_eigenfun, type="l", 
#         xlab = "timepoints (scaled)", ylab = "eigenfunctions")
```

```{r, message=F}
print(lapply(eigenlist, function(x) x$est_eigenval))
```

We notice that `face::face.sparse()` function is detecting three eigenfunction in most of the cases, for the same setup. However, the eigenvalue corresponding to the last eigenfunction seems to be small, and that is possibly the reason it was not detected by `fPASS::fpca.sc()` function. This provides us a more informed decision about the number of eigenfunctions to use. So, technically we should discard the last eigenfunction to compute the projection of mean difference onto the eigenfunctions, and use `npc_to_use = 2` in the `fPASS::PASS_Proj_Test_ufDA` function, so that the power function is not overinflated by detection of extra eigenfunctions. However, for better understanding of the reader, we will still use `npc_to_use = 3` for `fpca_method == 'face'` to find out how the distribution of minimum sample size required is allowing an extra eigenfunction (with possibly zero eigenvalues) for `fpca_method == 'face'`. The number of replication is kept at 10, to ensure that the vignette runs within a reasonable amount of time, please increase it for more comprehensive understanding of the effect of estimating eigenfunctions on range of the minimum sample size required. 

#####  Distribution of minimum sample size using `face::face.sparse()` function with pre-specified number of eigenfunction.

```{r}
npc_to_use <- 3
```

```{r, message=F, results='hide'}
fpca_method <- "face"

library(foreach)
required_sample_size <- foreach(i=1:5, .combine='c', .packages = "fPASS") %do% {
        mean_diff_fn <- mean_diff_fn
        PASS_Proj_Test_ufDA(sample_size = NULL, target.power = target.power,
                            sig.level = sig.level, nobs_per_subj = nobs_per_subj,
                            obs.design = obs.design, mean_diff_fnm = mean_diff_fnm,
                            cov.type = cov.type, cov.par = cov.par,
                            sigma2.e = sigma2.e, missing_type = missing_type,
                            missing_percent = missing_percent, eval_SS = 500,
                            alloc.ratio = alloc.ratio, fpca_method = fpca_method,
                            mean_diff_add_args = mean_diff_add_args,
                            fpca_optns = fpca_optns, npc_to_use = npc_to_use,
                            nsim = 1e3)$required_SS
}
```

```{r}
quantile(required_sample_size, probs = c(0.25,0.5,0.75), names = TRUE)
```

As expected, using a higher value of `npc_to_use` leads to over inflation of power, and consequently computed sample size is significantly less than which we obtained for `fpca_method = 'fpca.sc'`. But, this comes at cost of over inflation of type I error of the test. Therefore, for a specified correlation structure of the response, we must take extra care to NOT use unnecessary eigenfunctions, to avoid the risk of inflation of inflated type I error. 

#### Missing data

Under the same setup, if We want to compute the minimum sample size required when the percentage of missing observation at each time point is around $25\%$, then we can see the distribution of sample size required after running the function about 100 times as follows. 

```{r}
missing_type <- "constant"
missing_percent <- 0.25
```

```{r, message=FALSE, results='hide'}
fpca_method <- "fpca.sc"
eval_SS  <- 5000

library(foreach)
required_sample_size <- foreach(i=1:10, .combine='c', .packages = "fPASS") %do% {
        mean_diff_fn <- mean_diff_fn
        PASS_Proj_Test_ufDA(sample_size = NULL, target.power = target.power,
                            sig.level = sig.level, nobs_per_subj = nobs_per_subj,
                            obs.design = obs.design, mean_diff_fnm = "mean_diff_fn",
                            cov.type = cov.type, cov.par = cov.par,
                            sigma2.e = sigma2.e, missing_type = missing_type,
                            missing_percent = missing_percent, eval_SS = eval_SS,
                            alloc.ratio = alloc.ratio, fpca_method = fpca_method,
                            mean_diff_add_args = mean_diff_add_args,
                            fpca_optns = fpca_optns, npc_to_use = 2, nsim = 1e3)$required_SS
}
```

```{r, message=FALSE}
quantile(required_sample_size, probs = c(0.25,0.5,0.75), names = TRUE)
```

The IQR of the distribution of the minimum sample sizes when the percentage of missing is as significant as $25\%$ is still about the same as that when there is no missing, i.e. `missing_type == 'nomiss'` and `missing_percent = 0`. This reflects one of the biggest advantage of the projection-based test compared to traditional mixed model approaches, that as long as the missing percentage do not affect the estimation of the eigencomponents, the power of the test does not degrade with the increasing percentage of missing responses. 

Do try it out the same example in NCSS PASS software and check out the minimum sample size required for the test to achieve a power of $80\%$, and how does the minimum sample size increases when you increase the missing percentage to $25\%$, to get a fair comparison of the effectiveness of our procedure compared to traditional GEE type tests. To do this, go to NCSS PASS (2023) software, *GEE -> GEE Tests for the Slope of Two groups in a Repeated Measures Design (Continuous Outcome) -> Choose AR(1) proportional correlation structure and Number of measurement times as 6.* 

We report the total computation time it needed create the entire vignette is approximately `r round((proc.time()[3] - st[3])/60)` minutes.

