if (FALSE) {
##################################################################
######### DTI Data Example #########
##################################################################
##################################################################
# For more about this example, see Swihart et al. 2013
##################################################################
## load and reassign the data;
data(DTI2)
Y <- DTI2$pasat ## PASAT outcome
id <- DTI2$id ## subject id
W1 <- DTI2$cca ## Corpus Callosum
W2 <- DTI2$rcst ## Right corticospinal
V <- DTI2$visit ## visit
## prep scalar covariate
visit.1.rest <- matrix(as.numeric(V > 1), ncol=1)
covar.in <- visit.1.rest
## note there is missingness in the functional predictors
apply(is.na(W1), 2, mean)
apply(is.na(W2), 2, mean)
## fit two univariate models
pfr.obj.t1 <- pfr(Y = Y, covariates=covar.in, funcs = list(W1), subj = id, kz = 10, kb = 50)
pfr.obj.t2 <- pfr(Y = Y, covariates=covar.in, funcs = list(W2), subj = id, kz = 10, kb = 50)
### one model with two functional predictors using "smooth.face"
### for smoothing predictors
pfr.obj.t3 <- pfr(Y = Y, covariates=covar.in, funcs = list(W1, W2),
subj = id, kz = 10, kb = 50, nbasis=35,smooth.option="fpca.face")
## plot the coefficient function and bounds
dev.new()
par(mfrow=c(2,2))
ran <- c(-2,.5)
matplot(cbind(pfr.obj.t1$BetaHat[[1]], pfr.obj.t1$Bounds[[1]]),
type = 'l', lty = c(1,2,2), col = c(1,2,2), ylab = "BetaHat",
main = "CCA", xlab="Location", ylim=ran)
abline(h=0, col="blue")
matplot(cbind(pfr.obj.t2$BetaHat[[1]], pfr.obj.t2$Bounds[[1]]),
type = 'l', lty = c(1,2,2), col = c(1,2,2), ylab = "BetaHat",
main = "RCST", xlab="Location", ylim=ran)
abline(h=0, col="blue")
matplot(cbind(pfr.obj.t3$BetaHat[[1]], pfr.obj.t3$Bounds[[1]]),
type = 'l', lty = c(1,2,2), col = c(1,2,2), ylab = "BetaHat",
main = "CCA - mult.", xlab="Location", ylim=ran)
abline(h=0, col="blue")
matplot(cbind(pfr.obj.t3$BetaHat[[2]], pfr.obj.t3$Bounds[[2]]),
type = 'l', lty = c(1,2,2), col = c(1,2,2), ylab = "BetaHat",
main = "RCST - mult.", xlab="Location", ylim=ran)
abline(h=0, col="blue")
##################################################################
# use baseline data to regress continuous outcomes on functional
# predictors (continuous outcomes only recorded for case == 1)
##################################################################
data(DTI)
# subset data as needed for this example
cca = DTI$cca[which(DTI$visit ==1 & DTI$case == 1),]
rcst = DTI$rcst[which(DTI$visit ==1 & DTI$case == 1),]
DTI = DTI[which(DTI$visit ==1 & DTI$case == 1),]
# note there is missingness in the functional predictors
apply(is.na(cca), 2, mean)
apply(is.na(rcst), 2, mean)
# fit two models with single functional predictors and plot the results
fit.cca = pfr(Y=DTI$pasat, funcs = cca, kz=10, kb=50, nbasis=20)
fit.rcst = pfr(Y=DTI$pasat, funcs = rcst, kz=10, kb=50, nbasis=20)
par(mfrow = c(1,2))
matplot(cbind(fit.cca$BetaHat[[1]], fit.cca$Bounds[[1]]),
type = 'l', lty = c(1,2,2), col = c(1,2,2), ylab = "BetaHat",
main = "CCA")
matplot(cbind(fit.rcst$BetaHat[[1]], fit.rcst$Bounds[[1]]),
type = 'l', lty = c(1,2,2), col = c(1,2,2), ylab = "BetaHat",
main = "RCST")
# fit a model with two functional predictors and plot the results
fit.cca.rcst = pfr(Y=DTI$pasat, funcs = list(cca, rcst), kz=10, kb=30, nbasis=20)
par(mfrow = c(1,2))
matplot(cbind(fit.cca.rcst$BetaHat[[1]], fit.cca.rcst$Bounds[[1]]),
type = 'l', lty = c(1,2,2), col = c(1,2,2), ylab = "BetaHat",
main = "CCA")
matplot(cbind(fit.cca.rcst$BetaHat[[2]], fit.cca.rcst$Bounds[[2]]),
type = 'l', lty = c(1,2,2), col = c(1,2,2), ylab = "BetaHat",
main = "RCST")
##################################################################
# use baseline data to regress binary case-status outcomes on
# functional predictors
##################################################################
data(DTI)
# subset data as needed for this example
cca = DTI$cca[which(DTI$visit == 1),]
rcst = DTI$rcst[which(DTI$visit == 1),]
DTI = DTI[which(DTI$visit == 1),]
# fit two models with single functional predictors and plot the results
fit.cca = pfr(Y=DTI$case, funcs = cca, family = "binomial")
fit.rcst = pfr(Y=DTI$case, funcs = rcst, family = "binomial")
par(mfrow = c(1,2))
matplot(cbind(fit.cca$BetaHat[[1]], fit.cca$Bounds[[1]]),
type = 'l', lty = c(1,2,2), col = c(1,2,2), ylab = "BetaHat",
main = "CCA")
matplot(cbind(fit.rcst$BetaHat[[1]], fit.rcst$Bounds[[1]]),
type = 'l', lty = c(1,2,2), col = c(1,2,2), ylab = "BetaHat",
main = "RCST")
##################################################################
######### Octane Data Example #########
##################################################################
data(gasoline)
Y = gasoline$octane
funcs = gasoline$NIR
wavelengths = as.matrix(2*450:850)
# fit the model using pfr and the smoothing option "fpca.face"
fit = pfr(Y=Y, funcs=funcs, kz=15, kb=50,nbasis=35,smooth.option="fpca.face")
matplot(wavelengths, cbind(fit$BetaHat[[1]], fit$Bounds[[1]]),
type='l', lwd=c(2,1,1), lty=c(1,2,2), xlab = "Wavelengths",
ylab = "Coefficient Function", col=1)
}
Run the code above in your browser using DataLab