Nothing
expmGetSubsets <- function(m){
nr = nrow(m)
m <- matrix(suppressWarnings(as.numeric(m)),nr,nr)
m[is.na(m)] <- 999
subsets <- list()
indices <- 1:nrow(m)
indexmat <- 1:nrow(m)
while(length(indices) >0){
oldsubsetsize <- -1
subindices <- indices[1]
while(length(subindices) > oldsubsetsize){
oldsubsetsize <- length(subindices)
subindices <- unique(c(subindices,
which(m!=0 & col(m) %in% subindices,arr.ind = TRUE),
which(m!=0 & row(m) %in% subindices,arr.ind = TRUE)))
}
subsets[[length(subsets)+1]] <- subindices
indices <- indices[!indices %in% subindices]
}
subsets <- lapply(subsets, function(x) c(x,rep(0,nr-length(x))))
subsets <- t(array(as.integer(unlist(subsets)),dim =c(nr,length(subsets))))
return(subsets)
}
ctStanData <- function(ctm, datalong,optimize,sameInitialTimes=FALSE){
nsubjects <- length(unique(datalong[, ctm$subjectIDname]))
mats <- ctStanMatricesList()
#simply exponential?
driftdiagonly <- ifelse(all(!is.na(ctm$pars$value[ctm$pars$matrix == 'DRIFT' & ctm$pars$row != ctm$pars$col]) &
all(ctm$pars$value[ctm$pars$matrix == 'DRIFT' & ctm$pars$row != ctm$pars$col] == 0) ), 1, 0)
###data checks
if(any(!c(ctm$manifestNames,ctm$TDpredNames,ctm$TIpredNames) %in% colnames(datalong))) stop(paste0('
variables: ', paste0(c(ctm$manifestNames,ctm$TDpredNames,ctm$TIpredNames)[
which(!c(ctm$manifestNames,ctm$TDpredNames,ctm$TIpredNames) %in% colnames(datalong))], ', '),' not in data'))
if(!(ctm$subjectIDname %in% colnames(datalong))) stop(paste('id column', (ctm$subjectIDname), "not found in data"))
if(ctm$priors){
if(ctm$n.TIpred > 1 && any(abs(colMeans(datalong[,c(ctm$TIpredNames),drop=FALSE],na.rm=TRUE)) > .3)){
message('Uncentered TI predictors noted -- interpretability may be hindered and default priors may not be appropriate')
}
# #scale check
# if(naf(any(abs(colMeans(datalong[,c(ctm$manifestNames,ctm$TDpredNames),drop=FALSE],na.rm=TRUE)) > 5))){
# message('Uncentered data noted -- default priors *may* not be appropriate')
# }
#
#
# if(naf(any(abs(apply(datalong[,c(ctm$manifestNames,ctm$TDpredNames,ctm$TIpredNames),drop=FALSE],2,sd,na.rm=TRUE)) > 3))){
# message('Unscaled data noted -- default priors may not be appropriate')
# }
}
#id mapping
original <- unique(datalong[,ctm$subjectIDname])
datalong <- makeNumericIDs(datalong,ctm$subjectIDname,ctm$timeName)
new <- unique(datalong[,ctm$subjectIDname])
idmap <- data.frame(original, new)
if( (nrow(ctm$t0varstationary) + nrow(ctm$t0meansstationary)) >0 &&
length(c(ctm$modelmats$calcs$driftcint, ctm$modelmats$calcs$diffusion)) > 0) message('Stationarity assumptions based on initial states when using non-linear dynamics')
nindvarying <- max(ctm$modelmats$matsetup$indvarying)
nparams <- max(ctm$modelmats$matsetup$param[ctm$modelmats$matsetup$when %in% c(0,-1)])
nmatrices <- length(mats$base)
ctm$modelmats$matsetup[which(ctm$modelmats$matsetup$indvarying > 0),]
indvaryingindex <- ctm$modelmats$matsetup$param[which(ctm$modelmats$matsetup$indvarying > 0)]
indvaryingindex <- array(indvaryingindex[!duplicated(indvaryingindex)])
sdscale <- array(ctm$modelmats$matvalues$sdscale[match(indvaryingindex,ctm$modelmats$matsetup$param)])
#tipred data
if(ctm$n.TIpred > 0) {
tipreds <- datalong[match(unique(datalong[,ctm$subjectIDname]),datalong[,ctm$subjectIDname]),ctm$TIpredNames,drop=FALSE]
if(any(is.na(tipreds))) {
if(!optimize){
message(paste0("NA's in TIpreds - sampling ", sum(is.na(tipreds)),' values'))
tipreds[is.na(tipreds)] = 99999
}
if(optimize){
message(paste0("NA's in TIpreds - imputing ", sum(is.na(tipreds)),' NA\'s to allow optimization -- TIpred effect estimates may be overly confident.'))
# tipreds[is.na(tipreds)] = 0
timu <- apply(tipreds,2,mean,na.rm=TRUE)
tisd <- apply(tipreds,2,sd,na.rm=TRUE)
meandat <- data.table((datalong))[ , lapply(.SD, function(x)
sum(x,na.rm=TRUE)/sum(!is.na(x))) ,
# x[1]),
by=c(ctm$subjectIDname)]
# sddat <- data.table((datalong))[ , lapply(.SD, function(x)
# sd(x,na.rm=TRUE)) ,
# by=c(ctm$subjectIDname)]
# sddat<-sddat[,!colnames(sddat) %in% ctm$subjectIDname,with=FALSE]
meandat <- data.frame(scale(meandat[,apply(meandat,2,sd,na.rm=TRUE) > 1e-4,with=FALSE]))
meandat[is.na(meandat)] <- 0
# sddat <- sddat[,apply(sddat,2,sd,na.rm=TRUE) > 1e-4,with=FALSE]
# colnames(sddat) <- paste0('sd_',colnames(sddat))
# meandat <- cbind(meandat,sddat)
# cml <- covml(meandat,reg = TRUE)
for(i in 1:ctm$n.TIpred){
# cnames <- c(colnames(meandat)[!colnames(meandat) %in% ctm$TIpredNames[i]],ctm$TIpredNames[i])
#
#
# chol <- cml$cp$covm
# dimnames(chol) <- list(colnames(meandat),colnames(meandat))
# chol <- t(chol(chol[cnames,cnames]))
#
# tipreds[is.na(tipreds[,i]),i] <-
# apply(meandat[is.na(tipreds[,i]),head(cnames,length(cnames)-1)],1,function(x){
# sum(x * c(tail(chol,1))[-ncol(chol)])
# }) *tisd[i] + timu[i]
lmform = formula(paste0(ctm$TIpredNames[i],' ~ 1 + ',
paste0(colnames(meandat)[-which(colnames(meandat) %in% ctm$TIpredNames[i])],
collapse=' + ')))
# lmr <- MASS::lm.ridge(formula = lmform, data = meandat,lambda=0,na.action=na.exclude)
lmf <- lm(formula = lmform,data = meandat,na.action=na.exclude)
# lmf$coefficients[-1] <- lmr$coef
# lmf$coefficients[1] <- lmr$ym
# lmr$coef *
# apply(meandat[,-which(colnames(meandat) %in% ctm$TIpredNames[i]),with=FALSE],
# 2,sd,na.rm=TRUE) -
# lmf$coefficients[-1]
#
# plot(meandat[,ctm$TIpredNames[i]],predict(lmf),main=ctm$TIpredNames[i],col=as.numeric(is.na(tipreds[,i]))+1)
tipreds[is.na(tipreds[,i]),i] <- predict(lmf)[is.na(tipreds[,i])] * tisd[i] + timu[i]
}
}
}
}
if(sameInitialTimes){
datanew <- datalong[!duplicated(datalong[[ctm$subjectIDname]]),]
datanew[[ctm$timeName]] <- min(datanew[[ctm$timeName]])
datalong <- merge.data.frame(x = datalong,y = datanew[,c(ctm$subjectIDname,ctm$timeName)],
by = c(ctm$subjectIDname,ctm$timeName),all = TRUE)
}
datalong[,ctm$manifestNames][is.na(datalong[,ctm$manifestNames])]<-99999 #missing data
standata <- list(
Y=cbind(as.matrix(datalong[,ctm$manifestNames])),
subject=array(as.integer(datalong[,ctm$subjectIDname])),
time=array(datalong[,ctm$timeName]),
ndatapoints=as.integer(nrow(datalong)),
nobs_y=array(as.integer(apply(datalong[,ctm$manifestNames,drop=FALSE],1,function(x) length(x[x!=99999]))),dim=nrow(datalong)),
whichobs_y=matrix(as.integer(t(apply(datalong[,ctm$manifestNames,drop=FALSE],1,function(x) {
out<-as.numeric(which(x!=99999))
if(length(out)==0) out<-rep(0,ctm$n.manifest)
if(length(out)<ctm$n.manifest) out<-c(out,rep(0,ctm$n.manifest-length(out)))
out
}) )),nrow=c(nrow(datalong),ncol=ctm$n.manifest)),
nbinary_y=array(as.integer(apply(datalong[,ctm$manifestNames,drop=FALSE],1,function(x)
length(x[ctm$manifesttype==1 & x!=99999]))),dim=nrow(datalong)),
whichbinary_y=matrix(as.integer(t(apply(datalong[,ctm$manifestNames,drop=FALSE],1,function(x) {
out<-as.numeric(which(ctm$manifesttype==1 & x!=99999)) #conditional on whichobs_y
if(length(out)==0) out<-rep(0,ctm$n.manifest)
if(length(out)<ctm$n.manifest) out<-c(out,rep(0,ctm$n.manifest-length(out)))
out}) )),nrow=c(nrow(datalong),ncol=ctm$n.manifest)),
# nordinal_y=array(as.integer(apply(datalong[,ctm$manifestNames,drop=FALSE],1,function(x)
# length(x[ctm$manifesttype==2 & x!=99999]))),dim=nrow(datalong)),
# whichordinal_y=matrix(as.integer(t(apply(datalong[,ctm$manifestNames,drop=FALSE],1,function(x) {
# out<-as.numeric(which(ctm$manifesttype==2 & x!=99999)) #conditional on whichobs_y
# if(length(out)==0) out<-rep(0,ctm$n.manifest)
# if(length(out)<ctm$n.manifest) out<-c(out,rep(0,ctm$n.manifest-length(out)))
# out}) )),nrow=c(nrow(datalong),ncol=ctm$n.manifest)),
ncont_y=array(as.integer(apply(datalong[,ctm$manifestNames,drop=FALSE],1,function(x)
length(x[(ctm$manifesttype==0) & x!=99999]))),dim=nrow(datalong)),
whichcont_y=matrix(as.integer(t(apply(datalong[,ctm$manifestNames,drop=FALSE],1,function(x) {
out<-as.numeric(which( (ctm$manifesttype==0) & x!=99999)) #conditional on whichobs_y
if(length(out)==0) out<-rep(0,ctm$n.manifest)
if(length(out)<ctm$n.manifest) out<-c(out,rep(0,ctm$n.manifest-length(out)))
out
}) )),nrow=c(nrow(datalong),ncol=ctm$n.manifest))
)
if(ctm$n.TDpred > 0) {
tdpreds <- datalong[,ctm$TDpredNames,drop=FALSE]
if(any(is.na(tdpreds))) {
# if(NAtdpreds == 'error'
message("NA's in TDpreds! Replaced by zeroes, consider appropriateness...")
tdpreds[is.na(tdpreds)] <-0 ## rough fix for missingness
}
}
if(ctm$n.TDpred ==0) tdpreds <- matrix(0,standata$ndatapoints,0) #standata$ndatapoints,
standata$tdpreds=array(as.matrix(tdpreds),dim=c(nrow(tdpreds),ncol(tdpreds)))
#subset selection
if(is.null(ctm$dokalmanrows)) standata$dokalmanrows <-
array(rep(1L, standata$ndatapoints)) else standata$dokalmanrows <- array(as.integer(ctm$dokalmanrows))
standata$priormod = 1L#sum(standata$dokalmanrows)/standata$ndatapoints
standata<-c(standata,
list(
nsubjects=as.integer(nsubjects),
nmanifest=as.integer(ctm$n.manifest),
nlatentpop = ctm$nlatentpop,
Jstep = ctm$nlcontrol$Jstep,
maxtimestep = ctm$nlcontrol$maxtimestep,
dokalman=as.integer(is.null(ctm$dokalman)),
intoverstates=as.integer(ctm$intoverstates),
verbose=0L,
manifesttype=array(as.integer(ctm$manifesttype),dim=length(ctm$manifesttype)),
# nordinal=as.integer(sum(ctm$manifesttype==2)),
indvaryingindex=array(as.integer(indvaryingindex)),
intoverpopindvaryingindex=array(as.integer(ctm$intoverpopindvaryingindex)),
notindvaryingindex=array(as.integer(which(!(1:nparams) %in% indvaryingindex))),
continuoustime=as.integer(sum(ctm$continuoustime)),
nlatent=as.integer(ctm$n.latent),
ntipred=as.integer(ctm$n.TIpred),
ntdpred=as.integer(ctm$n.TDpred),
nparams=as.integer(nparams),
gendata=0L,
nindvarying=as.integer(nindvarying),
nindvaryingoffdiagonals=as.integer((nindvarying^2-nindvarying)/2),
nt0varstationary=as.integer(nrow(ctm$t0varstationary)),
nt0meansstationary=as.integer(nrow(ctm$t0meansstationary)),
t0varstationary=matrix(as.integer(ctm$t0varstationary),ncol=2),
t0meansstationary=matrix(as.integer(ctm$t0meansstationary),ncol=2),
driftdiagonly = as.integer(driftdiagonly),
intoverpop=as.integer(ctm$intoverpop),
# nlmeasurement=as.integer(nlmeasurement),
priors=as.integer(ctm$priors),
savescores=0L,
savesubjectmatrices=0L,
nJAxfinite = ifelse(ctm$recompile,0L,length(ctm$JAxfinite)),
nJyfinite = ifelse(ctm$recompile,0L,length(ctm$Jyfinite))
))
if(ctm$n.TIpred == 0) tipreds <- array(0,c(0,0))
standata$tipredsdata <- as.matrix(tipreds)
standata$nmissingtipreds <- as.integer(length(tipreds[tipreds== 99999]))
standata$ntipredeffects <- as.integer(ifelse(ctm$n.TIpred > 0, as.integer(max(ctm$modelmats$TIPREDEFFECTsetup)), 0))
standata$TIPREDEFFECTsetup <- apply(ctm$modelmats$TIPREDEFFECTsetup,c(1,2),as.integer,.drop=FALSE)
standata$tipredsimputedscale <- ctm$tipredsimputedscale
standata$tipredeffectscale <- ctm$tipredeffectscale
# ##ordinal categories
# # if(any(ctm$manifesttype==2)){
# standata$ncategories=array(as.integer(apply(datalong[,ctm$manifestNames,drop=FALSE][,ctm$manifesttype==2,drop=FALSE],2,function(x) length(unique(x)))))
# standata$nthresholdpars=as.integer(sum(standata$ncategories-1))
# # } else standata$ncategories=array(integer())
# if(as.integer(ctm$NOrdinalIntegrationPoints) %% 2 == 0) stop('NOrdinalIntegrationPoints must be odd')
# standata$nordinalintegrationpoints=as.integer(ctm$NOrdinalIntegrationPoints)
# midpoint = ceiling(standata$nordinalintegrationpoints/2) #ensure midpoint (mean) is included first, because catprobs is only recomputed for non mean points.
# standata$sigmapoints=array(round(statmod::gauss.quad.prob(n=ctm$NOrdinalIntegrationPoints,dist='normal')$nodes[
# c(midpoint,(1:standata$nordinalintegrationpoints)[-midpoint])],5))
# standata$sigweights=array(round(statmod::gauss.quad.prob(n=ctm$NOrdinalIntegrationPoints,dist='normal')$weights[
# c(midpoint,(1:standata$nordinalintegrationpoints)[-midpoint])],5))
# #drift, jacobian off diagonal check # seems to work ok internally with expm2
# mx=listOfMatrices(ctm$pars)
# driftcint <- rbind(cbind(mx$DRIFT[1:ctm$n.latent,1:ctm$n.latent,drop=FALSE],mx$CINT),0)
# z=c()
# for(ri in 1:nrow(driftcint)){
# if(all(driftcint[ri,-ri] %in% 0) && all(driftcint[-ri,ri] %in% 0)) z[ri]=0L else z[ri]=1L
# }
# standata$drcintoffdiag <- array(as.integer(c(z,1)),dim=nrow(driftcint))
#
# jac <- mx$JAx
# z=c()
# for(ri in 1:nrow(jac)){
# if(all(jac[ri,-ri] %in% 0) && all(jac[-ri,ri] %in% 0)) z[ri]=0L else z[ri]=1L
# }
# standata$jacoffdiag <- array(as.integer(c(z,1)),dim=nrow(jac))
# standata$jacoffdiagindex <- array(as.integer(sort(unique(c(1:ctm$n.latent,which(standata$jacoffdiag ==1))))))
# standata$njacoffdiagindex <- as.integer(length(standata$jacoffdiagindex))
standata$DRIFTsubsets <- expmGetSubsets(listOfMatrices(ctm$pars)$DRIFT)
standata$JAxsubsets <- expmGetSubsets(listOfMatrices(ctm$pars)$JAx)
standata$nDRIFTsubsets <- nrow(standata$DRIFTsubsets)
standata$nJAxsubsets <- nrow(standata$JAxsubsets)
tempdiff <- listOfMatrices(ctm$pars)$DIFFUSION
derrind <- unique(c(which(tempdiff!=0,arr.ind = TRUE))) #which latents have some non zero diffusion
derrind <- unique(c(standata$JAxsubsets[ #and which expm subsets are these part of
apply(standata$JAxsubsets,1,function(x) any(derrind %in% x)),]))
derrind <- array(sort(derrind[derrind <= standata$nlatent & derrind > 0])) #exclude any elements that are stable individual differences
standata$derrind <- derrind
standata$ndiffusion <- length(derrind)
standata$JAxfinite <- ctm$JAxfinite
standata$Jyfinite <- ctm$Jyfinite
standata$Jycolindexsize <- 1L
standata$Jycolindex <- array(1L)
standata$difftype <- 0L;
standata$dotipred <- 1L;
# #t0 index
# T0check<-rep(1,nrow(datalong))
# for(i in 2:nrow(datalong)){
# T0check[i]<- ifelse(datalong[i,ctm$subjectIDname] != datalong[i-1,ctm$subjectIDname], 1, 0)
# }
if (!(ctm$timeName %in% colnames(datalong))) stop(paste('time column', (ctm$timeName), "not found in data"))
if(any(is.na(datalong[,ctm$timeName]))) stop('Missings in time column!')
oldsubi<-datalong[1,ctm$subjectIDname]-1
dT<-rep(-1,length(datalong[,ctm$timeName]))
dt0rows <- c()
for(rowi in 1:length(datalong[,ctm$timeName])) {
subi<-datalong[rowi,ctm$subjectIDname]
# if(rowi==1 && subi!=1) stop('subject id column must ascend from 1 to total subjects without gaps')
# if(oldsubi!=subi && subi-oldsubi!=1) stop('subject id column must ascend from 1 to total subjects without gaps')
if(subi - oldsubi == 1) {
dT[rowi]<-0
subistartrow<-rowi
}
if(subi - oldsubi == 0) {
dT[rowi]<-datalong[rowi,ctm$timeName] - datalong[rowi-1,ctm$timeName]
if(dT[rowi] < 0) stop(paste0('A time interval of ', dT[rowi],' was found at row ',rowi))
if(dT[rowi] == 0) dt0rows <- c(dt0rows,rowi)
}
if(length(dt0rows)>0) warning(paste0('A time interval of 0 was found at row/s: ',paste0(dt0rows,collapse=', ')))
oldsubi<-subi
}
# if(mean(dT) > 3) message('Average time interval > 3 -- in typical settings this can be too large, consider time scaling...')
subindices <- lapply(mats$base,function(x) 0)
for(mati in mats$base){
if( (!ctm$intoverpop && any(ctm$pars$indvarying[ctm$pars$matrix %in% names(mats$base)[mati]])) ||
(ctm$n.TIpred >0 && (
any(unlist(ctm$pars[ctm$pars$matrix %in% names(mats$base)[mati],paste0(ctm$TIpredNames,'_effect')])) ||
any(ctm$pars$matrix %in% names(mats$base)[mati] & grepl('[',ctm$pars$param,fixed=TRUE)) )
)) subindex <- 1 else subindex <- 0
subindices[[names(mats$base)[mati]]] <- subindex
}
if(ctm$stationary || nrow(ctm$t0varstationary) > 0) subindices$T0VAR <-
max(c(subindices$T0VAR,subindices$DRIFT,subindices$DIFFUSION))
if(ctm$stationary || nrow(ctm$t0meansstationary) > 0) subindices$T0MEANS <-
max(c(subindices$T0MEANS,subindices$DRIFT,subindices$CINT))
standata$subindices <- as.integer(unlist(subindices))[order(mats$base)]
#state dependence
statedep <- rep(0L,as.integer(max(mats$all)))
lhscalcs <- sapply(unique(unlist(ctm$modelmats$calcs)),function(x) gsub('=.*','',x))
for(i in 1:length(statedep)){
matname <- try(names(mats$all[mats$all %in% i]),silent=TRUE)
if(length(matname)==0) next
statedep[i] <- ifelse(any(
ctm$modelmats$matsetup$when >0 &
ctm$modelmats$matsetup$matrix %in% i),1L,0L)
if(any(sapply(lhscalcs,function(calci)
grepl(matname,calci)))) statedep[i] <- 1L #if any calcs modify this matrix, make state dependent
}
#
# statedep=rep(0L,4)
# v = 0L
# if(any(ctm$modelmats$matsetup$when == 2) ||
# any(grepl('state[',ctm$jacobian$JAx,fixed=TRUE)) ||
# any(grepl('state[',ctm$modelmats$calcs$driftcint,fixed=TRUE)) ||
# any(grepl('state[',ctm$modelmats$calcs$diffusion,fixed=TRUE))
# ) statedep[2] = 1L
#
# if(any(ctm$modelmats$matsetup$when == 2 & ctm$modelmats$matsetup$matrix == 4) ||
# any(grepl('state[',ctm$modelmats$calcs$diffusion,fixed=TRUE))
# ) multiplicativenoise = 1L
standata$statedep <- statedep
# standata$nstatedep <- as.integer(length(statedep))
# standata$multiplicativenoise <- multiplicativenoise
standata$choleskymats<- 0L
if(!is.null(ctm$covmattransform)){
if(ctm$covmattransform=='rawcorr_indep') standata$choleskymats<- -1L
if(ctm$covmattransform=='cholesky') standata$choleskymats<- 1L
if(!ctm$covmattransform %in% c('rawcorr','rawcorr_indep','cholesky')) stop(
'covtransform must be either "rawcorr", "rawcorr_indep", or "cholesky"')
}
standata$matsetup <- apply(ctm$modelmats$matsetup[,-1],c(1,2),as.integer,.drop=FALSE) #remove parname and convert to int
standata$matvalues <- apply(ctm$modelmats$matvalues,c(1,2),as.numeric)
standata$nmatrices <- as.integer(nmatrices)
standata$matrixdims <- ctm$modelmats$matrixdims
standata$nrowmatsetup <- as.integer(nrow(ctm$modelmats$matsetup))
standata$sdscale <- array(as.numeric(sdscale),dim=length(sdscale))
standata$rawpopcovscale=1
standata$approxct <- 0L
if(!is.null(ctm$approxct)) standata$approxct <- as.integer(ctm$approxct)
standata$taylorheun <- 0L
if(!is.null(ctm$taylorheun)) standata$taylorheun <- as.integer(ctm$taylorheun)
#fixed hyper pars #ow disabled
# if(!is.null(ctm$fixedrawpopchol)) {
# standata$fixedrawpopmeans = array(ctm$fixedrawpopmeans)
# standata$fixedrawpopchol= ctm$fixedrawpopchol
# }
# standata$fixedhyper <- as.integer(ifelse(is.null(ctm$fixedrawpopchol),0,1))
# if(is.null(ctm$fixedrawpopchol)) {
# standata$fixedrawpopmeans = array(0,dim = 0)
# standata$fixedrawpopchol = matrix(0,0,0)
# }
# standata$fixedsubpars <- as.integer(!is.null(ctm$fixedsubpars))
# if(!is.null(ctm$fixedsubpars)) standata$fixedindparams <-
# ctm$fixedsubpars else standata$fixedindparams <-array(0,dim=c(0,0))
standata$idmap <- idmap
standata$popcovn=1000L
standata$llsinglerow=0L
# standata$doonesubject=0L
if(any(ctm$manifesttype > 0)){ #if any non gaussian observables, just use finite differences for measurement jacobian
standata$nJyfinite <- as.integer(ctm$n.latent)
standata$Jyfinite <- array(as.integer(1:ctm$n.latent))
}
# } else{
# standata$nJyfinite <- 0L
# standata$Jyfinite <- integer()
# }
if(!is.null(ctm$TIpredAuto) && ctm$TIpredAuto %in% c(1L,TRUE)) standata$TIpredAuto <- 1L else standata$TIpredAuto <- 0L
mc=c(ctStanMatricesList()$all)#base,ctStanMatricesList()$jacobian)
ms=data.frame(standata$matsetup)
ms=ms[order(ms$param),]
standata$whenmat <- array(0L,dim=c(max(mc),5)) #whenmat contains 0's when matrix isn't computed, 1's when it is. 'when 5' is indvaryig.
for(mi in mc){
mrows=which(ms$matrix==mi)
for(wheni in 1:4){
standata$whenmat[mi,wheni] = as.integer(any(
ms$when[mrows] %in% c(wheni,100))) #100 for PARS, find a better approach to work out exactly when PARS are needed
}
standata$whenmat[mi,5] = as.integer(any(
ms$param[mrows] > 0 & ms$when[mrows] %in% c(0,100) & (ms$indvarying[mrows] > 1 | ms$tipred[mrows] > 0) ))
}
rownames(standata$whenmat)[mc] <- names(mc)
if(optimize){ #indvarying elements need to be updated dynamically
for(ri in 1:nrow(standata$whenmat)){
if(standata$whenmat[ri,5]>0){ #if any of this matrices pars are indvarying
standata$whenmat[mi,
switch(rownames(standata$whenmat)[ri],PARS=5,T0MEANS=5,LAMBDA=4,DRIFT=2,DIFFUSION=2,MANIFESTVAR=4,MANIFESTMEANS=4,CINT=2,T0VAR=5,TDPREDEFFECT=3,
JAx=2,Jtd=3,Jy=4,asymCINT=2,asymDIFFUSIONcov=2,DIFFUSIONcov=2,MANIFESTcov=4,T0cov=5)] <- 1L
# standata$whenmat[ri,5] <- 0 #set indvarying to zero because updating dynamically -- don't do this not all are dynamically updated.
}
}
}
#this PARS when = 100 thing is annoyinh, improve it...
standata$whenvecp <- array(0L, c(2,standata$nparams)) #whenvecp contains 0's for unchanging pars
standata$whenvecp[1,] <- as.integer(1:standata$nparams) #base parameters
standata$whenvecp[2,ms$param[ms$when %in% c(0,100) & ms$copyrow <1 & (ms$tipred > 0 | ms$indvarying > 0) & ms$param > 0]] <-
as.integer(ms$param[ms$when %in% c(0,100) & ms$copyrow <1 & (ms$tipred > 0 | ms$indvarying > 0) & ms$param > 0])
# standata$whenvecp[3,] <- as.integer(1:ncol(standata$whenvecp))
standata$whenvecs <- array(0L,dim=c(6,standata$nlatentpop)) #when do we need to compute transformed states?
for(wheni in 1:4){ #whenvecs specifies array of when values for the state transform vector -- 1 to compute, 0 not
standata$whenvecs[wheni,ms$param[ms$when %in% c(wheni,100) & ms$copyrow <=0 & ms$param > 0]] <-
as.integer(ms$param[ms$when %in% c(wheni,100) & ms$copyrow <=0 & ms$param > 0]) #100 is for PARS - needed everywhere.
}
# why was this in the code? when do we need the below line... ind varying based on states?
# if(standata$intoverpop==1 && standata$nlatentpop > standata$nlatent){
# standata$whenvecs[5,ms$param[ms$when==0 & ms$param > 0 & ms$copyrow < 1 & (ms$indvarying > 0 | ms$tipred > 0)]] <-
# as.integer(ms$param[ms$when==0 & ms$param > 0 & ms$copyrow < 1 & (ms$indvarying > 0 | ms$tipred > 0)])
# }
# standata$whenvecs[6,] <- as.integer(1:ncol(standata$whenvecs))
#special matrix adjustments
standata$whenmat[mc[names(mc) == 'asymDIFFUSIONcov'],] <-
apply(standata$whenmat[ mc[names(mc) %in% c('DIFFUSION','DRIFT')],],2,max)
standata$whenmat[mc[names(mc) == 'asymCINT'],] <-
apply(standata$whenmat[mc[names(mc) %in% c('CINT','DRIFT')],],2,max)
standata$whenmat[mc[names(mc) == 'DIFFUSIONcov'],] <- standata$whenmat[mc[names(mc) == 'DIFFUSION'],]
standata$whenmat[mc[names(mc) == 'T0cov'],] <- standata$whenmat[mc[names(mc) == 'T0VAR'],]
standata$whenmat[mc[names(mc) == 'MANIFESTcov'],] <- standata$whenmat[mc[names(mc) == 'MANIFESTVAR'],]
standata$statedep[31:33] <- standata$statedep[c(4,5,8)]
#laplace priors
standata$laplaceprior <- array(rep(0L,standata$nparams))
standata$laplacetipreds <- 0L
if(!is.null(ctm$laplaceprior)){
if('tipreds' %in% ctm$laplaceprior) standata$laplacetipreds <- 1L
ms <- data.frame(standata$matsetup)
standata$laplaceprior[
ms$param[
ms$matrix %in% ctStanMatricesList()$all[names(ctStanMatricesList()$all) %in% ctm$laplaceprior] &
ms$param > 0 &
ms$row!=ms$col &
ms$when==0 &
ms$copyrow<1]
] <- 1L
}
standata$laplaceprioronly <- ifelse(is.null(ctm$laplaceprioronly),0L,as.integer(ctm$laplaceprioronly))
#CINT non zero
ms <- data.frame(standata$matsetup)
CINTnonzero <- c()#1:standata$nlatent
for(i in 1:standata$nlatent){
ri=which(ms$matrix %in% ctStanMatricesList()$all[names(ctStanMatricesList()$all) %in% 'CINT'] & ms$row %in% i)
if(ms$param[ri]==0 && standata$matvalues[ri,'value']==0) next else CINTnonzero <- c(CINTnonzero,i)
}
standata$CINTnonzero <- array(as.integer(CINTnonzero))
standata$CINTnonzerosize <- length(CINTnonzero)
#JAx different from DRIFT?
standata$JAxDRIFTequiv <- 1L
ml <- listOfMatrices(ctm$pars)
for(i in 1:nrow(ml$JAx)){
for(j in 1:ncol(ml$JAx)){
if(i <= nrow(ml$DRIFT) && j <= nrow(ml$DRIFT)){ #check drift equivalence
# if(i==j && !ctm$continuoustime && ml$JAx[i,j] %in% 1) ml$JAx[i,j] <-
if(ml$JAx[i,j] != ml$DRIFT[i,j]) standata$JAxDRIFTequiv <- 0L
}
if(i > nrow(ml$DRIFT) || j > nrow(ml$DRIFT)){ #check drift equivalence
if(i != j && !ml$JAx[i,j] %in% 0) standata$JAxDRIFTequiv <- 0L
if(i == j && !ml$JAx[i,j] %in% ifelse(ctm$continuoustime,0,1)) standata$JAxDRIFTequiv <- 0L
}
}
}
standata$recompile <- as.integer(ctm$recompile)
standata$nsubsets <- 1L
return(standata)
}
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.