So, I don’t know what’s wrong with your stan model but a) stan can probably fit this ok as far as I see, but b) the problem seems very sensitive to initial conditions. after rescaling the processes to match the observations (dividing by area), ctsem fit this ok with an extended kalman filter, using both optimization and stan’s HMC, and also fit ok when disabling the kalman filter and sampling the process states.
library(ctsem)
Simulated_Data$id=1 #needs an id column, designed for multiple subjects
Simulated_Data$C1s=Simulated_Data$C1/Area
Simulated_Data$C2s=Simulated_Data$C2/Area
m <- ctModel(time='YEAR',
type='standt', #discrete time model
manifestNames = c('OJ','OS','OA'),
latentNames = c('J','S','A'),
DRIFT=0, #set linear temporal dependencies to zero
CINT=c( #specify non linear change terms
'alpha*S*exp(-beta*S)',
'A * exp(-M/2)',
'S * exp(-M/2) + J*exp(-M)'),
LAMBDA=diag(0,3), #set linear observation model to zero,
MANIFESTMEANS=c( 'J', 'S', 'A'),
T0VAR=diag(.0001,3), #fix initial system noise to small value,
DIFFUSION=0, #set ongoing system noise to 0
T0MEANS=c( #estimate initial state values, specify prior based on standard normal
'j1|param*.01+.02',
's1|param*.01+.01',
'a1|param*.01+.02'),
TDpredNames =c('C1s','C2s'), #time dependent predictors
TDPREDEFFECT = c( #fix linear effect of time dependent (catch) predictors
0,0,
-1,0,
0,-1),
MANIFESTVAR=c( #specify observation noise matrix (sd's on diag)
'Jsd',0,0,
0,'Ssd',0,
0,0,'Asd'),
PARS=c('alpha|param+10','beta|param+100','M|param+1') #specify any parameters included in equations
)
ctModelLatex(m)
f <- ctStanFit(datalong = Simulated_Data,ctstanmodel = m,
cores=8,chains=8,optimize=F,iter=1000,nopriors = F,
intoverstates = T, #using extended kalman filter by default, turn off to sample latent states
plot=T)
ctKalman(f,plot=T,removeObs = 10)+coord_cartesian(ylim=c(0,.05)) #forward predictions keeping only every 10th obs
summary(f)
