if(identical(Sys.getenv("NOT_CRAN"), "true")& .Machine$sizeof.pointer != 4){ library(ctsem) library(testthat) cores = 2 context('misc') test_that("loo", { # for(cores in c(1,cores)){ #cores=6 data(AnomAuth) AnomAuthmodel<-ctModel(LAMBDA=matrix(c(1, 0, 0, 1), nrow=2, ncol=2), n.latent=2,n.manifest=2, # TIpredNames = 'Y1', # DRIFT=c('dr1||||Y1','dr12||||Y1','dr21||||Y1','dr22||||Y1'), MANIFESTVAR=diag(0,2), Tpoints=5) aa=ctDeintervalise(ctWideToLong(AnomAuth[1:500,], Tpoints = AnomAuthmodel$Tpoints,n.manifest = 2)) aa[4:20,AnomAuthmodel$manifestNames] <- NA sm <- ctStanModel(AnomAuthmodel,tipredDefault = FALSE) sm$pars$indvarying<- FALSE sf=ctStanFit(aa, ctstanmodel = sm, optimize=TRUE,verbose=0,savescores = FALSE,cores=cores) sdat <- sf$standata sdat$dokalmanrows[sdat$subject==1] <- 0L #remove 1 subject smf <- stan_reinitsf(sf$stanmodel,sdat) test_isclose( sf$stanfit$optimfit$value - rstan::log_prob(smf,sf$stanfit$rawest), #check ll equiv sum(sf$stanfit$transformedparsfull$llrow[sdat$subject==1])) loo=ctLOO(fit = sf,folds = 10,cores=cores,parallelFolds = T,subjectwise = T) loo2=ctLOO(fit = sf,folds = 10,cores=cores,parallelFolds = T,subjectwise = T) loo3=ctLOO(fit = sf,folds = sf$standata$nsubjects,cores=cores,parallelFolds = T,subjectwise = T,casewiseApproximation = T) if(F){ plot(loo3$outsampleLogLikRow-sf$stanfit$transformedparsfull$llrow[1,],loo$outsampleLogLikRow-sf$stanfit$transformedparsfull$llrow[1,]) abline(0,1) sd(apply(loo3$insampleLogLikRowFolds,2,sum,na.rm=T)) sd(apply(loo$insampleLogLikRowFolds,2,sum,na.rm=T)) l3=melt(data.table(ll=loo3$insampleLogLikeRowFolds,obs=1:nrow(loo3$insampleLogLikeRowFolds)),id.vars='obs') l3[,mean:=mean(value,na.rm=T),by=obs] l3[,sd:=sd(value,na.rm=T),by=obs] l=melt(data.table(ll=loo$insampleLogLikRowFolds,obs=1:nrow(loo$insampleLogLikRowFolds)),id.vars='obs') l[,mean:=mean(value,na.rm=T),by=obs] l[,sd:=sd(value,na.rm=T),by=obs] ggplot(l3[variable %in% variable[1],],aes(x=obs,y=mean-l[variable %in% variable[1],]$mean))+ geom_line()+ # geom_errorbar(aes(ymin=mean-sd,ymax=mean+sd))+ # geom_line(data=l[variable %in% variable[1],],aes(x=obs,y=mean),color='red')+ # geom_errorbar(data=l[variable %in% variable[1],],aes(ymin=mean-sd,ymax=mean+sd),color='red')+ theme_bw() cor(loo3$outsampleLogLikRow-sf$stanfit$transformedparsfull$llrow[1,],loo$outsampleLogLikRow-sf$stanfit$transformedparsfull$llrow[1,]) cor(loo2$outsampleLogLikRow-sf$stanfit$transformedparsfull$llrow[1,],loo$outsampleLogLikRow-sf$stanfit$transformedparsfull$llrow[1,]) pari = 6 range = range(c(loo$foldpars[pari,],loo3$foldpars[pari,])) plot(loo$foldpars[pari,],loo3$foldpars[pari,],xlim=range,ylim=range) } test_isclose( mean(loo2$outsampleLogLikRow-loo$outsampleLogLikRow), 0,tol=.01) test_isclose( mean(loo3$outsampleLogLikRow-loo$outsampleLogLikRow), 0,tol=.01) test_isclose( mean(abs(loo3$outsampleLogLikRow -loo$outsampleLogLikRow)), 0,tol=.05) # test_isclose( # mean(c(apply(loo$insampleLogLikRowFolds,1,sum,na.rm=TRUE)/9-sf$stanfit$transformedparsfull$llrow)), # 0, # tol=.1) # # test_isclose( # mean(c(apply(loo2$insampleLogLikRowFolds,1,sum,na.rm=TRUE)/9-sf$stanfit$transformedparsfull$llrow)), # 0, # tol=.1) # # # test_isclose( # mean(c(apply(loo3$insampleLogLikRowFolds,1,sum,na.rm=TRUE)/9-sf$stanfit$transformedparsfull$llrow)), # 0, # tol=.1) }) }