tests/testthat/test-bootHessian.R

In cdriveraus/ctsem: Continuous Time Structural Equation Modelling

if(FALSE){
  if(identical(Sys.getenv("NOT_CRAN"), "true")& .Machine$sizeof.pointer != 4){
    
    
    # library(future)
    # library(data.table)
    # plan(strategy = 'multisession',workers=10)
    
    testfunc <- function(i){
      message(i)
      nsubjects=200
      t0m <- 10+rnorm(nsubjects)
      cint <- t0m/2+rnorm(nsubjects)
      # cor(t0m,cint)
      for(subi in 1:nsubjects){
        gm <- suppressMessages(ctModel(Tpoints=10,
          LAMBDA=matrix(1), 
          DRIFT= -1,
          T0MEANS = t0m[subi],
          DIFFUSION=.5,
          MANIFESTVAR = 0.5,
          T0VAR = 0,
          MANIFESTMEANS = 0,
          CINT=cint[subi]))
        
        dd <- suppressMessages(data.frame(ctGenerate(ctmodelobj = gm,n.subjects = 1,
          burnin = 0,dtmean = 1,logdtsd = 0)))
        dd$id <- subi
        if(subi==1) d <- dd else d <- rbind(d,dd)
      }
      
      m <- ctModel(type='stanct',LAMBDA=matrix(1),CINT='cint',MANIFESTMEANS=0)
      # m$pars$indvarying=F
      
      f <- ctStanFit(datalong = d,ctstanmodel = m,cores=1,priors=TRUE,optimcontrol=list(carefulfit=F,stochastic=F,finishsamples=5000))
      s=summary(f)
      # print(s)
      
      scores=t(ctsem:::scorecalc(standata = f$standata,est = f$stanfit$rawest,stanmodel = f$stanmodel,subjectsonly = T,cores=1))
      
      fc <- ctStanFit(datalong = d,ctstanmodel = m,cores=1,priors=TRUE,optimcontrol=list(carefulfit=F,stochastic=F,finishsamples=5000,bootstrapUncertainty=F))
      sc=summary(fc)
      # print(sc)
      
      # print(sc$popmeans)
      # print(s$popmeans)
      colnames(sc$popmeans) <- paste0(colnames(sc$popmeans),'c')
      colnames(sc$rawpopcorr) <- paste0(colnames(sc$rawpopcorr),'c')
      colnames(sc$popsd) <- paste0(colnames(sc$popsd),'c')
      dt=data.frame(rbind(s$popmeans,s$popsd,s$rawpopcorr[,1:5]),rbind(sc$popmeans,sc$popsd,sc$rawpopcorr[,1:5]))
      rownames(dt)=c(rownames(s$popmeans),paste0('popsd_',rownames(s$popsd)),rownames(s$rawpopcorr))
      return(dt)
    }
    
    out <- list()
    for(i in 1:100){
      out[[i]] <- future(testfunc(i))
    }
    out <- value(out)
    
    
    truepars <- out[[1]][,'mean',drop=FALSE]
    truepars[] <- c(10,-1,2,.5,5,1,1.116,.45)
    
    out <- lapply(1:length(out),function(o) data.frame(run=o,par=rownames(out[[o]]),out[[o]]))
    outdt <- rbindlist(out)
    outdtb <- data.frame(outdt)[,colnames(outdt)[!grepl('c$',colnames(outdt))]]
    outdtc <- data.frame(outdt)[,c('run','par',colnames(outdt)[grepl('c$',colnames(outdt))])]
    colnames(outdtc)=gsub('c$','',colnames(outdtc))
    outdt <- rbind(data.table(type='bootstrap',outdtb),data.table(type='classic',outdtc))
    outdt <- melt(outdt,id.vars = c('type','par','run'))
    require(ggplot2)
    ggplot(outdt[variable %in% c('X97.5.','X2.5.','X50.'),],aes(x=variable,y=value,col=type,group=interaction(run,type)))+
      geom_jitter(alpha=.5,width=.2,height=0)+geom_line(alpha=.2)+theme_bw()+
      geom_hline(aes(yintercept=value),data=data.frame(par=rownames(truepars),value=truepars$mean))+
      facet_wrap(vars(par),scales = 'free')
    
    
    covered <- lapply(out,function(x) data.table(
      b= truepars > x$`X2.5.` & truepars < x$`X97.5.`, 
      c= truepars > x$`X2.5.c` & truepars < x$`X97.5.c`))
    
    coverage <- as.matrix(covered[[1]])
    for(i in 1:nrow(coverage)){
      for(j in 1:ncol(coverage)){
        coverage[i,j]<-sum(sapply(covered,function(x) as.matrix(x)[i,j]))/length(out)
      }
    }
    rownames(coverage) <- rownames(truepars)
    # print(coverage)
    
    
    
  }
}