In statKim/robfpca: Functional principal component analysis for partially observed elliptical process

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

library(robfpca)

Comparison between Proposed FPCA and Sparse FPCA

Generate partially observed functional data from heavy-tailed distribution

First, we generate 100 partially observed curves from heavy-tailed $t_3$ distribution.

# Generate partially observed curves from heavy-tailed distribution
set.seed(46)
X.list <- sim_delaigle(n = 100,
                       type = "partial",
                       dist = "tdist")
X <- list2matrix(X.list)   # transform list to matrix
gr <- seq(0, 1, length.out = ncol(X))   # observed timepoints
matplot(gr, t(X), 
        type = "l",
        xlab = "t", ylab = "X(t)")

Robust mean and covariance estimation

We estimate the mean and covariance functions using cov_ogk() in this package which implements proposed method via M-estimator and equation (3.4) in our paper.

# Robust mean and covariance functions
cov.obj <- cov_ogk(X,
                   type = "huber",
                   MM = TRUE,
                   smooth = TRUE,
                   bw = 0.3)
mu.ogk <- cov.obj$mean
cov.ogk <- cov.obj$cov

Competing methods

We consider the non-robust FPCA method for the incomplete functional data proposed by Yao et al. (2005).

### Yao et al. (2005)
library(fdapace)
bw <- 0.1   # fixed bandwidth
optns <- list(methodXi = "CE", dataType = "Sparse", kernel = "epan", 
              verbose = FALSE, error = FALSE,
              userBwMu = bw, userBwCov = bw)
mu.yao.obj <- GetMeanCurve(Ly = X.list$Ly, Lt = X.list$Lt, optns = optns)
cov.yao.obj <- GetCovSurface(Ly = X.list$Ly, Lt = X.list$Lt, optns = optns)
mu.yao <- mu.yao.obj$mu
cov.yao <- cov.yao.obj$cov

Following figures show the true and estimated correlation surfaces.

# Correlation surfaces
library(GA)
par(mfrow = c(1, 3))
cov.true <- get_delaigle_cov(gr, model = 2)   # True covariance function
persp3D(gr, gr, 
        cov2cor(cov.true),
        xlab = "s", ylab = "t",
        main = "True",
        theta = -70, phi = 30, expand = 1)
persp3D(gr, gr, 
        cov2cor(cov.yao),
        xlab = "s", ylab = "t", 
        main = "Yao et al. (2005)",
        theta = -70, phi = 30, expand = 1)
persp3D(gr, gr, 
        cov2cor(cov.ogk),
        xlab = "s", ylab = "t",
        main = "Proposed",
        theta = -70, phi = 30, expand = 1)

Functional principal component analysis

We perform FPCA based on the conditional expectation. (Yao et al. (2005))

# Yao et al.(2005)
pca.yao.obj <- funPCA(Lt = X.list$Lt, 
                      Ly = X.list$Ly,
                      mu = mu.yao, 
                      cov = cov.yao, 
                      work.grid = gr,
                      K = 4)

# Proposed
pca.ogk.obj <- funPCA(Lt = X.list$Lt, 
                      Ly = X.list$Ly,
                      mu = mu.ogk, 
                      cov = cov.ogk, 
                      work.grid = gr,
                      K = 4)

# Plot 1~3 eigenfunctions
eig.true <- get_delaigle_eigen(gr, model = 2)   # True eigenfunctions
eig.yao <- check_eigen_sign(pca.yao.obj$eig.fun, 
                            eig.true)   # match eigen directions
eig.ogk <- check_eigen_sign(pca.ogk.obj$eig.fun, 
                            eig.true)   # match eigen directions
par(mfrow = c(1, 3))
for (i in 1:3) {
    plot(gr, eig.true[, i],
         type = "l",
         lwd = 2,
         ylim = c(-3, 3),
         xlab = "t", 
         ylab = paste("PC", i))
    lines(gr, eig.yao[, i], col = 2, lwd = 2)
    lines(gr, eig.ogk[, i], col = 3, lwd = 2)

    if (i == 1) {
        legend("bottomright", 
               c("True","Yao et al.(2005)","Proposed"),
               lty = c(1, 1, 1),
               col = c(1, 2, 3))
    }
}

Completion and reconstruction for missing parts

# Curve reconstruction
pred_yao_mat <- predict(pca.yao.obj, K = 4)
pred_ogk_mat <- predict(pca.ogk.obj, K = 4)


### Completion and Reconstruction example
i <- 1   # 1st curve
NA_ind <- which(is.na(X[i, ]))   # index of missing periods
X_comp <- cbind(X[i, ],
                NA,
                NA)
X_comp[NA_ind, 2] <- pred_yao_mat[i, NA_ind]  # completion of Sparse FPCA 
X_comp[NA_ind, 3] <- pred_ogk_mat[i, NA_ind]  # completion of proposed method

par(mfrow = c(1, 2))
# Completion
matplot(gr, 
        X_comp,
        type = "l",
        lwd = 2,
        xlab = "t", ylab = "X(t)")
lines(gr[NA_ind], 
      X.list$x.full[i, NA_ind],
      col = "darkgray", lty = 2, lwd = 2)
grid()
legend("bottomleft", 
       c("True","Yao et al.(2005)","Proposed"),
       lty = 1:3,
       col = 1:3)

# Reconstruction
matplot(gr, 
        cbind(X[i, ], 
              pred_yao_mat[i, ],
              pred_ogk_mat[i, ]),
        type = "l",
        lwd = 2,
        xlab = "t", ylab = "X(t)")
lines(gr[NA_ind], 
      X.list$x.full[i, NA_ind],
      col = "darkgray", lty = 2, lwd = 2)
grid()

Reference

• Park, Y., Kim, H. & Lim, Y. (2023). Functional principal component analysis for partially observed elliptical process, Computational Statistics & Data Analysis, 184, 107745.
• Yao, F., Müller, H. G., & Wang, J. L. (2005). Functional data analysis for sparse longitudinal data. Journal of the American Statistical Association, 100(470), 577-590.

statKim/robfpca documentation built on April 15, 2023, 10:12 p.m.