gfpca_Bayes: gfpca_Bayes

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

Description

Implements a Bayesian approach to generalized functional principal components analysis for sparsely observed binary curves

Usage

gfpca_Bayes(data, npc = 3, output_index = NULL, nbasis = 10,
  iter = 1000, warmup = 400, method = c("vb", "sampling"))

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	Number of FPC basis functions to estimate; defaults to 3
output_index	Grid on which estimates should be computed. Defaults to NULL and returns estimates on the timepoints in the observed dataset
nbasis	Number of basis functions used in spline expansions; defaults to 10
iter	Number of sampler iterations; defaults to 1000
warmup	Number of iterations discarded as warmup; defaults to 400
method	Bayesian fitting method; vb (default) or sampling

Author(s)

Jeff Goldsmith [e-mail] and John Muschelli

References

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

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 Bayesian models
fit.Bayes = gfpca_Bayes(data = data.sparse)
plot(mu)
lines(fit.Bayes$mu, col=2)
plot_shiny(fit.Bayes)

## End(Not run)

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