Description Usage Arguments Details Value Author(s) See Also Examples
Estimation of functional linear mixed models (FLMMs) for functional data sampled on equal grids based on functional principal component analysis (FPCA). The implemented models are special cases of the general FLMM
Y_i(t_d) = μ(t_d,x_i) + z_i^\top U(t_d) + ε_i(t_d), i = 1, …,n, d = 1, …, D,
with Y_i(t_d) the value of the response of curve i at observation point
t_d, μ(t_d,x_i) is a mean function, which may depend on covariates
x_i and needs to be estimated beforehand. z_i is a covariate
vector, which is multiplied with the vector of functional random
effects U(t_d).
Usually, the functional random effects will additionally include a smooth error term which
is a functional random intercept with a special structure that captures deviations
from the mean which are correlated along the support of the functions.
In this case, the last block of z_i corresponds to an indicator vector of
indicators for each curve and the last block in U(t) consists of curve-specific
functional random effects.
ε_i(t_d) is independent and identically
distributed white noise measurement error with homoscedastic,
constant variance.
The vector-valued functional random effects can be subdivided into H
independent blocks of functional random effects
U(t_d) = (U_1(t_d)^\top, …, U_H(t_d)^\top)^\top,
with U_g(t_d) and U_h(t_d) independent
for g \neq h. Each block U_h(t_d) further contains L^{U_h} independent
copies U_{gl}(t_d), l=1, …, L^{U_h}, of a vector-valued stochastic process with
ρ^{U_h} vector components U_{gls}(t_d), s = 1,…, ρ^{U_h}.
The total number of functional random effects then amounts to q = ∑_{h=1}^H L^{U_h}ρ^{U_h}.
The code implements a very general functional linear mixed model for n
curves observed at D grid points. Grid points are assumed to be
equidistant and so far no missings are assumed.
The curves are assumed to be centered. The functional random effects for each grouping
factor are assumed to be correlated (e.g., random intercept and
slope curves). The code can handle group-specific functional random
effects including group-specific smooth errors. Covariates are assumed to be standardized.
1 2 3 |
Y |
n x D matrix of n curves observed on D grid points.
|
gridpoints |
vector of grid points along curves, corresponding to columns of |
Zlist |
list of length H of ρ^{U_g} design matrices
Z_{\cdot s^{U_g}},
g = 1,…, H, s = 1,…, ρ^{U_g}. Needed instead of |
G |
number of grouping factors not used for estimation of error variance.
Needed if |
Lvec |
vector of length H containing the number of levels for each grouping factor.
Needed if |
groups |
n \times G matrix with G grouping factors for the rows of |
Zvars |
list of length G with n \times ρ^{U_g} matrices of random variables
for grouping factor g, where G denotes the number of random grouping factors not
used for the estimation of the error variance. Random curves for each grouping
factor are assumed to be correlated (e.g., random intercept and slope).
Set to |
L |
pre-specified level of variance explained (between 0 and 1), determines number of functional principal components. |
NPC |
vector of length H with number of functional principal components to
keep for each functional random effect. Overrides |
smooth |
|
bf |
number of marginal basis functions used for all smooths.
Defaults to |
smoothalg |
smoothing algorithm used for covariance smoothing.
Available options are |
The model fit for centered curves Y_i(.) is
Y = ZU + ε,
with Y = [Y_i(t_d)]_{i = 1, …, n, d = 1, …, D},
Z consisting of
H blocks Z^{U_h} for H grouping factors, Z = [Z^{U_1}|…| Z^{U_H}],
and each Z^{U_h} = [Z_1^{U_h} |…| Z_{ρ^{U_h}}^{U_h}]. U is row-wise divided
into blocks U_1,…, U_H, corresponding to Z.
In case no group-specific functional random effects are specified, column j in Z_{s}^{U_g}, s=1,…,ρ^{U_g},
is assumed to be equal to the sth variable (column) in Zvars[[g]]
times an
indicator for the jth level of grouping factor g, g=1,…,G.
Note that G here denotes the number of random grouping factors not used for the estimation
of the error variance, i.e., all except the smooth error term(s).
The total number of grouping factors is denoted by H.
The estimated eigenvectors and eigenvalues are rescaled to ensure that the approximated eigenfunctions are
orthonormal with respect tot the L^2-inner product.
The estimation of the error variance takes place in two steps. In case of smoothing (smooth = TRUE
),
the error variance is first estimated as the average difference of the raw and the smoothed diagonal values.
In case no smoothing is applied, the estimated error variance is zero.
Subsequent to the eigen decomposition and selection of the eigenfunctions to keep for each grouping factor,
the estimated error variance is recalculated in order to capture the left out variability due to the truncation
of the infinite Karhunen-Loeve expansions.
The function returns a list containing the input arguments
Y
, gridpoints
, groups
, Zvars
, L
, smooth
, bf
,
and smoothalg
. It additionally contains:
Zlist
either the input argument Zlist
or if set to NA
,
the computed list of list of design matrices Z_{\cdot s}^{U_g},
g=1,…,H, s = 1,…, ρ^{U_g}.
G
either the input argument G
or if set to NA
,
the computed number of random grouping factors G not used for the estimation of the error variance.
Lvec
either the input argument Lvec
or if set to NA
,
the computed vector of length H containing the number of levels
for each grouping factor (including the smooth error(s)).
NPC
either the input argument NPC
or if set to NA
,
the number of functional principal components kept for each functional random effect (including the smooth error(s)).
rhovec
vector of length H of number of random effects for each grouping factor (including the smooth error(s)).
phi
list of length H of D x N^{U_g} matrices containing
the N^{U_g} functional principal components kept per grouping factor (including the smooth error(s)).
sigma2
estimated measurement error variance σ^2.
nu
list of length H of N^{U_g} x 1 vectors of estimated eigenvalues
ν_k^{U_g}.
xi
list of length H of L^{U_g} x N^{U_g} matrices containing
the predicted random basis weights. Within matrices, basis weights are ordered corresponding
to the ordered levels of the grouping factors, e.g., a grouping factor with levels 4, 2, 3, 1
(L^{U_g} = 4) will result in rows in xi[[g]]
corresponding to levels 1, 2, 3, 4.
totvar
total average variance of the curves.
exvar
level of variance explained by the selected functional principal components (+ error variance).
Sonja Greven, Jona Cederbaum
sparseFLMM
for the estimation of functional linear mixed models for irregularly
or sparsely sampled data based on functional principal component analysis.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 | # fit model with group-specific functional random intercepts for two groups
# and a non group-specific smooth error, i.e., G = 2, H = 1.
################
# load libraries
################
require(mvtnorm)
require(Matrix)
set.seed(123)
#########################
# specify data generation
#########################
nus <- list(c(0.5, 0.3), c(1, 0.4), c(2)) # eigenvalues
sigmasq <- 2.5e-05 # error variance
NPCs <- c(rep(2, 2), 1) # number of eigenfunctions
Lvec <- c(rep(2, 2), 480) # number of levels
H <- 3 # number of functional random effects (in total)
G <- 2 # number of functional random effects not used for the estimation of the error variance
gridpoints <- seq(from = 0, to = 1, length = 100) # grid points
class_nr <- 2 # number of groups
# define eigenfunctions
#######################
funB1<-function(k,t){
if(k == 1)
out <- sqrt(2) * sin(2 * pi * t)
if(k == 2)
out <- sqrt(2) * cos(2 * pi * t)
out
}
funB2<-function(k,t){
if(k == 1)
out <- sqrt(7) * (20 * t^3 - 30 * t^2 + 12 * t - 1)
if(k == 2)
out <- sqrt(3) * (2 * t - 1)
out
}
funE<-function(k,t){
if(k == 1)
out <- 1 + t - t
if(k == 2)
out <- sqrt(5) * (6 * t^2 - 6 * t + 1)
out
}
###############
# generate data
###############
D <- length(gridpoints) # number of grid points
n <- Lvec[3] # number of curves in total
class <- rep(1:class_nr, each = n / class_nr)
group1 <- rep(rep(1:Lvec[1], each = n / (Lvec[1] * class_nr)), class_nr)
group2 <- 1:n
data <- data.frame(class = class, group1 = group1, group2 = group2)
# get eigenfunction evaluations
###############################
phis <- list(sapply(1:NPCs[1], FUN = funB1, t = gridpoints),
sapply(1:NPCs[2], FUN = funB2, t = gridpoints),
sapply(1:NPCs[3], FUN = funE, t = gridpoints))
# draw basis weights
####################
xis <- vector("list", H)
for(i in 1:H){
if(NPCs[i] > 0){
xis[[i]] <- rmvnorm(Lvec[i], mean = rep(0, NPCs[i]), sigma = diag(NPCs[i]) * nus[[i]])
}
}
# construct functional random effects
#####################################
B1 <- xis[[1]] %*% t(phis[[1]])
B2 <- xis[[2]] %*% t(phis[[2]])
E <- xis[[3]] %*% t(phis[[3]])
B1_mat <- B2_mat <- E_mat <- matrix(0, nrow = n, ncol = D)
B1_mat[group1 == 1 & class == 1, ] <- t(replicate(n = n / (Lvec[1] * class_nr),
B1[1, ], simplify = "matrix"))
B1_mat[group1 == 2 & class == 1, ] <- t(replicate(n = n / (Lvec[1] * class_nr),
B1[2, ], simplify = "matrix"))
B2_mat[group1 == 1 & class == 2, ] <- t(replicate(n = n / (Lvec[1] * class_nr),
B2[1, ], simplify = "matrix"))
B2_mat[group1 == 2 & class == 2, ] <- t(replicate(n = n / (Lvec[1] * class_nr),
B2[2, ], simplify = "matrix"))
E_mat <- E
# draw white noise measurement error
####################################
eps <- matrix(rnorm(n * D, mean = 0, sd = sqrt(sigmasq)), nrow = n, ncol = D)
# construct curves
##################
Y <- B1_mat + B2_mat + E_mat + eps
#################
# construct Zlist
#################
Zlist <- list()
Zlist[[1]] <- Zlist[[2]] <- Zlist[[3]] <- list()
d1 <- data.frame(a = as.factor(data$group1[data$class == 1]))
Zlist[[1]][[1]] <- rBind(sparse.model.matrix(~ -1 + a, d1),
matrix(0, nrow = (1 / class_nr * n), ncol = (Lvec[1])))
d2 <- data.frame(a = as.factor(data$group1[data$class == 2]))
Zlist[[2]][[1]] <- rBind(matrix(0, nrow = (1 / class_nr * n),
ncol = (Lvec[2])), sparse.model.matrix(~ -1 + a, d2))
d3 <- data.frame(a = as.factor(data$group2))
Zlist[[3]][[1]] <- sparse.model.matrix(~ -1 + a, d3)
################
# run estimation
################
results <- denseFLMM(Y = Y, gridpoints = gridpoints, Zlist = Zlist, G = G, Lvec = Lvec,
groups = NA, Zvars = NA, L = 0.99999, NPC = NA,
smooth = FALSE)
###############################
# plot estimated eigenfunctions
###############################
with(results, matplot(gridpoints, phi[[1]], type = "l"))
with(results, matplot(gridpoints, phi[[2]], type = "l"))
with(results, matplot(gridpoints, phi[[3]], type = "l"))
|
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.