R/summary.ctsemFit.R

Defines functions ctSummaryMatrices ctParamsSummary omxSummary summary.ctsemFit summary.ctsemMultigroupFit

Documented in summary.ctsemFit summary.ctsemMultigroupFit

#' Summary function for ctsemMultigroupFit object
#'
#' Provides summary details for objects fitted with \code{\link{ctMultigroupFit}}.
#' 
#' @param object ctsemMultigroupFit object as generated by \code{\link{ctMultigroupFit}}
#' @param group character string of subgroup to display summary parameters for. Default of 'show chooser' displays list and lets you select.
#' @param ... additional parameters to pass to \code{\link{summary.ctsemFit}}.
#' @return Summary of ctsemMultigroupFit object
#' @method summary ctsemMultigroupFit
#' @export
summary.ctsemMultigroupFit<-function(object,group='show chooser',...){
  
  if(group=='show chooser')  {
    message(paste0(1:length(object$groups), '  ', object$groups, '\n'))
    
    selection<-as.numeric(readline(
      'Enter the number of the group you wish to display the summary matrices for (OpenMx summary portion is same across groups):  '))
    
    selection<-object$groups[[selection]]
  }
  
  if(group != 'show chooser') selection <- object$groups[[group]]
  
  #create ctFit object containing only single group
  tempctobj<-list(ctfitargs=object$ctfitargs, ctmodelobj=object$ctmodelobj, mxobj=object$mxobj[[selection]])
  class(tempctobj)<-'ctsemFit'
  
  out<-ctSummaryMatrices(tempctobj,...)
  # out<-c(out,ctParamsSummary(tempctobj, ctSummaryMatrices=out))
  
  #replace mxobj in single group ctFit with whole model object
  tempctobj$mxobj <- object$mxobj
  
  out$ctParameters<-lapply(object$groups,function(x){
    ctParamsSummary(object=tempctobj,ctSummaryMatrices=out,group=x)
  })
  out$ctParameters <- do.call(rbind,out$ctParameters)
  out$ctParameters <- out$ctParameters[!duplicated(rownames(out$ctParameters)),,drop=FALSE]
  out$ctparammessage<-'Note: Continuous time parameter estimates above are of the full variance-covariance matrices, not cholesky decompositions as used by ctModel.'
  # if(tempctobj$ctfitargs$transformedParams==TRUE) out$ctparammessage<- c(out$ctparammessage, 'Note: Covariance related standard errors are approximated with delta method..')
  
  out$omxsummary<-omxSummary(tempctobj,verbose=TRUE)
  return(out)
}

#' Summary function for ctsemFit object
#'
#' Provides summary details for ctsemFit objects.
#' 
#' @param object ctsemFit object as generated by ctFit.
#' @param ridging if TRUE, adds a small amount of variance to diagonals when calculating standardised (correlation) matrices,
#' should only be used if standardised matrices return NAN.
#' @param timeInterval positive numeric value specifying time interval to use for discrete parameter matrices, defaults to 1.
#' @param verbose Logical. If TRUE, displays the raw, internally transformed (when fitting with default arguments) OpenMx parameters and corresponding standard errors, as well as
#' additional summary matrices. Parameter transforms are described in the vignette, \code{vignette('ctsem')}. Additional summary matrices
#' include: 'discrete' matrices -- matrices representing
#' the effect for the given time interval (default of 1); 'asymptotic' matrices -- represents the effect as time interval
#' approaches infinity (therefore asymCINT describes mean level of processes at the asymptote, asymDIFFUSION describes total within-
#' subject variance at the asymptote, etc); 'standardised' matrices -- transforms covariance matrices to correlation matrices, and transforms 
#' discreteDRIFT based on DIFFUSION, to give effect sizes.
#' @param ... additional parameters to pass.
#' @return Summary of ctsemFit object
#' @method summary ctsemFit
#' @details Important: Although \code{\link{ctModel}} takes cholesky decomposed variance-covariance matrices as input,
#' the summary function displays the full variance-covariance matrices. These can be cholesky decomposed for comparison purposes using
#' \code{t(chol(summary(ctfitobject)$covariancematrix))}.
#' Standard errors are displayed in the $ctparameters section, however if \code{\link{ctFit}} was used with transformedParams=TRUE (the default, and recommended) 
#' covariance matrix standard errors will have been approximated using the delta method. For 
#' inferential purposes, maximum likelihood confidence intervals may be estimated using the \code{\link{ctCI}} function.
#' @examples 
#' ## Examples set to 'donttest' because they take longer than 5s. 
#' 
#' \donttest{
#' ### example from Driver, Oud, Voelkle (2015), 
#' ### simulated happiness and leisure time with unobserved heterogeneity.
#' 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)
#' summary(traitfit,timeInterval=1)
#' }
#' @export
summary.ctsemFit<-function(object,ridging=FALSE,timeInterval=1,verbose=FALSE,...){

  output<-ctSummaryMatrices(object=object,ridging=ridging,timeInterval=timeInterval,verbose=verbose)
  ctParameters<-ctParamsSummary(object=object,ctSummaryMatrices=output)
  output$ctparameters<-ctParameters
  # if(object$ctfitargs$transformedParams==TRUE) 
  output$ctparammessage<-'Note: Continuous time parameter estimates above are of the full variance-covariance matrices, not cholesky decompositions as used by ctModel.'
  if(object$ctfitargs$transformedParams==TRUE) output$ctparammessage<- c(output$ctparammessage, 'Note: Standard errors are approximated with delta method so are only rough approximations.')
  output$omxsummary<-omxSummary(object=object,verbose=verbose)
  if(verbose==FALSE) output$message<-'For additional summary matrices and raw OpenMx parameter output, use argument verbose=TRUE' 
  return(output)
  
}




omxSummary<-function(object,verbose=FALSE){
  omxver<-unlist(utils::packageVersion('OpenMx'))
  if(omxver[1] < 3 & omxver[2] < 5) {
    omxsummary<-methods::getMethod("summary","MxModel")(object$mxobj) #get openmx summary
  } else { 
    omxsummary<-utils::getS3method("summary","MxModel")(object$mxobj) #get openmx summary
  }
  
  output<-list()
  if(verbose==TRUE) output<-c(omxsummary['parameters'])
  if(verbose==TRUE & object$ctfitargs$transformedParams==TRUE) output$parammessage<-'Note: Parameters and std errors displayed in OpenMx section are based on internal transformations of the continuous time parameters, for optimization purposes. See the vignette for details.'
  output<-c(output,omxsummary['modelName'],omxsummary['wallTime'],omxsummary['mxVersion'],omxsummary['timestamp'],
    omxsummary['estimatedParameters'], omxsummary['CI'],
    omxsummary['AIC.Mx'],omxsummary['BIC.Mx'],omxsummary['observedStatistics'],
    omxsummary['npsolMessage'],
    omxsummary['degreesOfFreedom'],omxsummary['Minus2LogLikelihood'] )
  
  return(output)
}


ctParamsSummary<-function(object,ctSummaryMatrices,group=NA){
  
  # browser()
  if(is.na(group[1])) mxobj <- object$mxobj else mxobj <- object$mxobj[[group]]
  parnames<-names(omxGetParameters(mxobj))
  parvalues<-omxGetParameters(mxobj)
  newparvalues<-parvalues
  if(is.null(object$mxobj$output$standardErrors)) parsd<-matrix(NA,nrow=length(parvalues))
  else parsd<-object$mxobj$output$standardErrors
  parsd <- parsd[rownames(parsd) %in% parnames,,drop=FALSE]
  parmatrix<-rep(NA,length(parnames))
 
  # 
  for(parami in 1:length(parnames)){ #for every free param
    for(matrixi in names(ctSummaryMatrices)[names(ctSummaryMatrices) %in% names(object$ctmodelobj)]){ #check every matrix that is in both ctmodelobj and output
      if(parnames[parami] %in% object$ctmodelobj[[matrixi]]) { #if the free param is in the ctmodelobj matrix
        ind <- which(parnames[parami] == object$ctmodelobj[[matrixi]],arr.ind = TRUE)[1,,drop=FALSE]
        parmatrix[parami]<-matrixi
        newparvalues[parami]<-ctSummaryMatrices[[matrixi]][match(parnames[parami],object$ctmodelobj[[matrixi]])]
        parsd[parami]<- parsd[parami]
        #delta approx for sd pars
        if(newparvalues[parami] != parvalues[parami]) {
          # print(parami)
          parsd[parami]<- 
          suppressMessages(mxSE(paste0(ifelse(is.na(group[1]),'',paste0(group,'.')),matrixi,'[',ind[1,1],',',ind[1,2],']'),object$mxobj))
        }
         # (exp(parvalues[parami] + parsd[parami]) - exp(parvalues[parami] - parsd[parami]))/2
        # parsd[parami]<- abs(((newparvalues[parami]) / parvalues[parami]) * parsd[parami]) #old, bad? first order delta approximation of std error
        parvalues[parami]<- newparvalues[parami]
      }
    }
  }
  out<-data.frame(parvalues,parmatrix,parsd)
  colnames(out)<-c( 'Value', 'Matrix', 'StdError')
  out<-as.data.frame(out)
}





ctSummaryMatrices<-function(object,ridging=FALSE,timeInterval=1,verbose=FALSE,...){
  
  if(ridging==TRUE) ridging<- .0001 else ridging <- 0
  
  outlist<-c()
  
  asymptotes<-object$ctfitargs$asymptotes
  stationary<-object$ctfitargs$stationary
  discreteTime <- object$ctfitargs$discreteTime
  
  mxobj<-object$mxobj
  
  out<-list()
  
  n.latent<-nrow(OpenMx::mxEval(DRIFT,mxobj,compute=TRUE))
  n.TIpred<-object$ctmodelobj$n.TIpred
  n.TDpred<-object$ctmodelobj$n.TDpred
  n.manifest<-nrow(OpenMx::mxEval(LAMBDA,mxobj,compute=TRUE))
  
  
  try(latentNames<-object$ctmodelobj$latentNames)
  try(manifestNames<-object$ctmodelobj$manifestNames)
  output<-list()
  
  
  if(discreteTime==FALSE){
    
    LAMBDA<-tryCatch({ OpenMx::mxEval(LAMBDA, mxobj,compute=TRUE)}, error=function(e) e )
    tryCatch({  dimnames(LAMBDA)<-list(manifestNames,latentNames)}, error=function(e) e )
    outlist<-c(outlist,'LAMBDA')
    
    DRIFT<-tryCatch({ OpenMx::mxEval(DRIFT, mxobj,compute=TRUE)}, error=function(e) e )
    tryCatch({  dimnames(DRIFT)<-list(latentNames,latentNames)}, error=function(e) e )
    outlist<-c(outlist,'DRIFT')
    
    if(verbose==TRUE){
      discreteDRIFT<-tryCatch({ OpenMx::expm(DRIFT * timeInterval)}, error=function(e) e )
      tryCatch({  dimnames(discreteDRIFT)<-list(latentNames,latentNames)}, error=function(e) e )
      outlist<-c(outlist,'discreteDRIFT')
      
      outlist<-c(outlist,'discreteDRIFTstd')
    }
    
    DRIFTHATCH<-tryCatch({ (DRIFT %x% diag(n.latent)) + (diag(n.latent) %x% DRIFT)}, error=function(e) e )
    
    MANIFESTVAR<-tryCatch({ mxEval(MANIFESTVAR,mxobj, compute=TRUE)}, error=function(e) e )
    tryCatch({  dimnames(MANIFESTVAR)<-list(manifestNames,manifestNames)}, error=function(e) e )
    outlist<-c(outlist,'MANIFESTVAR')

    if(verbose==TRUE){
      MANIFESTVARdiag<-tryCatch({ diag(diag(MANIFESTVAR),n.manifest)+diag(c(ridging),n.manifest)}, error=function(e) e )
      MANIFESTVARstd<-tryCatch({ suppressWarnings(solve(sqrt(MANIFESTVARdiag)) %&% MANIFESTVAR)}, error=function(e) e )
      tryCatch({  dimnames(MANIFESTVARstd)<-list(latentNames,latentNames)}, error=function(e) e )
      outlist<-c(outlist,'MANIFESTVARstd')
    }
    
    MANIFESTMEANS<-tryCatch({ mxEval(MANIFESTMEANS,mxobj, compute=TRUE)}, error=function(e) e )
    tryCatch({  dimnames(MANIFESTMEANS)<-list(manifestNames,manifestNames)}, error=function(e) e )
    outlist<-c(outlist,'MANIFESTMEANS')
    
    CINT<-tryCatch({ mxobj$CINT$values}, error=function(e) e )
    if(asymptotes==TRUE) CINT <- tryCatch({ -DRIFT %*% CINT}, error=function(e) e )
    tryCatch({  rownames(CINT)<-latentNames}, error=function(e) e )
    outlist<-c(outlist,'CINT')
    
    if(verbose==TRUE){
      discreteCINT<-tryCatch({ solve(DRIFT) %*%(discreteDRIFT - diag(n.latent)) %*% CINT}, error=function(e) e )
      tryCatch({  rownames(discreteCINT)<-latentNames}, error=function(e) e )
      outlist<-c(outlist,'discreteCINT')
      
      asymCINT<-tryCatch({ -solve(DRIFT) %*% CINT}, error=function(e) e )
      tryCatch({  rownames(asymCINT)<-latentNames}, error=function(e) e )
      outlist<-c(outlist,'asymCINT')
    }
    
    
    
    DIFFUSION<-tryCatch({ mxEval(DIFFUSION, mxobj, compute=TRUE)}, error=function(e) e )
    if(asymptotes==TRUE) DIFFUSION <- tryCatch({ matrix(-DRIFTHATCH %*% OpenMx::rvectorize(DIFFUSION), nrow=n.latent, ncol=n.latent)}, error=function(e) e )
    tryCatch({  dimnames(DIFFUSION)<-list(latentNames,latentNames)}, error=function(e) e )
    outlist<-c(outlist,'DIFFUSION')
    
    if(verbose==TRUE){
      discreteDIFFUSION<-tryCatch({ matrix(solve(DRIFTHATCH) %*% ((OpenMx::expm(DRIFTHATCH * timeInterval)) - 
          (diag(n.latent) %x% diag(n.latent)) ) %*% OpenMx::rvectorize(DIFFUSION), nrow=n.latent) }, error=function(e) e )
      tryCatch({  dimnames(discreteDIFFUSION)<-list(latentNames,latentNames)}, error=function(e) e )
      outlist<-c(outlist,'discreteDIFFUSION')
      
      asymDIFFUSION<-tryCatch({ matrix(-solve(DRIFTHATCH) %*% OpenMx::cvectorize(DIFFUSION),nrow=n.latent)}, error=function(e) e )
      tryCatch({  dimnames(asymDIFFUSION)<-list(latentNames,latentNames)}, error=function(e) e )
      outlist<-c(outlist,'asymDIFFUSION')
      
      asymDIFFUSIONdiag<-tryCatch({ diag(diag(asymDIFFUSION),n.latent)+diag(c(ridging),n.latent)}, error=function(e) e )
      asymDIFFUSIONstd<-tryCatch({ suppressWarnings(solve(sqrt(asymDIFFUSIONdiag)) %&% asymDIFFUSION)}, error=function(e) e )
      tryCatch({  dimnames(asymDIFFUSIONstd)<-list(latentNames,latentNames)}, error=function(e) e )
      outlist<-c(outlist,'asymDIFFUSIONstd')
      
      #std dev of affecting latent divided by std dev of affected latent
      standardiser<-tryCatch({ suppressWarnings(rep(sqrt(diag(asymDIFFUSION)),each=n.latent) / rep(sqrt(diag(asymDIFFUSION)),times=n.latent))}, error=function(e) e )
      discreteDRIFTstd<-tryCatch({ discreteDRIFT * standardiser    }, error=function(e) e )
      tryCatch({  dimnames(discreteDRIFTstd)<-list(latentNames,latentNames)}, error=function(e) e )
    }
    #added to outlist above
    
    #             asymTOTALVAR<-tryCatch({ asymDIFFUSION }, error=function(e) e )#set totalvar to t0withinvar to begin
    #     tryCatch({  dimnames(asymTOTALVAR)<-list(latentNames,latentNames)}, error=function(e) e )
    #     outlist<-c(outlist,'asymTOTALVAR')
    #     
    #     
    #             asymTOTALVARdiag<-tryCatch({ diag(diag(asymTOTALVAR),n.latent)+diag(c(ridging),n.latent)}, error=function(e) e )
    #             asymTOTALVARstd<-tryCatch({ suppressWarnings(solve(sqrt(asymTOTALVARdiag)) %&% asymTOTALVAR)}, error=function(e) e )
    #     tryCatch({  dimnames(asymTOTALVARstd)<-list(latentNames,latentNames)}, error=function(e) e )
    #     outlist<-c(outlist,'asymTOTALVARstd')
    
    
    # if('T0VAR' %in% stationary == FALSE){ #then include base T0 matrices
      T0VAR<-tryCatch({ OpenMx::mxEval(T0VAR, mxobj,compute=TRUE)}, error=function(e) e )
      tryCatch({  dimnames(T0VAR)<-list(latentNames,latentNames)}, error=function(e) e )
      outlist<-c(outlist,'T0VAR')
      
      #       T0TOTALVAR<-tryCatch({ T0VAR }, error=function(e) e )#set t0totalvar to t0withinvar to begin
      #       tryCatch({  dimnames(T0TOTALVAR)<-list(latentNames,latentNames)}, error=function(e) e )
      #       outlist<-c(outlist,'T0TOTALVAR')
      
      if(verbose==TRUE){
        T0VARdiag<-tryCatch({ diag(diag(T0VAR),n.latent)+diag(c(ridging),n.latent)}, error=function(e) e )
        T0VARstd<-tryCatch({ suppressWarnings(solve(sqrt(T0VARdiag)) %&% T0VAR)}, error=function(e) e )
        tryCatch({  dimnames(T0VARstd)<-list(latentNames,latentNames)}, error=function(e) e )
        outlist<-c(outlist,'T0VARstd')
      }
      
    # } # end T0VAR matrices
    
    # if('T0MEANS' %in% stationary == FALSE){ #then include base T0 matrices
      T0MEANS<-tryCatch({ OpenMx::mxEval(T0MEANS, mxobj,compute=TRUE)}, error=function(e) e )
      tryCatch({  rownames(T0MEANS)<-latentNames}, error=function(e) e )
      outlist<-c(outlist,'T0MEANS')
    # }
    
    
    #trait matrices
    if(any(object$ctmodelobj$TRAITVAR != 0)) {
      
      TRAITVAR<-tryCatch({ OpenMx::mxEval(TRAITVAR, mxobj,compute=TRUE)}, error=function(e) e )
      # if(asymptotes==TRUE) TRAITVAR<-tryCatch({ DRIFT %&% TRAITVAR }, error=function(e) e )
      tryCatch({  dimnames(TRAITVAR)<-list(latentNames,latentNames)}, error=function(e) e )
      outlist<-c(outlist,'TRAITVAR')
      
            #   asymTRAITVAR<-tryCatch({ solve(DRIFT) %*% TRAITVAR %*% t(solve(DRIFT))}, error=function(e) e )
            # tryCatch({  dimnames(asymTRAITVAR)<-list(latentNames,latentNames)}, error=function(e) e )
            # outlist<-c(outlist,'asymTRAITVAR')


      #       discreteTRAITVAR<-tryCatch({ (diag(n.latent)-OpenMx::expm(DRIFT * timeInterval)) %*% asymTRAITVAR %*% t((diag(n.latent)-OpenMx::expm(DRIFT * timeInterval)))}, error=function(e) e )
      #       tryCatch({  dimnames(discreteTRAITVAR)<-list(latentNames,latentNames)}, error=function(e) e )
      #       outlist<-c(outlist,'discreteTRAITVAR')

      if(verbose==TRUE){
        TRAITVARdiag<-tryCatch({ diag(diag(TRAITVAR))}, error=function(e) e )
        TRAITVARstd<-tryCatch({ suppressWarnings(solve(sqrt(TRAITVARdiag)) %&% TRAITVAR)}, error=function(e) e )
        tryCatch({  dimnames(TRAITVARstd)<-list(latentNames,latentNames)}, error=function(e) e )
        outlist<-c(outlist,'TRAITVARstd')
      }
      #           asymTOTALVAR<-tryCatch({ asymDIFFUSION + asymTRAITVAR}, error=function(e) e )
      #       tryCatch({  dimnames(asymTOTALVAR)<-list(latentNames,latentNames)}, error=function(e) e )
      #       
      #           asymTOTALVARdiag<-tryCatch({ diag(diag(asymTOTALVAR),n.latent)+diag(c(ridging),n.latent)}, error=function(e) e )
      #           asymTOTALVARstd<-tryCatch({ suppressWarnings(solve(sqrt(asymTOTALVARdiag)) %&% asymTOTALVAR)}, error=function(e) e )
      #       tryCatch({  dimnames(asymTOTALVARstd)<-list(latentNames,latentNames)}, error=function(e) e )
      
      #       if('T0TRAITEFFECT' %in% stationary == FALSE) { #then include T0 trait matrices
      #         
      #          T0TRAITEFFECT<-tryCatch({mxEval(T0TRAITEFFECT, mxobj,compute=TRUE)}, error=function(e) e )
      #         tryCatch({  dimnames(T0TRAITEFFECT)<-list(latentNames,latentNames)}, error=function(e) e )
      #         outlist<-c(outlist,'T0TRAITEFFECT')
      #         
      #          T0TRAITVAR<-tryCatch({ T0TRAITEFFECT %*% TRAITVAR %*% t(T0TRAITEFFECT) }, error=function(e) e )#is this valid?
      #         if(asymptotes
      #         tryCatch({  dimnames(T0TRAITVAR)<-list(latentNames,latentNames)}, error=function(e) e )
      #         outlist<-c(outlist,'T0TRAITVAR')
      #         
      #         
      #          T0TRAITVARdiag<-tryCatch({ suppressWarnings(sqrt(diag(diag(T0TRAITVAR),n.latent)))}, error=function(e) e )
      #          T0TRAITVARstd<-tryCatch({ solve(T0TRAITVARdiag) %*% T0TRAITVAR %*% solve(T0TRAITVARdiag)}, error=function(e) e )
      #         tryCatch({  dimnames(T0TRAITVARstd)<-list(latentNames,latentNames)}, error=function(e) e )
      #         outlist<-c(outlist,'T0TRAITVARstd')
      #         
      #          T0TOTALVAR<-tryCatch({ T0TRAITVAR+T0VAR}, error=function(e) e )
      #         tryCatch({  dimnames(T0TOTALVAR)<-list(latentNames,latentNames)}, error=function(e) e )
      #         
      #          T0TOTALVARdiag<-tryCatch({ diag(diag(T0TOTALVAR),n.latent)}, error=function(e) e )
      #          T0TOTALVARstd<-tryCatch({ suppressWarnings(solve(sqrt(T0TOTALVARdiag)) %*% T0TOTALVAR %*% solve(sqrt(T0TOTALVARdiag)))}, error=function(e) e )
      #         tryCatch({  dimnames(T0TOTALVARstd)<-list(latentNames,latentNames)}, error=function(e) e )
      #         
      #       } #end T0 trait matrices
    }#end trait matrices
    
    
    
    
    
    
    
    
    if(any(object$ctmodelobj$MANIFESTTRAITVAR != 0)) {
      
      MANIFESTTRAITVAR<-tryCatch({ mxEval(MANIFESTTRAITVAR, mxobj,compute=TRUE)}, error=function(e) e )
      tryCatch({  dimnames(MANIFESTTRAITVAR)<-list(manifestNames,manifestNames)}, error=function(e) e )
      outlist<-c(outlist,'MANIFESTTRAITVAR')
      
      if(verbose==TRUE){
        MANIFESTTRAITVARdiag<-tryCatch({ diag(diag(MANIFESTTRAITVAR),n.manifest)+diag(c(ridging),n.manifest)}, error=function(e) e )
        MANIFESTTRAITVARstd<-tryCatch({ suppressWarnings(solve(sqrt(MANIFESTTRAITVARdiag)) %&% MANIFESTTRAITVAR)}, error=function(e) e )
        tryCatch({  dimnames(MANIFESTTRAITVARstd)<-list(manifestNames,manifestNames)}, error=function(e) e )
        outlist<-c(outlist,'MANIFESTTRAITVARstd')
      }
      
    }#end manifesttrait
    
    
    
    
    
    #TIpred matrices
    if(n.TIpred > 0) {
      TIpredNames<-object$ctmodelobj$TIpredNames      
      if(asymptotes==FALSE) TIPREDEFFECT<-tryCatch({mxobj$TIPREDEFFECT$values}, error=function(e) e )
      if(asymptotes==TRUE)   TIPREDEFFECT<-tryCatch({-DRIFT %*% mxobj$TIPREDEFFECT$values}, error=function(e) e )
      tryCatch({  dimnames(TIPREDEFFECT)<-list(latentNames,TIpredNames)}, error=function(e) e )
      outlist<-c(outlist,'TIPREDEFFECT')
      
      #       if (randomPredictors==TRUE) TIPREDVAR<-tryCatch({mxobj$S$values[TIpredNames,TIpredNames] }, error=function(e) e )
      #       if(randomPredictors==FALSE) 
      TIPREDVAR<-tryCatch({OpenMx::mxEval(TIPREDVAR,mxobj,compute=TRUE) }, error=function(e) e )
      outlist<-c(outlist,'TIPREDVAR')
      
      if(verbose==TRUE){
        TIPREDVARdiag<-tryCatch({diag(diag(TIPREDVAR),n.TIpred)+diag(c(ridging),n.TIpred)}, error=function(e) e )
        TIPREDVARstd<-tryCatch({suppressWarnings(solve(sqrt(TIPREDVARdiag)) %&% TIPREDVAR)}, error=function(e) e )
        tryCatch({  dimnames(TIPREDVARstd)<-list(TIpredNames,TIpredNames)}, error=function(e) e )
        outlist<-c(outlist,'TIPREDVARstd')
        
        discreteTIPREDEFFECT<-tryCatch({solve(DRIFT) %*%(discreteDRIFT - diag(n.latent)) %*% TIPREDEFFECT }, error=function(e) e )#disrete from continuous
        tryCatch({  dimnames(discreteTIPREDEFFECT)<-list(latentNames,TIpredNames) }, error=function(e) e )
        outlist<-c(outlist,'discreteTIPREDEFFECT')
        
        asymTIPREDEFFECT<-tryCatch({solve(diag(n.latent)-discreteDRIFT)  %*%  (discreteTIPREDEFFECT)}, error=function(e) e )
        tryCatch({  dimnames(asymTIPREDEFFECT)<-list(latentNames,TIpredNames)}, error=function(e) e )
        outlist<-c(outlist,'asymTIPREDEFFECT')
        
        #sqrt of affecting latent variance divided by sqrt of affected
        standardiser<-tryCatch({suppressWarnings(rep(sqrt(diag(TIPREDVAR)),each=n.latent) / rep(diag(sqrt(asymDIFFUSION)),times=n.TIpred))}, error=function(e) e )
        asymTIPREDEFFECTstd<-tryCatch({asymTIPREDEFFECT * standardiser}, error=function(e) e )
        tryCatch({  dimnames(asymTIPREDEFFECTstd)<-list(latentNames,TIpredNames)}, error=function(e) e )
        outlist<-c(outlist,'asymTIPREDEFFECTstd')
        
        addedTIPREDVAR <-tryCatch({asymTIPREDEFFECT %*% TIPREDVAR %*% t(asymTIPREDEFFECT) }, error=function(e) e )#valid?
        tryCatch({  dimnames(addedTIPREDVAR)<-list(latentNames,latentNames)}, error=function(e) e )
        outlist<-c(outlist,'addedTIPREDVAR')
        
        addedTIPREDVARdiag<-tryCatch({diag(diag(addedTIPREDVAR),n.latent)+diag(c(ridging),n.latent)}, error=function(e) e )
        addedTIPREDVARstd<-tryCatch({suppressWarnings(solve(sqrt(addedTIPREDVARdiag)) %&% addedTIPREDVAR)}, error=function(e) e )
        tryCatch({  dimnames(addedTIPREDVARstd)<-list(latentNames,latentNames)}, error=function(e) e )
        outlist<-c(outlist,'addedTIPREDVARstd')
      }
      
      #            asymTOTALVAR<-tryCatch({asymTOTALVAR + addedTIPREDVAR}, error=function(e) e )
      #             tryCatch({  dimnames(asymTOTALVAR)<-list(latentNames,latentNames)}, error=function(e) e )
      #       
      #               asymTOTALVARdiag<-tryCatch({diag(diag(asymTOTALVAR),n.latent)+diag(c(ridging),n.latent)}, error=function(e) e )
      #             asymTOTALVARstd<-tryCatch({suppressWarnings(solve(sqrt(asymTOTALVARdiag)) %&% asymTOTALVAR)}, error=function(e) e )
      #               tryCatch({  dimnames(asymTOTALVARstd)<-list(latentNames,latentNames)}, error=function(e) e )
      
      if('TOTIPRED' %in% stationary == FALSE) { #include TIPRED T0 matrices
        T0TIPREDEFFECT<-tryCatch({OpenMx::mxEval(T0TIPREDEFFECT, mxobj,compute=TRUE)}, error=function(e) e )
        tryCatch({  dimnames(T0TIPREDEFFECT)<-list(latentNames,TIpredNames)}, error=function(e) e )
        outlist<-c(outlist,'T0TIPREDEFFECT')
        
        #sqrt of affecting latent variance divided by sqrt of affected
        if(verbose==TRUE){
          standardiser<-tryCatch({suppressWarnings(rep(sqrt(diag(TIPREDVAR)),each=n.latent) / rep(diag(sqrt(T0VAR)),times=n.TIpred))}, error=function(e) e )
          T0TIPREDEFFECTstd<-tryCatch({T0TIPREDEFFECT * standardiser}, error=function(e) e )
          tryCatch({  dimnames(T0TIPREDEFFECTstd)<-list(latentNames,TIpredNames)}, error=function(e) e )
          outlist<-c(outlist,'T0TIPREDEFFECTstd')
          
          addedT0TIPREDVAR<-tryCatch({T0TIPREDEFFECT %*% TIPREDVAR %*% t(T0TIPREDEFFECT)}, error=function(e) e )#is this valid?
          tryCatch({  dimnames(addedT0TIPREDVAR)<-list(latentNames,latentNames)}, error=function(e) e )
          outlist<-c(outlist,'addedT0TIPREDVAR')
        }
        #       
        #       T0TOTALVAR<-tryCatch({addedT0TIPREDVAR+T0TOTALVAR }, error=function(e) e )#update totalvar
        #       tryCatch({  dimnames(T0TOTALVAR)<-list(latentNames,latentNames)}, error=function(e) e )
        #       
        #       T0TOTALVARdiag<-tryCatch({diag(diag(T0TOTALVAR),n.latent)+diag(c(ridging),n.latent)}, error=function(e) e )
        #       T0TOTALVARstd<-tryCatch({suppressWarnings(solve(sqrt(T0TOTALVARdiag)) %&% T0TOTALVAR)}, error=function(e) e )
        #       tryCatch({  dimnames(T0TOTALVARstd)<-list(latentNames,latentNames)}, error=function(e) e )
        
      } #end T0 TIPRED matrices
      
    } # end TIPRED matrices
    
    
    
    if(n.TDpred > 0) {
      
      TDpredNames<-object$ctmodelobj$TDpredNames      
      
      TDPREDEFFECT<-tryCatch({mxobj$TDPREDEFFECT$values}, error=function(e) e )
      tryCatch({  dimnames(TDPREDEFFECT)<-list(latentNames,TDpredNames)}, error=function(e) e )
      outlist<-c(outlist,'TDPREDEFFECT')
      
      # TDPREDVAR<-tryCatch({OpenMx::mxEval(TDPREDVAR,mxobj,compute=TRUE) }, error=function(e) e )
      # tryCatch({  dimnames(TDPREDVAR)<-list(TDpredNames,TDpredNames)}, error=function(e) e )
      # outlist<-c(outlist,'TDPREDVAR')
      # 
      # if(verbose==TRUE){
      #   TDPREDVARdiag<-tryCatch({diag((diag(TDPREDVAR))+diag(c(ridging),n.TDpred))}, error=function(e) e )
      #   TDPREDVARstd<-tryCatch({suppressWarnings(solve(sqrt(TDPREDVARdiag)) %&% TDPREDVAR)}, error=function(e) e )
      #   tryCatch({  dimnames(TDPREDVARstd)<-list(TDpredNames,TDpredNames)}, error=function(e) e )
      #   outlist<-c(outlist,'TDPREDVARstd')
      # }
      
    } # end TDpred section
    
    ctsummary<-mget(outlist)
    output<-ctsummary
  }  #end summary from continuous params
  return(output)
}

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.