Nothing
R/plot.ctsemFit.R
In ctsemOMX: Continuous Time SEM - 'OpenMx' Based Functions
Defines functions
ctPlot
plot.ctsemFit
plot.ctsemMultigroupFit
Documented in
ctPlot
plot.ctsemFit
plot.ctsemMultigroupFit
#' Plot function for ctsemMultigroupFit object
#'
#' Plots \code{\link{ctMultigroupFit}} objects.
#'
#' @param x ctsemMultigroupFit object as generated by \code{\link{ctMultigroupFit}}
#' @param group character string of subgroup to plot. Default of 'show chooser' displays list and lets you select.
#' @param ... additional parameters to pass to \code{\link{plot.ctsemFit}} function.
#' @return Nothing. Side-effect: plots graphs.
#' @method plot ctsemMultigroupFit
#' @export
plot.ctsemMultigroupFit<-function(x,group='show chooser',...){
if(group=='show chooser') {
message(paste0(1:length(x$groups), ' ', x$groups, '\n'))
selection<-as.numeric(readline(
'Enter the number of the group you wish to plot: '))
selection<-x$groups[[selection]]
}
if(group != 'show chooser') selection <- group
tempctobj<-list(ctfitargs=x$ctfitargs, ctmodelobj=x$ctmodelobj, mxobj=x$mxobj[[selection]])
class(tempctobj)<-'ctsemFit'
plot.ctsemFit(tempctobj,...)
}
#' Plotting function for object class ctsemFit
#'
#' Ouputs mean trajectories, autoregression, and crossregression plots.
#' For more customization possibilities, see \code{\link{ctPlot}}.
#'
#' @param x ctsemFit object as generated by \code{\link{ctFit}}.
#' @param resolution Numeric. Plot points between each unit of time. Default of 'auto' adapts to max.time and results in 500 in total.
#' @param wait If true, user is prompted to continue before plotting next graph. If false, graphs are plotted one after another without waiting.
#' @param max.time Time scale on which to plot parameters. If auto, parameters are plotted for full range of observed variables.
#' @param mean if TRUE, plot of means from 0 to max.time included in output.
#' @param withinVariance if TRUE, plot within subject variance / covariance.
#' @param AR if TRUE, plot of autoregressive values from 0 to max.time included in output.
#' @param CR if TRUE, plot of cross regressive values from 0 to max.time included in output.
#' @param standardiseCR if TRUE , cross regression values are standardised based on estimated within subject variance.
#' @param randomImpulse if TRUE (default), plots expected change in processes given a random fluctuation of +1 for each process --
#' plot is then a mixture of DIFFUSION and DRIFT characteristics.
#' @param experimentalImpulse if TRUE (default), plots expected change in processes given an exogenous input of +1 for each process --
#' alternate characterisation of autoregressive and cross regressive plots.
#' @param xlab X axis label.
#' @param ylab Y axis label.
#' @param meansylim Vector of min and max limits for mean trajectory plot. 'auto' calculates automatically.
#' @param ARylim Vector of min and max limits for autoregression plot. 'auto' is c(0,1), and expands if necessary.
#' @param CRylim Vector of min and max limits for cross regression plot. 'auto' is c(-1,1), and expands if necessary.
#' @param ... Other options passed to \code{plot()}.
#' @return Nothing. Side-effect: plots graphs.
#' @method plot ctsemFit
#' @examples
#' ## Examples set to 'donttest' because they take longer than 5s.
#'
#' ### example from Driver, Oud, Voelkle (2015),
#' ### simulated happiness and leisure time with unobserved heterogeneity.
#' \donttest{
#' data(ctExample1)
#' traitmodel <- ctModel(n.manifest=2, n.latent=2, Tpoints=6, LAMBDA=diag(2),
#' manifestNames=c('LeisureTime', 'Happiness'),
#' latentNames=c('LeisureTime', 'Happiness'), TRAITVAR="auto")
#' traitfit <- ctFit(dat=ctExample1, ctmodelobj=traitmodel)
#' plot(traitfit, wait=FALSE)
#' }
#' @export
plot.ctsemFit<-function(x,resolution=50,wait=TRUE,max.time="auto",mean=TRUE,
withinVariance=TRUE,AR=TRUE,CR=TRUE,standardiseCR=FALSE,randomImpulse=FALSE,
experimentalImpulse=FALSE,
xlab="Time",
meansylim='auto',ARylim='auto', CRylim='auto', ylab="Value",...){
message("Plotting fit")
ctfitobj<-x
mxobj<-ctfitobj$mxobj
ctmodelobj<-ctfitobj$ctmodelobj
ctsummary<-summary(ctfitobj,verbose=TRUE)
#read in values
for(i in 1:length(ctsummary)){ #this loop reads in the specified continuous time model so the objects are available
assign(names(ctsummary[i]),eval(parse(text = paste0("ctsummary","$",names(ctsummary[i])))))
}
asymptotes<-ctfitobj$ctfitargs$asymptotes
latentNames<-ctfitobj$ctmodelobj$latentNames
n.latent<-ctfitobj$ctmodelobj$n.latent
n.manifest<-ctfitobj$ctmodelobj$n.manifest
n.TIpred<-ctfitobj$ctmodelobj$n.TIpred
n.TDpred<-ctfitobj$ctmodelobj$n.TDpred
manifestNames <- ctfitobj$ctmodelobj$manifestNames
TDpredNames <- ctfitobj$ctmodelobj$TDpredNames
TIpredNames <- ctfitobj$ctmodelobj$TIpredNames
Tpoints<-ctfitobj$ctmodelobj$Tpoints
stationary<-ctfitobj$ctfitargs$stationary
T0MEANS<-OpenMx::mxEval(T0MEANS, ctfitobj$mxobj,compute=F)
T0VAR<-OpenMx::mxEval(T0VAR, ctfitobj$mxobj,compute=F)
paramlabels<-matrix(c(paste0(rep(latentNames,n.latent),'_',rep(latentNames,each=n.latent))),n.latent,n.latent)
if(max.time=='auto' & ctfitobj$ctfitargs$objective!='cov'){
if(max.time=="auto" & ctfitobj$ctfitargs$objective!='Kalman') max.time <- max(rowSums(as.matrix(mxobj$data$observed[,paste0('dT',1:(Tpoints-1)),drop=FALSE]),na.rm=T)) # max time of plot
if(max.time=="auto" & ctfitobj$ctfitargs$objective=='Kalman') max.time <- sum(mxobj$data$observed[mxobj$data$observed[,'id']==1,'dT1'],na.rm=T) # max time of plot
}
if(max.time=='auto' & ctfitobj$ctfitargs$objective=='cov') stop('max.time argument must be set when plotting covariance based data')
if(resolution=='auto') resolution <- ceiling(500/max.time)
colourvector <- grDevices::rainbow(ncol(DRIFT),v=.8) #set plot colours
plottimes<-matrix(seq(0,max.time,1/resolution)[-1],ncol=1) #time steps
if(mean==TRUE){# 2. plot mean trend
#TD predictors
tdpredeffect<-matrix(0,nrow=length(plottimes),ncol=n.latent) #default 0 matrix
#extract and round absolute times to resolution
# if(n.TDpred >0){
# if(ctfitobj$ctfitargs$objective!='Kalman') data<-mxobj$data$observed
# if(ctfitobj$ctfitargs$objective=='Kalman') {
# data<- ctLongToWide(mxobj$data$observed,id='id',time='dT1',
# manifestNames=manifestNames,TDpredNames=TDpredNames,TIpredNames=TIpredNames)
# data<-data[,-which(colnames(data)=='T0'),drop=F]
# colnames(data)[which(colnames(data) %in% paste0('T',1:(Tpoints-1)))]<-paste0('dT',1:(Tpoints-1))
# }
#
# obstimes<-cbind(
# matrix(rep(0,nrow(data)),ncol=1),
# matrix(
# apply(data[,paste0('dT',1:(Tpoints-1)),drop=F],1,function(x){
# obstimesi<-c()
# for(i in 1:(length(x))){
# obstimesi[i]<-round(sum(x[1:i])*resolution,digits=0) / resolution
# }
# return(obstimesi)
# }
# ),
# byrow=T,ncol=(Tpoints-1))
# )
#
# #extract tdpred observations and match with rounded times
# TDpredMeans<-matrix(NA,nrow=max(obstimes*resolution),ncol=n.TDpred)
# for(predi in 1:n.TDpred){
# TDpredObs<-matrix(0,nrow=max(obstimes*resolution),ncol=nrow(data))
# TDpredObs[cbind(c(obstimes * resolution),c(row(obstimes)))] <- data[,paste0(TDpredNames[predi],'_T',0:(Tpoints-2))]
# TDpredMeans[,predi]<-apply(TDpredObs,1,mean,na.rm=T)
# }
# TDpredMeans=TDpredMeans[1:length(plottimes),,drop=FALSE]
#
# TDPREDEFFECT <- OpenMx::mxEval(TDPREDEFFECT,mxobj,compute=T)
# ARforResolution <- OpenMx::expm(DRIFT*(1/resolution))
# for(i in 2:length(plottimes)){
# tdpredeffect[i,]<- ARforResolution %*% t(tdpredeffect[i-1, ,drop=F]) + TDPREDEFFECT %*% t(TDpredMeans[i-1, ,drop=F])
# }
# }
means<-matrix(apply(cbind(1:length(plottimes)),1,function(x) {
dT<-plottimes[x]
OpenMx::expm(DRIFT*dT) %*% T0MEANS +
solve(DRIFT) %*% (OpenMx::expm(DRIFT*dT)-diag(nrow(DRIFT))) %*% CINT
}), byrow=T,ncol=nrow(DRIFT))
if(meansylim=='auto') meansylim<- c(min(means),max(means))
graphics::plot(plottimes, means[,1], type = "l", xlab = xlab, ylab = ylab,
main="Process means",
xlim=c(0,max.time), ylim=meansylim, lwd=2,col=colourvector[1])
if(n.latent > 1) {
for(i in 2:n.latent){
graphics::points(plottimes, means[,i], type = "l", lwd=2,col=colourvector[i])
graphics::legend("topright",legend=latentNames,text.col=colourvector,bty="n")
}
}
if(wait==TRUE){
message("Press [enter] to display next graph")
readline()
}
} #end mean plot
if(withinVariance==TRUE){ #plot within subject variance
DRIFTHATCH<-(DRIFT %x% diag(n.latent) + diag(n.latent) %x% DRIFT)
withinvar<-matrix(apply(plottimes,1,function(x) matrix(
OpenMx::expm(DRIFT %x% x) %*% T0VAR %*% t(OpenMx::expm(DRIFT %x% x)) + #initial variance
matrix(solve(DRIFTHATCH) %*% ((OpenMx::expm(DRIFTHATCH %x% x)) - diag(n.latent^2) ) %*% #diffusion process
OpenMx::rvectorize(DIFFUSION),nrow=n.latent),nrow=n.latent)[row(diag(n.latent))>=col(diag(n.latent))]),
byrow=T,ncol=length(diag(n.latent)[row(diag(n.latent))>=col(diag(n.latent))]))
colnames(withinvar)<-paramlabels[lower.tri(paramlabels,diag=T)]
colourvector <- grDevices::rainbow(length(DRIFT[upper.tri(DRIFT,diag=T)==T]),v=.8)
graphics::plot(plottimes, withinvar[,1], type = "l", xlab = xlab, ylab = ylab,
main="Within subject variance / covariance",
xlim=c(0,max.time), ylim=c(min(withinvar),max(withinvar)), lwd=2,col=colourvector[1])
if(n.latent > 1) {
for(i in 2:ncol(withinvar)){
graphics::points(plottimes, withinvar[,i], type = "l", lwd=2,col=colourvector[i])
graphics::legend("topright",legend=colnames(withinvar),text.col=colourvector,bty="n")
}
}
if(wait==TRUE){
message("Press [enter] to display next graph")
readline()
}
} #end within variance plots
# if(betweenVariance==TRUE){ #plot between subject variance
# if(is.null(ctfitobj$ctmodelobj$TRAITVAR) & n.TIpred < 1) message('No between subject variance modelled - between subject variance trajectories not possible')
# if(!is.null(ctfitobj$ctmodelobj$TRAITVAR) | n.TIpred > 0) {
# traitvariance<-0
# tipredvariance<-0
#
#
# if(!is.null(ctfitobj$ctmodelobj$TRAITVAR)) {
# T0TRAITEFFECT<-OpenMx::mxEval(T0TRAITEFFECT, ctfitobj$mxobj,compute=T)
#
# if(asymptotes==FALSE) traitvariance<-matrix(apply(plottimes,1,function(x) matrix(
# (OpenMx::expm(DRIFT %x% x) %*% ( T0TRAITEFFECT) + (solve(DRIFT) %*% (OpenMx::expm(DRIFT %x% x) - diag(n.latent)))) %*%
# TRAITVAR %*% t(
# (OpenMx::expm(DRIFT %x% x) %*% ( T0TRAITEFFECT) + (solve(DRIFT) %*% (OpenMx::expm(DRIFT %x% x) - diag(n.latent)))) )
# ,nrow=n.latent)[row(diag(n.latent))>=col(diag(n.latent))]),
# byrow=T,ncol=length(diag(n.latent)[row(diag(n.latent))>=col(diag(n.latent))]))
#
# if(asymptotes==TRUE) traitvariance<-matrix(apply(plottimes,1,function(x) matrix(
# (OpenMx::expm(DRIFT %x% x) %*% ( T0TRAITEFFECT) + #remaining t0 effect
# (diag(n.latent)-OpenMx::expm(DRIFT*x))) %&% asymTRAITVAR, #plus discrete interval effect
# ,nrow=n.latent)[row(diag(n.latent))>=col(diag(n.latent))]),
# byrow=T,ncol=length(diag(n.latent)[row(diag(n.latent))>=col(diag(n.latent))]))
#
#
# }
# if(n.TIpred > 0) {
# T0TIPREDEFFECT<-OpenMx::mxEval(T0TIPREDEFFECT, ctfitobj$mxobj,compute=T)
# if(asymptotes==FALSE) tipredvariance<-matrix(apply(plottimes,1,function(x) matrix(
# ( (OpenMx::expm(DRIFT %x% x)) %*% (T0TIPREDEFFECT ) + #T0 loading
# (solve(DRIFT) %*%(OpenMx::expm(DRIFT %x% x) - diag(n.latent)) %*% TIPREDEFFECT ) ) %*% #discrete loading
# TIPREDVAR %*% t( #variance of ti predictor
# ( (OpenMx::expm(DRIFT %x% x)) %*% (T0TIPREDEFFECT ) + #T0 loading
# (solve(DRIFT) %*%(OpenMx::expm(DRIFT %x% x) - diag(n.latent)) %*% TIPREDEFFECT ) ) )
#
# ,nrow=n.latent)[row(diag(n.latent))>=col(diag(n.latent))]),
# byrow=T,ncol=length(diag(n.latent)[row(diag(n.latent))>=col(diag(n.latent))]))
#
# if(asymptotes==TRUE) tipredvariance<-matrix(apply(plottimes,1,function(x) matrix(
# ( (OpenMx::expm(DRIFT %x% x)) %*% (T0TIPREDEFFECT ) + #T0 loading
# (solve(DRIFT) %*%(OpenMx::expm(DRIFT %x% x) - diag(n.latent)) %*% TIPREDEFFECT ) ) %&% TIPREDVAR, #discrete loading
# nrow=n.latent)[row(diag(n.latent))>=col(diag(n.latent))]),
# byrow=T,ncol=length(diag(n.latent)[row(diag(n.latent))>=col(diag(n.latent))]))
#
#
# }
# betweenvariance<-traitvariance + tipredvariance
#
# colnames(betweenvariance)<-indexMatrix(dimension=n.latent,symmetrical=TRUE,unique=T,
# upper=FALSE,sep='_',indices=F,namesvector=latentNames)
#
# colourvector <- grDevices::rainbow(length(DRIFT[upper.tri(DRIFT,diag=T)==T]),v=.8)
# graphics::plot(plottimes, betweenvariance[,1], type = "l", xlab = xlab, ylab = ylab,
# main="Between subject variance / covariance",
# xlim=c(0,max.time), ylim=c(min(betweenvariance),max(betweenvariance)), lwd=2,col=colourvector[1])
# if(n.latent > 1) {
# for(i in 2:ncol(withinvar)){
# graphics::points(plottimes, betweenvariance[,i], type = "l", lwd=2,col=colourvector[i])
# graphics::legend("topright",legend=colnames(betweenvariance),text.col=colourvector,bty="n")
# }
# }
#
#
# if(wait==TRUE){
# message("Press [enter] to display next graph")
# readline()
# }
# }}#end between variance plots
# 3. plot autoregressive/cross-lagged coefficients as function of time interval
ar<- matrix(apply(plottimes,1,function(x) c(diag(OpenMx::expm(DRIFT*x)))),ncol=nrow(DRIFT),byrow=T)
if(standardiseCR==TRUE) {
standardiser<-try(suppressWarnings(rep(sqrt(diag(abs(asymDIFFUSION))),each=n.latent) / rep(diag(sqrt(abs(asymDIFFUSION))),times=n.latent)))
if('try-error' %in% class(standardiser)) {
message('Unable to standardardise cross regression - plotting unstandardised.')
standardiseCR<-FALSE
}
standardiser[is.nan(standardiser)]<-0
}
if(standardiseCR==FALSE) standardiser<-1
driftindices<-(suppressWarnings(is.na(as.numeric(ctmodelobj$DRIFT)))==TRUE | ctmodelobj$DRIFT != 0)
cl<- matrix(apply(plottimes,1,function(x) {
c((OpenMx::expm(DRIFT*x)*standardiser)[row(DRIFT)!=col(DRIFT) & driftindices==TRUE])
}
),ncol=length(DRIFT[row(DRIFT)!=col(DRIFT) & driftindices==TRUE]),byrow=T)
arvars<-c(diag(paramlabels))
if(AR==TRUE){
if(all(ARylim == 'auto')){
ARylim <- c(0,1)
if(max(ar) > 1) ARylim[2] <- max(ar)
if(min(ar) < 0) ARylim[1] <- min(ar)
}
#AR coefficient plot
graphics::plot(plottimes, ar[,1], type = "l", xlab = xlab, ylab = ylab,
ylim=ARylim, xlim=c(0,max.time), main="Autoregression",lwd=2,col=2,...)
if(ncol(DRIFT)>1){ #if there is more than one AR parameter to plot
for( i in 2:ncol(ar)){
graphics::points(plottimes, ar[,i], type = "l", lwd=2,col=colourvector[i],...)
}
}
arnames <- paste0(arvars)
graphics::legend("topright",legend=arnames,text.col=colourvector,bty="n")
}
if(CR==TRUE && ncol(cl)>0){ #CL coefficients
if(wait==TRUE){
message("Press [enter] to display next graph")
readline()
}
if(standardiseCR==TRUE) CRtitle<-"Standardised crossregression"
if(standardiseCR==FALSE) CRtitle<-"Unstandardised crossregression"
colourvector <- grDevices::rainbow(ncol(cl),v=.8) #set plot colours
clvars<-paramlabels[row(DRIFT)!=col(DRIFT) & driftindices==TRUE]
# if(length(clvars[is.na(as.numeric(clvars))]) > 0) {#if there is 1 or more estimated cross effects
# clvars <- clvars[is.na(as.numeric(clvars))] #remove fixed params from clvars list
if(all(CRylim =='auto')){
CRylim <- c(-1,1)
if(max(cl) > 1) CRylim[2] <- max(cl)
if(min(cl) < -1) CRylim[1] <- min(cl)
}
graphics::plot(plottimes, cl[,1], type = "l", xlab = xlab, ylab = ylab,
ylim=CRylim, lwd=2,col=colourvector[1], main=CRtitle,...)
if(ncol(cl)>1){
for(i in 2:(ncol(cl))){
graphics::points(plottimes, cl[,i], type = "l", lwd=2,col=colourvector[i],...)
}
}
graphics::legend("topright",legend=paste0(clvars),text.col=colourvector,bty="n")
}
if(randomImpulse==TRUE && ncol(DRIFT)>1){ #randomCR coefficients
colourvector <- grDevices::rainbow(ncol(DRIFT),v=.8) #set plot colours
for( coli in 1:ncol(DRIFT)){
if(wait==TRUE){
message("Press [enter] to display next graph")
readline()
}
impulse<-rep(0,ncol(DRIFT))
impulse[coli]<-1
DIFFUSIONdiag<-diag(diag(DIFFUSION),n.latent)
diag(DIFFUSIONdiag)[diag(DIFFUSIONdiag)==0]<- .0001
DIFFUSIONdiag<-solve(sqrt(DIFFUSIONdiag)) %*% DIFFUSION %*% t(solve(sqrt(DIFFUSIONdiag)))
impulse<- DIFFUSIONdiag %*% impulse
impulseresponse<-lapply(plottimes,function(x) expm(DRIFT *x) %*% impulse)
IRylim <- c(-1,1)
if(max(unlist(impulseresponse)) > 1) IRylim[2] <- max(unlist(impulseresponse))
if(min(unlist(impulseresponse)) < -1) IRylim[1] <- min(unlist(impulseresponse))
graphics::plot(plottimes, unlist(lapply(impulseresponse,function(x) x[1] )), type = "l", xlab = xlab, ylab = ylab,
ylim=IRylim, lwd=2,col=colourvector[1], main=paste0('Expec. given obs. change on ',latentNames[coli]),...)
graphics::legend("topright",legend=paste0(latentNames),text.col=colourvector,bty="n")
if(n.latent>1) {
for(rowi in 2:nrow(DRIFT)){
graphics::points(plottimes, unlist(lapply(impulseresponse,function(x) x[rowi] )), type = "l",
lwd=2,col=colourvector[rowi],...)
}
}
}
}
if(experimentalImpulse==TRUE && ncol(DRIFT)>1){ #randomCR coefficients
colourvector <- grDevices::rainbow(ncol(DRIFT),v=.8) #set plot colours
for( coli in 1:ncol(DRIFT)){
if(wait==TRUE){
message("Press [enter] to display next graph")
readline()
}
impulse<-rep(0,ncol(DRIFT))
impulse[coli]<-1
impulseresponse<-lapply(plottimes,function(x) expm(DRIFT *x) %*% impulse)
IRylim <- c(-1,1)
if(max(unlist(impulseresponse)) > 1) IRylim[2] <- max(unlist(impulseresponse))
if(min(unlist(impulseresponse)) < -1) IRylim[1] <- min(unlist(impulseresponse))
graphics::plot(plottimes, unlist(lapply(impulseresponse,function(x) x[1] )), type = "l", xlab = xlab, ylab = ylab,
ylim=IRylim, lwd=2,col=colourvector[1], main=paste0('Expec. given input on ',latentNames[coli]),...)
graphics::legend("topright",legend=paste0(latentNames),text.col=colourvector,bty="n")
if(n.latent>1) {
for(rowi in 2:nrow(DRIFT)){
graphics::points(plottimes, unlist(lapply(impulseresponse,function(x) x[rowi] )), type = "l",
lwd=2,col=colourvector[rowi],...)
}
}
}
}
}
#' ctPlot
#'
#' Plots mean trajectories, autoregression, and crossregression plots, for ctsemFit objects.
#' More customizeable than basic plot.ctsemFit function.
#'
#' @param x ctsemFit object as generated by \code{\link{ctFit}}.
#' @param subject numeric. Specifies the subject (row of data from the mxobj) to plot for factorScores type plot.
#' @param plotType string. "mean" for expectation independent of any data,
#' "AR" for autoregressions, "CR" for cross regressions,
#' "standardiseCR" for standardised cross regressions (standardised based on estimated within subject variance),
#' "withinVar" for within variance and covariance,
#' "randomImpulse" for expected change in processes given a random fluctuation of +1 for each process
#' (so a mixture of DIFFUSION and DRIFT characteristics),
#' "experimentalImpulse" for expected change in processes given an exogenous input of +1 for each process, provides
#' alternate characterisation of autoregressive and cross regressive plots.
#' @param resolution Numeric. Plot points between each unit of time. Default of 'auto' adapts to xlim and results in 500 points in total.
#' @param impulseIndex Numeric. Only required for impulse plot types, specifies which column of the DRIFT matrix the impulse relates to.
#' @param xlim vector. As per usual for plot(), but xlim may not be negative.
#' @param typeVector Vector of plot types to use for plotting.
#' @param colVector vector of colours to use for plotting.
#' @param ltyVector Vector of line types to use for plotting.
#' @param ... Other options passed to \code{plot()}. ylim is required.
#' @return Character vector of labels from the DRIFT matrix in order plotted - useful for legends. Side-effect: plots graphs.
#' @examples
#' ## Examples set to 'donttest' because they take longer than 5s.
#'
#' ### example from Driver, Oud, Voelkle (2016),
#' ### simulated happiness and leisure time with unobserved heterogeneity.
#' \donttest{
#' data(ctExample1)
#' traitmodel <- ctModel(n.manifest=2, n.latent=2, Tpoints=6, LAMBDA=diag(2),
#' manifestNames=c('LeisureTime', 'Happiness'),
#' latentNames=c('LeisureTime', 'Happiness'), TRAITVAR="auto")
#' traitfit <- ctFit(dat=ctExample1, ctmodelobj=traitmodel)
#' ctPlot(traitfit, plotType='CR', xlim=c(0,5),ylim=c(-1,1))
#' }
#' @export
ctPlot<-function(x,plotType,xlim,resolution=50,impulseIndex=NULL,subject=1,typeVector='auto',colVector='auto',ltyVector='auto',...){
ctfitobj<-x
mxobj<-ctfitobj$mxobj
ctsummary<-summary(ctfitobj,verbose=TRUE)
#read in values
for(i in 1:length(ctsummary)){ #this loop reads in the specified continuous time model so the objects are available
assign(names(ctsummary[i]),eval(parse(text = paste0("ctsummary","$",names(ctsummary[i])))))
}
asymptotes<-ctfitobj$ctfitargs$asymptotes
paramlabels<-paste0(rep(ctfitobj$ctmodelobj$latentNames,ctfitobj$ctmodelobj$n.latent),'_',
rep(ctfitobj$ctmodelobj$latentNames,each=ctfitobj$ctmodelobj$n.latent))
latentNames<-ctfitobj$ctmodelobj$latentNames
n.latent<-ctfitobj$ctmodelobj$n.latent
n.manifest<-ctfitobj$ctmodelobj$n.manifest
n.TIpred<-ctfitobj$ctmodelobj$n.TIpred
n.TDpred<-ctfitobj$ctmodelobj$n.TDpred
manifestNames <- ctfitobj$ctmodelobj$manifestNames
TDpredNames <- ctfitobj$ctmodelobj$TDpredNames
TIpredNames <- ctfitobj$ctmodelobj$TIpredNames
Tpoints<-ctfitobj$ctmodelobj$Tpoints
stationary<-ctfitobj$ctfitargs$stationary
T0MEANS<-OpenMx::mxEval(T0MEANS, ctfitobj$mxobj,compute=F)
T0VAR<-OpenMx::mxEval(T0VAR, ctfitobj$mxobj,compute=F)
if(resolution=='auto') resolution <- ceiling(500/xlim[2])
plottimes<-matrix(seq(0,xlim[2],1/resolution)[-1],ncol=1) #time steps
if(plotType=='mean'){# 2. plot mean trend
#TD predictors
tdpredeffect<-matrix(0,nrow=length(plottimes),ncol=n.latent) #default 0 matrix
#extract and round absolute times to resolution
# if(n.TDpred >0){
# if(ctfitobj$ctfitargs$objective!='Kalman') data<-mxobj$data$observed
# if(ctfitobj$ctfitargs$objective=='Kalman') {
# data<- ctLongToWide(mxobj$data$observed,id='id',time='dT1',
# manifestNames=manifestNames,TDpredNames=TDpredNames,TIpredNames=TIpredNames)
# data<-data[,-which(colnames(data)=='T0'),drop=F]
# colnames(data)[which(colnames(data) %in% paste0('T',1:(Tpoints-1)))]<-paste0('dT',1:(Tpoints-1))
# }
#
# obstimes<-cbind(
# matrix(rep(0,nrow(data)),ncol=1),
# matrix(
# apply(data[,paste0('dT',1:(Tpoints-1)),drop=F],1,function(x){
# obstimesi<-c()
# for(i in 1:(length(x))){
# obstimesi[i]<-round(sum(x[1:i])*resolution,digits=0) / resolution
# }
# return(obstimesi)
# }
# ),
# byrow=T,ncol=(Tpoints-1))
# )
# # obstimes=obstimes[obstimes <=xlim]
#
# #extract tdpred observations and match with rounded times
# TDpredMeans<-matrix(NA,nrow=max(obstimes*resolution),ncol=n.TDpred)
# for(predi in 1:n.TDpred){
# TDpredObs<-matrix(0,nrow=max(obstimes*resolution),ncol=nrow(data))
# TDpredObs[cbind(c(obstimes * resolution),c(row(obstimes)))] <- data[,paste0(TDpredNames[predi],'_T',0:(Tpoints-2))]
# TDpredMeans[,predi]<-apply(TDpredObs,1,mean,na.rm=T)
# }
# TDpredMeans=TDpredMeans[1:length(plottimes),,drop=FALSE]
#
# TDPREDEFFECT <- OpenMx::mxEval(TDPREDEFFECT,mxobj,compute=T)
# ARforResolution <- OpenMx::expm(DRIFT*(1/resolution))
# for(i in 2:length(plottimes)){
# tdpredeffect[i,]<- ARforResolution %*% t(tdpredeffect[i-1, ,drop=F]) + TDPREDEFFECT %*% t(TDpredMeans[i-1, ,drop=F])
# }
# }
means<-matrix(apply(cbind(1:length(plottimes)),1,function(x) {
dT<-plottimes[x]
OpenMx::expm(DRIFT*dT) %*% T0MEANS +
solve(DRIFT) %*% (OpenMx::expm(DRIFT*dT)-diag(nrow(DRIFT))) %*% CINT
}), byrow=T,ncol=nrow(DRIFT))
plotdata<-means
legend<-latentNames
} #end mean plot
if(plotType=='withinVariance'){ #plot within subject variance
DRIFTHATCH<-(DRIFT %x% diag(n.latent) + diag(n.latent) %x% DRIFT)
withinvar<-matrix(apply(plottimes,1,function(x) matrix(
OpenMx::expm(DRIFT %x% x) %*% T0VAR %*% t(OpenMx::expm(DRIFT %x% x)) + #initial variance
matrix(solve(DRIFTHATCH) %*% ((OpenMx::expm(DRIFTHATCH %x% x)) - diag(n.latent^2) ) %*% #diffusion process
OpenMx::rvectorize(DIFFUSION),nrow=n.latent),nrow=n.latent)[row(diag(n.latent))>=col(diag(n.latent))]),
byrow=T,ncol=length(diag(n.latent)[row(diag(n.latent))>=col(diag(n.latent))]))
colnames(withinvar)<-paramlabels[lower.tri(paramlabels,diag=T)]
plotdata<-withinvar
legend<-colnames(withinvar)
} #end within variance plots
if(plotType=='AR'){
ar<- matrix(apply(plottimes,1,function(x) c(diag(OpenMx::expm(DRIFT*x)))),ncol=nrow(DRIFT),byrow=T)
plotdata<-ar
legend<-paste0(c(diag(paramlabels)))
}
if((plotType=='CR' | plotType=='standardiseCR') & ncol(DRIFT)<2) stop('No cross effects to plot!')
if((plotType=='CR' | plotType=='standardiseCR') & ncol(DRIFT)>1){
if(plotType=='standardiseCR') {
standardiser<-try(suppressWarnings(rep(sqrt(diag(abs(asymDIFFUSION))),each=n.latent) / rep(diag(sqrt(abs(asymDIFFUSION))),times=n.latent)))
if('try-error' %in% class(standardiser)) stop('Unable to standardardise cross regression - try unstandardised.')
standardiser[is.nan(standardiser)]<-0
}
if(plotType!='standardiseCR') standardiser<-1
if(ctfitobj$ctfitargs$transformedParams==TRUE) driftindices<-(mxobj$DRIFT$free==TRUE | mxobj$DRIFT$values != 0)
if(ctfitobj$ctfitargs$transformedParams==FALSE) driftindices<-(mxobj$DRIFT$free==TRUE | mxobj$DRIFT$values != 0)
cl<- matrix(apply(plottimes,1,function(x) {
c((OpenMx::expm(DRIFT*x)*standardiser)[row(DRIFT)!=col(DRIFT) & driftindices==TRUE])
}
),ncol=length(DRIFT[row(DRIFT)!=col(DRIFT) & driftindices==TRUE]),byrow=T)
clvars<-paramlabels[row(DRIFT)!=col(DRIFT) & driftindices==TRUE]
plotdata<-cl
legend<-paste0(clvars)
}
if(plotType=='randomImpulse'){
# if(is.null(impulseIndex)) stop('impulseIndex must be specified - integer from 1 to nrow(DRIFT)')
# if(as.integer(impulseIndex) != impulseIndex | impulseIndex < n.latent | impulseIndex > n.latent) stop('impulseIndex must be an integer between 1 and n.latent')
impulse<-rep(0,ncol(DRIFT))
impulse[impulseIndex]<-1
DIFFUSIONdiag<-diag(diag(DIFFUSION),n.latent)
diag(DIFFUSIONdiag)[diag(DIFFUSIONdiag)<=0]<-.0001
# if(standardiseCR==TRUE)
DIFFUSIONdiag<-solve(sqrt(DIFFUSIONdiag)) %*% DIFFUSION %*% t(solve(sqrt(DIFFUSIONdiag)))
impulse<- DIFFUSIONdiag %*% impulse
impulseresponse<-lapply(plottimes,function(x) expm(DRIFT *x) %*% impulse)
IRylim <- c(-1,1)
if(max(unlist(impulseresponse)) > 1) IRylim[2] <- max(unlist(impulseresponse))
if(min(unlist(impulseresponse)) < -1) IRylim[1] <- min(unlist(impulseresponse))
plotdata<-matrix(unlist(lapply(1:length(impulseresponse[[1]]),function(z) unlist(lapply(impulseresponse,function(x) x[z] )))),ncol=nrow(DRIFT))
legend<-paste0(latentNames)
}
if(plotType=='experimentalImpulse'){
if(is.null(impulseIndex)) stop('impulseIndex must be specified - integer from 1 to nrow(DRIFT)')
impulse<-rep(0,ncol(DRIFT))
impulse[impulseIndex]<-1
impulseresponse<-lapply(plottimes,function(x) expm(DRIFT *x) %*% impulse)
IRylim <- c(-1,1)
if(max(unlist(impulseresponse)) > 1) IRylim[2] <- max(unlist(impulseresponse))
if(min(unlist(impulseresponse)) < -1) IRylim[1] <- min(unlist(impulseresponse))
plotdata<-matrix(unlist(lapply(1:length(impulseresponse[[1]]),function(z) unlist(lapply(impulseresponse,function(x) x[z] )))),ncol=nrow(DRIFT))
legend<-paste0(latentNames)
}
if(plotType=='factorScores'){
newdata<-mxobj@data@observed[subject,,drop=F]
smallmxobj<-mxModel(mxobj,mxData(newdata,type='raw'))
plottimes <- newdata[colnames(newdata) %in% paste0('dT',1:(Tpoints-1))]
for(i in 2:length(plottimes)){
plottimes[i]<-plottimes[i-1]+plottimes[i]
}
plottimes<-c(0,plottimes)
newdata<-newdata[1:(n.manifest*Tpoints)]
scores<-mxFactorScores(smallmxobj,type='Regression')
amat <- mxEval(A, mxobj, defvar.row=subject, compute=TRUE)
mmat <- mxEval(M, mxobj, defvar.row=subject, compute=TRUE)
dataRow1 <- testd[1,]
fullscores <- solve(diag(1, n.latent*Tpoints) - amat[1:(n.latent*Tpoints), 1:(n.latent*Tpoints)]) %*% (scores[1,,1] - matrix(mmat)[1:(n.latent*Tpoints),])
plotdata<-cbind(matrix(newdata,ncol=n.manifest,byrow=T),matrix(fullscores,ncol=n.latent,byrow=T))
errorbars<-matrix(scores[,,2],ncol=n.latent,byrow=T)
if(ltyVector[1]=='auto') ltyVector <- c(rep(1,n.manifest),rep(2,n.latent))
if(typeVector[1]=='auto') typeVector <- rep('b',ncol(plotdata))
}
if(colVector[1]=='auto') colVector <- grDevices::rainbow(ncol(plotdata),v=.8) #set plot colours
if(ltyVector[1]=='auto') ltyVector <- rep(1,ncol(plotdata))
if(typeVector[1]=='auto') typeVector <- rep('l',ncol(plotdata))
graphics::plot(plottimes,rep(NA,length(plottimes)),...)
graphics::grid()
graphics::points(plottimes,plotdata[,1],type=typeVector[1], col=colVector[1],lty=ltyVector[1],...)
if(ncol(plotdata) > 1) {
for(rowi in 2:ncol(plotdata)){
graphics::points(plottimes, plotdata[,rowi],type=typeVector[rowi], col=colVector[rowi],lty=ltyVector[rowi],...)
if(plotType=='factorScores' & rowi > n.manifest){
graphics::arrows(plottimes, plotdata[,rowi] - errorbars[,rowi-n.manifest], plottimes, plotdata[,rowi] + errorbars[,rowi-n.manifest],
col=colVector[rowi],length=0.05, angle=90, code=3)
}
}
}
return(legend)
}
Try the ctsemOMX package in your browser
Any scripts or data that you put into this service are public.
ctsemOMX documentation built on Oct. 5, 2023, 5:06 p.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 ctPlot plot.ctsemFit plot.ctsemMultigroupFit
Documented in ctPlot plot.ctsemFit plot.ctsemMultigroupFit
#' Plot function for ctsemMultigroupFit object
#'
#' Plots \code{\link{ctMultigroupFit}} objects.
#'
#' @param x ctsemMultigroupFit object as generated by \code{\link{ctMultigroupFit}}
#' @param group character string of subgroup to plot. Default of 'show chooser' displays list and lets you select.
#' @param ... additional parameters to pass to \code{\link{plot.ctsemFit}} function.
#' @return Nothing. Side-effect: plots graphs.
#' @method plot ctsemMultigroupFit
#' @export
plot.ctsemMultigroupFit<-function(x,group='show chooser',...){
if(group=='show chooser') {
message(paste0(1:length(x$groups), ' ', x$groups, '\n'))
selection<-as.numeric(readline(
'Enter the number of the group you wish to plot: '))
selection<-x$groups[[selection]]
}
if(group != 'show chooser') selection <- group
tempctobj<-list(ctfitargs=x$ctfitargs, ctmodelobj=x$ctmodelobj, mxobj=x$mxobj[[selection]])
class(tempctobj)<-'ctsemFit'
plot.ctsemFit(tempctobj,...)
}
#' Plotting function for object class ctsemFit
#'
#' Ouputs mean trajectories, autoregression, and crossregression plots.
#' For more customization possibilities, see \code{\link{ctPlot}}.
#'
#' @param x ctsemFit object as generated by \code{\link{ctFit}}.
#' @param resolution Numeric. Plot points between each unit of time. Default of 'auto' adapts to max.time and results in 500 in total.
#' @param wait If true, user is prompted to continue before plotting next graph. If false, graphs are plotted one after another without waiting.
#' @param max.time Time scale on which to plot parameters. If auto, parameters are plotted for full range of observed variables.
#' @param mean if TRUE, plot of means from 0 to max.time included in output.
#' @param withinVariance if TRUE, plot within subject variance / covariance.
#' @param AR if TRUE, plot of autoregressive values from 0 to max.time included in output.
#' @param CR if TRUE, plot of cross regressive values from 0 to max.time included in output.
#' @param standardiseCR if TRUE , cross regression values are standardised based on estimated within subject variance.
#' @param randomImpulse if TRUE (default), plots expected change in processes given a random fluctuation of +1 for each process --
#' plot is then a mixture of DIFFUSION and DRIFT characteristics.
#' @param experimentalImpulse if TRUE (default), plots expected change in processes given an exogenous input of +1 for each process --
#' alternate characterisation of autoregressive and cross regressive plots.
#' @param xlab X axis label.
#' @param ylab Y axis label.
#' @param meansylim Vector of min and max limits for mean trajectory plot. 'auto' calculates automatically.
#' @param ARylim Vector of min and max limits for autoregression plot. 'auto' is c(0,1), and expands if necessary.
#' @param CRylim Vector of min and max limits for cross regression plot. 'auto' is c(-1,1), and expands if necessary.
#' @param ... Other options passed to \code{plot()}.
#' @return Nothing. Side-effect: plots graphs.
#' @method plot ctsemFit
#' @examples
#' ## Examples set to 'donttest' because they take longer than 5s.
#'
#' ### example from Driver, Oud, Voelkle (2015),
#' ### simulated happiness and leisure time with unobserved heterogeneity.
#' \donttest{
#' data(ctExample1)
#' traitmodel <- ctModel(n.manifest=2, n.latent=2, Tpoints=6, LAMBDA=diag(2),
#' manifestNames=c('LeisureTime', 'Happiness'),
#' latentNames=c('LeisureTime', 'Happiness'), TRAITVAR="auto")
#' traitfit <- ctFit(dat=ctExample1, ctmodelobj=traitmodel)
#' plot(traitfit, wait=FALSE)
#' }
#' @export
plot.ctsemFit<-function(x,resolution=50,wait=TRUE,max.time="auto",mean=TRUE,
withinVariance=TRUE,AR=TRUE,CR=TRUE,standardiseCR=FALSE,randomImpulse=FALSE,
experimentalImpulse=FALSE,
xlab="Time",
meansylim='auto',ARylim='auto', CRylim='auto', ylab="Value",...){
message("Plotting fit")
ctfitobj<-x
mxobj<-ctfitobj$mxobj
ctmodelobj<-ctfitobj$ctmodelobj
ctsummary<-summary(ctfitobj,verbose=TRUE)
#read in values
for(i in 1:length(ctsummary)){ #this loop reads in the specified continuous time model so the objects are available
assign(names(ctsummary[i]),eval(parse(text = paste0("ctsummary","$",names(ctsummary[i])))))
}
asymptotes<-ctfitobj$ctfitargs$asymptotes
latentNames<-ctfitobj$ctmodelobj$latentNames
n.latent<-ctfitobj$ctmodelobj$n.latent
n.manifest<-ctfitobj$ctmodelobj$n.manifest
n.TIpred<-ctfitobj$ctmodelobj$n.TIpred
n.TDpred<-ctfitobj$ctmodelobj$n.TDpred
manifestNames <- ctfitobj$ctmodelobj$manifestNames
TDpredNames <- ctfitobj$ctmodelobj$TDpredNames
TIpredNames <- ctfitobj$ctmodelobj$TIpredNames
Tpoints<-ctfitobj$ctmodelobj$Tpoints
stationary<-ctfitobj$ctfitargs$stationary
T0MEANS<-OpenMx::mxEval(T0MEANS, ctfitobj$mxobj,compute=F)
T0VAR<-OpenMx::mxEval(T0VAR, ctfitobj$mxobj,compute=F)
paramlabels<-matrix(c(paste0(rep(latentNames,n.latent),'_',rep(latentNames,each=n.latent))),n.latent,n.latent)
if(max.time=='auto' & ctfitobj$ctfitargs$objective!='cov'){
if(max.time=="auto" & ctfitobj$ctfitargs$objective!='Kalman') max.time <- max(rowSums(as.matrix(mxobj$data$observed[,paste0('dT',1:(Tpoints-1)),drop=FALSE]),na.rm=T)) # max time of plot
if(max.time=="auto" & ctfitobj$ctfitargs$objective=='Kalman') max.time <- sum(mxobj$data$observed[mxobj$data$observed[,'id']==1,'dT1'],na.rm=T) # max time of plot
}
if(max.time=='auto' & ctfitobj$ctfitargs$objective=='cov') stop('max.time argument must be set when plotting covariance based data')
if(resolution=='auto') resolution <- ceiling(500/max.time)
colourvector <- grDevices::rainbow(ncol(DRIFT),v=.8) #set plot colours
plottimes<-matrix(seq(0,max.time,1/resolution)[-1],ncol=1) #time steps
if(mean==TRUE){# 2. plot mean trend
#TD predictors
tdpredeffect<-matrix(0,nrow=length(plottimes),ncol=n.latent) #default 0 matrix
#extract and round absolute times to resolution
# if(n.TDpred >0){
# if(ctfitobj$ctfitargs$objective!='Kalman') data<-mxobj$data$observed
# if(ctfitobj$ctfitargs$objective=='Kalman') {
# data<- ctLongToWide(mxobj$data$observed,id='id',time='dT1',
# manifestNames=manifestNames,TDpredNames=TDpredNames,TIpredNames=TIpredNames)
# data<-data[,-which(colnames(data)=='T0'),drop=F]
# colnames(data)[which(colnames(data) %in% paste0('T',1:(Tpoints-1)))]<-paste0('dT',1:(Tpoints-1))
# }
#
# obstimes<-cbind(
# matrix(rep(0,nrow(data)),ncol=1),
# matrix(
# apply(data[,paste0('dT',1:(Tpoints-1)),drop=F],1,function(x){
# obstimesi<-c()
# for(i in 1:(length(x))){
# obstimesi[i]<-round(sum(x[1:i])*resolution,digits=0) / resolution
# }
# return(obstimesi)
# }
# ),
# byrow=T,ncol=(Tpoints-1))
# )
#
# #extract tdpred observations and match with rounded times
# TDpredMeans<-matrix(NA,nrow=max(obstimes*resolution),ncol=n.TDpred)
# for(predi in 1:n.TDpred){
# TDpredObs<-matrix(0,nrow=max(obstimes*resolution),ncol=nrow(data))
# TDpredObs[cbind(c(obstimes * resolution),c(row(obstimes)))] <- data[,paste0(TDpredNames[predi],'_T',0:(Tpoints-2))]
# TDpredMeans[,predi]<-apply(TDpredObs,1,mean,na.rm=T)
# }
# TDpredMeans=TDpredMeans[1:length(plottimes),,drop=FALSE]
#
# TDPREDEFFECT <- OpenMx::mxEval(TDPREDEFFECT,mxobj,compute=T)
# ARforResolution <- OpenMx::expm(DRIFT*(1/resolution))
# for(i in 2:length(plottimes)){
# tdpredeffect[i,]<- ARforResolution %*% t(tdpredeffect[i-1, ,drop=F]) + TDPREDEFFECT %*% t(TDpredMeans[i-1, ,drop=F])
# }
# }
means<-matrix(apply(cbind(1:length(plottimes)),1,function(x) {
dT<-plottimes[x]
OpenMx::expm(DRIFT*dT) %*% T0MEANS +
solve(DRIFT) %*% (OpenMx::expm(DRIFT*dT)-diag(nrow(DRIFT))) %*% CINT
}), byrow=T,ncol=nrow(DRIFT))
if(meansylim=='auto') meansylim<- c(min(means),max(means))
graphics::plot(plottimes, means[,1], type = "l", xlab = xlab, ylab = ylab,
main="Process means",
xlim=c(0,max.time), ylim=meansylim, lwd=2,col=colourvector[1])
if(n.latent > 1) {
for(i in 2:n.latent){
graphics::points(plottimes, means[,i], type = "l", lwd=2,col=colourvector[i])
graphics::legend("topright",legend=latentNames,text.col=colourvector,bty="n")
}
}
if(wait==TRUE){
message("Press [enter] to display next graph")
readline()
}
} #end mean plot
if(withinVariance==TRUE){ #plot within subject variance
DRIFTHATCH<-(DRIFT %x% diag(n.latent) + diag(n.latent) %x% DRIFT)
withinvar<-matrix(apply(plottimes,1,function(x) matrix(
OpenMx::expm(DRIFT %x% x) %*% T0VAR %*% t(OpenMx::expm(DRIFT %x% x)) + #initial variance
matrix(solve(DRIFTHATCH) %*% ((OpenMx::expm(DRIFTHATCH %x% x)) - diag(n.latent^2) ) %*% #diffusion process
OpenMx::rvectorize(DIFFUSION),nrow=n.latent),nrow=n.latent)[row(diag(n.latent))>=col(diag(n.latent))]),
byrow=T,ncol=length(diag(n.latent)[row(diag(n.latent))>=col(diag(n.latent))]))
colnames(withinvar)<-paramlabels[lower.tri(paramlabels,diag=T)]
colourvector <- grDevices::rainbow(length(DRIFT[upper.tri(DRIFT,diag=T)==T]),v=.8)
graphics::plot(plottimes, withinvar[,1], type = "l", xlab = xlab, ylab = ylab,
main="Within subject variance / covariance",
xlim=c(0,max.time), ylim=c(min(withinvar),max(withinvar)), lwd=2,col=colourvector[1])
if(n.latent > 1) {
for(i in 2:ncol(withinvar)){
graphics::points(plottimes, withinvar[,i], type = "l", lwd=2,col=colourvector[i])
graphics::legend("topright",legend=colnames(withinvar),text.col=colourvector,bty="n")
}
}
if(wait==TRUE){
message("Press [enter] to display next graph")
readline()
}
} #end within variance plots
# if(betweenVariance==TRUE){ #plot between subject variance
# if(is.null(ctfitobj$ctmodelobj$TRAITVAR) & n.TIpred < 1) message('No between subject variance modelled - between subject variance trajectories not possible')
# if(!is.null(ctfitobj$ctmodelobj$TRAITVAR) | n.TIpred > 0) {
# traitvariance<-0
# tipredvariance<-0
#
#
# if(!is.null(ctfitobj$ctmodelobj$TRAITVAR)) {
# T0TRAITEFFECT<-OpenMx::mxEval(T0TRAITEFFECT, ctfitobj$mxobj,compute=T)
#
# if(asymptotes==FALSE) traitvariance<-matrix(apply(plottimes,1,function(x) matrix(
# (OpenMx::expm(DRIFT %x% x) %*% ( T0TRAITEFFECT) + (solve(DRIFT) %*% (OpenMx::expm(DRIFT %x% x) - diag(n.latent)))) %*%
# TRAITVAR %*% t(
# (OpenMx::expm(DRIFT %x% x) %*% ( T0TRAITEFFECT) + (solve(DRIFT) %*% (OpenMx::expm(DRIFT %x% x) - diag(n.latent)))) )
# ,nrow=n.latent)[row(diag(n.latent))>=col(diag(n.latent))]),
# byrow=T,ncol=length(diag(n.latent)[row(diag(n.latent))>=col(diag(n.latent))]))
#
# if(asymptotes==TRUE) traitvariance<-matrix(apply(plottimes,1,function(x) matrix(
# (OpenMx::expm(DRIFT %x% x) %*% ( T0TRAITEFFECT) + #remaining t0 effect
# (diag(n.latent)-OpenMx::expm(DRIFT*x))) %&% asymTRAITVAR, #plus discrete interval effect
# ,nrow=n.latent)[row(diag(n.latent))>=col(diag(n.latent))]),
# byrow=T,ncol=length(diag(n.latent)[row(diag(n.latent))>=col(diag(n.latent))]))
#
#
# }
# if(n.TIpred > 0) {
# T0TIPREDEFFECT<-OpenMx::mxEval(T0TIPREDEFFECT, ctfitobj$mxobj,compute=T)
# if(asymptotes==FALSE) tipredvariance<-matrix(apply(plottimes,1,function(x) matrix(
# ( (OpenMx::expm(DRIFT %x% x)) %*% (T0TIPREDEFFECT ) + #T0 loading
# (solve(DRIFT) %*%(OpenMx::expm(DRIFT %x% x) - diag(n.latent)) %*% TIPREDEFFECT ) ) %*% #discrete loading
# TIPREDVAR %*% t( #variance of ti predictor
# ( (OpenMx::expm(DRIFT %x% x)) %*% (T0TIPREDEFFECT ) + #T0 loading
# (solve(DRIFT) %*%(OpenMx::expm(DRIFT %x% x) - diag(n.latent)) %*% TIPREDEFFECT ) ) )
#
# ,nrow=n.latent)[row(diag(n.latent))>=col(diag(n.latent))]),
# byrow=T,ncol=length(diag(n.latent)[row(diag(n.latent))>=col(diag(n.latent))]))
#
# if(asymptotes==TRUE) tipredvariance<-matrix(apply(plottimes,1,function(x) matrix(
# ( (OpenMx::expm(DRIFT %x% x)) %*% (T0TIPREDEFFECT ) + #T0 loading
# (solve(DRIFT) %*%(OpenMx::expm(DRIFT %x% x) - diag(n.latent)) %*% TIPREDEFFECT ) ) %&% TIPREDVAR, #discrete loading
# nrow=n.latent)[row(diag(n.latent))>=col(diag(n.latent))]),
# byrow=T,ncol=length(diag(n.latent)[row(diag(n.latent))>=col(diag(n.latent))]))
#
#
# }
# betweenvariance<-traitvariance + tipredvariance
#
# colnames(betweenvariance)<-indexMatrix(dimension=n.latent,symmetrical=TRUE,unique=T,
# upper=FALSE,sep='_',indices=F,namesvector=latentNames)
#
# colourvector <- grDevices::rainbow(length(DRIFT[upper.tri(DRIFT,diag=T)==T]),v=.8)
# graphics::plot(plottimes, betweenvariance[,1], type = "l", xlab = xlab, ylab = ylab,
# main="Between subject variance / covariance",
# xlim=c(0,max.time), ylim=c(min(betweenvariance),max(betweenvariance)), lwd=2,col=colourvector[1])
# if(n.latent > 1) {
# for(i in 2:ncol(withinvar)){
# graphics::points(plottimes, betweenvariance[,i], type = "l", lwd=2,col=colourvector[i])
# graphics::legend("topright",legend=colnames(betweenvariance),text.col=colourvector,bty="n")
# }
# }
#
#
# if(wait==TRUE){
# message("Press [enter] to display next graph")
# readline()
# }
# }}#end between variance plots
# 3. plot autoregressive/cross-lagged coefficients as function of time interval
ar<- matrix(apply(plottimes,1,function(x) c(diag(OpenMx::expm(DRIFT*x)))),ncol=nrow(DRIFT),byrow=T)
if(standardiseCR==TRUE) {
standardiser<-try(suppressWarnings(rep(sqrt(diag(abs(asymDIFFUSION))),each=n.latent) / rep(diag(sqrt(abs(asymDIFFUSION))),times=n.latent)))
if('try-error' %in% class(standardiser)) {
message('Unable to standardardise cross regression - plotting unstandardised.')
standardiseCR<-FALSE
}
standardiser[is.nan(standardiser)]<-0
}
if(standardiseCR==FALSE) standardiser<-1
driftindices<-(suppressWarnings(is.na(as.numeric(ctmodelobj$DRIFT)))==TRUE | ctmodelobj$DRIFT != 0)
cl<- matrix(apply(plottimes,1,function(x) {
c((OpenMx::expm(DRIFT*x)*standardiser)[row(DRIFT)!=col(DRIFT) & driftindices==TRUE])
}
),ncol=length(DRIFT[row(DRIFT)!=col(DRIFT) & driftindices==TRUE]),byrow=T)
arvars<-c(diag(paramlabels))
if(AR==TRUE){
if(all(ARylim == 'auto')){
ARylim <- c(0,1)
if(max(ar) > 1) ARylim[2] <- max(ar)
if(min(ar) < 0) ARylim[1] <- min(ar)
}
#AR coefficient plot
graphics::plot(plottimes, ar[,1], type = "l", xlab = xlab, ylab = ylab,
ylim=ARylim, xlim=c(0,max.time), main="Autoregression",lwd=2,col=2,...)
if(ncol(DRIFT)>1){ #if there is more than one AR parameter to plot
for( i in 2:ncol(ar)){
graphics::points(plottimes, ar[,i], type = "l", lwd=2,col=colourvector[i],...)
}
}
arnames <- paste0(arvars)
graphics::legend("topright",legend=arnames,text.col=colourvector,bty="n")
}
if(CR==TRUE && ncol(cl)>0){ #CL coefficients
if(wait==TRUE){
message("Press [enter] to display next graph")
readline()
}
if(standardiseCR==TRUE) CRtitle<-"Standardised crossregression"
if(standardiseCR==FALSE) CRtitle<-"Unstandardised crossregression"
colourvector <- grDevices::rainbow(ncol(cl),v=.8) #set plot colours
clvars<-paramlabels[row(DRIFT)!=col(DRIFT) & driftindices==TRUE]
# if(length(clvars[is.na(as.numeric(clvars))]) > 0) {#if there is 1 or more estimated cross effects
# clvars <- clvars[is.na(as.numeric(clvars))] #remove fixed params from clvars list
if(all(CRylim =='auto')){
CRylim <- c(-1,1)
if(max(cl) > 1) CRylim[2] <- max(cl)
if(min(cl) < -1) CRylim[1] <- min(cl)
}
graphics::plot(plottimes, cl[,1], type = "l", xlab = xlab, ylab = ylab,
ylim=CRylim, lwd=2,col=colourvector[1], main=CRtitle,...)
if(ncol(cl)>1){
for(i in 2:(ncol(cl))){
graphics::points(plottimes, cl[,i], type = "l", lwd=2,col=colourvector[i],...)
}
}
graphics::legend("topright",legend=paste0(clvars),text.col=colourvector,bty="n")
}
if(randomImpulse==TRUE && ncol(DRIFT)>1){ #randomCR coefficients
colourvector <- grDevices::rainbow(ncol(DRIFT),v=.8) #set plot colours
for( coli in 1:ncol(DRIFT)){
if(wait==TRUE){
message("Press [enter] to display next graph")
readline()
}
impulse<-rep(0,ncol(DRIFT))
impulse[coli]<-1
DIFFUSIONdiag<-diag(diag(DIFFUSION),n.latent)
diag(DIFFUSIONdiag)[diag(DIFFUSIONdiag)==0]<- .0001
DIFFUSIONdiag<-solve(sqrt(DIFFUSIONdiag)) %*% DIFFUSION %*% t(solve(sqrt(DIFFUSIONdiag)))
impulse<- DIFFUSIONdiag %*% impulse
impulseresponse<-lapply(plottimes,function(x) expm(DRIFT *x) %*% impulse)
IRylim <- c(-1,1)
if(max(unlist(impulseresponse)) > 1) IRylim[2] <- max(unlist(impulseresponse))
if(min(unlist(impulseresponse)) < -1) IRylim[1] <- min(unlist(impulseresponse))
graphics::plot(plottimes, unlist(lapply(impulseresponse,function(x) x[1] )), type = "l", xlab = xlab, ylab = ylab,
ylim=IRylim, lwd=2,col=colourvector[1], main=paste0('Expec. given obs. change on ',latentNames[coli]),...)
graphics::legend("topright",legend=paste0(latentNames),text.col=colourvector,bty="n")
if(n.latent>1) {
for(rowi in 2:nrow(DRIFT)){
graphics::points(plottimes, unlist(lapply(impulseresponse,function(x) x[rowi] )), type = "l",
lwd=2,col=colourvector[rowi],...)
}
}
}
}
if(experimentalImpulse==TRUE && ncol(DRIFT)>1){ #randomCR coefficients
colourvector <- grDevices::rainbow(ncol(DRIFT),v=.8) #set plot colours
for( coli in 1:ncol(DRIFT)){
if(wait==TRUE){
message("Press [enter] to display next graph")
readline()
}
impulse<-rep(0,ncol(DRIFT))
impulse[coli]<-1
impulseresponse<-lapply(plottimes,function(x) expm(DRIFT *x) %*% impulse)
IRylim <- c(-1,1)
if(max(unlist(impulseresponse)) > 1) IRylim[2] <- max(unlist(impulseresponse))
if(min(unlist(impulseresponse)) < -1) IRylim[1] <- min(unlist(impulseresponse))
graphics::plot(plottimes, unlist(lapply(impulseresponse,function(x) x[1] )), type = "l", xlab = xlab, ylab = ylab,
ylim=IRylim, lwd=2,col=colourvector[1], main=paste0('Expec. given input on ',latentNames[coli]),...)
graphics::legend("topright",legend=paste0(latentNames),text.col=colourvector,bty="n")
if(n.latent>1) {
for(rowi in 2:nrow(DRIFT)){
graphics::points(plottimes, unlist(lapply(impulseresponse,function(x) x[rowi] )), type = "l",
lwd=2,col=colourvector[rowi],...)
}
}
}
}
}
#' ctPlot
#'
#' Plots mean trajectories, autoregression, and crossregression plots, for ctsemFit objects.
#' More customizeable than basic plot.ctsemFit function.
#'
#' @param x ctsemFit object as generated by \code{\link{ctFit}}.
#' @param subject numeric. Specifies the subject (row of data from the mxobj) to plot for factorScores type plot.
#' @param plotType string. "mean" for expectation independent of any data,
#' "AR" for autoregressions, "CR" for cross regressions,
#' "standardiseCR" for standardised cross regressions (standardised based on estimated within subject variance),
#' "withinVar" for within variance and covariance,
#' "randomImpulse" for expected change in processes given a random fluctuation of +1 for each process
#' (so a mixture of DIFFUSION and DRIFT characteristics),
#' "experimentalImpulse" for expected change in processes given an exogenous input of +1 for each process, provides
#' alternate characterisation of autoregressive and cross regressive plots.
#' @param resolution Numeric. Plot points between each unit of time. Default of 'auto' adapts to xlim and results in 500 points in total.
#' @param impulseIndex Numeric. Only required for impulse plot types, specifies which column of the DRIFT matrix the impulse relates to.
#' @param xlim vector. As per usual for plot(), but xlim may not be negative.
#' @param typeVector Vector of plot types to use for plotting.
#' @param colVector vector of colours to use for plotting.
#' @param ltyVector Vector of line types to use for plotting.
#' @param ... Other options passed to \code{plot()}. ylim is required.
#' @return Character vector of labels from the DRIFT matrix in order plotted - useful for legends. Side-effect: plots graphs.
#' @examples
#' ## Examples set to 'donttest' because they take longer than 5s.
#'
#' ### example from Driver, Oud, Voelkle (2016),
#' ### simulated happiness and leisure time with unobserved heterogeneity.
#' \donttest{
#' data(ctExample1)
#' traitmodel <- ctModel(n.manifest=2, n.latent=2, Tpoints=6, LAMBDA=diag(2),
#' manifestNames=c('LeisureTime', 'Happiness'),
#' latentNames=c('LeisureTime', 'Happiness'), TRAITVAR="auto")
#' traitfit <- ctFit(dat=ctExample1, ctmodelobj=traitmodel)
#' ctPlot(traitfit, plotType='CR', xlim=c(0,5),ylim=c(-1,1))
#' }
#' @export
ctPlot<-function(x,plotType,xlim,resolution=50,impulseIndex=NULL,subject=1,typeVector='auto',colVector='auto',ltyVector='auto',...){
ctfitobj<-x
mxobj<-ctfitobj$mxobj
ctsummary<-summary(ctfitobj,verbose=TRUE)
#read in values
for(i in 1:length(ctsummary)){ #this loop reads in the specified continuous time model so the objects are available
assign(names(ctsummary[i]),eval(parse(text = paste0("ctsummary","$",names(ctsummary[i])))))
}
asymptotes<-ctfitobj$ctfitargs$asymptotes
paramlabels<-paste0(rep(ctfitobj$ctmodelobj$latentNames,ctfitobj$ctmodelobj$n.latent),'_',
rep(ctfitobj$ctmodelobj$latentNames,each=ctfitobj$ctmodelobj$n.latent))
latentNames<-ctfitobj$ctmodelobj$latentNames
n.latent<-ctfitobj$ctmodelobj$n.latent
n.manifest<-ctfitobj$ctmodelobj$n.manifest
n.TIpred<-ctfitobj$ctmodelobj$n.TIpred
n.TDpred<-ctfitobj$ctmodelobj$n.TDpred
manifestNames <- ctfitobj$ctmodelobj$manifestNames
TDpredNames <- ctfitobj$ctmodelobj$TDpredNames
TIpredNames <- ctfitobj$ctmodelobj$TIpredNames
Tpoints<-ctfitobj$ctmodelobj$Tpoints
stationary<-ctfitobj$ctfitargs$stationary
T0MEANS<-OpenMx::mxEval(T0MEANS, ctfitobj$mxobj,compute=F)
T0VAR<-OpenMx::mxEval(T0VAR, ctfitobj$mxobj,compute=F)
if(resolution=='auto') resolution <- ceiling(500/xlim[2])
plottimes<-matrix(seq(0,xlim[2],1/resolution)[-1],ncol=1) #time steps
if(plotType=='mean'){# 2. plot mean trend
#TD predictors
tdpredeffect<-matrix(0,nrow=length(plottimes),ncol=n.latent) #default 0 matrix
#extract and round absolute times to resolution
# if(n.TDpred >0){
# if(ctfitobj$ctfitargs$objective!='Kalman') data<-mxobj$data$observed
# if(ctfitobj$ctfitargs$objective=='Kalman') {
# data<- ctLongToWide(mxobj$data$observed,id='id',time='dT1',
# manifestNames=manifestNames,TDpredNames=TDpredNames,TIpredNames=TIpredNames)
# data<-data[,-which(colnames(data)=='T0'),drop=F]
# colnames(data)[which(colnames(data) %in% paste0('T',1:(Tpoints-1)))]<-paste0('dT',1:(Tpoints-1))
# }
#
# obstimes<-cbind(
# matrix(rep(0,nrow(data)),ncol=1),
# matrix(
# apply(data[,paste0('dT',1:(Tpoints-1)),drop=F],1,function(x){
# obstimesi<-c()
# for(i in 1:(length(x))){
# obstimesi[i]<-round(sum(x[1:i])*resolution,digits=0) / resolution
# }
# return(obstimesi)
# }
# ),
# byrow=T,ncol=(Tpoints-1))
# )
# # obstimes=obstimes[obstimes <=xlim]
#
# #extract tdpred observations and match with rounded times
# TDpredMeans<-matrix(NA,nrow=max(obstimes*resolution),ncol=n.TDpred)
# for(predi in 1:n.TDpred){
# TDpredObs<-matrix(0,nrow=max(obstimes*resolution),ncol=nrow(data))
# TDpredObs[cbind(c(obstimes * resolution),c(row(obstimes)))] <- data[,paste0(TDpredNames[predi],'_T',0:(Tpoints-2))]
# TDpredMeans[,predi]<-apply(TDpredObs,1,mean,na.rm=T)
# }
# TDpredMeans=TDpredMeans[1:length(plottimes),,drop=FALSE]
#
# TDPREDEFFECT <- OpenMx::mxEval(TDPREDEFFECT,mxobj,compute=T)
# ARforResolution <- OpenMx::expm(DRIFT*(1/resolution))
# for(i in 2:length(plottimes)){
# tdpredeffect[i,]<- ARforResolution %*% t(tdpredeffect[i-1, ,drop=F]) + TDPREDEFFECT %*% t(TDpredMeans[i-1, ,drop=F])
# }
# }
means<-matrix(apply(cbind(1:length(plottimes)),1,function(x) {
dT<-plottimes[x]
OpenMx::expm(DRIFT*dT) %*% T0MEANS +
solve(DRIFT) %*% (OpenMx::expm(DRIFT*dT)-diag(nrow(DRIFT))) %*% CINT
}), byrow=T,ncol=nrow(DRIFT))
plotdata<-means
legend<-latentNames
} #end mean plot
if(plotType=='withinVariance'){ #plot within subject variance
DRIFTHATCH<-(DRIFT %x% diag(n.latent) + diag(n.latent) %x% DRIFT)
withinvar<-matrix(apply(plottimes,1,function(x) matrix(
OpenMx::expm(DRIFT %x% x) %*% T0VAR %*% t(OpenMx::expm(DRIFT %x% x)) + #initial variance
matrix(solve(DRIFTHATCH) %*% ((OpenMx::expm(DRIFTHATCH %x% x)) - diag(n.latent^2) ) %*% #diffusion process
OpenMx::rvectorize(DIFFUSION),nrow=n.latent),nrow=n.latent)[row(diag(n.latent))>=col(diag(n.latent))]),
byrow=T,ncol=length(diag(n.latent)[row(diag(n.latent))>=col(diag(n.latent))]))
colnames(withinvar)<-paramlabels[lower.tri(paramlabels,diag=T)]
plotdata<-withinvar
legend<-colnames(withinvar)
} #end within variance plots
if(plotType=='AR'){
ar<- matrix(apply(plottimes,1,function(x) c(diag(OpenMx::expm(DRIFT*x)))),ncol=nrow(DRIFT),byrow=T)
plotdata<-ar
legend<-paste0(c(diag(paramlabels)))
}
if((plotType=='CR' | plotType=='standardiseCR') & ncol(DRIFT)<2) stop('No cross effects to plot!')
if((plotType=='CR' | plotType=='standardiseCR') & ncol(DRIFT)>1){
if(plotType=='standardiseCR') {
standardiser<-try(suppressWarnings(rep(sqrt(diag(abs(asymDIFFUSION))),each=n.latent) / rep(diag(sqrt(abs(asymDIFFUSION))),times=n.latent)))
if('try-error' %in% class(standardiser)) stop('Unable to standardardise cross regression - try unstandardised.')
standardiser[is.nan(standardiser)]<-0
}
if(plotType!='standardiseCR') standardiser<-1
if(ctfitobj$ctfitargs$transformedParams==TRUE) driftindices<-(mxobj$DRIFT$free==TRUE | mxobj$DRIFT$values != 0)
if(ctfitobj$ctfitargs$transformedParams==FALSE) driftindices<-(mxobj$DRIFT$free==TRUE | mxobj$DRIFT$values != 0)
cl<- matrix(apply(plottimes,1,function(x) {
c((OpenMx::expm(DRIFT*x)*standardiser)[row(DRIFT)!=col(DRIFT) & driftindices==TRUE])
}
),ncol=length(DRIFT[row(DRIFT)!=col(DRIFT) & driftindices==TRUE]),byrow=T)
clvars<-paramlabels[row(DRIFT)!=col(DRIFT) & driftindices==TRUE]
plotdata<-cl
legend<-paste0(clvars)
}
if(plotType=='randomImpulse'){
# if(is.null(impulseIndex)) stop('impulseIndex must be specified - integer from 1 to nrow(DRIFT)')
# if(as.integer(impulseIndex) != impulseIndex | impulseIndex < n.latent | impulseIndex > n.latent) stop('impulseIndex must be an integer between 1 and n.latent')
impulse<-rep(0,ncol(DRIFT))
impulse[impulseIndex]<-1
DIFFUSIONdiag<-diag(diag(DIFFUSION),n.latent)
diag(DIFFUSIONdiag)[diag(DIFFUSIONdiag)<=0]<-.0001
# if(standardiseCR==TRUE)
DIFFUSIONdiag<-solve(sqrt(DIFFUSIONdiag)) %*% DIFFUSION %*% t(solve(sqrt(DIFFUSIONdiag)))
impulse<- DIFFUSIONdiag %*% impulse
impulseresponse<-lapply(plottimes,function(x) expm(DRIFT *x) %*% impulse)
IRylim <- c(-1,1)
if(max(unlist(impulseresponse)) > 1) IRylim[2] <- max(unlist(impulseresponse))
if(min(unlist(impulseresponse)) < -1) IRylim[1] <- min(unlist(impulseresponse))
plotdata<-matrix(unlist(lapply(1:length(impulseresponse[[1]]),function(z) unlist(lapply(impulseresponse,function(x) x[z] )))),ncol=nrow(DRIFT))
legend<-paste0(latentNames)
}
if(plotType=='experimentalImpulse'){
if(is.null(impulseIndex)) stop('impulseIndex must be specified - integer from 1 to nrow(DRIFT)')
impulse<-rep(0,ncol(DRIFT))
impulse[impulseIndex]<-1
impulseresponse<-lapply(plottimes,function(x) expm(DRIFT *x) %*% impulse)
IRylim <- c(-1,1)
if(max(unlist(impulseresponse)) > 1) IRylim[2] <- max(unlist(impulseresponse))
if(min(unlist(impulseresponse)) < -1) IRylim[1] <- min(unlist(impulseresponse))
plotdata<-matrix(unlist(lapply(1:length(impulseresponse[[1]]),function(z) unlist(lapply(impulseresponse,function(x) x[z] )))),ncol=nrow(DRIFT))
legend<-paste0(latentNames)
}
if(plotType=='factorScores'){
newdata<-mxobj@data@observed[subject,,drop=F]
smallmxobj<-mxModel(mxobj,mxData(newdata,type='raw'))
plottimes <- newdata[colnames(newdata) %in% paste0('dT',1:(Tpoints-1))]
for(i in 2:length(plottimes)){
plottimes[i]<-plottimes[i-1]+plottimes[i]
}
plottimes<-c(0,plottimes)
newdata<-newdata[1:(n.manifest*Tpoints)]
scores<-mxFactorScores(smallmxobj,type='Regression')
amat <- mxEval(A, mxobj, defvar.row=subject, compute=TRUE)
mmat <- mxEval(M, mxobj, defvar.row=subject, compute=TRUE)
dataRow1 <- testd[1,]
fullscores <- solve(diag(1, n.latent*Tpoints) - amat[1:(n.latent*Tpoints), 1:(n.latent*Tpoints)]) %*% (scores[1,,1] - matrix(mmat)[1:(n.latent*Tpoints),])
plotdata<-cbind(matrix(newdata,ncol=n.manifest,byrow=T),matrix(fullscores,ncol=n.latent,byrow=T))
errorbars<-matrix(scores[,,2],ncol=n.latent,byrow=T)
if(ltyVector[1]=='auto') ltyVector <- c(rep(1,n.manifest),rep(2,n.latent))
if(typeVector[1]=='auto') typeVector <- rep('b',ncol(plotdata))
}
if(colVector[1]=='auto') colVector <- grDevices::rainbow(ncol(plotdata),v=.8) #set plot colours
if(ltyVector[1]=='auto') ltyVector <- rep(1,ncol(plotdata))
if(typeVector[1]=='auto') typeVector <- rep('l',ncol(plotdata))
graphics::plot(plottimes,rep(NA,length(plottimes)),...)
graphics::grid()
graphics::points(plottimes,plotdata[,1],type=typeVector[1], col=colVector[1],lty=ltyVector[1],...)
if(ncol(plotdata) > 1) {
for(rowi in 2:ncol(plotdata)){
graphics::points(plottimes, plotdata[,rowi],type=typeVector[rowi], col=colVector[rowi],lty=ltyVector[rowi],...)
if(plotType=='factorScores' & rowi > n.manifest){
graphics::arrows(plottimes, plotdata[,rowi] - errorbars[,rowi-n.manifest], plottimes, plotdata[,rowi] + errorbars[,rowi-n.manifest],
col=colVector[rowi],length=0.05, angle=90, code=3)
}
}
}
return(legend)
}
Try the ctsemOMX 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.