hd_lfpca | R Documentation |
For V-by-J dimensioanl matrix of Y, HDLPCA first reduce the dimension of p to J if V>J: Y = UDV. Then, longitudinal PCA is applied to the projected scores of Y on the p-dimensional subspace: Y*=DY. The dimension of the subspace is selected based on the variance explained by first p-components. The default value is projectthresh=1.
hd_lfpca(
Y,
T,
J,
I,
visit,
verbose = FALSE,
prefix = date(),
Nx = NA,
Nw = NA,
varthresh = 0.95,
projectthresh = 1,
timeadjust = FALSE
)
Y |
p-by-J matrix, |
T |
Time of theimage collection |
J |
Total number of observations |
I |
Total number of subjects |
visit |
Vector of number of visits per subjects |
verbose |
(default=FALSE) |
Nx |
Dimension of the subject-specific components |
Nw |
Dimension of the subject-visit specific components |
varthresh |
(default=0.99) Threshold for variance explained for both subject-specific and subject-visit specific compoents for dimension selection |
projectthresh |
Threshold for variance explain in the first step of SVD |
timeadjust |
(default=FALSE) Scale time per subject |
xi : Subject-specific principal component scores.
phix0: Subject-specific principal comonent loadings corresponding to intercept.
phix1: Subject-specific principal comonent loadings corresponding to slope.
zeta : Subject-visit-specific principal component scores.
phiw : Subject-visit-specific principal comonent loadings
Seonjoo Lee, sl3670@cumc.columbia.edu
Zipunnikov, V., Greven, S., Shou, H., Caffo, B., Reich, D. S., & Crainiceanu, C. (2014). Longitudinal High-Dimensional Principal Components Analysis with Application to Diffusion Tensor Imaging of Multiple Sclerosis. The Annals of Applied Statistics, 8(4), 2175–2202.
set.seed(12345678)
I=100
visit=rpois(I,1)+3
time = unlist(lapply(visit, function(x) scale(c(0,cumsum(rpois(x-1,1)+1)))))
J = sum(visit)
V=2500
phix0 = matrix(0,V,3);phix0[1:50,1]<-.1;phix0[1:50 + 50,2]<-.1;phix0[1:50 + 100,3]<-.1
phix1 = matrix(0,V,3);phix1[1:50+150,1]<-.1;phix1[1:50 + 200,2]<-.1;phix1[1:50 + 250,3]<-.1
phiw = matrix(0,V,3);phiw[1:50+300,1]<-.1;phiw[1:50 + 350,2]<-.1;phiw[1:50 + 400,3]<-.1
xi = t(matrix(rnorm(I*3),ncol=I)*c(8,4,2))*3
zeta = t(matrix(rnorm(J*3),ncol=J)*c(8,4,2))*2
Y = phix0%*% t(xi[rep(1:I, visit),]) + phix1%*% t(time * xi[rep(1:I, visit),]) + phiw %*% t(zeta) + matrix(rnorm(V*J,0,.1),V,J)
re<-hd_lfpca(Y,T=scale(time,center=TRUE,scale=TRUE),J=J,I=I,visit=visit, varthresh=0.95, projectthresh=1,timeadjust=FALSE)
cor(phix0, re$phix0)
cor(phix1, re$phix1)
library(gplots)
par(mfrow=c(2,2),mar=rep(0.5,4),bg="gray")
bs=c(-100:100)/1000*1.5
image(phix0, axes=F,col=bluered(200),breaks=bs)
image(re$phix0[,1:3], axes=F,col=bluered(200),breaks=bs)
image(phix1, axes=F,col=bluered(200),breaks=bs)
image(re$phix1[,1:3], axes=F,col=bluered(200),breaks=bs)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.