Prediction of the Period of Psychotic Episode in Individual Schizophrenics by Simulation-Data Construction Approach

Abstract

Although schizophrenia can be treated, most patients still experience inevitable psychotic episodes from time to time. Precautious actions can be taken if the next onset can be predicted. However, sufficient information is always lacking in the clinical scenario. A possible solution is to use the virtual data generated from limited of original data. Data construction method (DCM) has been shown to generate the virtual felt earthquake data effectively and used in the prediction of further events. Here we investigated the performance of DCM in deriving the membership functions and discrete-event simulations (DES) in predicting the period embracing the initiation and termination time-points of the next psychotic episode of 35 individual schizophrenic patients. The results showed that 21 subjects had a success of simulations (RSS) ≥70%. Further analysis demonstrated that the co-morbidity of coronary heart diseases (CHD), risks of CHD, and the frequency of previous psychotic episodes increased the RSS. (150 words)

Appendix

Appendix 1. The algorithm of DCM

• Step 0: Given k>1 and set t = 1;

• Step 1: Input a sample $A_{o} = {\left\{ {\left. {{\left( {a_{{i_{o} }} ,\lambda _{i} } \right)}} \right|a_{{i_{o} }} \in \mathbb{R},\;\lambda _{i} \in \mathbb{N},\;a_{{{\left( {i + 1} \right)}_{o} }} >a_{{i_{o} }} ,\;1 \leqslant i,i_{o} \leqslant n} \right\}}$ and determine an arbitrary set C = {(1, 1), (c, 1)} with c > 1;

• Step 2: Take mode(A _o ), i.e., Eq. A.1, and let $A_{o} = {\left\{ {\left. {{\left( {a_{{i_{{}} }} ,\lambda _{i} } \right)}} \right|a_{{i_{{}} }} \in \mathbb{R},\;\lambda _{i} \in \mathbb{N},\; \leqslant i \leqslant n} \right\}}$ where a _i =$a_{{i_{{^{o} }} }}$ − mode(A _o );

  Note that

  $$\left\{ {\begin{array}{*{20}c}{{{\text{Initiation: mode}}{\left( {A_{o} } \right)} = \left\{ {\begin{array}{*{20}c}{{\max _{{\lambda _{i} }} \{ {\left( {a_{{i_{o} }} ,\lambda _{i} } \right)}\} ,{\text{ if there is an unique mode}}}} \\{{{\text{Pctl}}[\frac{{{\text{age of a patient}}}}{{67}}],{\text{ if there are multiple modes}}}} \\\end{array} } \right.}} \\{{{\text{Termination: mode}}{\left( {A_{o} } \right)} = \left\{ {\begin{array}{*{20}c}{{\max _{{\lambda _{i} }} \{ {\left( {a_{{i_{o} }} ,\lambda _{i} } \right)}\} ,{\text{ if there is an unique mode}}}} \\{{{\text{Pctl}}[1 - \frac{{{\text{age of a patient}}}}{{67}}],{\text{ if there are multiple modes}}}} \\\end{array} } \right.}} \\\end{array} } \right.$$

  (A.1)

  where $A_{o}$ is a sample formulating by the multiset form, $a_{{i_{o} }}$ is an element value, $\lambda _{i}$ is the frequency of $a_{{i_{o} }}$, Pctl[.] is the percentile of A _o . Presuming younger age a critical factor of remission (termination) from psychotic symptoms, the rules of mode (value that occurs the most frequently in a data set or a probability distribution) selection of periods with respect to the IRTs and TRTs of psychotic episodes are denoted by Eq. A.1. Statistics in 2006 showed that the average life-span among patients with schizophrenia in Taiwan was approximately 67 years while that of general population was 77 years [22]. If multiple modes of the response times of initiation and termination of psychotic episodes exist, then Eq. A.1 says that the older patients have longer response times for jumping from initiation to termination; on the other hand, their response times for jumping from termination to initiation are smaller.

• Step 3: $Z^{t} \equiv A{\left( : \right)}C^{t} = {\left\{ {\left. {{\left( {z_{{q^{{_{t} }} }} ,f_{{q^{{_{t} }} }} } \right)}} \right|z_{{q^{{_{t} }} }} \in {\left[ {a_{1} ,\;a_{n} } \right]},\;f_{{q^{{_{t} }} }} \in \mathbb{N},\;1 \leqslant q \leqslant Q} \right\}},\;Q \leqslant 2^{t} *{\left| A \right|},\;\forall t \in \mathbb{N};$;

• Step 4: Set a constant w with $w \geqslant 1 - {\left( {{\left( {{\left| {Z^{t} } \right|} + 1} \right)}{\left( {k^{2} + {\left| {Z^{t} } \right|} - 1} \right)}} \right)}^{{ - 1}} {\left( {{\left| {Z^{t} } \right|}k^{2} } \right)},$, $\overline{x} ^{t} = {\left( {{{\left( {{\sum\nolimits_q {z_{{q^{{_{t} }} }} \times f_{{q^{{_{t} }} }} } }} \right)}} \mathord{\left/ {\vphantom {{{\left( {{\sum\nolimits_q {z_{{q^{{_{t} }} }} \times f_{{q^{{_{t} }} }} } }} \right)}} {{\left| {Z^{t} } \right|}}}} \right. \kern-\nulldelimiterspace} {{\left| {Z^{t} } \right|}}} \right)},{\text{ and}}\;s^{t} = {\left( {{\left( {{\left| {Z^{t} } \right|}^{{ - 1}} {\sum\nolimits_q {z^{2}_{{q^{{_{t} }} }} \times f_{{q^{{_{t} }} }} } }} \right)} - {\left( {\overline{x} ^{t} } \right)}^{2} } \right)}^{{0.5}} .$.

  If $P{\left( {{\left| {z_{{q^{{_{t} }} }} - \overline{x} ^{t} } \right|} \leqslant ks^{t} } \right)} \geqslant w$ has been satisfied, go to step 5; Otherwise, t = t + 1 and go to step 3;

• Step 5: Assume t = T and calculate $z^{\prime }_{{q^{{_{T} }} }} = z_{{q^{T} }} + {\text{mode}}{\left( {A_{o} } \right)}$to output$Z^{{T\prime }} = {\left\{ {{\left( {z^{\prime }_{{q^{{_{T} }} }} ,f_{{q^{{_{T} }} }} } \right)}} \right\}}$.

• Step 6: Because the data distribution of Z _T is multimodality with M modals, let each modal be denoted by ${\left\{ {\left. {z_{m} {\left( {f_{m} } \right)}} \right|z_{m} \in Z_{T} ,f_{m} \in \mathbb{N},1 \leqslant m \leqslant M} \right\}}$.

• Step 7: Let $H = {\left\{ {{\left( {x_{1} ,y_{1} } \right)},{\left( {x_{2} ,y_{2} } \right)}, \ldots ,{\left( {x_{{M + 2}} ,y_{{M + 2}} } \right)}} \right\}} = {\left\{ {{\left( {{z_{{1_{T} }} ,f_{{1_{T} }} } \mathord{\left/ {\vphantom {{z_{{1_{T} }} ,f_{{1_{T} }} } {{\sum\nolimits_m {{\left( {f_{m} + f_{{1_{T} }} + f_{{Q_{T} }} } \right)}} }}}} \right. \kern-\nulldelimiterspace} {{\sum\nolimits_m {{\left( {f_{m} + f_{{1_{T} }} + f_{{Q_{T} }} } \right)}} }}} \right)},\;{\left\{ {{\left( {{Z_{m} ,\;f_{m} ,\;} \mathord{\left/ {\vphantom {{Z_{m} ,\;f_{m} ,\;} {{\sum\nolimits_m {{\left( {f_{m} + f_{{1_{T} }} + f_{{Q_{T} }} } \right)}} }}}} \right. \kern-\nulldelimiterspace} {{\sum\nolimits_m {{\left( {f_{m} + f_{{1_{T} }} + f_{{Q_{T} }} } \right)}} }}} \right)}\left| {1 \leqslant m \leqslant M} \right.} \right\}},{\left( {{z_{{Q_{T} }} ,f_{{Q_{T} }} } \mathord{\left/ {\vphantom {{z_{{Q_{T} }} ,f_{{Q_{T} }} } {{\sum\nolimits_m {{\left( {f_{m} + f_{{1_{T} }} + f_{{Q_{T} }} } \right)}} }}}} \right. \kern-\nulldelimiterspace} {{\sum\nolimits_m {{\left( {f_{m} + f_{{1_{T} }} + f_{{Q_{T} }} } \right)}} }}} \right)}} \right\}}.$

If the distribution of H is quasi-concave, fit {(x ₁, log(y ₁)), (x ₂, log(y ₂)), …, (x _{M  + 2}, log(y _{M  + 2}))} by using least square method; otherwise, fit {(x ₁, y ₁), (x ₂, y ₂), …, (x _{M  + 2}, y _{M  + 2})} directly by using Least square method (see Appendix 2). Then, an estimation density function can be obtained to represent A _o .

Appendix 2. Least square method to function fitting

Suppose to use the (M + 2) pair values of (x ₁, y ₁), (x ₂, y ₂), …, (x _{M  + 2}, y _{M  + 2}) to fit a three-power function

$$y = f{\left( x \right)} = b_{0} + b_{{1x}} + b_{2} x^{2} + b_{3} x^{3} $$

(A.2)

where b ₀, b ₁, b ₂, and b ₃ are unknown.

For minimizing $\Omega = {\sum\limits_{j = 1}^{M + 2} {{\left[ {y_{j} - {\left( {b_{0} + b_{1} x_{j} + b_{2} x^{2}_{j} + b_{3} x^{3}_{j} } \right)}} \right]}^{2} } }$, the following conditions should be met:

$$\frac{{\partial \Omega }}{{\partial b_{0} }} = - 2{\sum\limits_{j = 1}^{M + 2} {{\left[ {y_{j} - {\left( {b_{0} + b_{1} x_{j} + b_{2} x^{2}_{j} + b_{3} x^{3}_{j} } \right)}} \right]}} } = 0$$

(A.3)

$$\frac{{\partial \Omega }}{{\partial b_{1} }} = - 2{\sum\limits_{j = 1}^{M + 2} {x_{j} {\left[ {y_{j} - {\left( {b_{0} + b_{1} x_{j} + b_{2} x^{2}_{j} + b_{3} x^{3}_{j} } \right)}} \right]}} } = 0$$

(A.4)

$$\frac{{\partial \Omega }}{{\partial b_{2} }} = - 2{\sum\limits_{j = 1}^{M + 2} {x^{2}_{j} {\left[ {y_{j} - {\left( {b_{0} + b_{1} x_{j} + b_{2} x^{2}_{j} + b_{3} x^{3}_{j} } \right)}} \right]}} } = 0$$

(A.5)

$$\frac{{\partial \Omega }}{{\partial b_{3} }} = - 2{\sum\limits_{j = 1}^{M + 2} {x^{3}_{j} {\left[ {y_{j} - {\left( {b_{0} + b_{1} x_{j} + b_{2} x^{2}_{j} + b_{3} x^{3}_{j} } \right)}} \right]}} } = 0$$

(A.6)

Then, Eqs. A.3, A.4, A.5 and A.6 can be rearranged as

$${\sum\limits_{j = 1}^{M + 2} {y_{j} } } = (M + 2)b_{0} + b_{1} {\sum\limits_{j = 1}^{M + 2} {x_{j} } } + b_{2} {\sum\limits_{j = 1}^{M + 2} {x^{2}_{j} } } + b_{3} {\sum\limits_{j = 1}^{M + 2} {x^{3}_{j} } }$$

(A.7)

$${\sum\limits_{j = 1}^{M + 2} {x_{j} y_{j} } } = b_{0} {\sum\limits_{j = 1}^{M + 2} {x_{j} } } + b_{1} {\sum\limits_{j = 1}^{M + 2} {x^{2}_{j} } } + b_{2} {\sum\limits_{j = 1}^{M + 2} {x^{3}_{j} } } + b_{3} {\sum\limits_{j = 1}^{M + 2} {x^{4}_{j} } }$$

(A.8)

$${\sum\limits_{j = 1}^{M + 2} {x^{2}_{j} y_{j} } } = b_{0} {\sum\limits_{j = 1}^{M + 2} {x^{2}_{j} } } + b_{1} {\sum\limits_{j = 1}^{M + 2} {x^{3}_{j} } } + b_{2} {\sum\limits_{j = 1}^{M + 2} {x^{4}_{j} } } + b_{3} {\sum\limits_{j = 1}^{M + 2} {x^{5}_{j} } }$$

(A.9)

$${\sum\limits_{j = 1}^{M + 2} {x^{3}_{j} y_{j} } } = b_{0} {\sum\limits_{j = 1}^{M + 2} {x^{3}_{j} } } + b_{1} {\sum\limits_{j = 1}^{M + 2} {x^{4}_{j} } } + b_{2} {\sum\limits_{j = 1}^{M + 2} {x^{5}_{j} } } + b_{3} {\sum\limits_{j = 1}^{M + 2} {x^{6}_{j} } }$$

(A.10)

Let $Y = {\left[ {\begin{array}{*{20}c} {{y_{1} }} \\ { \vdots } \\ {{y_{{M + 2}} }} \end{array} } \right]},\;X = {\left[ {\begin{array}{*{20}c} {1} {{x_{1} }} {{x^{2}_{1} }} {{x^{3}_{1} }} \\ { \vdots } { \vdots } { \vdots } { \vdots } \\ {1} {{x_{{M + 2}} }} {{x^{2}_{{M + 2}} }} {{x^{3}_{{M + 2}} }} \end{array} } \right]}$, then

$$X^{T} X = {\left[ {\begin{array}{*{20}c} {{M + 2}} {{{\sum\nolimits_j {x_{j} } }}} {{{\sum\nolimits_j {x^{2}_{j} } }}} {{{\sum\nolimits_j {x^{3}_{j} } }}} \\ {{{\sum\nolimits_j {x_{j} } }}} {{{\sum\nolimits_j {x^{2}_{j} } }}} {{{\sum\nolimits_j {x^{3}_{j} } }}} {{{\sum\nolimits_j {x^{4}_{j} } }}} \\ {{{\sum\nolimits_j {x^{2}_{j} } }}} {{{\sum\nolimits_j {x^{3}_{j} } }}} {{{\sum\nolimits_j {x^{4}_{j} } }}} {{{\sum\nolimits_j {x^{5}_{j} } }}} \\ {{{\sum\nolimits_j {x^{3}_{j} } }}} {{{\sum\nolimits_j {x^{4}_{j} } }}} {{{\sum\nolimits_j {x^{5}_{j} } }}} {{{\sum\nolimits_j {x^{6}_{j} } }}} \end{array} } \right]},\;X^{T} Y = {\left[ {\begin{array}{*{20}c} {{{\sum\nolimits_j {y_{j} } }}} \\ {{{\sum\nolimits_j {x_{j} y_{j} } }}} \\ {{{\sum\nolimits_j {x^{2}_{j} y_{j} } }}} \\ {{{\sum\nolimits_j {x^{3}_{j} y_{j} } }}} \end{array} } \right]}$$

Therefore,

$$X^{T} Y = X^{T} X{\left[ {\begin{array}{*{20}c}{{b_{0} }} \\{{b_{1} }} \\{{b_{2} }} \\{{b_{3} }} \\\end{array} } \right]}$$

(A.11)

Finally, b ₀, b ₁, b ₂, and b ₃ can be solved by

$${\left[ {\begin{array}{*{20}c}{{b_{0} }} \\{{b_{1} }} \\{{b_{2} }} \\{{b_{3} }} \\\end{array} } \right]} = {\left( {X^{T} X} \right)}^{{ - 1}} {\left( {X^{T} Y} \right)}$$

(A.12)

Appendix 3. The computational criterion of the discrete-event simulation model

Let the initiation (I) and termination (T) of the psychotic episode be the events in the simulation model where the initiation of the psychotic episode is the initial event. Assume that there are 4 IRTs and 4 TRTs initially. According to historical time-points I = {i ₁, i ₂, …, i ₅} and T = {t ₁, t ₂, …, t ₄}, two functions, f _ΔI and f _ΔT , can be obtained to predict t ₅. Then, when using one trial of discrete-event simulation, a time sequence can be estimated as

$$E = {\left\{ {\widehat{i}_{1} ,\;\widehat{t}_{1} ,\;\widehat{i}_{2} ,\;\widehat{t}_{2} ,\;\widehat{i}_{3} ,\widehat{t}_{3} ,\;\widehat{i}_{4} ,\widehat{t}_{4} ,\;\widehat{i}_{5} ,\;\widehat{t}_{5} } \right\}}$$

where $\widehat{t}_{5}$ is the predictive time-point of the oncoming termination. Let e (days) be an acceptable bias, then if

$${{\left( {{\sum\limits_{r = 1}^5 {|\widehat{i}_{r} - i_{r} |} } + {\sum\limits_{r = 1}^4 {{\left| {\widehat{i}_{r} - t_{r} } \right|}} }} \right)}} \mathord{\left/ {\vphantom {{{\left( {{\sum\limits_{r = 1}^5 {|\widehat{i}_{r} - i_{r} |} } + {\sum\limits_{r = 1}^4 {{\left| {\widehat{i}_{r} - t_{r} } \right|}} }} \right)}} 9}} \right. \kern-\nulldelimiterspace} 9 \leqslant e$$

(A.13)

the prediction result will be recorded till J trails matched Eq. A.13 as

$$E_{j} = {\left\{ {\widehat{i}^{j}_{1} ,\;\widehat{t}^{j}_{1} ,\;\widehat{i}^{j}_{2} ,\;\widehat{t}^{j}_{2} ,\;\widehat{i}^{j}_{3} ,\;\widehat{t}^{j}_{3} ,\;\widehat{i}^{j}_{4} ,\;\widehat{t}^{j}_{4} ,\;\widehat{i}^{j}_{5} ,\;\widehat{t}^{j}_{5} } \right\}},\,j = 1,\,2,\, \ldots ,\,J.$$

(A.14)

Reversely, when setting termination (T) as the initial event, the estimated time sequence can be estimated as

$$E = {\left\{ {\widehat{t}_{1} ,\;\widehat{i}_{1} ,\;\widehat{t}_{2} ,\;\widehat{i}_{2} ,\;\widehat{t}_{3} ,\;\widehat{i}_{3} ,\;\widehat{t}_{4} ,\;\widehat{i}_{4} ,\;\widehat{t}_{5} ,\;\widehat{i}_{5} } \right\}}$$

where $\widehat{i}_{5}$ is the predictive time-point of the oncoming initiation. Then, Eqs. A.13 and A.14 can be rewritten into

$${{\left( {{\sum\limits_{r = 1}^5 {|\widehat{t}_{r} - t_{r} |} }} \right)} + {\sum\limits_{r = 1}^4 {{\left| {\widehat{i}_{r} - i_{r} } \right|}} }} \mathord{\left/ {\vphantom {{{\left( {{\sum\limits_{r = 1}^5 {|\widehat{t}_{r} - t_{r} |} }} \right)} + {\sum\limits_{r = 1}^4 {{\left| {\widehat{i}_{r} - i_{r} } \right|}} }} 9}} \right. \kern-\nulldelimiterspace} 9 \leqslant e$$

(A.15)

$$E_{j} = {\left\{ {\widehat{t}^{j}_{1} ,\;\widehat{i}^{j}_{1} ,\;\widehat{t}^{j}_{2} ,\;\widehat{i}^{j}_{2} ,\widehat{t}^{j}_{3} ,\;\widehat{i}^{j}_{3} ,\;\widehat{t}^{j}_{4} ,\;\widehat{i}^{j}_{4} ,\;\widehat{t}^{j}_{5} ,\;\widehat{i}^{j}_{5} } \right\}},\,j = 1,\,2,\, \ldots ,\,J.$$

(A.16)

Depending on Eqs. A.13 and A.15, the discrete-event simulation model can be validated by controlling the value of e. When a smaller value of e was selected, the prediction results would be more precise yet with higher simulation costs. This property can help decision makers to judge a suitable trade-off. Then, the time-point of the oncoming event, initiation or termination, can be denoted by $\widehat{i}_{5}$ or $\widehat{t}_{5}$. The corresponding interval estimators can be defined in the following.

Using Eq. A.13, ${\left\{ {\widehat{t}^{1}_{5} ,\;\widehat{t}^{2}_{5} ,\; \cdots ,\;\widehat{t}^{J}_{5} } \right\}}$ is obtained. Then, when setting a confidence level α, i.e., 90%, 95%, or 99%, a confidence interval can be computed by

$${\left( {\overline{{\widehat{t}}} _{5} - Z_{{\alpha \mathord{\left/ {\vphantom {\alpha 2}} \right. \kern-\nulldelimiterspace} 2}} \times \frac{{s_{{\widehat{t}_{5} }} }}{{{\sqrt J }}},\overline{{\widehat{{t_{5} }}}} + Z_{{\alpha \mathord{\left/ {\vphantom {\alpha 2}} \right. \kern-\nulldelimiterspace} 2}} \times \frac{{s_{{\widehat{t}_{5} }} }}{{{\sqrt J }}}} \right)}$$

(A.17)

where $\overline{{\widehat{t}}} _{5}$ and $s_{{\widehat{t}_{5} }}$ are the mean and standard deviation of $\widehat{t}_{5}$.

When setting termination (T) as the initial event, a confidence interval can be computed by

$${\left( {\overline{{\widehat{i}}} _{5} - Z_{{\alpha /2}} \times \frac{{s_{{\widehat{i}_{5} }} }}{{{\sqrt J }}},\overline{{\widehat{{i_{5} }}}} + Z_{{\alpha /2}} \times \frac{{s_{{\widehat{i}_{5} }} }}{{{\sqrt J }}}} \right)}$$

(A.18)

where $\overline{{\widehat{i}}} _{5}$ and $s_{{\widehat{i}_{5} }}$ are the mean and standard deviation of $\widehat{i}_{5}$.If $t_{5}$ or $i_{5}$ exactly falls into the confidence interval resulted from Eq. A.17 or A.18, this is a simulation with success.