Nothing
R/dynrPlot.R
In dynr: Dynamic Models with Regime-Switching
Defines functions
autoplot.dynrTaste
autoplot.dynrCook
dynr.ggplot
plotFormula
dynr.plotFreq
plotdf
plot.dynrCook
Mode
findRegime
multiplot
Documented in
autoplot.dynrCook
autoplot.dynrTaste
dynr.ggplot
dynr.plotFreq
plot.dynrCook
plotFormula
#------------------------------------------------------------------------------
# Plotting methods
multiplot <- function(..., plotlist=NULL, file, cols=1, layout=NULL) {
# Make a list from the ... arguments and plotlist
plots <- c(list(...), plotlist)
numPlots = length(plots)
# If layout is NULL, then use 'cols' to determine layout
if (is.null(layout)) {
# Make the panel
# ncol: Number of columns of plots
# nrow: Number of rows needed, calculated from # of cols
layout <- matrix(seq(1, cols * ceiling(numPlots/cols)),
ncol = cols, nrow = ceiling(numPlots/cols))
}
if (numPlots==1) {
print(plots[[1]])
} else {
# Set up the page
grid::grid.newpage()
grid::pushViewport(grid::viewport(layout = grid::grid.layout(nrow(layout), ncol(layout))))
# Make each plot, in the correct location
for (i in 1:numPlots) {
# Get the i,j matrix positions of the regions that contain this subplot
matchidx <- as.data.frame(which(layout == i, arr.ind = TRUE))
print(plots[[i]], vp = grid::viewport(layout.pos.row = matchidx$row,
layout.pos.col = matchidx$col))
}
}
}
findRegime = function(prop){
RegimeIndices = 1:length(prop)
MostLikelyRegime = RegimeIndices[prop==max(prop)]
return(MostLikelyRegime[1])
#If the probabilities are equal, return the first regime
}
Mode <- function(y) {
uy <- unique(y)
uy[which.max(tabulate(match(y, uy)))]
}
#Old plot function for dynrCook object
#setMethod("plot", "dynrCook",
# function(x, y=NULL, data.dynr, graphingPar=par(no.readonly = TRUE), ylab = #"Smoothed state values", xlab = "Time", numSubjDemo=2, legend.cex=1.2){
# opar = par(no.readonly = TRUE)
# par(graphingPar)
# thesmooth = data.frame(t(res@eta_smooth_final))
#
# dim_latent_var=dim(res@eta_smooth_final)[1]
# num_regime=dim(res@pr_t_given_T)[1]
#
# colnames(thesmooth) = paste0("state",1:dim_latent_var)
#
# ID = data.dynr[["id"]]
# rowIndex = 1:length(ID)
# uniID <- sort(unique(ID))
# thes = sort(sample(1:length(uniID), numSubjDemo))
#
# par(mfrow=c(ifelse(numSubjDemo%%2 > 0,
# numSubjDemo,numSubjDemo/2),
# ifelse(numSubjDemo%%2 > 0,
# 1,2)))
# if (num_regime > 1){
# thePr = t(x@pr_t_given_T)
# mostLikelyRegime = apply(thePr,1,findRegime)
# theR = Mode(mostLikelyRegime)
# }
# for (s in 1:numSubjDemo){
# T <- length(ID[ID %in% uniID[thes[s]]])
# therow = rowIndex[ID %in% uniID[thes[s]]]
# plot(
# c(1,T),
# c(min(thesmooth)-1, max(thesmooth)+quantile(unlist(thesmooth),.1)),
# ylab=ifelse(exists("ylab"),ylab,"State values"),
# xlab=ifelse(exists("xlab"),xlab,"Time"),
# main=ifelse(exists("main"),main,""),
# type='n')
# if (num_regime > 1){
# times=1:T
# thepri = mostLikelyRegime[therow]
# rect(times[thepri==theR]-.5,quantile(unlist(thesmooth),.01),times[th#epri==theR]+.5,quantile(unlist(thesmooth),.99),col="yellow",density=30)
# #legend('topleft',paste0("Regime",theR),
# #bty="n",cex=1.4,col=c(NA),fill=c("yellow"),
# #density=c(100))
# }
# time2 = if(T>500){
# seq(1,T,9)
# }else {1:T}
#
# for (j in 1:dim_latent_var){
# lines(1:T,thesmooth[therow,j], lty=1, lwd=1, col=j)
# points(time2, thesmooth[therow,j][time2], pch=as.character(j), col=j#)
# }
# stateNames = NULL
# for(j in 1:dim_latent_var){
# stateNames = c(stateNames, paste0('State ',j))
# }
# legend('topright',
# paste0("State",as.character(1:dim_latent_var)),
# lty=1:dim_latent_var,
# lwd=2, col=1:dim_latent_var, bty="n",
# pch=as.character(1:dim_latent_var), cex=legend.cex)
# }
# on.exit(par(opar))
# }
#)
##' Plot method for dynrCook objects
##'
##' @param x dynrCook object
##' @param dynrModel model object
##' @param style The style of the plot in the first panel. If style is 1 (default), user-selected smoothed state variables are plotted. If style is 2, user-selected observed-versus-predicted values are plotted.
##' @param names.state (optional) The names of the states to be plotted, which should be a subset of the state.names slot of the measurement slot of dynrModel.
##' @param names.observed (optional) The names of the observed variables to be plotted, which should be a subset of the obs.names slot of the measurement slot of dynrModel.
##' @param printDyn A logical value indicating whether or not to plot the formulas for the dynamic model
##' @param printMeas A logical value indicating whether or not to plot the formulas for the measurement model
##' @param textsize numeric. Font size used in the plot.
##' @param ... Further named arguments
##'
##' @return ggplot object.
##'
##' @details
##' This is a wrapper around \code{\link{dynr.ggplot}}. A great benefit of it is that it shows the model equations in a plot.
plot.dynrCook <- function(x, dynrModel, style = 1, names.state, names.observed, printDyn=TRUE, printMeas=TRUE, textsize=4, ...) {
#The first panel is the ggplot
if(missing(names.observed)){names.observed=dynrModel@measurement@obs.names}
if(missing(names.state)){names.state=dynrModel@measurement@state.names}
p1 <- dynr.ggplot(x, dynrModel, style, numSubjDemo=1, names.state=names.state, names.observed=names.observed, ...)
#If there are more than 2 regimes, the second panel shows the histogram of the most probable regimes across time and subjects.
if (dynrModel$num_regime>1){
regime <- NULL
highProbR <- data.frame(regime=apply(x@pr_t_given_T,2,which.max))
p2 <- ggplot2::ggplot(data = highProbR, ggplot2::aes(factor(regime))) +
ggplot2::geom_bar() +
ggplot2::scale_fill_brewer(palette = 3) +
ggplot2::xlab('Regime') + ggplot2::ylab('Counts') + ggplot2::labs(fill = '') +
ggplot2::ggtitle('Counts by most probable regime') +
ggplot2::theme(axis.title = ggplot2::element_text(size=14),
axis.text=ggplot2::element_text(size=12),
plot.title = ggplot2::element_text(size = 12, colour = "black", face = "bold", hjust = 0.5, vjust=0.01))
}
#The third panel plots the model formulae
p3 <- plotFormula(dynrModel, ParameterAs=sprintf("%.2f", x@transformed.parameters), printDyn=printDyn, printMeas=printMeas, textsize=textsize)
#Organize the panels using the multiplot function
if (dynrModel$num_regime > 1){
multiplot(p1, p2, p3, cols = 1, layout=matrix(c(1,1,2,3), nrow=2, byrow=TRUE))
}else{
multiplot(p1, p3, layout=matrix(c(1,1,2,2), nrow=2, byrow=TRUE))
}
}
plotdf <- function(vec_tex){
dataframe=data.frame(text=sapply(paste0("$",vec_tex,"$"),function(x){as.character(latex2exp::TeX(x))}))
dataframe$x<-0
return(dataframe)
}
##' Plot of the estimated frequencies of the regimes across all individuals and time points
##' based on their smoothed regime probabilities
##'
##' @param res The dynr object returned by dynr.cook().
##' @param dynrModel The model object to plot.
##' @param names.regime (optional) Names of the regimes (must match the length of the number of regimes)
##' @param title (optional) Title of the plot.
##' @param xlab (optional) Label of the x-axis.
##' @param ylab (optional) Label of the y-axis.
##' @param textsize (default = 12) Text size for the axis labels and title (= textsize + 2).
##' @param print (default = TRUE) A flag for whether the plot should be printed.
##'
##' @return ggplot object
dynr.plotFreq <- function(res, dynrModel, names.regime, title, xlab, ylab, textsize=12,print=TRUE) {
if(missing(names.regime)){names.regime<-paste0("Regime",1:dynrModel@num_regime)}
if(missing(title)){title<-'Counts by most probable regime'}
if(missing(xlab)){xlab<-'Regime'}
if(missing(ylab)){ylab<-'Counts'}
regime <- NULL
highProbR <- data.frame(regime=factor(apply(res@pr_t_given_T,2,which.max),
levels=1:dynrModel@num_regime,
labels=names.regime))
p2 <- ggplot2::ggplot(data = highProbR, ggplot2::aes(factor(regime))) +
ggplot2::geom_bar() +
ggplot2::scale_fill_brewer(palette = 3) +
ggplot2::xlab(xlab) + ggplot2::ylab(ylab) + ggplot2::labs(fill = '') +
ggplot2::ggtitle(title) +
ggplot2::theme(axis.title = ggplot2::element_text(size=textsize+2),
axis.text=ggplot2::element_text(size=textsize),
plot.title = ggplot2::element_text(size = textsize, colour = "black", face = "bold", hjust = 0.5, vjust=0.01))
if (print) print(p2)
return(p2)
}
##' Plot the formula from a model
##'
##' @param dynrModel The model object to plot.
##' @param ParameterAs The parameter values or names to plot. The underscores in parameter names are saved for
##' use of subscripts. Greek letters can be specified as corresponding LaTeX symbols without backslashes (e.g., "lambda")
##' and printed as greek letters.
##' @param printDyn A logical value indicating whether or not to plot the formulas for the dynamic model.
##' @param printMeas A logical value indicating whether or not to plot the formulas for the measurement model
##' @param printRS logical. Whether or not to print the regime-switching model. The default is FALSE.
##' @param textsize The text size use in the plot.
##'
##' @details
##' This function typesets a set of formulas that represent the model. Typical inputs to the \code{ParameterAs} argument are (1) the starting values for a model, (2) the final estimated values for a model, and (3) the parameter names. These are accessible with (1) \code{model$xstart}, (2) \code{coef(cook)}, and (3) \code{model$param.names} or \code{names(coef(cook))}, respectively.
##'
##' @return ggplot object
plotFormula <- function(dynrModel, ParameterAs, printDyn=TRUE, printMeas=TRUE,
printRS=FALSE, textsize=4){
dynrModel <- PopBackModel(dynrModel, LaTeXnames(ParameterAs, latex = FALSE))
#Dynamic model
if (printDyn) {
if ("dynrDynamicsMatrix" %in% class(dynrModel$dynamics)){
state.names <- (dynrModel$measurement)$state.names
exo.names <- (dynrModel$dynamics)$covariates
dynrModel@dynamics@values.dyn <- lapply((dynrModel$dynamics)$values.dyn, preProcessNames,state.names,state.names)
dynrModel@dynamics@values.exo <- lapply((dynrModel$dynamics)$values.exo, preProcessNames,state.names,exo.names)
dynrModel@dynamics@values.int <- lapply((dynrModel$dynamics)$values.int, preProcessNames,state.names)
}
dyn.df <- data.frame(text="bold('Dynamic Model')",x=0)
nRegime=ifelse(class(dynrModel@dynamics)=="dynrDynamicsFormula",
length(dynrModel@dynamics@formula),
length(dynrModel@dynamics@values.dyn))
dyn_tex=printex(dynrModel@dynamics,AsMatrix=FALSE)
nEq <- length(dyn_tex[[1]])
# noise_tex <- paste0(" + ", ifelse(dynrModel@dynamics@isContinuousTime, "d", ""), "w_{",1:nEq, "}(t)")
noise_tex <- paste0(ifelse(dynrModel@dynamics@isContinuousTime, "d", ""),
"w_{", 1:nEq, "}(t)")
for (i in 1:nRegime){
if (nRegime>1){
dyn.df <- rbind(dyn.df,data.frame(text=paste0("'Regime ",i,":'"),x=0))
}
if (length(dynrModel@noise@values.latent)==1){
values.latent.mat <- dynrModel@noise@values.latent[[1]]
}else if (length(dynrModel@noise@values.latent)==nRegime){
values.latent.mat <- dynrModel@noise@values.latent[[i]]
}else{stop("The number of regimes implied by the dynamic noise structure does not match the number of regimes in the dynamic model.")}
pos.zero <- diag(values.latent.mat) == 0
# noise_tex[diag(values.latent.mat)==0] <- ""
# dyn.df <- rbind(dyn.df,plotdf(LaTeXnames(paste0(dyn_tex[[i]], noise_tex))))
dyn.df <- rbind(dyn.df,
plotdf(LaTeXnames(
paste0(dyn_tex[[i]],
ifelse(pos.zero, "",
paste0(" + ", noise_tex))))),
plotdf(LaTeXnames(
paste0(ifelse(pos.zero, "",
paste0(noise_tex, " \\sim ", "N(0,\\,",
diag(values.latent.mat), ")"))))) )
}
} else {
dyn.df <- plotdf("")
}
#Measurement model
if (printMeas) {
meas.df <- data.frame(text="bold('Measurement Model')", x=0)
nRegime <- length(dynrModel$measurement$values.load)
meas_tex <- printex(dynrModel$measurement, AsMatrix=FALSE)
nEq <- length(meas_tex[[1]])
# noise_tex <- paste0(" + ", "epsilon_{", 1:nEq, "}")
noise_tex <- paste0("epsilon_{", 1:nEq, "}")
for (i in 1:nRegime){
if (nRegime > 1){
meas.df <- rbind(meas.df, data.frame(text=paste0("'Regime ",i,":'"), x=0))
}
if (length(dynrModel$noise$values.observed) == 1){
values.observed.mat <- dynrModel$noise$values.observed[[1]]
}else if (length(dynrModel$noise$values.observed) == nRegime){
values.observed.mat <- dynrModel$noise$values.observed[[i]]
}else{stop("The number of regimes implied by the measurement noise structure does not match the number of regimes in the measurement model.")}
pos.zero <- diag(values.observed.mat) == 0
# noise_tex[diag(values.observed.mat)==0] <- ""
# meas.df<-rbind(meas.df,plotdf(LaTeXnames(paste0(meas_tex[[i]], noise_tex))))
meas.df <- rbind(meas.df,
plotdf(LaTeXnames(
paste0(meas_tex[[i]],
ifelse(pos.zero, "", paste0(" + ", noise_tex, ",")),
"\\;\\,",
ifelse(pos.zero,
"",
paste0(noise_tex, " \\sim ", "N(0,\\,", diag(values.observed.mat), ")"))))))
}
} else {
meas.df <- plotdf("")
}
# Regime-switching model
if (printRS) {
rs.df <- data.frame(text="bold('Regime-switching Model')", x=0)
if (any(dim(dynrModel@regimes@params) == 0)) {
warning("No Regime-switching Model. Please turn off 'printRS'.")
} else {
# helper function to make sub. e.g., a_111 to a_{111}
to_sub <- function(parnames) {
idx_sub <- grepl('_', parnames)
if ( sum(idx_sub)==0 ) { return(parnames) }
subs <- parnames[idx_sub]
subs_spl <- strsplit(subs, '_')
subbed <- sapply(subs_spl, function(sub) {
paste0(sub[1], '_', paste0('{', sub[2], '}') )
})
parnames[idx_sub] <- subbed
parnames
}
covars <- dynrModel@regimes@covariates
if ( length(covars) == 0 ) {
covars <- c("")
} else {
covars <- c("", paste0("*", covars) )
}
covars_n <- length(covars)
# TODO: check nregime to fit with others
nregime <- nrow(dynrModel@regimes@params)
# column index for the intercepts of regimes
inter_col <- cumsum(c(1, rep(covars_n, nregime-2)))
# i:(i+ncov-1)
parnames <- to_sub( dynrModel@regimes@paramnames )
param_rs <- matrix(parnames, nregime)
param_inter <- param_rs[, inter_col, drop=FALSE]# intercepts
refrow <- dynrModel@regimes@refRow
if (length(refrow) != 0) {
for ( co in inter_col ) {
param_rs[-refrow, co] <-
paste0(param_rs[refrow, co], ' + ', param_rs[-refrow, co])
}
}
param_rs <- matrix(sweep(param_rs, 2, covars, FUN=paste0),
nrow=nregime)
ll <- vector('list', length=length(inter_col))
for ( i in 1:length(inter_col) ) {
ll[[i]] <-
param_rs[, inter_col[i]:(inter_col[i]+covars_n-1),
drop=FALSE]
}
sl <- sapply(ll, function(l) {
apply(l, 1, paste, collapse=" + ")
})
# transition probability subs
transub <- array(c(row(sl), col(sl)), c(nrow(sl), ncol(sl), 2))
transubed <- apply(transub, c(1,2), paste, collapse="")
rs.df <- rbind(rs.df,
plotdf(paste0("Log-odds(p_{",
transubed, "}) = ",
sl)))
}
} else {
rs.df <- plotdf("")
}
plot.df <- rbind(meas.df, dyn.df, rs.df)
plot.df$y=seq((dim(plot.df)[1]-1)*5+1, 1, by=-5)
x <- NULL
y <- NULL
fig<-ggplot2::ggplot(plot.df, ggplot2::aes(x=x, y=y, label=text))+
ggplot2::geom_text(parse=TRUE,size=textsize)+
ggplot2::theme_void() + ggplot2::ylim(min(plot.df$y)-5, max(plot.df$y)+5)
return(fig)
}
##' The ggplot of the smoothed state estimates and the most likely regimes
##'
##' @aliases autoplot.dynrCook
##'
##' @param res The dynr object returned by \code{dynr.cook()}.
##' @param dynrModel The model object to plot.
##' @param style The style of the plot. If style is 1 (default), user-selected smoothed state variables are plotted. If style is 2, user-selected observed-versus-predicted values are plotted.
##' @param numSubjDemo The number of subjects to be randomly selected for plotting.
##' @param idtoPlot Values of the ID variable to plot.
##' @param names.state (optional) The names of the states to be plotted, which should be a subset of the state.names slot of the measurement slot of dynrModel.
##' @param names.observed (optional) The names of the observed variables to be plotted, which should be a subset of the obs.names slot of the measurement slot of dynrModel.
##' @param names.regime (optional) The names of the regimes to be plotted, which can be missing.
##' @param shape.values (optional) A vector of values that correspond to the shapes of the points, which can be missing. See the R documentation on pch for details on possible shapes.
##' @param title (optional) A title of the plot.
##' @param ylab (optional) The label of the y axis.
##' @param is.bw Is plot in black and white? The default is FALSE.
##' @param colorPalette A color palette for lines and dots. It is a value passed to the palette argument of the \code{ggplot2::scale_colour_brewer()} function. These palettes are in the R package \pkg{RColorBrewer}. One can find them by attaching the package with \code{library(RColorBrewer)} and run \code{display.brewer.all()}.
##' @param fillPalette A color palette for blocks. It is a value passed to the palette argument of the \code{ggplot2::scale_fill_brewer()} function. These palettes are in the package \pkg{RColorBrewer}. One can find them by attaching the package with \code{library(RColorBrewer)} and run \code{display.brewer.all()}.
##' @param mancolorPalette (optional) A color palette for manually scaling the colors of lines and dots. It is a vector passed to the values argument of the \code{ggplot2::scale_colour_manual} function.
##' @param manfillPalette (optional) A color palette for manually scaling the colors of filled blocks. It is a vector passed to the values argument of the \code{ggplot2::scale_fill_manual} function.
##' @param ... A list of elements that modify the existing ggplot theme. Consult the \code{ggplot2::theme()} function in the R package \pkg{ggplot2} for more options.
##'
##' @details
##' This function outputs a ggplot layer that can be modified using functions in the package \pkg{ggplot2}. That is, one can add layers, scales, coords and facets with the "+" sign. In an example below, the \code{ggplot2::ylim()} function is used to modify the limits of the y axis of the graph. More details can be found on \url{https://ggplot2.tidyverse.org/} and \url{https://ggplot2.tidyverse.org/reference/}.
##'
##' The two functions \code{dynr.ggplot()} and \code{autoplot()} as identical aliases of one another. The \code{autoplot()} function is an S3 method from the package \pkg{ggplot2} that allows many objects to be plotted and works like the base \code{plot()} function.
##'
##' @return ggplot object
##'
##' @examples
##' # The following code is part of a demo example in dynr
##' \dontrun{
##' demo(RSLinearDiscreteYang, package='dynr')
##' p <- dynr.ggplot(yum, dynrModel = rsmod, style = 1,
##' names.regime = c("Deactivated", "Activated"),
##' title = "(B) Results from RS-AR model", numSubjDemo = 1,
##' shape.values = c(1),
##' text = element_text(size = 16),
##' is.bw = TRUE)
##' # One can modify the limits on the y axis by using '+'
##' p + ggplot2::ylim(-2, 4)
##'
##' autoplot(yum, dynrModel = rsmod, style = 1,
##' names.regime = c("Deactivated", "Activated"),
##' title = "(B) Results from RS-AR model", numSubjDemo = 1,
##' shape.values = c(1),
##' text = element_text(size = 16),
##' is.bw = TRUE)
##' }
dynr.ggplot <- function(res, dynrModel, style = 1,
numSubjDemo=2, idtoPlot=c(),
names.state,
names.observed,
names.regime,
shape.values,
title,
ylab,
is.bw=FALSE,
colorPalette="Set2",
fillPalette="Set2",
mancolorPalette, manfillPalette, ...){
data.dynr=dynrModel@data
dim_latent_var=dim(res@eta_smooth_final)[1]
num_regime=dim(res@pr_t_given_T)[1]
if(missing(names.observed)){names.observed=dynrModel@measurement@obs.names}
observed=which(dynrModel@measurement@obs.names%in%names.observed)
if ((length(observed)>8)&(missing(mancolorPalette))&(style==2)){stop("You provided too many variables than the default color palette can handle.\nPlease consider specify the mancolorPalette argument.")}
if(missing(names.state)){names.state=dynrModel@measurement@state.names}
states=which(dynrModel@measurement@state.names%in%names.state)
if(missing(shape.values)){
if (style==1) shape.values=states#48+states
if (style==2) shape.values=c(rep(32,length(observed)),observed)#pre-obs
}else{
if (style==2) shape.values=c(rep(32,length(shape.values)), shape.values)#pre-obs
}
if (style==1) line.values=rep(1,length(states))
if (style==2) line.values=rep(c(1,0),each=length(observed))#pre-obs
num_sbj=length(unique(data.dynr$id))
if(length(idtoPlot)<1){
randid =sample(unique(data.dynr$id),numSubjDemo)
}else {randid=idtoPlot}
if(missing(title)){
if (style==1) title=" "
if (style==2) title="Observed vs Predicted"
}
if(missing(ylab)){
if (style==1) ylab="Smoothed State Values"
if (style==2) ylab="Values"
}
data.plot<-data.frame(id=as.factor(data.dynr$id),time=data.dynr$time)
if(num_regime==1){
names.id.vars <- c("id","time")
if (style==1){
names.measure.vars <- names.state
data.plot <- cbind(data.plot, matrix(res@eta_smooth_final[states,],ncol=length(states),byrow=TRUE))
lines.var <- names.state
points.var <- names.state
}#end of style 1
if (style==2){
dynrModel=PopBackModel(dynrModel, LaTeXnames(res@transformed.parameters, latex = FALSE))
dynrModel@measurement@values.load=lapply(dynrModel@measurement@values.load, function(x){matrix(as.numeric(x),ncol=ncol(x))})
num_regime_meas <- length(dynrModel@measurement@values.load)
if (length(dynrModel@measurement@values.exo)!=0){
dynrModel@measurement@values.exo=lapply(dynrModel@measurement@values.exo, function(x){matrix(as.numeric(x),ncol=ncol(x))})
}
if (length(dynrModel@measurement@values.int)!=0){
dynrModel@measurement@values.int=lapply(dynrModel@measurement@values.int, function(x){matrix(as.numeric(x),ncol=ncol(x))})
}
#predicted scores
if (num_regime_meas==1){#one-regime measurement model
predicted=dynrModel@measurement@values.load[[1]]%*%res@eta_smooth_final
if (length(dynrModel@measurement@values.int)!=0){
predicted=predicted+dynrModel@measurement@values.int[[1]]%*%matrix(rep(1,ncol(res@eta_smooth_final)), nrow=1)
}
if (length(dynrModel@measurement@values.exo)!=0){
exo.mat=t(as.matrix(dynrModel@data$covariates[,paste0("covar",sapply(dynrModel@measurement@exo.names, function(x){which(dynrModel@data$covariate.names==x)}))]))
predicted=predicted+dynrModel@measurement@values.exo[[1]]%*%exo.mat
}
}#end of one-regime meas
names.measure.vars <- c(paste0(names.observed, ".predicted"), paste0(names.observed, ".observed"))
data.plot <- cbind(data.plot, matrix(predicted[observed,],ncol=length(observed),byrow=TRUE, dimnames=list(NULL, paste0(names.observed, ".predicted"))))
data.plot <- cbind(data.plot, data.dynr$observed[,observed])
lines.var <- paste0(names.observed, ".predicted")
points.var <- paste0(names.observed, ".observed")
if (length(observed)<8){
fancy8 <- c("#edaac1","#a7dab0","#dfb3e2","#8bc5ed","#bcb8ec","#e6b296","#d2d39d","#82d9d6")
mancolorPalette=rep(fancy8[observed],2)
}
}#end of style 2
names(data.plot) <- c(names.id.vars, names.measure.vars)
data_long <- reshape2::melt(data.plot[data.plot$id%in%randid,],
id.vars=names.id.vars,
measure.vars=names.measure.vars,
value.name="value")
data_long$variable<-factor(data_long$variable, levels=sort(names.measure.vars))
value <- NULL
variable <- NULL
partial.plot<-ggplot2::ggplot(data_long, ggplot2::aes(x=time, y=value, colour = variable, shape = variable, linetype = variable)) +
ggplot2::geom_line(data=data_long[data_long$variable%in%lines.var,], size=1) +
ggplot2::geom_point(data=data_long[data_long$variable%in%points.var,], size=4) +
ggplot2::scale_linetype_manual(labels=sort(names.measure.vars), values=line.values[order(names.measure.vars)])+
ggplot2::scale_shape_manual(labels=sort(names.measure.vars), values=shape.values[order(names.measure.vars)])+
ggplot2::facet_wrap(~id)+
ggplot2::labs(title = title, y=ylab, ggtitle="")+
ggplot2::theme(plot.title=ggplot2::element_text(lineheight=.8, face="bold", hjust = 0.5, vjust=0.01), legend.position="bottom",legend.text=ggplot2::element_text(lineheight=0.8, face="bold"),...)
}else{#more than two regimes
if(missing(names.regime)){names.regime = 1:num_regime}
findRegime = function(prop){
MostLikelyRegime = names.regime[prop==max(prop)]
return(MostLikelyRegime[1])
#If the probabilities are equal, return the first regime
}
addendtime <- function(data_frame){
data_frame$endtime <- c(data_frame$time[2:nrow(data_frame)], data_frame$time[nrow(data_frame)]+min(diff(data_frame$time)))
return(data_frame)
}
data.plot <- plyr::ddply(data.plot,"id",addendtime)
data.plot$regime<-as.factor(apply(res@pr_t_given_T,2,findRegime))
names.id.vars<-c("id","time","endtime","regime")
if (style==1){
names.measure.vars<-names.state
data.plot<-cbind(data.plot,matrix((res@eta_smooth_final[states,]),ncol=length(states),byrow=TRUE))
lines.var <- names.state
points.var <- names.state
}#end of style 1
if (style==2){
dynrModel=PopBackModel(dynrModel, LaTeXnames(res@transformed.parameters, latex = FALSE))
dynrModel@measurement@values.load=lapply(dynrModel@measurement@values.load, function(x){matrix(as.numeric(x),ncol=ncol(x))})
num_regime_meas <- length(dynrModel@measurement@values.load)
if (length(dynrModel@measurement@values.exo)!=0){
dynrModel@measurement@values.exo=lapply(dynrModel@measurement@values.exo, function(x){matrix(as.numeric(x),ncol=ncol(x))})
}
if (length(dynrModel@measurement@values.int)!=0){
dynrModel@measurement@values.int=lapply(dynrModel@measurement@values.int, function(x){matrix(as.numeric(x),ncol=ncol(x))})
}
#predicted scores
if (num_regime_meas==1){#one-regime measurement model
predicted=matrix(dynrModel@measurement@values.load[[1]][observed,],nrow=length(observed))%*%res@eta_smooth_final
if (length(dynrModel@measurement@values.int)!=0){
predicted=predicted+matrix(dynrModel@measurement@values.int[[1]][observed,],nrow=length(observed))%*%matrix(rep(1,ncol(res@eta_smooth_final)), nrow=1)
}
if (length(dynrModel@measurement@values.exo)!=0){
exo.mat=t(as.matrix(dynrModel@data$covariates[,paste0("covar",sapply(dynrModel@measurement@exo.names, function(x){which(dynrModel@data$covariate.names==x)}))]))
predicted=predicted+matrix(dynrModel@measurement@values.exo[[1]][observed,],nrow=length(observed))%*%exo.mat
}
#end of one-regime meas
}else if (num_regime_meas==num_regime){
ntotaltime=ncol(res@eta_smooth_final)
predicted=matrix(rep(0,length(observed)*ntotaltime),ncol=ntotaltime)
for (index in 1:ntotaltime){
predicted[,index]=matrix(dynrModel@measurement@values.load[[data.plot$regime[index]]][observed,],nrow=length(observed))%*%matrix(res@eta_smooth_final[,index], ncol=1)
if (length(dynrModel@measurement@values.int)!=0){
predicted[,index]=matrix(predicted[,index]+dynrModel@measurement@values.int[[data.plot$regime[index]]][observed,],nrow=length(observed))
}
if (length(dynrModel@measurement@values.exo)!=0){
exo.mat=t(as.matrix(dynrModel@data$covariates[index,paste0("covar",sapply(dynrModel@measurement@exo.names, function(x){which(dynrModel@data$covariate.names==x)}))]))
predicted[,index]=matrix(predicted[,index]+dynrModel@measurement@values.exo[[data.plot$regime[index]]][observed,],nrow=length(observed))%*%exo.mat
}
}
#end of regime specific measurement
}
names.measure.vars <- c(paste0(names.observed, ".predicted"), paste0(names.observed, ".observed"))
data.plot <- cbind(data.plot, matrix(predicted,ncol=length(observed), byrow=TRUE, dimnames=list(NULL, paste0(names.observed, ".predicted"))))
data.plot <- cbind(data.plot, data.dynr$observed[,observed])
lines.var <- paste0(names.observed, ".predicted")
points.var <- paste0(names.observed, ".observed")
if (length(observed)<=8){
fancy8 <- c("#edaac1","#a7dab0","#dfb3e2","#8bc5ed","#bcb8ec","#e6b296","#d2d39d","#82d9d6")
manfillPalette=fancy8[9-(1:num_regime)]
mancolorPalette=rep(fancy8[observed],2)
}
}#end of style 2
names(data.plot)<-c(names.id.vars,names.measure.vars)
data_long<- reshape2::melt(data.plot[data.plot$id%in%randid,],
id.vars=names.id.vars,
measure.vars=names.measure.vars,
value.name="value")
data_long$variable<-factor(data_long$variable, levels=sort(names.measure.vars))
#data_long$statenumber<-as.factor(sub("state","",data_long$variable))
endtime <- NULL
regime <- NULL
partial.plot<-ggplot2::ggplot(data_long,ggplot2::aes(x=time, y=value, group=variable)) +
ggplot2::geom_rect(ggplot2::aes(xmin=time, xmax=endtime, ymin=-Inf, ymax=Inf, fill=regime), alpha=.15) +
ggplot2::geom_line(data=data_long[data_long$variable%in%lines.var,], size=1, ggplot2::aes(color=variable, linetype=variable)) +
ggplot2::geom_point(data=data_long[data_long$variable%in%points.var,], size=3, ggplot2::aes(color=variable, shape=variable)) +
ggplot2::scale_linetype_manual(labels=sort(names.measure.vars), values=line.values[order(names.measure.vars)])+
ggplot2::scale_shape_manual(labels=sort(names.measure.vars), values=shape.values[order(names.measure.vars)])+
#geom_text(size=1, ggplot2::aes(label=statenumber,color=variable))+
ggplot2::facet_wrap(~id)+
ggplot2::labs(title = title,y=ylab,ggtitle="")+
ggplot2::theme(plot.title=ggplot2::element_text(lineheight=.8, face="bold", hjust = 0.5, vjust=0.01), legend.position="bottom",legend.text=ggplot2::element_text(lineheight=0.8, face="bold"),...)
}
if(is.bw){
bw.plot<-partial.plot+
ggplot2::scale_fill_grey(start = 0.1, end = 0.9)+
ggplot2::scale_color_grey(labels=sort(names.measure.vars), start = 0.2, end = 0.7)
return(bw.plot)
}else{
default.plot<-partial.plot+
ggplot2::scale_colour_brewer(labels=sort(names.measure.vars), palette=colorPalette)+
ggplot2::scale_fill_brewer(palette=fillPalette)
if(!missing(mancolorPalette)){
manual.plot<-partial.plot+ggplot2::scale_colour_manual(labels=sort(names.measure.vars), values=mancolorPalette[order(names.measure.vars)])
}
if(!missing(manfillPalette)){
if(!exists("manual.plot")){manual.plot<-partial.plot}
manual.plot<-manual.plot+ggplot2::scale_fill_manual(values=manfillPalette)
}
if (exists("manual.plot")){
return(manual.plot) #print(manual.plot)
}else{
return(default.plot)#print(default.plot)
}
}
}
##' @rdname dynr.ggplot
##'
##' @param object The same as res. The dynr object returned by dynr.cook().
##'
##' @return ggplot object
autoplot.dynrCook <- function(object, dynrModel, style = 1,
numSubjDemo=2, idtoPlot=c(),
names.state,
names.observed,
names.regime,
shape.values,
title,
ylab,
is.bw=FALSE,
colorPalette="Set2",
fillPalette="Set2",
mancolorPalette, manfillPalette, ...){
dynr.ggplot(object, dynrModel=dynrModel, style = style,
numSubjDemo=numSubjDemo, idtoPlot=idtoPlot,
names.state=names.state,
names.observed=names.observed,
names.regime=names.regime,
shape.values=shape.values,
title=title,
ylab=ylab,
is.bw=is.bw,
colorPalette=colorPalette,
fillPalette=fillPalette,
mancolorPalette=mancolorPalette, manfillPalette=manfillPalette, ...)
}
##' The ggplot of the outliers estimates.
##'
##' @param object A dynrTaste object.
##' @param numSubjDemo The number of subjects, who have
##' largest joint chi-square statistic, to be selected for plotting.
##' @param idtoPlot Values of the ID variable to plot.
##' @param names.state (optional) The names of the states to be plotted, which should be a subset of the state.names slot of the measurement slot of dynrModel. If NULL, the t statistic plots for all state variables will be included.
##' @param names.observed (optional) The names of the observed variables to be plotted, which should be a subset of the obs.names slot of the measurement slot of dynrModel. If NULL, the t statistic plots for all observed variables will be included.
##' @param ... Place holder for other arguments. Please do not use.
##'
##' @return a list of ggplot objects for each ID.
##' The plots of chi-square statistics (joint and independent),
##' and the plots of t statistic for \code{names.state} and \code{names.observed} will be included.
##' Users can modify the ggplot objects using ggplot grammar.
##' If a \code{filename} is provided, a pdf of plots will be saved additionally.
autoplot.dynrTaste <- function(object,
numSubjDemo=2, idtoPlot=NULL,
names.state=NULL, names.observed=NULL, ...) {
if ( !inherits(object, "dynrTaste") ) {
stop("dynrTaste object is required.")
}
stopifnot(numSubjDemo >= 1)
numSubjDemo <- as.integer(numSubjDemo)
conf_level <- object$conf.level
tstart <- object$tstart
chi_jnt <- object$chi.jnt
id <- object$id
id_unq <- unique(id)
id_n <- length(id_unq)
lat_name <- rownames(object$t.inn)
obs_name <- rownames(object$t.add)
lat_n <- length(lat_name)
obs_n <- length(obs_name)
if( is.null(names.state) ) {
lat_toplot <- lat_name
} else {
if ( !all(names.state %in% lat_name) ) {
stop("'names.state' should be a subset of the latent variables.")
}
lat_toplot <- names.state
}
if( is.null(names.observed) ) {
obs_toplot <- obs_name
} else {
if ( !all(names.observed %in% obs_name) ) {
stop("'names.observed' should be a subset of the observed variables.")
}
obs_toplot <- names.observed
}
if ( is.null(idtoPlot) ) {
chi_max <- vector("numeric", id_n)
for(i in 1:id_n){
begT <- tstart[i] + 1
endT <- tstart[i+1]
chi_max[i] <- max( chi_jnt[begT:endT] )
}
id_toplot <- id_unq[ order(chi_max, decreasing=TRUE) ][
1:numSubjDemo ]
} else {
id_toplot <- idtoPlot
if ( !all(id_toplot %in% id_unq) ) {
stop("Not all ID are in the data.")
}
}
### create ggplot objects
ggadd_hline <- function(yintercept) {
geom_hline(yintercept=yintercept,
colour="black", linetype="dashed", alpha=0.5)
}
ggadd_point <- function(data, aes_y) {
geom_point(data=data, aes_string(x="time", y=aes_y),
colour="red", size=2, alpha=0.5)
}
ggadd_theme <- function() {
theme(plot.title=element_text(size = rel(0.85)),
axis.title.x=element_blank(),
#axis.text.x=element_blank(),
axis.ticks.x=element_blank(),
panel.grid.major=element_blank(),
panel.grid.minor=element_blank(),
panel.background=element_blank(),
plot.margin = margin(0, 0, 0, 0),
axis.line=element_line(colour = "black"))
}
## Joint chi-square
id_toplot_n <- length(id_toplot)
time <- object$time
chi_jnt_shk <- NULL
chi_df_jnt <- data.frame(id = id, time = time,
chi_jnt = object$chi.jnt,
chi_jnt_shk = object$chi.jnt.shock)
# subset data for id to plot
chi_df_jnt_plot <- subset(chi_df_jnt, is.element(id, id_toplot))
# subset data only for shocks
chi_df_jnt_shk <- subset(chi_df_jnt_plot, chi_jnt_shk)
qchi_jnt <- qchisq(conf_level, lat_n + obs_n)
plots_jnt <- lapply(id_toplot, function(id_i) {
chi_df_jnt_shk_i <- subset(chi_df_jnt_shk, id==id_i)
ggplot2::ggplot(subset(chi_df_jnt_plot, id==id_i),
aes(x=time, y=chi_jnt)) +
geom_line(alpha=0.5) +
labs(title="") +
ylab(expression(paste("Joint ", Chi^2))) +
list(ggadd_point(chi_df_jnt_shk_i, "chi_jnt"),
ggadd_hline(qchi_jnt), ggadd_theme())
})
#### Independend chi-square for innovative ####
chi_inn <- chi_inn_shk <- NULL
chi_df_inn <- data.frame(id = id, time = time,
chi_inn = object$chi.inn,
chi_inn_shk = object$chi.inn.shock)
chi_df_inn_plot <- subset(chi_df_inn, is.element(id, id_toplot))
chi_df_inn_shk <- subset(chi_df_inn_plot, chi_inn_shk)
qchi_inn <- qchisq(conf_level, lat_n)
plots_inn <- lapply(id_toplot, function(id_i) {
chi_df_inn_shk_i <- subset(chi_df_inn_shk, id==id_i)
ggplot2::ggplot(subset(chi_df_inn_plot, id==id_i),
aes(x=time, y=chi_inn)) +
geom_line(alpha=0.5) +
labs(title="") +
ylab(expression(paste("Indep. ", Chi^2, " - innovative"))) +
list(ggadd_point(chi_df_inn_shk_i, "chi_inn"),
ggadd_hline(qchi_inn), ggadd_theme())
})
#### Independent chi-square for additive ####
chi_add <- chi_add_shk <- NULL
chi_df_add <- data.frame(id = id, time = time,
chi_add = object$chi.add,
chi_add_shk = object$chi.add.shock)
chi_df_add_plot <- subset(chi_df_add, is.element(id, id_toplot))
chi_df_add_shk <- subset(chi_df_add_plot, chi_add_shk)
plots_add <- lapply(id_toplot, function(id_i) {
chi_df_add_shk_i <- subset(chi_df_add_shk, id==id_i)
ggplot2::ggplot(subset(chi_df_add_plot, id==id_i),
aes(x=time, y=chi_add)) +
geom_line(alpha=0.5) +
labs(title="") +
ylab(expression(paste("Indep. ", Chi^2, " - additive"))) +
list(ggadd_point(chi_df_add_shk_i, "chi_add"),
ggadd_hline(qchisq(conf_level, lat_n)), ggadd_theme())
})
#### t-statistic for innovative ####
t_inn <- object$t.inn
t_inn_name <- row.names(t_inn)
row.names(t_inn) <- paste0(t_inn_name, "_t")
t_df_inn <- data.frame(id = id, time = time,
t(t_inn), t(object$t.inn.shock))
t_df_inn_plot <- subset(t_df_inn, is.element(id, id_toplot))
plots_t_inn <- lapply(id_toplot, function(id_i) {
t_df_inn_plot_i <- subset(t_df_inn_plot, id==id_i)
lapply(lat_toplot, function(lat_i) {
t_df_inn_shk_i <- subset(t_df_inn_plot_i,
t_df_inn_plot_i[[lat_i]])
qt_i <- qt(conf_level, nrow(t_df_inn_plot_i) - obs_n)
ggplot2::ggplot(t_df_inn_plot_i,
aes_string(x="time", y=paste0(lat_i, "_t"))) +
geom_line(alpha=0.5) +
labs(title="") +
ylab(paste("t_[", lat_i, "]")) +
list(ggadd_point(t_df_inn_shk_i, paste0(lat_i, "_t")),
ggadd_hline(qt_i), ggadd_hline(-qt_i), ggadd_theme())
})
})
#### t-statistic for additive ####
t_add <- object$t.add
t_add_name <- row.names(t_add)
row.names(t_add) <- paste0(t_add_name, "_t")
t_df_add <- data.frame(id = id, time = time,
t(t_add), t(object$t.add.shock))
t_df_add_plot <- subset(t_df_add, is.element(id, id_toplot))
plots_t_add <- lapply(id_toplot, function(id_i) {
t_df_add_plot_i <- subset(t_df_add_plot, id==id_i)
t_add_df_i <- nrow(t_df_add_plot_i)
lapply(obs_toplot, function(obs_i) {
t_df_add_shk_i <- subset(t_df_add_plot_i,
t_df_add_plot_i[[obs_i]])
qt_i <- qt(conf_level, t_add_df_i - obs_n)
ggplot2::ggplot(t_df_add_plot_i,
aes_string(x="time", y=paste0(obs_i, "_t"))) +
geom_line(alpha=0.5) +
labs(title="") +
ylab(paste("t_[", obs_i, "]")) +
list(ggadd_point(t_df_add_shk_i, paste0(obs_i, "_t")),
ggadd_hline(qt_i), ggadd_hline(-qt_i), ggadd_theme())
})
})
gg_objects <- mapply(function(jnt_i, inn_i, add_i,
t_inn_i, t_add_i, id_toplot_i) {
gg_i <- c(list(jnt_i, inn_i, add_i), t_inn_i, t_add_i)
names(gg_i) <- c("chi_jnt", "chi_inn", "chi_add",
paste0("t_[", lat_toplot, "]"),
paste0("t_[", obs_toplot, "]"))
gg_i
},
plots_jnt, plots_inn, plots_add, plots_t_inn, plots_t_add, id_toplot,
SIMPLIFY=FALSE)
names(gg_objects) <- id_toplot
# if ( !is.null(filename) && is.character(filename) ) {
# gg_plots <- mapply(function(plots_i, id_toplot_i) {
# plots_arr <- ggpubr::ggarrange(
# plotlist=plots_i, nrow=length(plots_i), ncol=1, align="v")
# ggpubr::annotate_figure(
# plots_arr,
# top=ggpubr::text_grob(paste0("ID: ", id_toplot_i),
# color="black", face="bold", size=14),
# bottom=ggpubr::text_grob("time", color="black", size=11))
# },
# gg_objects, id_toplot,
# SIMPLIFY=FALSE)
#
# id_v <- chi_df_jnt_plot[["id"]]
# id_time_length <- vector("numeric", id_toplot_n)
# # time points for each id
# for (i in 1:id_toplot_n) {
# id_time_length[i] <- sum(is.element(id_v, id_toplot[i]))
# }
# # width: time/6, height: n of gg object * 2
# ggpubr::ggexport(gg_plots, filename=filename,
# width=mean(id_time_length)/6,
# height=length(gg_objects[[1]])*2)
# }
invisible(gg_objects)
}
Try the dynr package in your browser
Any scripts or data that you put into this service are public.
dynr documentation built on May 29, 2024, 2:49 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions autoplot.dynrTaste autoplot.dynrCook dynr.ggplot plotFormula dynr.plotFreq plotdf plot.dynrCook Mode findRegime multiplot
Documented in autoplot.dynrCook autoplot.dynrTaste dynr.ggplot dynr.plotFreq plot.dynrCook plotFormula
#------------------------------------------------------------------------------
# Plotting methods
multiplot <- function(..., plotlist=NULL, file, cols=1, layout=NULL) {
# Make a list from the ... arguments and plotlist
plots <- c(list(...), plotlist)
numPlots = length(plots)
# If layout is NULL, then use 'cols' to determine layout
if (is.null(layout)) {
# Make the panel
# ncol: Number of columns of plots
# nrow: Number of rows needed, calculated from # of cols
layout <- matrix(seq(1, cols * ceiling(numPlots/cols)),
ncol = cols, nrow = ceiling(numPlots/cols))
}
if (numPlots==1) {
print(plots[[1]])
} else {
# Set up the page
grid::grid.newpage()
grid::pushViewport(grid::viewport(layout = grid::grid.layout(nrow(layout), ncol(layout))))
# Make each plot, in the correct location
for (i in 1:numPlots) {
# Get the i,j matrix positions of the regions that contain this subplot
matchidx <- as.data.frame(which(layout == i, arr.ind = TRUE))
print(plots[[i]], vp = grid::viewport(layout.pos.row = matchidx$row,
layout.pos.col = matchidx$col))
}
}
}
findRegime = function(prop){
RegimeIndices = 1:length(prop)
MostLikelyRegime = RegimeIndices[prop==max(prop)]
return(MostLikelyRegime[1])
#If the probabilities are equal, return the first regime
}
Mode <- function(y) {
uy <- unique(y)
uy[which.max(tabulate(match(y, uy)))]
}
#Old plot function for dynrCook object
#setMethod("plot", "dynrCook",
# function(x, y=NULL, data.dynr, graphingPar=par(no.readonly = TRUE), ylab = #"Smoothed state values", xlab = "Time", numSubjDemo=2, legend.cex=1.2){
# opar = par(no.readonly = TRUE)
# par(graphingPar)
# thesmooth = data.frame(t(res@eta_smooth_final))
#
# dim_latent_var=dim(res@eta_smooth_final)[1]
# num_regime=dim(res@pr_t_given_T)[1]
#
# colnames(thesmooth) = paste0("state",1:dim_latent_var)
#
# ID = data.dynr[["id"]]
# rowIndex = 1:length(ID)
# uniID <- sort(unique(ID))
# thes = sort(sample(1:length(uniID), numSubjDemo))
#
# par(mfrow=c(ifelse(numSubjDemo%%2 > 0,
# numSubjDemo,numSubjDemo/2),
# ifelse(numSubjDemo%%2 > 0,
# 1,2)))
# if (num_regime > 1){
# thePr = t(x@pr_t_given_T)
# mostLikelyRegime = apply(thePr,1,findRegime)
# theR = Mode(mostLikelyRegime)
# }
# for (s in 1:numSubjDemo){
# T <- length(ID[ID %in% uniID[thes[s]]])
# therow = rowIndex[ID %in% uniID[thes[s]]]
# plot(
# c(1,T),
# c(min(thesmooth)-1, max(thesmooth)+quantile(unlist(thesmooth),.1)),
# ylab=ifelse(exists("ylab"),ylab,"State values"),
# xlab=ifelse(exists("xlab"),xlab,"Time"),
# main=ifelse(exists("main"),main,""),
# type='n')
# if (num_regime > 1){
# times=1:T
# thepri = mostLikelyRegime[therow]
# rect(times[thepri==theR]-.5,quantile(unlist(thesmooth),.01),times[th#epri==theR]+.5,quantile(unlist(thesmooth),.99),col="yellow",density=30)
# #legend('topleft',paste0("Regime",theR),
# #bty="n",cex=1.4,col=c(NA),fill=c("yellow"),
# #density=c(100))
# }
# time2 = if(T>500){
# seq(1,T,9)
# }else {1:T}
#
# for (j in 1:dim_latent_var){
# lines(1:T,thesmooth[therow,j], lty=1, lwd=1, col=j)
# points(time2, thesmooth[therow,j][time2], pch=as.character(j), col=j#)
# }
# stateNames = NULL
# for(j in 1:dim_latent_var){
# stateNames = c(stateNames, paste0('State ',j))
# }
# legend('topright',
# paste0("State",as.character(1:dim_latent_var)),
# lty=1:dim_latent_var,
# lwd=2, col=1:dim_latent_var, bty="n",
# pch=as.character(1:dim_latent_var), cex=legend.cex)
# }
# on.exit(par(opar))
# }
#)
##' Plot method for dynrCook objects
##'
##' @param x dynrCook object
##' @param dynrModel model object
##' @param style The style of the plot in the first panel. If style is 1 (default), user-selected smoothed state variables are plotted. If style is 2, user-selected observed-versus-predicted values are plotted.
##' @param names.state (optional) The names of the states to be plotted, which should be a subset of the state.names slot of the measurement slot of dynrModel.
##' @param names.observed (optional) The names of the observed variables to be plotted, which should be a subset of the obs.names slot of the measurement slot of dynrModel.
##' @param printDyn A logical value indicating whether or not to plot the formulas for the dynamic model
##' @param printMeas A logical value indicating whether or not to plot the formulas for the measurement model
##' @param textsize numeric. Font size used in the plot.
##' @param ... Further named arguments
##'
##' @return ggplot object.
##'
##' @details
##' This is a wrapper around \code{\link{dynr.ggplot}}. A great benefit of it is that it shows the model equations in a plot.
plot.dynrCook <- function(x, dynrModel, style = 1, names.state, names.observed, printDyn=TRUE, printMeas=TRUE, textsize=4, ...) {
#The first panel is the ggplot
if(missing(names.observed)){names.observed=dynrModel@measurement@obs.names}
if(missing(names.state)){names.state=dynrModel@measurement@state.names}
p1 <- dynr.ggplot(x, dynrModel, style, numSubjDemo=1, names.state=names.state, names.observed=names.observed, ...)
#If there are more than 2 regimes, the second panel shows the histogram of the most probable regimes across time and subjects.
if (dynrModel$num_regime>1){
regime <- NULL
highProbR <- data.frame(regime=apply(x@pr_t_given_T,2,which.max))
p2 <- ggplot2::ggplot(data = highProbR, ggplot2::aes(factor(regime))) +
ggplot2::geom_bar() +
ggplot2::scale_fill_brewer(palette = 3) +
ggplot2::xlab('Regime') + ggplot2::ylab('Counts') + ggplot2::labs(fill = '') +
ggplot2::ggtitle('Counts by most probable regime') +
ggplot2::theme(axis.title = ggplot2::element_text(size=14),
axis.text=ggplot2::element_text(size=12),
plot.title = ggplot2::element_text(size = 12, colour = "black", face = "bold", hjust = 0.5, vjust=0.01))
}
#The third panel plots the model formulae
p3 <- plotFormula(dynrModel, ParameterAs=sprintf("%.2f", x@transformed.parameters), printDyn=printDyn, printMeas=printMeas, textsize=textsize)
#Organize the panels using the multiplot function
if (dynrModel$num_regime > 1){
multiplot(p1, p2, p3, cols = 1, layout=matrix(c(1,1,2,3), nrow=2, byrow=TRUE))
}else{
multiplot(p1, p3, layout=matrix(c(1,1,2,2), nrow=2, byrow=TRUE))
}
}
plotdf <- function(vec_tex){
dataframe=data.frame(text=sapply(paste0("$",vec_tex,"$"),function(x){as.character(latex2exp::TeX(x))}))
dataframe$x<-0
return(dataframe)
}
##' Plot of the estimated frequencies of the regimes across all individuals and time points
##' based on their smoothed regime probabilities
##'
##' @param res The dynr object returned by dynr.cook().
##' @param dynrModel The model object to plot.
##' @param names.regime (optional) Names of the regimes (must match the length of the number of regimes)
##' @param title (optional) Title of the plot.
##' @param xlab (optional) Label of the x-axis.
##' @param ylab (optional) Label of the y-axis.
##' @param textsize (default = 12) Text size for the axis labels and title (= textsize + 2).
##' @param print (default = TRUE) A flag for whether the plot should be printed.
##'
##' @return ggplot object
dynr.plotFreq <- function(res, dynrModel, names.regime, title, xlab, ylab, textsize=12,print=TRUE) {
if(missing(names.regime)){names.regime<-paste0("Regime",1:dynrModel@num_regime)}
if(missing(title)){title<-'Counts by most probable regime'}
if(missing(xlab)){xlab<-'Regime'}
if(missing(ylab)){ylab<-'Counts'}
regime <- NULL
highProbR <- data.frame(regime=factor(apply(res@pr_t_given_T,2,which.max),
levels=1:dynrModel@num_regime,
labels=names.regime))
p2 <- ggplot2::ggplot(data = highProbR, ggplot2::aes(factor(regime))) +
ggplot2::geom_bar() +
ggplot2::scale_fill_brewer(palette = 3) +
ggplot2::xlab(xlab) + ggplot2::ylab(ylab) + ggplot2::labs(fill = '') +
ggplot2::ggtitle(title) +
ggplot2::theme(axis.title = ggplot2::element_text(size=textsize+2),
axis.text=ggplot2::element_text(size=textsize),
plot.title = ggplot2::element_text(size = textsize, colour = "black", face = "bold", hjust = 0.5, vjust=0.01))
if (print) print(p2)
return(p2)
}
##' Plot the formula from a model
##'
##' @param dynrModel The model object to plot.
##' @param ParameterAs The parameter values or names to plot. The underscores in parameter names are saved for
##' use of subscripts. Greek letters can be specified as corresponding LaTeX symbols without backslashes (e.g., "lambda")
##' and printed as greek letters.
##' @param printDyn A logical value indicating whether or not to plot the formulas for the dynamic model.
##' @param printMeas A logical value indicating whether or not to plot the formulas for the measurement model
##' @param printRS logical. Whether or not to print the regime-switching model. The default is FALSE.
##' @param textsize The text size use in the plot.
##'
##' @details
##' This function typesets a set of formulas that represent the model. Typical inputs to the \code{ParameterAs} argument are (1) the starting values for a model, (2) the final estimated values for a model, and (3) the parameter names. These are accessible with (1) \code{model$xstart}, (2) \code{coef(cook)}, and (3) \code{model$param.names} or \code{names(coef(cook))}, respectively.
##'
##' @return ggplot object
plotFormula <- function(dynrModel, ParameterAs, printDyn=TRUE, printMeas=TRUE,
printRS=FALSE, textsize=4){
dynrModel <- PopBackModel(dynrModel, LaTeXnames(ParameterAs, latex = FALSE))
#Dynamic model
if (printDyn) {
if ("dynrDynamicsMatrix" %in% class(dynrModel$dynamics)){
state.names <- (dynrModel$measurement)$state.names
exo.names <- (dynrModel$dynamics)$covariates
dynrModel@dynamics@values.dyn <- lapply((dynrModel$dynamics)$values.dyn, preProcessNames,state.names,state.names)
dynrModel@dynamics@values.exo <- lapply((dynrModel$dynamics)$values.exo, preProcessNames,state.names,exo.names)
dynrModel@dynamics@values.int <- lapply((dynrModel$dynamics)$values.int, preProcessNames,state.names)
}
dyn.df <- data.frame(text="bold('Dynamic Model')",x=0)
nRegime=ifelse(class(dynrModel@dynamics)=="dynrDynamicsFormula",
length(dynrModel@dynamics@formula),
length(dynrModel@dynamics@values.dyn))
dyn_tex=printex(dynrModel@dynamics,AsMatrix=FALSE)
nEq <- length(dyn_tex[[1]])
# noise_tex <- paste0(" + ", ifelse(dynrModel@dynamics@isContinuousTime, "d", ""), "w_{",1:nEq, "}(t)")
noise_tex <- paste0(ifelse(dynrModel@dynamics@isContinuousTime, "d", ""),
"w_{", 1:nEq, "}(t)")
for (i in 1:nRegime){
if (nRegime>1){
dyn.df <- rbind(dyn.df,data.frame(text=paste0("'Regime ",i,":'"),x=0))
}
if (length(dynrModel@noise@values.latent)==1){
values.latent.mat <- dynrModel@noise@values.latent[[1]]
}else if (length(dynrModel@noise@values.latent)==nRegime){
values.latent.mat <- dynrModel@noise@values.latent[[i]]
}else{stop("The number of regimes implied by the dynamic noise structure does not match the number of regimes in the dynamic model.")}
pos.zero <- diag(values.latent.mat) == 0
# noise_tex[diag(values.latent.mat)==0] <- ""
# dyn.df <- rbind(dyn.df,plotdf(LaTeXnames(paste0(dyn_tex[[i]], noise_tex))))
dyn.df <- rbind(dyn.df,
plotdf(LaTeXnames(
paste0(dyn_tex[[i]],
ifelse(pos.zero, "",
paste0(" + ", noise_tex))))),
plotdf(LaTeXnames(
paste0(ifelse(pos.zero, "",
paste0(noise_tex, " \\sim ", "N(0,\\,",
diag(values.latent.mat), ")"))))) )
}
} else {
dyn.df <- plotdf("")
}
#Measurement model
if (printMeas) {
meas.df <- data.frame(text="bold('Measurement Model')", x=0)
nRegime <- length(dynrModel$measurement$values.load)
meas_tex <- printex(dynrModel$measurement, AsMatrix=FALSE)
nEq <- length(meas_tex[[1]])
# noise_tex <- paste0(" + ", "epsilon_{", 1:nEq, "}")
noise_tex <- paste0("epsilon_{", 1:nEq, "}")
for (i in 1:nRegime){
if (nRegime > 1){
meas.df <- rbind(meas.df, data.frame(text=paste0("'Regime ",i,":'"), x=0))
}
if (length(dynrModel$noise$values.observed) == 1){
values.observed.mat <- dynrModel$noise$values.observed[[1]]
}else if (length(dynrModel$noise$values.observed) == nRegime){
values.observed.mat <- dynrModel$noise$values.observed[[i]]
}else{stop("The number of regimes implied by the measurement noise structure does not match the number of regimes in the measurement model.")}
pos.zero <- diag(values.observed.mat) == 0
# noise_tex[diag(values.observed.mat)==0] <- ""
# meas.df<-rbind(meas.df,plotdf(LaTeXnames(paste0(meas_tex[[i]], noise_tex))))
meas.df <- rbind(meas.df,
plotdf(LaTeXnames(
paste0(meas_tex[[i]],
ifelse(pos.zero, "", paste0(" + ", noise_tex, ",")),
"\\;\\,",
ifelse(pos.zero,
"",
paste0(noise_tex, " \\sim ", "N(0,\\,", diag(values.observed.mat), ")"))))))
}
} else {
meas.df <- plotdf("")
}
# Regime-switching model
if (printRS) {
rs.df <- data.frame(text="bold('Regime-switching Model')", x=0)
if (any(dim(dynrModel@regimes@params) == 0)) {
warning("No Regime-switching Model. Please turn off 'printRS'.")
} else {
# helper function to make sub. e.g., a_111 to a_{111}
to_sub <- function(parnames) {
idx_sub <- grepl('_', parnames)
if ( sum(idx_sub)==0 ) { return(parnames) }
subs <- parnames[idx_sub]
subs_spl <- strsplit(subs, '_')
subbed <- sapply(subs_spl, function(sub) {
paste0(sub[1], '_', paste0('{', sub[2], '}') )
})
parnames[idx_sub] <- subbed
parnames
}
covars <- dynrModel@regimes@covariates
if ( length(covars) == 0 ) {
covars <- c("")
} else {
covars <- c("", paste0("*", covars) )
}
covars_n <- length(covars)
# TODO: check nregime to fit with others
nregime <- nrow(dynrModel@regimes@params)
# column index for the intercepts of regimes
inter_col <- cumsum(c(1, rep(covars_n, nregime-2)))
# i:(i+ncov-1)
parnames <- to_sub( dynrModel@regimes@paramnames )
param_rs <- matrix(parnames, nregime)
param_inter <- param_rs[, inter_col, drop=FALSE]# intercepts
refrow <- dynrModel@regimes@refRow
if (length(refrow) != 0) {
for ( co in inter_col ) {
param_rs[-refrow, co] <-
paste0(param_rs[refrow, co], ' + ', param_rs[-refrow, co])
}
}
param_rs <- matrix(sweep(param_rs, 2, covars, FUN=paste0),
nrow=nregime)
ll <- vector('list', length=length(inter_col))
for ( i in 1:length(inter_col) ) {
ll[[i]] <-
param_rs[, inter_col[i]:(inter_col[i]+covars_n-1),
drop=FALSE]
}
sl <- sapply(ll, function(l) {
apply(l, 1, paste, collapse=" + ")
})
# transition probability subs
transub <- array(c(row(sl), col(sl)), c(nrow(sl), ncol(sl), 2))
transubed <- apply(transub, c(1,2), paste, collapse="")
rs.df <- rbind(rs.df,
plotdf(paste0("Log-odds(p_{",
transubed, "}) = ",
sl)))
}
} else {
rs.df <- plotdf("")
}
plot.df <- rbind(meas.df, dyn.df, rs.df)
plot.df$y=seq((dim(plot.df)[1]-1)*5+1, 1, by=-5)
x <- NULL
y <- NULL
fig<-ggplot2::ggplot(plot.df, ggplot2::aes(x=x, y=y, label=text))+
ggplot2::geom_text(parse=TRUE,size=textsize)+
ggplot2::theme_void() + ggplot2::ylim(min(plot.df$y)-5, max(plot.df$y)+5)
return(fig)
}
##' The ggplot of the smoothed state estimates and the most likely regimes
##'
##' @aliases autoplot.dynrCook
##'
##' @param res The dynr object returned by \code{dynr.cook()}.
##' @param dynrModel The model object to plot.
##' @param style The style of the plot. If style is 1 (default), user-selected smoothed state variables are plotted. If style is 2, user-selected observed-versus-predicted values are plotted.
##' @param numSubjDemo The number of subjects to be randomly selected for plotting.
##' @param idtoPlot Values of the ID variable to plot.
##' @param names.state (optional) The names of the states to be plotted, which should be a subset of the state.names slot of the measurement slot of dynrModel.
##' @param names.observed (optional) The names of the observed variables to be plotted, which should be a subset of the obs.names slot of the measurement slot of dynrModel.
##' @param names.regime (optional) The names of the regimes to be plotted, which can be missing.
##' @param shape.values (optional) A vector of values that correspond to the shapes of the points, which can be missing. See the R documentation on pch for details on possible shapes.
##' @param title (optional) A title of the plot.
##' @param ylab (optional) The label of the y axis.
##' @param is.bw Is plot in black and white? The default is FALSE.
##' @param colorPalette A color palette for lines and dots. It is a value passed to the palette argument of the \code{ggplot2::scale_colour_brewer()} function. These palettes are in the R package \pkg{RColorBrewer}. One can find them by attaching the package with \code{library(RColorBrewer)} and run \code{display.brewer.all()}.
##' @param fillPalette A color palette for blocks. It is a value passed to the palette argument of the \code{ggplot2::scale_fill_brewer()} function. These palettes are in the package \pkg{RColorBrewer}. One can find them by attaching the package with \code{library(RColorBrewer)} and run \code{display.brewer.all()}.
##' @param mancolorPalette (optional) A color palette for manually scaling the colors of lines and dots. It is a vector passed to the values argument of the \code{ggplot2::scale_colour_manual} function.
##' @param manfillPalette (optional) A color palette for manually scaling the colors of filled blocks. It is a vector passed to the values argument of the \code{ggplot2::scale_fill_manual} function.
##' @param ... A list of elements that modify the existing ggplot theme. Consult the \code{ggplot2::theme()} function in the R package \pkg{ggplot2} for more options.
##'
##' @details
##' This function outputs a ggplot layer that can be modified using functions in the package \pkg{ggplot2}. That is, one can add layers, scales, coords and facets with the "+" sign. In an example below, the \code{ggplot2::ylim()} function is used to modify the limits of the y axis of the graph. More details can be found on \url{https://ggplot2.tidyverse.org/} and \url{https://ggplot2.tidyverse.org/reference/}.
##'
##' The two functions \code{dynr.ggplot()} and \code{autoplot()} as identical aliases of one another. The \code{autoplot()} function is an S3 method from the package \pkg{ggplot2} that allows many objects to be plotted and works like the base \code{plot()} function.
##'
##' @return ggplot object
##'
##' @examples
##' # The following code is part of a demo example in dynr
##' \dontrun{
##' demo(RSLinearDiscreteYang, package='dynr')
##' p <- dynr.ggplot(yum, dynrModel = rsmod, style = 1,
##' names.regime = c("Deactivated", "Activated"),
##' title = "(B) Results from RS-AR model", numSubjDemo = 1,
##' shape.values = c(1),
##' text = element_text(size = 16),
##' is.bw = TRUE)
##' # One can modify the limits on the y axis by using '+'
##' p + ggplot2::ylim(-2, 4)
##'
##' autoplot(yum, dynrModel = rsmod, style = 1,
##' names.regime = c("Deactivated", "Activated"),
##' title = "(B) Results from RS-AR model", numSubjDemo = 1,
##' shape.values = c(1),
##' text = element_text(size = 16),
##' is.bw = TRUE)
##' }
dynr.ggplot <- function(res, dynrModel, style = 1,
numSubjDemo=2, idtoPlot=c(),
names.state,
names.observed,
names.regime,
shape.values,
title,
ylab,
is.bw=FALSE,
colorPalette="Set2",
fillPalette="Set2",
mancolorPalette, manfillPalette, ...){
data.dynr=dynrModel@data
dim_latent_var=dim(res@eta_smooth_final)[1]
num_regime=dim(res@pr_t_given_T)[1]
if(missing(names.observed)){names.observed=dynrModel@measurement@obs.names}
observed=which(dynrModel@measurement@obs.names%in%names.observed)
if ((length(observed)>8)&(missing(mancolorPalette))&(style==2)){stop("You provided too many variables than the default color palette can handle.\nPlease consider specify the mancolorPalette argument.")}
if(missing(names.state)){names.state=dynrModel@measurement@state.names}
states=which(dynrModel@measurement@state.names%in%names.state)
if(missing(shape.values)){
if (style==1) shape.values=states#48+states
if (style==2) shape.values=c(rep(32,length(observed)),observed)#pre-obs
}else{
if (style==2) shape.values=c(rep(32,length(shape.values)), shape.values)#pre-obs
}
if (style==1) line.values=rep(1,length(states))
if (style==2) line.values=rep(c(1,0),each=length(observed))#pre-obs
num_sbj=length(unique(data.dynr$id))
if(length(idtoPlot)<1){
randid =sample(unique(data.dynr$id),numSubjDemo)
}else {randid=idtoPlot}
if(missing(title)){
if (style==1) title=" "
if (style==2) title="Observed vs Predicted"
}
if(missing(ylab)){
if (style==1) ylab="Smoothed State Values"
if (style==2) ylab="Values"
}
data.plot<-data.frame(id=as.factor(data.dynr$id),time=data.dynr$time)
if(num_regime==1){
names.id.vars <- c("id","time")
if (style==1){
names.measure.vars <- names.state
data.plot <- cbind(data.plot, matrix(res@eta_smooth_final[states,],ncol=length(states),byrow=TRUE))
lines.var <- names.state
points.var <- names.state
}#end of style 1
if (style==2){
dynrModel=PopBackModel(dynrModel, LaTeXnames(res@transformed.parameters, latex = FALSE))
dynrModel@measurement@values.load=lapply(dynrModel@measurement@values.load, function(x){matrix(as.numeric(x),ncol=ncol(x))})
num_regime_meas <- length(dynrModel@measurement@values.load)
if (length(dynrModel@measurement@values.exo)!=0){
dynrModel@measurement@values.exo=lapply(dynrModel@measurement@values.exo, function(x){matrix(as.numeric(x),ncol=ncol(x))})
}
if (length(dynrModel@measurement@values.int)!=0){
dynrModel@measurement@values.int=lapply(dynrModel@measurement@values.int, function(x){matrix(as.numeric(x),ncol=ncol(x))})
}
#predicted scores
if (num_regime_meas==1){#one-regime measurement model
predicted=dynrModel@measurement@values.load[[1]]%*%res@eta_smooth_final
if (length(dynrModel@measurement@values.int)!=0){
predicted=predicted+dynrModel@measurement@values.int[[1]]%*%matrix(rep(1,ncol(res@eta_smooth_final)), nrow=1)
}
if (length(dynrModel@measurement@values.exo)!=0){
exo.mat=t(as.matrix(dynrModel@data$covariates[,paste0("covar",sapply(dynrModel@measurement@exo.names, function(x){which(dynrModel@data$covariate.names==x)}))]))
predicted=predicted+dynrModel@measurement@values.exo[[1]]%*%exo.mat
}
}#end of one-regime meas
names.measure.vars <- c(paste0(names.observed, ".predicted"), paste0(names.observed, ".observed"))
data.plot <- cbind(data.plot, matrix(predicted[observed,],ncol=length(observed),byrow=TRUE, dimnames=list(NULL, paste0(names.observed, ".predicted"))))
data.plot <- cbind(data.plot, data.dynr$observed[,observed])
lines.var <- paste0(names.observed, ".predicted")
points.var <- paste0(names.observed, ".observed")
if (length(observed)<8){
fancy8 <- c("#edaac1","#a7dab0","#dfb3e2","#8bc5ed","#bcb8ec","#e6b296","#d2d39d","#82d9d6")
mancolorPalette=rep(fancy8[observed],2)
}
}#end of style 2
names(data.plot) <- c(names.id.vars, names.measure.vars)
data_long <- reshape2::melt(data.plot[data.plot$id%in%randid,],
id.vars=names.id.vars,
measure.vars=names.measure.vars,
value.name="value")
data_long$variable<-factor(data_long$variable, levels=sort(names.measure.vars))
value <- NULL
variable <- NULL
partial.plot<-ggplot2::ggplot(data_long, ggplot2::aes(x=time, y=value, colour = variable, shape = variable, linetype = variable)) +
ggplot2::geom_line(data=data_long[data_long$variable%in%lines.var,], size=1) +
ggplot2::geom_point(data=data_long[data_long$variable%in%points.var,], size=4) +
ggplot2::scale_linetype_manual(labels=sort(names.measure.vars), values=line.values[order(names.measure.vars)])+
ggplot2::scale_shape_manual(labels=sort(names.measure.vars), values=shape.values[order(names.measure.vars)])+
ggplot2::facet_wrap(~id)+
ggplot2::labs(title = title, y=ylab, ggtitle="")+
ggplot2::theme(plot.title=ggplot2::element_text(lineheight=.8, face="bold", hjust = 0.5, vjust=0.01), legend.position="bottom",legend.text=ggplot2::element_text(lineheight=0.8, face="bold"),...)
}else{#more than two regimes
if(missing(names.regime)){names.regime = 1:num_regime}
findRegime = function(prop){
MostLikelyRegime = names.regime[prop==max(prop)]
return(MostLikelyRegime[1])
#If the probabilities are equal, return the first regime
}
addendtime <- function(data_frame){
data_frame$endtime <- c(data_frame$time[2:nrow(data_frame)], data_frame$time[nrow(data_frame)]+min(diff(data_frame$time)))
return(data_frame)
}
data.plot <- plyr::ddply(data.plot,"id",addendtime)
data.plot$regime<-as.factor(apply(res@pr_t_given_T,2,findRegime))
names.id.vars<-c("id","time","endtime","regime")
if (style==1){
names.measure.vars<-names.state
data.plot<-cbind(data.plot,matrix((res@eta_smooth_final[states,]),ncol=length(states),byrow=TRUE))
lines.var <- names.state
points.var <- names.state
}#end of style 1
if (style==2){
dynrModel=PopBackModel(dynrModel, LaTeXnames(res@transformed.parameters, latex = FALSE))
dynrModel@measurement@values.load=lapply(dynrModel@measurement@values.load, function(x){matrix(as.numeric(x),ncol=ncol(x))})
num_regime_meas <- length(dynrModel@measurement@values.load)
if (length(dynrModel@measurement@values.exo)!=0){
dynrModel@measurement@values.exo=lapply(dynrModel@measurement@values.exo, function(x){matrix(as.numeric(x),ncol=ncol(x))})
}
if (length(dynrModel@measurement@values.int)!=0){
dynrModel@measurement@values.int=lapply(dynrModel@measurement@values.int, function(x){matrix(as.numeric(x),ncol=ncol(x))})
}
#predicted scores
if (num_regime_meas==1){#one-regime measurement model
predicted=matrix(dynrModel@measurement@values.load[[1]][observed,],nrow=length(observed))%*%res@eta_smooth_final
if (length(dynrModel@measurement@values.int)!=0){
predicted=predicted+matrix(dynrModel@measurement@values.int[[1]][observed,],nrow=length(observed))%*%matrix(rep(1,ncol(res@eta_smooth_final)), nrow=1)
}
if (length(dynrModel@measurement@values.exo)!=0){
exo.mat=t(as.matrix(dynrModel@data$covariates[,paste0("covar",sapply(dynrModel@measurement@exo.names, function(x){which(dynrModel@data$covariate.names==x)}))]))
predicted=predicted+matrix(dynrModel@measurement@values.exo[[1]][observed,],nrow=length(observed))%*%exo.mat
}
#end of one-regime meas
}else if (num_regime_meas==num_regime){
ntotaltime=ncol(res@eta_smooth_final)
predicted=matrix(rep(0,length(observed)*ntotaltime),ncol=ntotaltime)
for (index in 1:ntotaltime){
predicted[,index]=matrix(dynrModel@measurement@values.load[[data.plot$regime[index]]][observed,],nrow=length(observed))%*%matrix(res@eta_smooth_final[,index], ncol=1)
if (length(dynrModel@measurement@values.int)!=0){
predicted[,index]=matrix(predicted[,index]+dynrModel@measurement@values.int[[data.plot$regime[index]]][observed,],nrow=length(observed))
}
if (length(dynrModel@measurement@values.exo)!=0){
exo.mat=t(as.matrix(dynrModel@data$covariates[index,paste0("covar",sapply(dynrModel@measurement@exo.names, function(x){which(dynrModel@data$covariate.names==x)}))]))
predicted[,index]=matrix(predicted[,index]+dynrModel@measurement@values.exo[[data.plot$regime[index]]][observed,],nrow=length(observed))%*%exo.mat
}
}
#end of regime specific measurement
}
names.measure.vars <- c(paste0(names.observed, ".predicted"), paste0(names.observed, ".observed"))
data.plot <- cbind(data.plot, matrix(predicted,ncol=length(observed), byrow=TRUE, dimnames=list(NULL, paste0(names.observed, ".predicted"))))
data.plot <- cbind(data.plot, data.dynr$observed[,observed])
lines.var <- paste0(names.observed, ".predicted")
points.var <- paste0(names.observed, ".observed")
if (length(observed)<=8){
fancy8 <- c("#edaac1","#a7dab0","#dfb3e2","#8bc5ed","#bcb8ec","#e6b296","#d2d39d","#82d9d6")
manfillPalette=fancy8[9-(1:num_regime)]
mancolorPalette=rep(fancy8[observed],2)
}
}#end of style 2
names(data.plot)<-c(names.id.vars,names.measure.vars)
data_long<- reshape2::melt(data.plot[data.plot$id%in%randid,],
id.vars=names.id.vars,
measure.vars=names.measure.vars,
value.name="value")
data_long$variable<-factor(data_long$variable, levels=sort(names.measure.vars))
#data_long$statenumber<-as.factor(sub("state","",data_long$variable))
endtime <- NULL
regime <- NULL
partial.plot<-ggplot2::ggplot(data_long,ggplot2::aes(x=time, y=value, group=variable)) +
ggplot2::geom_rect(ggplot2::aes(xmin=time, xmax=endtime, ymin=-Inf, ymax=Inf, fill=regime), alpha=.15) +
ggplot2::geom_line(data=data_long[data_long$variable%in%lines.var,], size=1, ggplot2::aes(color=variable, linetype=variable)) +
ggplot2::geom_point(data=data_long[data_long$variable%in%points.var,], size=3, ggplot2::aes(color=variable, shape=variable)) +
ggplot2::scale_linetype_manual(labels=sort(names.measure.vars), values=line.values[order(names.measure.vars)])+
ggplot2::scale_shape_manual(labels=sort(names.measure.vars), values=shape.values[order(names.measure.vars)])+
#geom_text(size=1, ggplot2::aes(label=statenumber,color=variable))+
ggplot2::facet_wrap(~id)+
ggplot2::labs(title = title,y=ylab,ggtitle="")+
ggplot2::theme(plot.title=ggplot2::element_text(lineheight=.8, face="bold", hjust = 0.5, vjust=0.01), legend.position="bottom",legend.text=ggplot2::element_text(lineheight=0.8, face="bold"),...)
}
if(is.bw){
bw.plot<-partial.plot+
ggplot2::scale_fill_grey(start = 0.1, end = 0.9)+
ggplot2::scale_color_grey(labels=sort(names.measure.vars), start = 0.2, end = 0.7)
return(bw.plot)
}else{
default.plot<-partial.plot+
ggplot2::scale_colour_brewer(labels=sort(names.measure.vars), palette=colorPalette)+
ggplot2::scale_fill_brewer(palette=fillPalette)
if(!missing(mancolorPalette)){
manual.plot<-partial.plot+ggplot2::scale_colour_manual(labels=sort(names.measure.vars), values=mancolorPalette[order(names.measure.vars)])
}
if(!missing(manfillPalette)){
if(!exists("manual.plot")){manual.plot<-partial.plot}
manual.plot<-manual.plot+ggplot2::scale_fill_manual(values=manfillPalette)
}
if (exists("manual.plot")){
return(manual.plot) #print(manual.plot)
}else{
return(default.plot)#print(default.plot)
}
}
}
##' @rdname dynr.ggplot
##'
##' @param object The same as res. The dynr object returned by dynr.cook().
##'
##' @return ggplot object
autoplot.dynrCook <- function(object, dynrModel, style = 1,
numSubjDemo=2, idtoPlot=c(),
names.state,
names.observed,
names.regime,
shape.values,
title,
ylab,
is.bw=FALSE,
colorPalette="Set2",
fillPalette="Set2",
mancolorPalette, manfillPalette, ...){
dynr.ggplot(object, dynrModel=dynrModel, style = style,
numSubjDemo=numSubjDemo, idtoPlot=idtoPlot,
names.state=names.state,
names.observed=names.observed,
names.regime=names.regime,
shape.values=shape.values,
title=title,
ylab=ylab,
is.bw=is.bw,
colorPalette=colorPalette,
fillPalette=fillPalette,
mancolorPalette=mancolorPalette, manfillPalette=manfillPalette, ...)
}
##' The ggplot of the outliers estimates.
##'
##' @param object A dynrTaste object.
##' @param numSubjDemo The number of subjects, who have
##' largest joint chi-square statistic, to be selected for plotting.
##' @param idtoPlot Values of the ID variable to plot.
##' @param names.state (optional) The names of the states to be plotted, which should be a subset of the state.names slot of the measurement slot of dynrModel. If NULL, the t statistic plots for all state variables will be included.
##' @param names.observed (optional) The names of the observed variables to be plotted, which should be a subset of the obs.names slot of the measurement slot of dynrModel. If NULL, the t statistic plots for all observed variables will be included.
##' @param ... Place holder for other arguments. Please do not use.
##'
##' @return a list of ggplot objects for each ID.
##' The plots of chi-square statistics (joint and independent),
##' and the plots of t statistic for \code{names.state} and \code{names.observed} will be included.
##' Users can modify the ggplot objects using ggplot grammar.
##' If a \code{filename} is provided, a pdf of plots will be saved additionally.
autoplot.dynrTaste <- function(object,
numSubjDemo=2, idtoPlot=NULL,
names.state=NULL, names.observed=NULL, ...) {
if ( !inherits(object, "dynrTaste") ) {
stop("dynrTaste object is required.")
}
stopifnot(numSubjDemo >= 1)
numSubjDemo <- as.integer(numSubjDemo)
conf_level <- object$conf.level
tstart <- object$tstart
chi_jnt <- object$chi.jnt
id <- object$id
id_unq <- unique(id)
id_n <- length(id_unq)
lat_name <- rownames(object$t.inn)
obs_name <- rownames(object$t.add)
lat_n <- length(lat_name)
obs_n <- length(obs_name)
if( is.null(names.state) ) {
lat_toplot <- lat_name
} else {
if ( !all(names.state %in% lat_name) ) {
stop("'names.state' should be a subset of the latent variables.")
}
lat_toplot <- names.state
}
if( is.null(names.observed) ) {
obs_toplot <- obs_name
} else {
if ( !all(names.observed %in% obs_name) ) {
stop("'names.observed' should be a subset of the observed variables.")
}
obs_toplot <- names.observed
}
if ( is.null(idtoPlot) ) {
chi_max <- vector("numeric", id_n)
for(i in 1:id_n){
begT <- tstart[i] + 1
endT <- tstart[i+1]
chi_max[i] <- max( chi_jnt[begT:endT] )
}
id_toplot <- id_unq[ order(chi_max, decreasing=TRUE) ][
1:numSubjDemo ]
} else {
id_toplot <- idtoPlot
if ( !all(id_toplot %in% id_unq) ) {
stop("Not all ID are in the data.")
}
}
### create ggplot objects
ggadd_hline <- function(yintercept) {
geom_hline(yintercept=yintercept,
colour="black", linetype="dashed", alpha=0.5)
}
ggadd_point <- function(data, aes_y) {
geom_point(data=data, aes_string(x="time", y=aes_y),
colour="red", size=2, alpha=0.5)
}
ggadd_theme <- function() {
theme(plot.title=element_text(size = rel(0.85)),
axis.title.x=element_blank(),
#axis.text.x=element_blank(),
axis.ticks.x=element_blank(),
panel.grid.major=element_blank(),
panel.grid.minor=element_blank(),
panel.background=element_blank(),
plot.margin = margin(0, 0, 0, 0),
axis.line=element_line(colour = "black"))
}
## Joint chi-square
id_toplot_n <- length(id_toplot)
time <- object$time
chi_jnt_shk <- NULL
chi_df_jnt <- data.frame(id = id, time = time,
chi_jnt = object$chi.jnt,
chi_jnt_shk = object$chi.jnt.shock)
# subset data for id to plot
chi_df_jnt_plot <- subset(chi_df_jnt, is.element(id, id_toplot))
# subset data only for shocks
chi_df_jnt_shk <- subset(chi_df_jnt_plot, chi_jnt_shk)
qchi_jnt <- qchisq(conf_level, lat_n + obs_n)
plots_jnt <- lapply(id_toplot, function(id_i) {
chi_df_jnt_shk_i <- subset(chi_df_jnt_shk, id==id_i)
ggplot2::ggplot(subset(chi_df_jnt_plot, id==id_i),
aes(x=time, y=chi_jnt)) +
geom_line(alpha=0.5) +
labs(title="") +
ylab(expression(paste("Joint ", Chi^2))) +
list(ggadd_point(chi_df_jnt_shk_i, "chi_jnt"),
ggadd_hline(qchi_jnt), ggadd_theme())
})
#### Independend chi-square for innovative ####
chi_inn <- chi_inn_shk <- NULL
chi_df_inn <- data.frame(id = id, time = time,
chi_inn = object$chi.inn,
chi_inn_shk = object$chi.inn.shock)
chi_df_inn_plot <- subset(chi_df_inn, is.element(id, id_toplot))
chi_df_inn_shk <- subset(chi_df_inn_plot, chi_inn_shk)
qchi_inn <- qchisq(conf_level, lat_n)
plots_inn <- lapply(id_toplot, function(id_i) {
chi_df_inn_shk_i <- subset(chi_df_inn_shk, id==id_i)
ggplot2::ggplot(subset(chi_df_inn_plot, id==id_i),
aes(x=time, y=chi_inn)) +
geom_line(alpha=0.5) +
labs(title="") +
ylab(expression(paste("Indep. ", Chi^2, " - innovative"))) +
list(ggadd_point(chi_df_inn_shk_i, "chi_inn"),
ggadd_hline(qchi_inn), ggadd_theme())
})
#### Independent chi-square for additive ####
chi_add <- chi_add_shk <- NULL
chi_df_add <- data.frame(id = id, time = time,
chi_add = object$chi.add,
chi_add_shk = object$chi.add.shock)
chi_df_add_plot <- subset(chi_df_add, is.element(id, id_toplot))
chi_df_add_shk <- subset(chi_df_add_plot, chi_add_shk)
plots_add <- lapply(id_toplot, function(id_i) {
chi_df_add_shk_i <- subset(chi_df_add_shk, id==id_i)
ggplot2::ggplot(subset(chi_df_add_plot, id==id_i),
aes(x=time, y=chi_add)) +
geom_line(alpha=0.5) +
labs(title="") +
ylab(expression(paste("Indep. ", Chi^2, " - additive"))) +
list(ggadd_point(chi_df_add_shk_i, "chi_add"),
ggadd_hline(qchisq(conf_level, lat_n)), ggadd_theme())
})
#### t-statistic for innovative ####
t_inn <- object$t.inn
t_inn_name <- row.names(t_inn)
row.names(t_inn) <- paste0(t_inn_name, "_t")
t_df_inn <- data.frame(id = id, time = time,
t(t_inn), t(object$t.inn.shock))
t_df_inn_plot <- subset(t_df_inn, is.element(id, id_toplot))
plots_t_inn <- lapply(id_toplot, function(id_i) {
t_df_inn_plot_i <- subset(t_df_inn_plot, id==id_i)
lapply(lat_toplot, function(lat_i) {
t_df_inn_shk_i <- subset(t_df_inn_plot_i,
t_df_inn_plot_i[[lat_i]])
qt_i <- qt(conf_level, nrow(t_df_inn_plot_i) - obs_n)
ggplot2::ggplot(t_df_inn_plot_i,
aes_string(x="time", y=paste0(lat_i, "_t"))) +
geom_line(alpha=0.5) +
labs(title="") +
ylab(paste("t_[", lat_i, "]")) +
list(ggadd_point(t_df_inn_shk_i, paste0(lat_i, "_t")),
ggadd_hline(qt_i), ggadd_hline(-qt_i), ggadd_theme())
})
})
#### t-statistic for additive ####
t_add <- object$t.add
t_add_name <- row.names(t_add)
row.names(t_add) <- paste0(t_add_name, "_t")
t_df_add <- data.frame(id = id, time = time,
t(t_add), t(object$t.add.shock))
t_df_add_plot <- subset(t_df_add, is.element(id, id_toplot))
plots_t_add <- lapply(id_toplot, function(id_i) {
t_df_add_plot_i <- subset(t_df_add_plot, id==id_i)
t_add_df_i <- nrow(t_df_add_plot_i)
lapply(obs_toplot, function(obs_i) {
t_df_add_shk_i <- subset(t_df_add_plot_i,
t_df_add_plot_i[[obs_i]])
qt_i <- qt(conf_level, t_add_df_i - obs_n)
ggplot2::ggplot(t_df_add_plot_i,
aes_string(x="time", y=paste0(obs_i, "_t"))) +
geom_line(alpha=0.5) +
labs(title="") +
ylab(paste("t_[", obs_i, "]")) +
list(ggadd_point(t_df_add_shk_i, paste0(obs_i, "_t")),
ggadd_hline(qt_i), ggadd_hline(-qt_i), ggadd_theme())
})
})
gg_objects <- mapply(function(jnt_i, inn_i, add_i,
t_inn_i, t_add_i, id_toplot_i) {
gg_i <- c(list(jnt_i, inn_i, add_i), t_inn_i, t_add_i)
names(gg_i) <- c("chi_jnt", "chi_inn", "chi_add",
paste0("t_[", lat_toplot, "]"),
paste0("t_[", obs_toplot, "]"))
gg_i
},
plots_jnt, plots_inn, plots_add, plots_t_inn, plots_t_add, id_toplot,
SIMPLIFY=FALSE)
names(gg_objects) <- id_toplot
# if ( !is.null(filename) && is.character(filename) ) {
# gg_plots <- mapply(function(plots_i, id_toplot_i) {
# plots_arr <- ggpubr::ggarrange(
# plotlist=plots_i, nrow=length(plots_i), ncol=1, align="v")
# ggpubr::annotate_figure(
# plots_arr,
# top=ggpubr::text_grob(paste0("ID: ", id_toplot_i),
# color="black", face="bold", size=14),
# bottom=ggpubr::text_grob("time", color="black", size=11))
# },
# gg_objects, id_toplot,
# SIMPLIFY=FALSE)
#
# id_v <- chi_df_jnt_plot[["id"]]
# id_time_length <- vector("numeric", id_toplot_n)
# # time points for each id
# for (i in 1:id_toplot_n) {
# id_time_length[i] <- sum(is.element(id_v, id_toplot[i]))
# }
# # width: time/6, height: n of gg object * 2
# ggpubr::ggexport(gg_plots, filename=filename,
# width=mean(id_time_length)/6,
# height=length(gg_objects[[1]])*2)
# }
invisible(gg_objects)
}
Try the dynr package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.