#' Interactive Plotting for Multilevel Functional Principal Component Analysis
#'
#' Produces an interactive plot illustrating a multilevel functional principal component
#' analysis.
#'
#' @param obj mfpca object to be plotted.
#' @param xlab x axis label
#' @param ylab y axis label
#' @param title plot title
#' @param ... additional arguments passed to plotting functions
#'
#' @author Julia Wrobel \email{julia.wrobel@@cuanschutz.edu},
#' Jeff Goldsmith \email{jeff.goldsmith@@columbia.edu}
#'
#' @seealso \code{\link{plot_shiny}}
#' @import shiny
#' @import ggplot2
#' @importFrom gridExtra grid.arrange
#' @importFrom dplyr mutate
#' @importFrom reshape2 melt
#'
#' @export
#' @return No object is returned. This function takes in objects of class 'mfpca' and outputs a shiny application for that object.
#'
plot_shiny.mfpca = function(obj, xlab = "", ylab="", title = "", ...) {
mfpca.obj <- obj
### NULLify global values called in ggplot
Y = id = k = lambda = value = mu_visit = mu_subj = visit = subj = time = grid = muPC_call1 = muPC_call2 = muPCtext = NULL
################################
## code for processing tabs
################################
npc = mfpca.obj$npc
mu = data.frame(grid = 1:length(mfpca.obj$mu), value = mfpca.obj$mu)
efunctions = mfpca.obj$efunctions;
sqrt.evalues = lapply(mfpca.obj$evalues, function(i) diag(sqrt(i)))
scaled_efunctions = lapply(1:2, function(i) efunctions[[i]] %*% sqrt.evalues[[i]])
names(scaled_efunctions) <- names(sqrt.evalues) <- names(efunctions) <- c("level1", "level2")
#################################
## Tab 3: mean +/- FPC plot
#################################
muPC.help = "Solid black line indicates population mean. For the FPC selected below, blue and red lines
indicate the population mean +/- the FPC times 2 SDs of the associated score distribution."
muPC.call = as.list(rep(NA, 2))
muPC.call[[1]] = selectInput(inputId = "PCchoice1", label = ("Select Level 1 FPC"), choices = 1:npc$level1, selected = 1)
muPC.call[[2]] = selectInput(inputId = "PCchoice2", label = ("Select Level 2 FPC"), choices = 1:npc$level2, selected = 1)
#################################
## Tab 2: scree plot
#################################
scree <- lapply(names(npc), function(level) {
data.frame(k = rep(1:npc[[level]], 2),
lambda = c(mfpca.obj$evalues[[level]], cumsum(mfpca.obj$evalues[[level]])/ sum(mfpca.obj$evalues[[level]])),
type = rep(c("Score Variance", "Percent Variance Explained"), each = npc[[level]]))
})
names(scree) <- c("level1", "level2")
levelVariance = lapply(1:2, function(i) round(sum(mfpca.obj$evalues[[i]])/(sum(mfpca.obj$evalues[[1]])+sum(mfpca.obj$evalues[[2]])), 3)*100)
scree.help1 = paste0("Scree plots for level 1; the left panel shows the plot of eigenvalues and the right panel shows the
cumulative percent variance explained. Level 1 accounts for ", levelVariance[[1]], "% of total variance." )
scree.help2 = paste0("Scree plots for level 2; the left panel shows the plot of eigenvalues and the right panel shows the
cumulative percent variance explained. Level 2 accounts for ", levelVariance[[2]], "% of total variance." )
#################################
## Tab 3:Linear combination of PCs
#################################
varpercents = lapply(c(1, 2, 12), function(i) varPercent(i, mfpca.obj))
plot.npc = list(min(3, mfpca.obj$npc[[1]]), min(3, mfpca.obj$npc[[2]]))
mfpca.calls = mfpcaCalls(plot.npc = plot.npc, mfpca.obj, varpercents)
calls <- mfpca.calls$calls
PCs <- mfpca.calls$PCs
#################################
## Tab 4: subject fits
#################################
ids = unique(mfpca.obj$Y.df$id)
Y.df = mfpca.obj$Y.df
Yhat.subj = mfpca.obj$Yhat.subject
Yhat = mfpca.obj$Yhat
rownames(Y.df) = rownames(Yhat) = rownames(Yhat.subj)
#################################
## Tab 5: score plots
#################################
scoreTextA = "Use the drop down menus to select FPCs for the X and Y axis. Plot shows observed
score scatterplot for selected FPCs; click and drag on the scatterplot to select subjects."
## repeat level 1 scores so they are the same dimension as input dataset
scores_new = list(mfpca.obj$scores$level1[rep(seq(nrow(mfpca.obj$scores$level1)), data.frame(table(Y.df$id))$Freq),], mfpca.obj$scores$level2)
Yhat.level = lapply(1:2, function(i) scores_new[[i]] %*% t(efunctions[[i]]))
scoredata = as.list(rep(NA,2))
for(i in 1:2){
scoredata[[i]] = data.frame(scores_new[[i]])
colnames(scoredata[[i]]) = c(paste0("PC", 1:npc[[i]]))
scoredata[[i]]$Yhat.level = Yhat.level[[i]]
scoredata[[i]]$Yhat = Yhat
}
Yhat.all.m = melt(mfpca.obj$Yhat)
Yhat.level.all.m = lapply(1:2, function(i) melt(scoredata[[i]]$Yhat.level) )
colnames(Yhat.level.all.m[[1]]) = colnames(Yhat.level.all.m[[2]]) = colnames(Yhat.all.m) = c("subj", "time", "value")
####### set some defaults for ggplot
plotDefaults = list(theme_bw(), xlab(xlab), ylab(ylab), ylim(c(range(Yhat)[1], range(Yhat)[2])),
scale_x_continuous(breaks = seq(0, length(mfpca.obj$mu)-1, length=6),
labels = paste0(c(0, 0.2, 0.4, 0.6, 0.8, 1))) )
#################################
## App
#################################
shinyApp(
#################################
## UI
#################################
ui = navbarPage(title = strong(style = "color: #ACD6FF; padding: 0px 0px 10px 10px; opacity: 0.95; ", "MFPCA Plot"),
windowTitle = "refund.shiny", collapsible = FALSE, id = "nav", inverse = TRUE, header = NULL,
tabPanel("Mean +/- FPCs", icon = icon("stats", lib = "glyphicon"),
tabsetPanel(
tabPanelModuleUI("muPC1", tabTitle = "Level 1", calls = muPC.call[[1]], helperText = muPC.help ),
tabPanelModuleUI("muPC2", tabTitle = "Level 2", calls = muPC.call[[2]], helperText = muPC.help )
) ),
tabPanel("Scree Plot", icon = icon("medkit"),
tabsetPanel(
tabPanelModuleUI("scree1", tabTitle = "Level 1", helperText = scree.help1 ),
tabPanelModuleUI("scree2", tabTitle = "Level 2", helperText = scree.help2 )
) ),
tabPanel("Linear Combinations", icon = icon("chart-line"), withMathJax(),
column(3, h4("Sliders for Levels 1 and 2"),
helpText("Plot shows the linear combination of mean and FPCs with the scores
specified using the sliders below. Black curve is population mean; blue curve
is subject-specific mean; red curve is subject-visit specific mean."), hr(),
tabsetPanel(
tabPanel("Level 1", eval(calls[[1]]) ),
tabPanel("Level 2", eval(calls[[2]]) )
)
),
column(9, h4("Linear Combination of Mean and FPCs"), plotOutput('LinCom') )
),
tabPanel("Subject Fits", icon = icon("user"),
column(3,
helpText("Plot shows observed data and fitted values for the subject selected below. Blue curve is
subject-specific mean, red curves are subject-visit specific fitted values, and red
points are observed data."), hr(),
selectInput("subject", label = ("Select Subject"), choices = ids, selected =ids[1]),
uiOutput("visitnum"),
checkboxInput("colorVisit", label="Color by Visit", value =FALSE),
helpText("If 'Color by Visit' is selected, observed values and subject-visit specific fitted values
are colored by visit number.") ),
column(9, h4("Fitted and Observed Values for Selected Subject"), plotOutput("Subject"))
),
tabPanel("Score Scatterplot",icon = icon("binoculars"),
tabsetPanel(
tabPanel("Level 1",
fluidRow(
column(3,
helpText(scoreTextA), hr(),
selectInput("PCX1", label = ("Select X-axis FPC"), choices = 1:npc[[1]], selected = 1),
selectInput("PCY1", label = ("Select Y-axis FPC"), choices = 1:npc[[1]], selected = 2), hr()
), ## end column3
column(9, h4("Score Scatterplot for Selected FPCs"),
plotOutput("ScorePlot1",
brush=brushOpts(
id = "ScorePlotL1_brush",
resetOnNew = TRUE)
)
) ## end column9
),
fluidRow(
column(3, helpText("The left panel shows all fitted curves, for all subjects and visits; the
right panel shows fitted curves using only Level 1 FPCs -- the mean and
Level 2 effects are omitted. Blue curves correspond to subjects selected
in the graph above.")
),
column(9, plotOutput("YhatPlot1"))
)
), ## end tabPanel
tabPanel("Level 2",
fluidRow(
column(3,
helpText(scoreTextA), hr(),
selectInput("PCX2", label = ("Select X-axis FPC"), choices = 1:npc[[2]], selected = 1),
selectInput("PCY2", label = ("Select Y-axis FPC"), choices = 1:npc[[2]], selected = 2), hr()
), ## end column3
column(9, h4("Score Scatterplot for Selected FPCs"),
plotOutput("ScorePlot2",
brush=brushOpts(
id = "ScorePlotL2_brush",
resetOnNew = TRUE)
)
) ## end column9
),
fluidRow(
column(3, helpText("The left panel shows all fitted curves, for all subjects and visits; the
right panel shows fitted curves using only Level 2 FPCs -- the mean and
Level 1 effects are omitted. Blue curves correspond to subjects selected
in the graph above.")
),
column(9, plotOutput("YhatPlot2"))
)
) ## end tabPanel
) ## end tabsetPanel
) # end big tabPanel
), ## end UI
#################################
## Server
#################################
server = function(input, output){
#################################
## Code for mu PC plot
#################################
plotInputMuPC <- reactive({
PCchoice = list(as.numeric(input$PCchoice1), as.numeric(input$PCchoice2))
names(PCchoice) <- c("level1", "level2")
scaled_efuncs = lapply(1:2, function(i) scaled_efunctions[[i]][,PCchoice[[i]]])
p1 <- lapply(1:2, function(i){
ggplot(mu, aes(x = grid, y = value)) + geom_line(lwd=1) + plotDefaults +
geom_point(data = data.frame(grid =1:length(mfpca.obj$mu),value = mfpca.obj$mu + 2*scaled_efuncs[[i]]), color = "blue", size = 4, shape = '+') +
geom_point(data = data.frame(grid =1:length(mfpca.obj$mu), value = mfpca.obj$mu - 2*scaled_efuncs[[i]]), color = "indianred", size = 4, shape = "-") +
ggtitle(bquote(psi[.(PCchoice[[i]])]~(t) ~ "," ~.(100*round(mfpca.obj$evalues[[i]][PCchoice[[i]]]/sum(mfpca.obj$evalues[[i]]),3)) ~ "% Variance"))
})
})
plotInputMuPC1 <- reactive({ plotInputMuPC()[[1]] }) ; plotInputMuPC2 <- reactive({ plotInputMuPC()[[2]] })
callModule(tabPanelModule, "muPC1", plotObject = plotInputMuPC1, plotName = "muPC1")
callModule(tabPanelModule, "muPC2", plotObject = plotInputMuPC2, plotName = "muPC2")
#################################
## Code for scree plot
#################################
plotInputScree <- reactive({
p2 <-screeplots <- lapply(scree, function(i){
ggplot(i, aes(x=k, y=lambda))+geom_line(linetype=1, lwd=1.5, color="black")+
geom_point(size = 4, color = "black")+ theme_bw() + xlab("Principal Component") + ylab("") +
facet_wrap(~type, scales = "free_y") + ylim(0, NA)
})
})
plotInputScree1 <- reactive({ plotInputScree()[[1]] }) ; plotInputScree2 <- reactive({ plotInputScree()[[2]] })
callModule(tabPanelModule, "scree1", plotObject = plotInputScree1, plotName = "scree1")
callModule(tabPanelModule, "scree2", plotObject = plotInputScree2, plotName = "scree2")
#################################
## Code for linear combinations
#################################
plotInputLinCom <- reactive({
PCweights = lapply(1:2, function(i) rep(NA, plot.npc[[i]])) ;
names(PCweights) <- c("level1", "level2")
for(i in 1:plot.npc[[1]]){PCweights$level1[i] = input[[PCs$level1[i]]]}
for(i in 1:plot.npc[[2]]){PCweights$level2[i] = input[[PCs$level2[i]]]}
df = data.frame(1:length(mfpca.obj$mu),
as.matrix(mfpca.obj$mu)+efunctions$level1[,1:plot.npc[[1]]] %*% sqrt.evalues$level1[1:plot.npc[[1]], 1:plot.npc[[1]]] %*% PCweights$level1 +
efunctions$level2[,1:plot.npc[[2]]] %*% sqrt.evalues$level2[1:plot.npc[[2]], 1:plot.npc[[2]]] %*% PCweights$level2,
as.matrix(mfpca.obj$mu)+efunctions$level1[,1:plot.npc[[1]]] %*% sqrt.evalues$level1[1:plot.npc[[1]], 1:plot.npc[[1]]] %*% PCweights$level1 )
names(df) = c("grid", "mu_visit", "mu_subj")
p3 <- ggplot(mu, aes(x=grid, y=value))+geom_line(lwd=0.75, aes(color= "mu"))+ plotDefaults + theme(legend.key = element_blank())+
geom_line(data = df, lwd = 1, aes(x=grid, y = mu_visit, color = "visit")) +
geom_line(data = df, lwd = 1.5, aes(x=grid, y = mu_subj, color = "subject")) +
scale_color_manual("Line Legend", values = c(mu = "gray", visit = "indianred", subject = "cornflowerblue"), guide = FALSE)
})
output$LinCom <- renderPlot( print(plotInputLinCom()) )
#################################
## Code for subject plots
#################################
dataInputSubject <- reactive({
id.cur = as.numeric(input$subject)
df.obs = mutate(melt(as.matrix(subset(Y.df, id == id.cur, select = Y))), grid=rep(1:ncol(Y.df$Y), each = length(which(Y.df$id == id.cur))))
df.Yhat.subj = mutate(melt(as.matrix(subset(Yhat.subj, Y.df$id == id.cur))), grid=rep(1:ncol(Y.df$Y), each = length(which(Y.df$id == id.cur))))
df.Yhat = mutate(melt(as.matrix(subset(Yhat, Y.df$id == id.cur))), grid=rep(1:ncol(Y.df$Y), each = length(which(Y.df$id == id.cur))))
names(df.obs) = names(df.Yhat.subj) = names(df.Yhat) = c("visit", "Ynames", "value", "time"); df.obs$visit <- df.Yhat$visit
numVisits = length(unique(df.Yhat$visit))
return(list(df.obs, df.Yhat, df.Yhat.subj, numVisits))
})
output$visitnum <- renderUI({ selectInput("visit","Select Visit", c(Choose='', 1:dataInputSubject()[[4]]), selectize = TRUE, multiple = TRUE)})
plotInputSubject <- reactive({
p4 <- ggplot(dataInputSubject()[[1]], aes(x = time, y = value, group = visit)) + plotDefaults +
geom_line(data = mu, aes(x = grid, y = value, group=NULL), col="gray")
visits <- input[["visit"]]
df.Yhat.visits <- subset(dataInputSubject()[[2]], visit %in% visits)
df.obs.visits <- subset(dataInputSubject()[[1]], visit %in% visits)
if(length(visits)>=1){p4 = p4 + geom_path(data = dataInputSubject()[[3]], col="cornflowerblue") +
geom_point(data = df.obs.visits, col = "indianred") +
geom_path(data = df.Yhat.visits, col="indianred")
}
else if(input$colorVisit) {p4 = p4 + geom_point(aes(col = factor(visit))) + geom_path(data=dataInputSubject()[[2]], aes(col = factor(visit)))+
theme(legend.position="none") + geom_path(data = dataInputSubject()[[3]], col="cornflowerblue", lwd = 1.25)
}
else{p4 = p4 + geom_path(data=dataInputSubject()[[2]], col = "indianred") + geom_point(col = "indianred", alpha = 1/5) +
geom_path(data = dataInputSubject()[[3]], col="cornflowerblue", lwd = 1.25)
}
})
output$Subject <- renderPlot( print(plotInputSubject()) )
#################################
## Code for score plots
#################################
PCX <- reactive({ list(level1 = paste0("PC", input$PCX1), level2 = paste0("PC", input$PCX2) ) })
PCY <- reactive({ list(level1 = paste0("PC", input$PCY1), level2 = paste0("PC", input$PCY2) ) })
## Level 1 scatter plot of scores
output$ScorePlot1 <- renderPlot({
ggplot(scoredata[[1]], aes_string(x = PCX()[[1]], y = PCY()[[1]]))+geom_point(color = "blue", alpha = 1/5, size = 3)+
theme_bw()+xlab(paste("Scores for FPC", input$PCX1))+ylab(paste("Scores for FPC", input$PCY1))
})
#########################################
# Store Plots that use Brush Values
#######################################
YhatPlots <- reactive({
brushes <- list(input$ScorePlotL1_brush, input$ScorePlotL2_brush)
plots = list()
for (i in 1:2){
if(!is.null(brushes[[i]])){
points = brushedPoints(scoredata[[i]], brushes[[i]], xvar=PCX()[[i]], yvar = PCY()[[i]])
Yhat.m = melt(as.matrix(points$Yhat))
Yhat.level.m = melt(as.matrix(points$Yhat.level))
colnames(Yhat.level.m) <- c("subj", "time", "value")
}else{ Yhat.m = data.frame(1, 1:length(mfpca.obj$mu), mfpca.obj$mu) }
colnames(Yhat.m) <- c("subj", "time", "value")
plots[[i]] <- ggplot(Yhat.all.m, aes(x=time, y=value, group = subj)) + geom_line(alpha = 1/5, color="black") + plotDefaults+
geom_line(data = Yhat.m, aes(x=as.numeric(time), y=value, group = subj), color="cornflowerblue")
plots[[i+2]] <- ggplot(Yhat.level.all.m[[i]], aes(x=time, y=value, group = subj)) + geom_line(alpha = 1/5, color="black") +
plotDefaults[-c(4)] + ylim(c(range(scoredata[[i]]$Yhat.level)[1], range(scoredata[[1]]$Yhat.level)[2]))
if(!is.null(brushes[[i]])){
plots[[i+2]] = plots[[i+2]] + geom_line(data = Yhat.level.m, aes(x=as.numeric(time), y=value, group = subj), color="cornflowerblue")
} else{ plots[[i+2]] = plots[[i+2]] }
}
return(plots)
})
output$YhatPlot1 <- renderPlot({
grid.arrange(YhatPlots()[[1]],YhatPlots()[[3]], ncol = 2)
})
## Level 2 Plot 1
output$ScorePlot2 <- renderPlot({
ggplot(scoredata[[2]], aes_string(x = PCX()[[2]], y = PCY()[[2]]))+geom_point(color = "blue", alpha = 1/5, size = 3)+theme_bw()+
xlab(paste("Scores for FPC", input$PCX2))+ylab(paste("Scores for FPC", input$PCY2))
})
output$YhatPlot2 <- renderPlot({
grid.arrange(YhatPlots()[[2]],YhatPlots()[[4]], ncol = 2)
})
} ## end server
)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.