Nothing
if(identical(Sys.getenv("NOT_CRAN"), "true")& .Machine$sizeof.pointer != 4){
library(ctsem)
library(testthat)
cores=2
context("ukfpopcheck")
test_that("ukfpopcheck1", {
set.seed(2)
Tpoints<-10
n.latent=2
n.manifest=3
nsubjects=160
burnin=3
dtmat = matrix(exp(rnorm(burnin+Tpoints-1,-.3,.5)),1)
par1<-rnorm(nsubjects,-.3,.8)
par2 <- par1*.4 + rnorm(nsubjects,0,.3)
for(i in 1:nsubjects){
gm<-suppressMessages(ctModel(type='omx',n.latent=2,n.manifest=n.manifest,Tpoints=Tpoints,LAMBDA=matrix(c(1,.4,0,0,0,1),3,ncol=2),
DRIFT=diag(-.4,2),
CINT=matrix(c(par1[i],par2[i]),2),
T0VAR=diag(1,2),
T0MEANS=matrix(c(30,50),n.latent),
MANIFESTVAR=t(chol(diag(.5,n.manifest))),
DIFFUSION=t(chol(diag(3,2)))))
if(i==1) cd<-suppressMessages(ctGenerate(gm,n.subjects=1,burnin=burnin,wide=FALSE,dtmat = dtmat)) else {
newdat <- suppressMessages(ctGenerate(gm,n.subjects=1,burnin=burnin,wide=FALSE,dtmat = dtmat))
newdat[,'id'] <- i
cd<-rbind(cd,newdat)
}
}
m1<-ctModel(type='omx',n.latent=4,n.manifest=n.manifest,Tpoints=Tpoints,
LAMBDA=cbind(gm$LAMBDA,0,0),
MANIFESTMEANS=matrix(0,nrow=n.manifest)
)
m1$DRIFT[3:4,] <- 0
m1$DRIFT[,3:4] <- 0
diag(m1$DRIFT)[3:4] <- -1e-5
m1$DRIFT[1,3] = 1
m1$DRIFT[2,4] = 1
m1$DIFFUSION[3:4,] <- 0
m1$DIFFUSION[,3:4] <- 0
# m1$T0VAR[3:4,1:2] <- 0
#model for ukf ctsem
m2<-ctModel(type='omx',n.latent=2,n.manifest=n.manifest,Tpoints=Tpoints,
MANIFESTMEANS=matrix(0,n.manifest),
CINT=matrix(paste0('cint',1:gm$n.latent)),
LAMBDA=gm$LAMBDA
)
#fit 1
sm1 <- ctStanModel(m1)
sm1$pars$indvarying <- FALSE
sf1 <- ctStanFit(cd,sm1,cores=cores,verbose=0)
sf1d=sf1$stanfit$transformedparsfull$pop_DRIFT[1,1:2,1:2]
sf1ll=sf1$stanfit$optimfit$value
#fit 2
sm2 <- ctStanModel(m2)
sm2$pars$sdscale <- .2
sf2 <- ctStanFit(cd,sm2,cores=cores)
sf2d=sf2$stanfit$transformedparsfull$pop_DRIFT[1,1:2,1:2]
sf2ll=sf2$stanfit$optimfit$value
test_isclose(sf1ll,sf2ll,tol=1e-3)
test_isclose(sf1$stanfit$transformedparsfull$pop_T0cov[1,,]/
sf2$stanfit$transformedparsfull$popcov[1,,],
matrix(1,4,4),tol=1)
test_isclose(cov2cor(sf1$stanfit$transformedparsfull$pop_T0cov[1,,]),
cov2cor(sf2$stanfit$transformedparsfull$rawpopcov[1,,]),tol=1e-2)
test_isclose(sf1d,sf2d,tol=1e-2)
})
test_that("ukfpopcheck2_randomdrift", {
set.seed(1)
nsubjects=60
Tpoints=20
driftsd=.2
driftmu= log(.5)
dt=1
drift= -2*log1p(exp(2*(rnorm(nsubjects,driftmu,driftsd))))
for(si in 1:nsubjects){
m=suppressMessages(ctModel(LAMBDA=diag(1), Tpoints=Tpoints, DRIFT=matrix(drift[si]),T0MEANS = matrix(3),
T0VAR=matrix(sqrt(.1)), DIFFUSION=diag(1,1), CINT=matrix(-2),MANIFESTVAR=matrix(sqrt(.1))))
d=suppressMessages(ctGenerate(m,n.subjects = 1,burnin = 0,wide = FALSE,dtmean = dt))[,-1]
if(si==1) dat=cbind(si,d) else dat=rbind(dat,cbind(si,d))
}
colnames(dat)[1]='id'
dtm <- ctModel(LAMBDA=diag(1), type='stanct',
DRIFT=matrix('dr11| -2*log1p(exp(-2*param))'),
CINT=matrix('cint'),
MANIFESTMEANS = matrix(0)
)
dtm$pars$indvarying <- FALSE
dtm$pars$indvarying[dtm$pars$matrix %in% 'DRIFT'] <- TRUE
dtf = ctStanFit(datalong = dat,ctstanmodel = dtm,optimize=TRUE,
verbose=0,optimcontrol=list(estonly=F),savescores = F)
s1=summary(dtf,parmatrices = F,priorcheck = F,residualcov = F)
dtm2 <- ctModel(LAMBDA=matrix(c(1,0),1,2), type='stanct',
DIFFUSION=matrix(c('diff',0,0,0),2,2),DRIFT=matrix(c('-2*log1p(exp(-2*state[2]))',0,0,-.00001),2,2),
CINT=matrix(c('cint',0),2,1),T0MEANS=matrix(c('t0m1','t0m2|param'),2,1),
T0VAR=matrix(c('t0v11',0,0,paste0('t0var22 | ',gsub(".* sdscale","",gsub('rawpopsdbase','param',dtm$rawpopsdtransform),fixed=TRUE))),2,2),
MANIFESTMEANS = matrix(0))
dtm2$pars$indvarying <- FALSE
dtf2=ctStanFit(datalong = dat,ctstanmodel = dtm2,optimize = TRUE)
s2=summary(dtf2,parmatrices = F,priorcheck = F,residualcov = F)
test_isclose(s1$ll,s2$ll,tol=1e-3)
test_isclose(sort(dtf2$stanfit$rawest),sort(dtf$stanfit$rawest),tol=1e-3) #sorting is an ugly hack! could improve...
})
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.