gfpca_TwoStep: gfpca_TwoStep

In jeff-goldsmith/gfpca: Generalized Functional Principal Components Analysis

Description

Implements a two-step approach to generalized functional principal components analysis for sparsely observed binary curves

Usage

gfpca_TwoStep(data, npc = NULL, pve = 0.9, output_index = NULL,
  type = c("approx", "naive", "fixed"), basis = NULL, nbasis = 10)

Arguments

data	A dataframe containing observed data. Should have column names index for observation times, value for observed responses, and id for curve indicators.
npc	prespecified value for the number of principal components (if given, this overrides pve).
pve	proportion of variance explained; used to choose the number of principal components.
output_index	Grid on which estimates should be computed. Defaults to NULL and returns estimates on the timepoints in the observed dataset
type	Type of estimate for the FPCs; either approx or naive
basis	Basis used as FPCs; in simulations, allows the use of the true FPCs
nbasis	Number of basis functions used in spline expansions

Author(s)

Jan Gertheiss [e-mail]

References

Gertheiss, J., Goldsmith, J., and Staicu, A.-M. (2016). A note on modeling sparse exponential-family functional response curves. Under Review.

See Also

gfpca_TwoStep, gfpca_Bayes.

Examples

## Not run: 
library(mvtnorm)
library(boot)
library(refund.shiny)

## set simulation design elements

bf = 10                           ## number of bspline fns used in smoothing the cov
D = 101                           ## size of grid for observations
Kp = 2                            ## number of true FPC basis functions
grid = seq(0, 1, length = D)

## sample size and sparsity
I <- 300
mobs <- 7:10

## mean structure
mu <- 8*(grid - 0.4)^2 - 3

## Eigenfunctions /Eigenvalues for cov:
psi.true = matrix(NA, 2, D)
psi.true[1,] = sqrt(2)*cos(2*pi*grid)
psi.true[2,] = sqrt(2)*sin(2*pi*grid)

lambda.true = c(1, 0.5)

## generate data

set.seed(1)

## pca effects: xi_i1 phi1(t)+ xi_i2 phi2(t)
c.true = rmvnorm(I, mean = rep(0, Kp), sigma = diag(lambda.true))
Zi = c.true %*% psi.true

Wi = matrix(rep(mu, I), nrow=I, byrow=T) + Zi
pi.true = inv.logit(Wi)  # inverse logit is defined by g(x)=exp(x)/(1+exp(x))
Yi.obs = matrix(NA, I, D)
for(i in 1:I){
  for(j in 1:D){
    Yi.obs[i,j] = rbinom(1, 1, pi.true[i,j])
  }
}

## "sparsify" data
for (i in 1:I)
{
  mobsi <- sample(mobs, 1)
  obsi <- sample(1:D, mobsi)
  Yi.obs[i,-obsi] <- NA
}

Y.vec = as.vector(t(Yi.obs))
subject <- rep(1:I, rep(D,I))
t.vec = rep(grid, I)

data.sparse = data.frame(
  index = t.vec,
  value = Y.vec,
  id = subject
)

data.sparse = data.sparse[!is.na(data.sparse$value),]

## fit models

# 2-step with eigenfunctions from FPCA
fit.2step = gfpca_TwoStep(data = data.sparse, type="naive")
plot(mu)
lines(fit.2step$mu, col=2)

# 2-step with eigenfunctions estimated following Hall et al. (2008)
fit.2step = gfpca_TwoStep(data = data.sparse, type="approx")
plot(mu)
lines(fit.2step$mu, col=2)
plot_shiny(fit.2step)


## End(Not run)

jeff-goldsmith/gfpca documentation built on May 19, 2019, 1:45 a.m.