DynamicalCompPsych.R

#generate data from python script

library('reticulate')
# reticulate::install_python()

a=py_run_string("
import matplotlib.pyplot as plt
import numpy as np
from random import random
from random import seed
import matplotlib.pyplot as plt
plt.rcParams.update({'font.size': 28})



#Bipolar cycle rapide:
Par={'Smax':10, 'Rs':1, 'lambdas':0.1, 'taux':14, 'P': 10, 'Rb':1.04, 'lambdab': 0.05, 'L':1.01, 'tauy':14, 'S':10, 'alpha':0.5, 'beta':0.5, 'tauz':1 }





Li=[o for o in range(1200000)]
dt=0.01
Lt=[dt*i for i in Li]
y=0.1
x=0.0 
z=0.0
f=0.0
Ly=[]
Lx=[]
Lz=[]
Lf=[]
seed(10)


Smax = Par['Smax']
Rs = Par['Rs']
lambdas = Par['lambdas']
taux = Par['taux']
P = Par['P']
Rb = Par['Rb']
lambdab = Par['lambdab']
L = Par['L']
tauy = Par['tauy']
S = Par['S']
alpha = Par['alpha']
beta = Par['beta']
tauz = Par['tauz']


for i in Li:

    #b=random()
    a=(random()-0.5)
   # La.append(3.4+a)
    #print(a)
    dx=(Smax/(1+np.exp((Rs-y)/lambdas))-x)/taux
    dy=(P/(1+np.exp((Rb-y)/lambdab))+L*f-y*x-z)/tauy
    dz=(S*(alpha*x+beta*y)*(a)-z)/tauz
    df=(y-1.0*f)/720
    y=y+dy*dt#+a
    x=x+dx*dt
    z=z+dz*dt
    f=f+df*dt

   # changes in parameters during simulation
   # pharmacological intervention
    if i>547500:
       P = 7.



    Ly.append(y)
    Lx.append(x)
    Lz.append(z)
    Lf.append(f)


cutT=0
")

#setup data in R
dat <- data.frame(id=1,time=a$Lt,x=unlist(a$Lx),y=unlist(a$Ly),z=unlist(a$Lz),f=unlist(a$Lf))
dat <- dat[1:547500,]

#model
library(ctsem)

# dx=(Smax/(1+np.exp((Rs-y)/lambdas))-x)/taux
# dy=(P/(1+np.exp((Rb-y)/lambdab))+L*f-y*x-z)/tauy
# dz=(S*(alpha*x+beta*y)*(a)-z)/tauz
# df=(y-1.0*f)/720

m <- ctModel(LAMBDA = cbind(diag(4),0),
  latentNames = c('x','y','z','f','a'),
  manifestNames = c('x','y','z','f'),
  type='stanct',
  DRIFT=c(
    0,0,0,0,0,
    0,0,0,0,0,
    0,0,0,0,0,
    0,0,0,0,0,
    0,0,0,0,-10), 
  CINT=c('((Smax/(1+exp((Rs-y)/lambdas))-x) / taux)',
    '((P/(1+exp((Rb-y)/lambdab))+L*f-y*x-z) / tauy)',
    '((S*(alpha*x+beta*y)*(a)-z)/tauz)',
    '((y-1.0*f)/720)',
    0),
  DIFFUSION=c(
    0,0,0,0,0,
    0,0,0,0,0,
    0,0,0,0,0,
    0,0,0,0,0,
    0,0,0,0,1),
  T0MEANS=c('t0x','t0y','t0z','t0f',0),
  MANIFESTMEANS = 0,
  # PARS=c('Smax|param+10','Rs|param+1','lambdas|param+.1','taux|param+14','tauy|param+14',
  #   'S|param+10','P|param+10','Rb|param+1.04','lambdab|param+.05','L|param+1.01','alpha|param+.5','beta|param+.5','tauz|param+1')
  PARS=c('Smax','Rs','lambdas','taux','tauy',
    'S','P','Rb','lambdab','L','alpha','beta','tauz')
)

ctModelLatex(m)


m$pars$transform[m$pars$matrix %in% 'PARS'] <- 'log1p_exp(param)'

subsamplestep <- 10
fit <- ctStanFit(datalong = dat[seq(1,floor(nrow(dat)/8),subsamplestep),],ctstanmodel = m,
  cores=6,plot=10,nlcontrol=list(maxtimestep=.1),verbose=1,optimcontrol=list(stochastic=F))
summary(fit)  

#show predictions leaving every nth observation
ctKalman(fit,plot=T,removeObs = 1000,kalmanvec=c('yprior'))+coord_cartesian(ylim=c(0,1))