# NOT RUN {
data(simdat)
# add hypothetical correlated term:
simdat$predictor <- (simdat$Trial+10)^.75 + rnorm(nrow(simdat))
# principal components analysis:
pca <- prcomp(simdat[, c('Trial', 'predictor')])
# only first PC term contributes:
summary(pca)
# get rotation (weights of predictors in PC):
pcar <- pca$rotation
# add PC1 to data:
simdat$PC1 <- pca$x[,1]
# }
# NOT RUN {
# model:
m1 <- bam(Y ~ Group + te(Time, PC1, by=Group)
+ s(Time, Subject, bs='fs', m=1, k=5), data=simdat)
# inspect surface:
fvisgam(m1, view=c('Time', 'PC1'), cond=list(Group='Children'),
rm.ranef=TRUE)
# how does Trial contribute?
p <- get_pca_predictions(m1, pca.term='PC1', weights=pcar[,'PC1'],
view=c('Time', 'Trial'), cond=list(Group='Children'),
rm.ranef=TRUE, partial=FALSE)
# Note that the range of Trial is estimated based on the values of PC1.
# A better solution is to specify the range:
p <- get_pca_predictions(m1, pca.term='PC1', weights=pcar[,'PC1'],
view=list(Time=range(simdat$Time), Trial=range(simdat$Trial)),
cond=list(Group='Children'),rm.ranef=TRUE, partial=FALSE)
# plotting of the surface:
plot_pca_surface(m1, pca.term='PC1', weights=pcar[,'PC1'],
view=c('Time', 'Trial'), cond=list(Group='Children'),rm.ranef=TRUE)
# }
Run the code above in your browser using DataLab