國
立 政 治 大 學
‧
N a tio na
l C h engchi U ni ve rs it y
Chapter 4
Numerical Experiments and Analysis
In this chapter we consider Figure 2.1 again. Since we need to analyze the natural history of chronic hepatitis B virus infection only, the state “curative therapy” is combined with the state “death/transplantation” in Figure 2.1. In addition, after discussing with Dr. Chen in Chang Gung Medical Foundation Keelung Branch, we think that it is reasonable and acceptable to combine those two states “remission” and
“inactive carrier” since patients take a long time to be in state ”remission” after seroconversion. Moreover, patients with remission are hardly differentiated from normal people, that is, inactive carriers. Therefore, we have the new transition diagram shown as Figure 4.1.
As in Figure 2.1, we do not know exactly about the probabilities from “HBeAg(+) hepatitis HBV-DNA>2 10 6~7IU/ml” to “HBeAg(+) hepatitis HBV-DNA>2 10 4~5 IU/ml” and “HBeAg seroconversion”. The current imformation is that the probabilities are bounded by 2% and 15%. Hence, we know that the real life expectancy will also be bounded by the results of defining the probabilities as 2% and 15%. Therefore, we will assume that the sum of probabilities from “HBeAg(+) hepatitis HBV-DNA>2 10 6~7IU/ml” to “HBeAg(+) hepatitis HBV-DNA>2 10 4~5 IU/ml” and “HBeAg(+) hepatitis HBV-DNA> 2 10 6~7 IU/ml” to “HBeAg seroconversion” will be 2% and 15% respectively.
‧ 國
立 政 治 大 學
‧
N a tio na
l C h engchi U ni ve rs it y
In addition, two other parts in Figure 2.1 needs to be specified. One is 90%~95%
without indicating ”/year” from “HBeAg seroconversion” to “Remission”; the other is 5% without indicating “/year” from “HBeAg seroconversion” to “HBeAg(-) hepatitis HBV-DNA> 2 10 3~ 4 IU/ml”. These numbers seem denoting two proportions associated with “HBeAg seroconversion” instead of annual transition rates with respect to “HBeAg(-) hepatitis HBV-DNA>2 10 3~ 4IU/ml” and “Remission”. In order to make use of the method with transition probabilities, we take two scenarios for computation:
(1) Assume that 5% and 90~95% are annual rates, like other numbers given in Figure 2.1. Thus it is written and drawn in Figure 4.1.
(2) Assume that 5% and 90~95% are two proportions associated with “HBeAg seroconversion”. Accordingly, we draw similar progression flow for description in Figure 4.10.
Figure 4.1: The new transition diagram after rearrangement of states.
‧ 國
立 政 治 大 學
‧
N a tio na
l C h engchi U ni ve rs it y
One approach taken in Figure 4.1 is to consider the ratio of the probabilities from
“HBeAg(+) hepatitis HBV-DNA> 2 10 6~7 IU/ml” to “HBeAg(+) hepatitis HBV-DNA>2 10 4~5IU/ml” and from “HBeAg(+) hepatitis HBV-DNA>2 10 6~7 IU/ml” to “HBeAg seroconversion”. According to Chu et al. [9], we assume that the ratio approximately is 2:1. Therefore, suppose that by adopting 2% and 15% for investigation there will be about 1.3% or 10% to “HBeAg(+) hepatitis HBV-DNA>
2 10 4~5IU/ml” or 0.7% and 5% to “HBeAg seroconversion” respectively.
Two examples which illustrate the use of the revised matrix are given below. In short, we use Case 1, Case 2, Case 3, and Case 4 to represent 4 different cases. Case 1 and Case 2 are to be discussed in Example 4.1, and Case 3 and Case 4 are to be discussed in Example 4.2, correspondingly referred to Figure 4.1 and Figure 4.10:
Case 1:
In the transition model of HBV infection in Figure 4.1, assume the probability from
“HBeAg(+) hepatitis HBV-DNA> 2 10 6~7 IU/ml” to “HBeAg(+) hepatitis HBV-DNA>2 10 4~5IU/ml” and “HBeAg seroconversion” is 2%.
Case 2:
In the transition model of HBV infection in Figure 4.1, assume the probability from
“HBeAg(+) hepatitis HBV-DNA> 2 10 6~7 IU/ml” to “HBeAg(+) hepatitis HBV-DNA>2 10 4~5IU/ml” and “HBeAg seroconversion” is 15%.
Case 3:
In the transition model of HBV infection in Figure 4.10, assume the probability from
“HBeAg(+) hepatitis HBV-DNA> 2 10 6~7 IU/ml” to “HBeAg(+) hepatitis HBV-DNA>2 10 4~5IU/ml” and “HBeAg seroconversion” is 2%.
‧
In the transition model of HBV infection in Figure 4.10, assume the probability from
“HBeAg(+) hepatitis HBV-DNA> 2 10 6~7 IU/ml” to “HBeAg(+) hepatitis HBV-DNA>2 10 4~5IU/ml” and “HBeAg seroconversion” is 15%.
Example 4.1 Consider scenario (1) with the transition diagram shown in Figure 4.1.
First, we compute the probability that a person with initial age s at state i reaches the state j for the first time after n years without applying the life table. Our computation method is applied as mentioned in Chapter 3. In this example, we assume that there is a person who had the chronic hepatitis B virus infection when he is 25 years old. The original Markov transition matrix is written as follows.
0.98 0.013 0.007 0 0 0 0 0 0 0
Hence, P is a homogeneous matrix since it is time independent.
Then we compute the total probabilities of reaching death and “HBsAg loss” by Theorem 2.3 and equation (2.2.4). More precisely, we can also estimate the probabilities for a person had the chronic hepatitis B virus infection at age 25 at the initially healthy state hits “death” for the first time after n years by equation (2.3.5) and the average transition time from the initially healthy state to death. We show the simulation results as the following Figures 4.2, 4.3, 4.4, 4.5, and Table 4.1.
‧ 國
立 政 治 大 學
‧
N a tio na
l C h engchi U ni ve rs it y
Computing fi m( ),n of Case 1 and Case 2.
Let fi m( ),n(25), n1, 2,, be the probability of death within one year assuming a patient was aged 25 starting HBV infection. For example, according to equation (3.3.3), we know that fi m,( )n (25) is the probability for an individual who had the HBV infection at age 25 dies when he is exactly at age 25n. An individual will die eventually, so the sum of the probability of death at different ages is 100%. At first, one individual will not die easily, so the probability of death is small and then increasing year by year because of diseases or accidents. Consider Case 1. The highest mortality is at about age 70, and then the probability decreases since the sum of mortality is 100%. Note from Figure 4.2 that the corresponding probability is about 0.7% when he is aged 75, under the conditions that P1,21.3% and P1,30.7%.
Figure 4.2: The mortality of a person who had the HBV infection at age 25 without applying the population mortality rate (Case 1).
0 50 100 150 200 250 300 350
0 1 2 3 4 5 6 7x 10-3
age
probability of death within one year
‧ 國
立 政 治 大 學
‧
N a tio na
l C h engchi U ni ve rs it y
Figure 4.3: The mortality of a person who had the HBV infection at age 25 without applying the population mortality rate (Case 2).
Consider Case 2, for example, we know that fi m( ),n (25) is the probability for an individual who had the HBV infection at age 25 dies when he is exactly at age 25n. Note from Figure 4.3 that the corresponding probability for an individual who had the HBV infection at age 25 dies when he is exactly at the age 52 is about 1.6%, under the conditions that P1,210% and P1,3 5%.
0 50 100 150 200 250
0 0.002 0.004 0.006 0.008 0.01 0.012 0.014 0.016
age
probability of death within one year
‧
Computing the survival curves for Case 1 and Case 2.
Naturally assume at initial age for the observed patient the probability of survival is 100%. Then the survival probability decreases year by year eventually approaching to 0%. Note that from equation (3.1.13), the survival function is
( ) The corresponding survival functions FX( )x are shown graphically in Figure 4.4 and Figure 4.5, for Case 1 and Case 2 respectively.
Figure 4.4: The survival probability of a person who had the HBV infection at age 25 (Case 1).
For example, Figure 4.4 shows that the probability for an individual who had the HBV infection at age 25 to survive up to more than 100 years old is about 18%, under the conditions that P1,21.3% and P1,30.7%.
‧ 國
立 政 治 大 學
‧
N a tio na
l C h engchi U ni ve rs it y
Figure 4.5: The survival probability of a person who had the HBV infection at age 25 (Case 2).
Similarly, note from Figure 4.5 that the probability to survive up to more than 100 years old is about 10%, under the conditions that P1,210% and P1,35%.
Table 4.1 summarizes the mean remaining life expectancy, the standard deviation, probability of death after 100 years old, and the maximal life year.
Table 4.1: Summarized results of Case 1 and Case 2 without the life table.
Case 1 (2%) Case 2 (15%)
Mean remaining life expectancy 86.5/year 42.2/year
The standard deviation 59.3/year 32.3/year
Probability of death after 100 years old 34.3% 7.9%
Maximal life year 371 223
0 50 100 150 200 250 300 350 400 450 500
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
age
survival probablilty
‧ 國
立 政 治 大 學
‧
N a tio na
l C h engchi U ni ve rs it y
Since the real probability for the transition from “HBeAg(+) hepatitis HBV-DNA>2 10 6~7IU/ml” to “HBeAg(+) hepatitis HBV-DNA>2 10 4~5IU/ml” and
“HBeAg(+) hepatitis HBV-DNA>2 10 6~7IU/ml” to “HBeAg seroconversion” are bounded by 2% and 15%. Then, clearly, the real average life expectancy, the standard deviation, the probability of death after 100 years old, and the maximal life year are all bounded by the results of Case 1 and Case 2.
Consider Figure 4.2, Figure 4.3, Figure 4.4, Figure 4.5, and Table 4.1 presented above, it is not difficult to find that the mortality is unreasonable. The probability of death after 100 years old is 7.9%-34.3% in Case 1 and Case 2 respectively. We found that it is illegitimate for a man or woman staying alive up to almost for 223-371 years old. Obviously, it violates the common human expectancy we have known nowadays.
According to the data from ministry of interior of Taiwan [30], people in Taiwan usually die before 100 years old. There are only few people in Taiwan who have been alive more than one hundred years. Moreover, for a person who had the HBV infection the mean life expectancy is about 42-86 years, which is ever longer than the one who does not get the HBV infection.
‧
Now, we conduct the computation by the model proposed in this thesis in Chapter 3.
For instance, suppose that there is a person who had the chronic hepatitis B virus infection when he is 25 years old. We calculate the probability that a person with an initial age t at the state i reaches the state j for the first time after n years.
Since the initial age is 25, we know that px p250.9827 in Table 3.1. For example, we have the one-step transition matrix P25,26 as follows. By Algorithm 3.1, we know that the transition probability matrix from disease symptoms to death, 1M9 1 is a matrix whose elements are all 1,0M1 9 is a zero matrix, and IM1 1 is the identity matrix.
Hence, we have
25,26
0.9631 0.0128 0.0069 0 0 0 0 0 0 0.0173
0 0.9307 0 0 0 0.0393 0 0.0049 0.0079 0.0173
0.0491 0 0 0.0491 0.8845 0 0 0 0 0.0173
‧
Computing the survival curves for Case 1 and Case 2 with the life table.
By equation (3.3.1), (3.3.4) and (4.1.1), we may conduct the computation. The outcomes are displayed in the following Figure 4.6, 4.7 and Table 4.2. Note that equation (4.1.1) is time independent, and the following equation
( ) is time dependent. By Theorem 3.3 and equation (4.1.2), we have the survival function shown in Figure 4.6 and Figure 4.7 by applying the life table, respectively.
Figure 4.6: The survival probability of different ages starting at age 25 by applying the life table (Case 1).
‧
For example, after combining the life table with the Markov model, the probability for an individual who had the HBV infection at age 25 to survive up to 50 years old is about 40% in natural course of HBV infection, under the conditions that
1,2 1.3%
P and P1,30.7%.
Figure 4.7: The survival probability of different ages starting at age 25 by applying the life table (Case 2).
Similarly observe that the probability for an individual who had the HBV infection at age 25 to survive up to 50 years old is about 30% in natural course of HBV infection, under the conditions that P1,210% and P1,35%.
Table 4.2 summarizes the mean remaining life expectancy, standard deviation, probability of death after 100 years old, and the maximal life year by applying the life table in Case 1 and Case 2.
‧ 國
立 政 治 大 學
‧
N a tio na
l C h engchi U ni ve rs it y
Table 4.2 : Summarized results of Case 1 and Case 2 with the life table.
Case 1 (2%) Case 2 (15%)
Mean remaining life expectancy 22.2/year 18.5/year
The standard deviation 11.6/year 11.3/year
Probability of death after 100 years old 0% 0%
Maximal life year 86 80
Computing the confidence interval of life expectancy.
Consider Figure 4.6, Figure 4.7, and Table 4.2 shown above, the mean remaining life expectancy will be 18.5-22.2 years, and the standard deviation will be 11.3-11.6 years. The probability of death after 100 years is approximately 0%. These simulation data indicate that conducting mathematical computation by applying the life table is much more reasonable than before. Therefore, our proposed model is more accurate and authentic.
Next we construct the 95% confidence interval for the mean remaining life expectancy. The 95% confidence interval for the mean remaining life expectancy is shown in Table 4.3.
Table 4.3: The 95% confidence interval for life expectancy in Case 1 and Case 2.
Case 1 (2%) Case 2 (15%) Mean remaining life
expectancy
22.2/year 18.5/year
The standard deviation 11.6/year 11.3/year
Confidence interval for
life expectancy [ 19.29, 25.11 ] [ 15.51, 21.49 ]
‧
Regression models of survival curves for Case 1 and Case 2.
Now we provide the regression curve of survival functions for Case1 and Case 2 as shown in Figure 4.8 and Figure 4.9, respectively.
Figure 4.8 shows the regression model for the survival curve. Note that x is the age of the patient, the x axis represents the natural log value of age x , which is denoted by ln x, and the y axis is the survival probability corresponding to age x.
Figure 4.8: A survival curve by applying the life table (Case 1).
4.4
(24.5) (29.9) (36.6) (44.7) (54.6) (66.7) (81.5)
We have the cubic regression model y -73.57 60.89( ln ) -16.32( ln ) +1.426( ln ) x x 2 x3
in Case 1. R-square is the determination coefficient, where R-square=99.9% and adjusted R-square=99.9%. Hence, the regression curve can explain 99.9% variation.
Only 0.1% variation cannot be explained by this regression model.
‧
Figure 4.9: A survival curve by applying the life table (Case 2).
4.4
(24.5) (29.9) (36.6) (44.7) (54.6) (66.7) (81.5)
We give the cubic regression model y -63.67 53.83( ln ) -14.67( ln ) x x 21.302( ln )x 3
for Case 2. Since R-square=100% and adjusted R-square=100%, the regression curve can explain 100% variation. Only 0% variation cannot be explained by this regression model.
‧ 國
立 政 治 大 學
‧
N a tio na
l C h engchi U ni ve rs it y
ANOVA tables for Case 1 and Case 2.
The ANOVA tables for Case 1 and Case 2 are constructed with Minitab [21].
Table 4.4 and 4.5 give the analysis of variance and the sequential analysis of variance.
Table 4.4 : The ANOVA Table for Case 1.
Analysis of Variance
Source DF SS MS F P
Regression 3 6.37245 2.12415 15386.26 0.000
Error 53 0.00732 0.00014
Total 56 6.37977
Note that DF means the degree of freedom, SS is the sum of squares, MS is the mean sum of squares. SSE is the sum of squares error. Since p-value=0, we know that this regression model is considerably significant. Therefore, we can use this regression model to do further prediction.
Table 4.5 : The sequential analysis of variance for Case 1.
Sequential Analysis of Variance
Source DF SS F P
Linear 1 6.25198 2690.79 0.000
Quadratic 1 0.01877 9.30 0.004
Cubic 1 0.10171 736.71 0.000
We know that this regression model in Case 1 is significant, and the p-values of the linear term, quadratic term, and cubic term are 0, 0.004, and 0 respectively. Hence,
‧ 國
立 政 治 大 學
‧
N a tio na
l C h engchi U ni ve rs it y
we know that the linear term, quadratic term, and cubic term are significant as well.
Thus, we accept the regression model y -73.57 60.89( ln ) -16.32 ( ln ) x x 21.426( ln )x3
as our prediction for the survival curve function.
Table 4.6 : The ANOVA table for Case 2.
Analysis of Variance
Source DF SS MS F P
Regression 3 5.95127 1.98376 48373.95 0.000
Error 52 0.00213 0.00004
Total 55 5.95340
Since p-value=0, we know that this regression model is considerably significant.
Therefore, we can use this regression model to do further prediction.
Table 4.7 : The sequential analysis of variance for Case 2.
Sequential Analysis of Variance
Source DF SS F P
Linear 1 5.80386 2095.75 0.000
Quadratic 1 0.06911 45.54 0.000
Cubic 1 0.07830 1909.30 0.000
We know that this regression model in Case 2 is significant, and the p-values of the linear term, quadratic term, and cubic term are 0 respectively. Hence, we can tell that the linear term, quadratic term, and cubic term are significant as well. Hence, we accept the regression model y -63.67 53.83( ln ) -14.67( ln ) x x21.302( ln )x3 as our
‧ 國
立 政 治 大 學
‧
N a tio na
l C h engchi U ni ve rs it y
prediction for the survival curve function. According the simulation results of the proposed model, the mean remaining life expectancy is estimated about 18.5-22.2 years.
Example 4.2 Considering Case 3 and Case 4 of scenario (2) with transition rates and respectively as shown in Figure 4.10, we take the same approach described in Example 4.1. In brief, we consider Figure 4.10 and classify it into 11 states. They are summarized in Table 4.8.
Figure 4.10: An alternative interpretation for the transition diagram of Figure 4.1.
Let Case 3 be the case it takes “2%” into consideration and Case 4 for “15%”, since the real probability for the transition from “HBeAg(+) hepatitis HBV-DNA>
‧
hepatitis HBV-DNA>2 10 6~7IU/ml” to “HBeAg seroconversion” are bounded by between 2% and 15%. It runs the same analytic procedure just as what we have done in Example 4.1.Table 4.8: The symbols of states for Case 3 and 4.
Symptoms State symbol
HBeAg(+) hepatitis HBD-DNA> 2*10^6~7 IU/mL s1'
HBeAg(+) hepatitis HBD-DNA> 2*10^4~5 IU/mL s2'
HBeAg seroconversion
HBeAg(-) hepatitis HBD-DNA> 2*10^3~4 IU/mL 3
' s
HBeAg seroconversion
HBeAg(-) hepatitis HBD-DNA> 2*10^3~4 IU/mL Remission
HBeAg seroconversion s6'
Decompensation s7'
Liver cirrhosis s8'
HBsAg loss s9'
HCC s10'
End of HBV infection s11'
Let represent the annual probability from “HBeAg(+) hepatitis HBV-DNA” to
“HBeAg(+) hepatitis HBV-DNA> 2 10 4~5 IU/ml” and represent the annual probability from “HBeAg(+) hepatitis HBV-DNA” to “HBeAg seroconversion”.
Furthermore, set 1.3 and 0.7 in Case 3; 10 and 5 in Case 4. Since we assume the number “90%” is a proportion, we classify those 10 states in Figure 4.1 into 11 states as shown in Figure 4.10 and Table 4.8. The reason is that we rearrange the original states s3, s4, and s5 in Figure 4.1 into three states s3', s4,
s5, and s6 as shown in Figure 4.11. Thus, for example, the probability from
‧
“HBeAg(+) hepatitis HBV-DNA>2 10 6~7IU/ml” to “Seroconversion Remission” is
0.9 % as shown in Figure 4.12. Hence, we have new transition diagram shown in Figure 4.10.
Figure 4.11: The original transition diagram.
Figure 4.12: The transition diagram after separation.
HBeAg(+) hepatitis
‧ 國
立 政 治 大 學
‧
N a tio na
l C h engchi U ni ve rs it y
Computing fi m( ),n of Case 3 and Case 4.
Here we start our simulation. We also assume that the individual who got the HBV infection is at age 25. The following are the results of simulation.
Figure 4.13: The mortality of a person who had the HBV infection at age 25 (Case 3).
0 50 100 150 200 250 300 350 400
0 1 2 3 4 5 6 7x 10-3
age
mortality
‧ 國
立 政 治 大 學
‧
N a tio na
l C h engchi U ni ve rs it y
Figure 4.14: The mortality of a person who had the HBV infection at age 25 (Case 4).
0 50 100 150 200 250
0 0.002 0.004 0.006 0.008 0.01 0.012 0.014 0.016
age
mortality
‧ 國
立 政 治 大 學
‧
N a tio na
l C h engchi U ni ve rs it y
Computing the survival curves for Case 3 and Case 4.
The survival curves without the life table are shown in Figure 4.15 and 4.16.
Figure 4.15: The survival probability of a person who had the HBV infection at age 25 (Case 3).
0 50 100 150 200 250 300 350 400
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
age
survival probability
‧ 國
立 政 治 大 學
‧
N a tio na
l C h engchi U ni ve rs it y
Figure 4.16: The survival probability of a person who had the HBV infection at age 25 (Case 4).
Table 4.8 summarizes the mean remaining life expectancy, the standard deviation, the probability of death after 100 years old, and the maximal life year by applying the life table in Case 3 and Case 4.
Table 4.9: Summarized results of Case 3 and Case 4 without the life table.
Case 3 (2%) Case 4 (15%)
Mean remaining life expectancy 87.5/year 43.1/year
The standard deviation 58.6/year 30.8/year
Probability of death after 100 years old 20.8% 9.2%
Maximal life year 366 233
0 50 100 150 200 250
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
age
survival probability
‧ 國
立 政 治 大 學
‧
N a tio na
l C h engchi U ni ve rs it y
Similarly, since the real probabilities for the transitions are bounded by 2% and 15%. Clearly, the real average life expectancy, the standard deviation, the probability of death after 100 years old, and the maximal life year are all bounded by the results of Case 3 and Case 4 for the same reason. Obviously we find something unreasonable about Figure 4.13, Figure 4.14, Figure 4.15, Figure 4.16, and Table 4.8, too. The probability of death after 100 years old is 9.2%-20.8% in Case 3 and Case 4. We found that it is illegitimate for an individual staying alive up to almost 233-366 years old. Hence, this model has the same problem just like the cases in example 4.1.
Computing the survival curves for Case 3 and Case 4 with the life table.
Now, according to the model and method proposed in Chapter 3, the results follow.
Figure 4.17: The survival probability of different ages starting at age 25 by applying the life table (Case 3).
20 30 40 50 60 70 80 90
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
age
survival probability
‧
Figure 4.18: The survival probability of different ages starting at age 25 by applying the life table (Case 4).
Table 4.10: Summarized results of Case 3 and Case 4 with the life table.
Case 3 (2%) Case 4 (15%)
Mean remaining life expectancy 21.6/year 19.4/year
The standard deviation 11.8/year 11.1/year
Probability of dying after 100 years old 0% 0%
Maximal life year 81 80
Consider Figure 4.17, Figure 4.18, and Table 4.9, the mean remaining life expectancy is 19.4-21.6 years, and the standard deviation is 11.1-11.8 years. The probability of death after 100 years is approximately 0%. Those simulation data indicates that conducting mathematical computation by applying the life table is much more reasonable than before, too.
20 30 40 50 60 70 80
‧
Computing the confidence interval of life expectancy.
We also construct the 95% confidence interval for the mean remaining life expectancy as Table 4.11.
Table 4.11: The 95% confidence interval for life expectancy in Case 3 and Case 4.
Case 3 (2%) Case 4 (15%) Mean remaining life
expectancy 21.6/year 19.4/year
The standard deviation 11.8/year 11.1/year Confidence interval for
life expectancy [ 18.51 , 24.69 ] [ 16.47 , 22.33 ]
Regression models of survival curves for Case 3 and Case 4.
Now we have the regression models shown in Figure 4.19 and Figure 4.20.
Figure 4.19: A survival curve by applying the life table (Case 3).
Figure 4.19: A survival curve by applying the life table (Case 3).