if (FALSE) {
data(cd4)
# obtain a subsample of the data with 25 subjects
set.seed(1236)
sample = sample(1:dim(cd4)[1], 25)
Y.sub = cd4[sample,]
# obtain a mixed-model based FPCA decomposition
Fit.MM = fpca.sc(Y.sub, var = TRUE, simul = TRUE)
# use iterated variance to obtain curve estimates and variances
Fit.IV = ccb.fpc(Y.sub, n.boot = 25, simul = TRUE)
# for one subject, examine curve estimates, pointwise and simultaneous itervals
EX = 2
EX.IV = cbind(Fit.IV$Yhat[EX,],
Fit.IV$Yhat[EX,] + 1.96 * sqrt(Fit.IV$diag.var[EX,]),
Fit.IV$Yhat[EX,] - 1.96 * sqrt(Fit.IV$diag.var[EX,]),
Fit.IV$Yhat[EX,] + Fit.IV$crit.val[EX] * sqrt(Fit.IV$diag.var[EX,]),
Fit.IV$Yhat[EX,] - Fit.IV$crit.val[EX] * sqrt(Fit.IV$diag.var[EX,]))
EX.MM = cbind(Fit.MM$Yhat[EX,],
Fit.MM$Yhat[EX,] + 1.96 * sqrt(Fit.MM$diag.var[EX,]),
Fit.MM$Yhat[EX,] - 1.96 * sqrt(Fit.MM$diag.var[EX,]),
Fit.MM$Yhat[EX,] + Fit.MM$crit.val[EX] * sqrt(Fit.MM$diag.var[EX,]),
Fit.MM$Yhat[EX,] - Fit.MM$crit.val[EX] * sqrt(Fit.MM$diag.var[EX,]))
# plot data for one subject, with curve and interval estimates
d = as.numeric(colnames(cd4))
plot(d[which(!is.na(Y.sub[EX,]))], Y.sub[EX,which(!is.na(Y.sub[EX,]))], type = 'o',
pch = 19, cex=.75, ylim = range(0, 3400), xlim = range(d),
xlab = "Months since seroconversion", lwd = 1.2, ylab = "Total CD4 Cell Count",
main = "Est. & CI - Sampled Data")
matpoints(d, EX.IV, col = 2, type = 'l', lwd = c(2, 1, 1, 1, 1), lty = c(1,1,1,2,2))
matpoints(d, EX.MM, col = 4, type = 'l', lwd = c(2, 1, 1, 1, 1), lty = c(1,1,1,2,2))
legend("topright", c("IV Est", "IV PW Int", "IV Simul Int",
expression(paste("MM - ", hat(theta), " Est", sep = "")),
expression(paste("MM - ", hat(theta), " PW Int", sep = "")),
expression(paste("MM - ", hat(theta), " Simul Int", sep = ""))),
lty=c(1,1,2,1,1,2), lwd = c(2.5,.75,.75,2.5,.75,.75),
col = c("red","red","red","blue","blue","blue"))
}
Run the code above in your browser using DataLab