program VALVES; {$M 65520,0,655360} (****************************************************************************) (* 27.05.93 *) (* LIFE EXPECTANCY AFTER PRIMARY AORTIC VALVE REPLACEMENT. *) (* This program calculates with the aid of a Markov model the expected *) (* number of life years of a patient (characterized by age and sex) since *) (* primary aortic valve replacement. *) (* *) (* Rob de Kruyk, CKB/EUR. *) (****************************************************************************) (****************************************************************************) (* 05.05.94 *) (* UNCERTAINTY ANALYSIS OF LIFE EXPECTANCY (ROBUSTNESS OF THE PROGRAM) *) (* The point estimates of the parameters of the Markov model are varied in *) (* order to quantify the influence of parameter variations on the life *) (* expectancies of four valve types. The point estimates are regarded as *) (* the medians of the individual probability distributions representing *) (* parameter uncertainty. To each parameter that is varied a 0.975th *) (* quantile is specified. Both the medians and quantiles are stored in the *) (* file "simul.dat". Five possible correlations between different valve *) (* types are defined in the source program: 0, 0.25, 0.50, 0.75 and 1. *) (* The simulation results are stored in different files "res1.out", ..., *) (* "res5.out" for each of these correlations. *) (* B.v.Vliet *) (****************************************************************************) {$N+} uses Crt, Dos; const NORM_95 = 1.96; (* 0.975th quantile of standard normal distr. *) minlft = 19; (* minimum age of patient *) maxvry = 40; (* maximum #years without valve replacement *) maxop = 10; (* maximum #operations per patient *) aantkleppen = 4; (* #valves *) aantvar = 20; (* #variables *) plus = 0.5; (* correlation used between variables for the *) NC = -plus; (* same valve type, see above *) simulatie = 2; (* sample size (>1) of uncertainty distribution *) type real = single; int = integer; ARR_KKV = array[1..aantkleppen+1,1..aantkleppen+1,1..aantvar] of real; rang = record eerste : int; tweede : int; derde : int; vierde : int; end; kans = record Pe, {incidence late endocarditis} ve, {odds ratio incidence early endocarditis} riscper, {length period early endocarditis in months} Pes, {mortality endocarditis} Pei, {morbidity endocarditis} Pk, {incidence valve thrombosis} Pks, {mortality valve thrombosis} Pki, {morbidity valve thrombosis} Pt, {incidence thromboembolism} Pts, {mortality thromboembolism} Pti, {morbidity thromboembolism} Pb, {incidence anticoagulantia-related haemorrhage} Pbs, {mortality anticoagulantia-related haemorrhage} Pbi, {morbidity anticoagulantia-related haemorrhage} Pf, {incidence non-structural dysfunction} Pok, {percentage operations at non-structural dysfunction} Pns, {Pok*Pf} beta, {shape parameter Weibull distribution} sigmac, {constant part scale parameter Weibull distribution} sigmal, {age dependent part scale parameter Weibull distr.} Phgs, {mortality thromboembolism or haemorrhage over some period} Ps, {probability of structural deterioration over some period} Pen, {probability of endocarditis over some period} eps {error in constant term of operation mortality, see fcn Oper} : real; Pos, {mortality after surgery} Poi {morbidity after surgery} : array[minlft..109] of real; end; pkans = ^kans; (* pointer to type record above *) arcbs = array[minlft..109] of real; (* array for age- and sex-specific mortality *) var h, vrschl_lft : real; rho,aa : ARR_KKV; out : array[1..5] of text; {output-files for different cor's} f : text; {input-file with incidences, mortalities and morbidities} age : INT; {age of patient at primary valve replacement} sex : char; rangklep : array[1..aantkleppen] of rang; i,z,j,k : integer; Po, {probability that a patient dies in some year} Pi, {probability to be 'alive with morbidity' in some year} Ph, {probability to be 'alive without morbidity' in a year} Pih, {probability to be 'alive with morbidity, valve } {replaced by homograft} Phh {probability to be 'alive without morbidity', valve} {replaced by homograft} : array[1..aantkleppen,minlft..109] of real; gemov, {life expectancy} gemgez {life expectancy without morbidity} : array[1..simulatie,1..aantkleppen] of real; median : array[1..6] of real; {medians life expectancies} dif_lifexp : array[1..simulatie,1..6] of real; {differences life expectancies between valve types} (****************************************************************************) Function OPER(age, opera: integer; x_i : real): real; (* returns probability to die after i-th re-operation for some age *) (* x_i is a deviation from mean value *) (* a belongs to (0.00628 - 0.005 , 0.00628 + 0.005) *) (* b --"-- (1.022 - 0.008 , 1.022 + 0.008) *) var a,b,p,x,pi,sx: real; begin {logistic equation for operation mortality after primary valve replacement} {for a Bjork-Shiley convexoconcave valve, see Van der Meulen [1993]} a := 0.015 / exp( 40 * ln(1.022) ) + x_i; b := 1.022 + x_i * 8 / 5; x:= a * exp( age * ln(b) ); p:= x/(1+x); {odds ratio = 1.7 for each re-operation} Oper:= p*exp(opera*ln(1.7))/(1-p+p*exp(opera*ln(1.7))); end; (* Oper *) (****************************************************************************) Function KLEPFALEN(age,jso, beta,sigmac,sigmal: real): real; (* returns probability of structural deterioration for age and #years since *) (* operation *) var sigma: real; begin {scale parameter consists of a constant and an age dependent part} sigma:= sigmac + ((age-jso+0.5)*sigmal); if sigma <= 0 then Klepfalen:= 0 else if jso = 0 then Klepfalen:= exp(beta*ln(1/sigma)) else Klepfalen:= 1-exp(exp(beta*ln(jso/sigma))-exp(beta*ln((jso+1)/sigma))); end; (* Klepfalen *) (****************************************************************************) Function ENDO(jso, Pe,ve,riscper: real): real; (* returns probability of endocarditis for #years since operation [=jso] *) var Pve: real; begin Pve:= ve * Pe; if (jso = 0) and (riscper > 0) then Endo:= 1 - exp(-Pe*(12-riscper)/12-Pve*riscper/12) else Endo:= 1 - exp(-Pe); end; (* Endo *) (****************************************************************************) Function logit(x: real): real; begin logit := ln(x/(1-x)); end; (* logit *) (****************************************************************************) Procedure INLEES(var CBS: arcbs; age: integer; sex: char); (* reads age- and sex-specific mortalities until age 109, with steps of *) (* one year *) var x : real; i : integer; f : text; begin if sex = 'f' then assign(f,'g:\data\chessa\mrfem90.dat') (* file for females *) else assign(f,'g:\data\chessa\mrmal90.dat'); (* file for males *) reset(f); read(f,x); readln(f,x); {birth} for i:= 0 to age-1 do begin read(f,x); readln(f,x); end; for i:= age to 109 do begin read(f,x); readln(f,x); CBS[i]:= x; end; close(f); end; (* Inlees *) (****************************************************************************) Procedure BEREKEN; (* main procedure of program *) Type rec = record mu : Real; bov : Real; sigma: Real; end; ARR_KV = array[1..aantkleppen,1..aantvar] of real; ARR_KK = array[1..aantkleppen,1..aantkleppen] of real; var tiep, {record with probabilities for patients with original valve type} allo : pkans; {similar, but now for patients with valve replaced by homograft} CBS : arcbs; {array with age- and sex-specific mort.} i,j,k,r,n, lft, {auxiliary variable for age} g : integer; {valve type: 1:homograft; 2:bioprosthesis; 3:mechanical; 4:autograft} ch : char; {the event structural deterioration is is split up into a part causing} {death (ss) and a part merely leading to morbidity or complete cure (si)} Pss, {probability of "ss" for patients with original valve type} Pssh, {probability for patients with original valve replaced by homograft} Psi, {probability of "si" for patients with original valve type} Psih, {probability for patients with original valve replaced by homograft} cum,cumi,p1,p2, perhh : real; {percentage of patients with homograft} som,som1 : real; b : 0..1; {b=1 for autograft and 0 for other types} a : ARR_KV ; {array for storing probabilities and incidences from input-file} Klep, {probability of structural deterioration for some age and #years since} {last operation for patients with original valve type} En, {same for endocarditis, independent of age} Kleph, {same as Klep, but original valve replaced by homograft} Enh, {same as En, but original valve replaced by homograft} Opovl : real; {probability of surviving operation at some age and #previous operations} qj, {qj[i,k] is the probability that in year j a patient has undergone k ops} {and that i years have passed since last operation, given that patient} {is still alive with original valve} qj_1, {same for year j-1} hqj, {same for valve replaced by homograft} hqj_1 {same for year j-1} : array[0..maxvry,1..maxop] of real; P : array[1..aantkleppen,1..aantvar] of rec; gem_verschil, sigma_verschil, verschil : ARR_KK; rangtussenkleppen, rang_plusdelta, rang_mindelta : array[1..aantkleppen,1..aantkleppen] of int; m,s,xn, PreviousNormalDeviate, C : real; {correlation between valves} ThereIsAPreviousNormal : BOOLEAN; AcceptSim : BOOLEAN; {variable for retaining or rejecting} {a realization of the uncert. distr.} InvN : array[1..16] of real; {values of inverse of N(0,1)} Help : array[1..16] of real; E_logit,E_log : array[1..aantkleppen,1..aantvar] of real; (****************************************************************************) (* Function RandNormal1 returns a random deviate from N(Mu,Sigma). *) (* Simulation method: the Box-Muller method. *) (* (source: Press et al,'Numerical recipes in fortran' p202-203). *) (* As this method always returns two standard-Normal Random deviates *) (* every time the function is called it will be checked whether a *) (* deviate is left from the previous call. If so, this value will be *) (* used, which will be assigned to the global variable *) (* PreviousNormalDeviate. *) (****************************************************************************) function RandNormal(Mu,Sigma:real):real; var v1,v2,R,Factor, StandardOrdinate : real; begin if not ThereIsAPreviousNormal then begin repeat v1:=2*Random - 1; v2:=2*Random -1; R:=v1*v1+v2*v2; until R<=1; Factor:=Sqrt(-2*Ln(R)/R); PreviousNormalDeviate:=v1*Factor; StandardOrdinate:=v2*Factor; ThereIsAPreviousNormal:=true; end else begin StandardOrdinate:=PreviousNormalDeviate; ThereIsAPreviousNormal := false; end; RandNormal:=StandardOrdinate*Sigma+Mu; end; PROCEDURE INLEZEN_CORR (VAR cor: real; VAR rho1 : ARR_KKV); VAR i1,g1,g2,N : INT; begin for i1 := 1 to aantvar do begin if NOT ( ( i1 = 15 ) OR ( i1 = 12 ) ) then N := aantkleppen else N := aantkleppen + 1; for g1 := 1 to N do begin rho1[g1,g1,i1] := 1; for g2 := 1 to g1-1 do begin if ( ( i1 = 2 ) OR ( i1 = 16 ) ) then begin if ( g1 = 3 ) and ( g2 = 2 ) then rho1[g1,g2,i1] := cor else rho1[g1,g2,i1] := 0; end else if ( ( i1 = 3 ) OR ( i1 = 5 ) OR ( i1 = 7 ) OR ( i1 = 8 ) OR ( i1 = 11 ) OR ( i1 = 14 ) OR ( i1 = 17 ) OR ( i1 = 19 ) ) then rho1[g1,g2,i1] := 0 else if ( i1 = 6 ) then begin if ( ( g1 = 3 ) OR ( g1 = 2 ) ) then begin if ( g1 = 2 ) then rho1[g1,g2,i1] := NC else begin if ( g2 = 1 ) then rho1[g1,g2,i1] := -NC else rho1[g1,g2,i1] := NC; end; end else rho1[g1,g2,i1] := 0; end else if ( i1 = 18 ) then begin if ( ( g1 = 3 ) OR ( g2 = 3 ) ) then rho1[g1,g2,i1] := 0 else rho1[g1,g2,i1] := cor; end else if ( i1 = 12 ) then begin if ( ( g1 <> 1 ) AND ( g2 = 1 ) ) then rho1[g1,g2,i1] := NC else rho1[g1,g2,i1] := cor; end else if ( i1 = 15 ) then begin if ( ( g1 <> 1 ) AND ( g2 = 1 ) ) then rho1[g1,g2,i1] := plus else rho1[g1,g2,i1] := cor; end else rho1[g1,g2,i1] := cor; end; end; end; end; PROCEDURE CHOLESKY_MATRIX (cor: real; VAR rho2,aa1 : ARR_KKV); VAR i1,g1,g2,h1,h2,k1,N,teller : INT; som,t1,t2 : REAL; begin INLEZEN_CORR(cor,rho2); for i1 := 1 to aantvar do begin if ( ( i1 = 1 ) OR ( i1 = 4 ) OR ( i1 = 9 ) OR ( i1 = 10 ) OR ( i1 = 12 ) OR ( i1 = 13 ) OR ( i1 = 15 ) OR ( i1 = 20 ) ) then begin if NOT ( ( i1 = 15 ) OR ( i1 = 12 ) ) then begin N := aantkleppen; aa1[1,1,i1] := P[1,i1].sigma; end else begin N := aantkleppen + 1; if ( i1 = 12 ) then aa1[1,1,i1] := P[g1,9].sigma else aa1[1,1,i1] := P[g1,1].sigma; end; for g1 := 2 to N do begin if NOT ( ( i1 =12 ) OR ( i1 = 15 ) ) then h1 := g1 else h1 := g1 - 1; aa1[g1,1,i1] := rho2[g1,1,i1] * P[h1,i1].sigma; if ( g1 > 2 ) then begin for g2 := 2 to (g1-1) do begin if NOT ( ( i1 = 12 ) OR ( i1 = 15 ) ) then h2 := g2 else h2 := g2 - 1; som := 0; teller := 0; for k1 := 1 to (g2-1) do som := som + aa1[g1,k1,i1] * aa1[g2,k1,i1]; aa1[g1,g2,i1] := ( rho2[g1,g2,i1] * P[h1,i1].sigma * P[h2,i1].sigma - som ) / aa1[g2,g2,i1]; end; end; if ( g1 > 1 ) then begin som := 0; for k1 := 1 to (g1-1) do som := som + SQR(aa1[g1,k1,i1]); aa1[g1,g1,i1] := SQRT( SQR(P[h1,i1].sigma) - som ); end; end; end; if ( i1 = 2 ) OR ( i1 = 16 ) then begin aa1[1,1,i1] := P[2,i1].sigma; aa1[2,1,i1] := rho2[3,2,i1] * P[3,i1].sigma; aa1[2,2,i1] := P[3,i1].sigma * SQRT( 1-SQR(rho2[3,2,i1]) ); for g1 := 3 to 4 do for g2 := 1 to g1 do aa1[g1,g2,i1] := 0; end; if ( i1 = 18 ) then begin aa1[1,1,i1] := P[1,i1].sigma; aa1[2,1,i1] := rho2[2,1,i1] * P[2,i1].sigma; aa1[2,2,i1] := P[2,i1].sigma * SQRT( 1-SQR(rho2[2,1,i1]) ); aa1[3,1,i1] := rho2[4,1,i1] * P[4,i1].sigma; aa1[3,2,i1] := P[4,i1].sigma * ( rho2[4,2,i1] - rho2[4,1,i1] * rho2[2,1,i1] ) / SQRT( 1 - SQR(rho2[2,1,i1]) ); aa1[3,3,i1] := SQRT( SQR( P[4,i1].sigma ) - SQR ( aa1[3,1,i1] ) - SQR ( aa1[3,2,i1] ) ); for g2 := 1 to 4 do aa1[4,g2,i1] := 0; end; if ( i1 = 6 ) then begin aa1[1,1,i1] := P[3,9].sigma; aa1[2,1,i1] := rho2[2,1,i1] * P[3,12].sigma; aa1[2,2,i1] := P[3,12].sigma * SQRT( 1-SQR(rho2[2,1,i1]) ); aa1[3,1,i1] := rho2[3,1,i1] * P[3,i1].sigma; aa1[3,2,i1] := P[3,i1].sigma * ( rho2[3,2,i1] - rho2[3,1,i1] * rho2[2,1,i1] ) / SQRT( 1 - SQR(rho2[2,1,i1]) ); aa1[3,3,i1] := SQRT( SQR( P[3,i1].sigma ) - SQR ( aa1[3,2,i1] ) - SQR ( aa1[3,1,i1] ) ); for g2 := 1 to 4 do aa1[4,g2,i1] := 0; end; if ( ( i1 = 3 ) OR ( i1 = 5 ) OR ( i1 = 7 ) OR ( i1 = 8 ) OR ( i1 = 11 ) OR ( i1 = 14 ) OR ( i1 = 17 ) OR ( i1 = 19 ) ) then begin for g1 := 1 to aantkleppen do for g2 := 1 to g1 do aa1[g1,g2,i1] := 0; end; end; end; PROCEDURE REALISATIES (aa2 : ARR_KKV); VAR j1,h1,h2,k1,k2 : INT; som1,t1,t2 : REAL; zzz : array[1..aantkleppen,1..aantvar] of REAL; begin {assign probabilities read to variables} for j1 := 1 to aantvar do begin for h1 := 1 to aantkleppen do zzz[h1,j1] := RandNormal(0,1); if ( ( j1 = 1 ) OR ( j1 = 4 ) OR ( j1 = 9 ) OR ( j1 = 10 ) OR ( j1 = 13 ) OR ( j1 = 20 ) ) then begin for h1 := 1 to aantkleppen do begin som1 := 0; for h2 := 1 to h1 do som1 := som1 + aa2[h1,h2,j1] * zzz[h2,j1]; a[h1,j1] := som1 + P[h1,j1].mu; end; end; if ( ( j1 = 12 ) OR ( j1 = 15 ) ) then begin for h1 := 1 to aantkleppen do begin som1 := 0; for h2 := 1 to h1+1 do begin if ( h2 > 1 ) then begin k2 := j1; k1 := h2 - 1; end else if ( j1 = 12 ) then begin k1 := h1; k2 := 9; end else begin k1 := h1; k2 := 1; end; som1 := som1 + aa2[h1+1,h2,j1] * zzz[k1,k2]; end; a[h1,j1] := som1 + P[h1,j1].mu; end; end; if ( ( j1 = 2 ) OR ( j1 = 16 ) ) then begin a[1,j1] := P[1,j1].mu; a[4,j1] := P[4,j1].mu; zzz[3,j1] := RandNormal(0,1); zzz[2,j1] := RandNormal(0,1); a[2,j1] := aa2[1,1,j1] * zzz[2,j1] + P[2,j1].mu; a[3,j1] := aa2[2,1,j1] * zzz[2,j1] + aa2[2,2,j1] * zzz[3,j1] + P[3,j1].mu; end; if ( j1 = 18 ) then begin a[3,j1] := P[3,j1].mu; zzz[1,j1] := RandNormal(0,1); zzz[2,j1] := RandNormal(0,1); zzz[4,j1] := RandNormal(0,1); a[1,j1] := aa2[1,1,j1] * zzz[1,j1] + P[1,j1].mu; a[2,j1] := aa2[2,1,j1] * zzz[1,j1] + aa2[2,2,j1] * zzz[2,j1] + P[2,j1].mu; a[4,j1] := aa2[3,1,j1] * zzz[1,j1] + aa2[3,2,j1] * zzz[2,j1] + aa2[3,3,j1] * zzz[3,j1] + P[4,j1].mu; end; if ( ( j1 = 6 ) OR ( j1 = 7 ) ) then begin zzz[3,j1] := RandNormal(0,1); for h1 := 1 to aantkleppen do if ( h1 = 3 ) then a[h1,j1] := zzz[3,j1] * P[h1,j1].sigma + P[h1,j1].mu else a[h1,j1] := 0; end; if ( ( j1 = 3 ) OR ( j1 = 5 ) OR ( j1 = 8 ) OR ( j1 = 11 ) OR ( j1 = 14 ) OR ( j1 = 17 ) OR ( j1 = 19 ) ) then begin for h1 := 1 to aantkleppen do a[h1,j1] := P[h1,j1].mu; end; end; {transform the realizations a[.,.] drawn} for h1:=1 to aantkleppen do for j1:=1 to aantvar-1 do if P[h1,j1].mu <> P[h1,j1].bov then if (j1<>2) and (j1<>18) then a[h1,j1] := exp(a[h1,j1])/(1 + exp(a[h1,j1])) else a[h1,j1] := exp(a[h1,j1]); end; begin (* bereken *) Help[1] := 0.002; Help[2] := 0.003; Help[3] := 0.004; Help[4] := 0.005; Help[5] := 0.007; Help[6] := 0.012; Help[7] := 0.016; Help[8] := 0.02; Help[9] := 0.025; Help[10]:= 0.07; Help[11]:= 0.1; Help[12]:= 0.13; Help[13]:= 0.14; Help[14]:= 0.25; Help[15]:= 0.5; Help[16]:= 0.66; InvN[1] := -2.88; InvN[2] := -2.75; InvN[3] := -2.65; InvN[4] := -2.575; InvN[5] := -2.455;InvN[6] := -2.26; InvN[7] := -2.145;InvN[8] := -2.05; InvN[9] := -1.96; InvN[10]:= -1.48; InvN[11]:= -1.28; InvN[12]:= -1.13; InvN[13]:= -1.08; InvN[14]:= -0.675;InvN[15]:= 0.0; InvN[16]:= 0.41; ThereIsAPreviousNormal := false; vrschl_lft := 1.225; write('Sex of patient(m/f): '); readln(sex); write('Age of patient(in years): '); readln(age); lft := age; assign(f,'g:\data\chessa\simul.dat'); reset(f); for i := minlft to 109 do CBS[i]:= 0; Inlees(CBS,age,sex); readln(f); for g := 1 to aantkleppen do begin for i := 1 to (aantvar-1) do read(f,P[g,i].mu,P[g,i].bov); readln(f); readln(f); P[g,20].mu := 0; (* added *) P[g,20].bov := 0.00001; P[g,20].sigma := P[g,20].bov/NORM_95; end; close(f); { for g:=1 to aantkleppen do for i:=1 to aantvar-1 do begin m := P[g,i].mu; s := P[g,i].bov; if m <> s then begin if (i<>2) and (i<>18) then begin if ((m > 0.499) and (m < 0.501)) then P[g,i].mu := logit(m) else begin P[g,i].mu := 2*logit(s)*(1-(2*m)) + sqr(NORM_95); P[g,i].mu := P[g,i].mu - sqrt(4*sqr(NORM_95)*(1-(2*m))*(logit(s)-logit(m)) + exp(4*ln(NORM_95))); P[g,i].mu := P[g,i].mu/(2*(1-(2*m))); end; P[g,i].sigma := (logit(s) - P[g,i].mu)/NORM_95; end else begin P[g,i].sigma := 0.5*(sqrt(sqr(NORM_95) + (4*ln(s/m))) - NORM_95); P[g,i].mu := sqr(P[g,i].sigma) + ln(m); end; end; end; for g:=1 to aantkleppen do begin for i:=1 to aantvar-1 do begin m := P[g,i].mu; s := P[g,i].bov; if m <> s then begin if (i<>2) and (i<>18) then begin if ((m > 0.49999) and (m < 0.50001)) then P[g,i].mu := 0 else if ((m > 0.024999) and (m < 0.025001)) or (m = 0.975) then P[g,i].mu := (sqr(logit(m)) + sqr(logit(s)))/(2*logit(s)) else begin r:=1; while Help[r] < m do r:=r+1; P[g,i].mu := (-NORM_95/InvN[r])* sqrt(sqr(logit(s)) - sqr(logit(m)) + sqr(logit(m)*NORM_95/InvN[r])); P[g,i].mu := (P[g,i].mu + logit(s))/(1 - sqr(NORM_95/InvN[r])); end; P[g,i].sigma := (logit(s) - P[g,i].mu)/NORM_95; end else begin P[g,i].sigma := NORM_95 - sqrt(sqr(NORM_95) - ln(s/m)); P[g,i].mu := ln(m) - (0.5*sqr(P[g,i].sigma)); end; end; write(P[g,i].mu:6:4,' '); end; writeln; end; } for g:=1 to aantkleppen do begin for i:=1 to aantvar-1 do begin m := P[g,i].mu; s := P[g,i].bov; if m <> s then begin if (i<>2) and (i<>18) then begin P[g,i].mu := logit(m); P[g,i].sigma := (logit(s) - P[g,i].mu)/NORM_95; { E_logit[g,i] := 0; for n:=1 to 9999 do begin xn := n/10000; E_logit[g,i] := E_logit[g,i] + exp(sqr(logit(xn) - P[g,i].mu)/(-2*sqr(P[g,i].sigma)))/(1-xn); end; E_logit[g,i] := E_logit[g,i]/(10000*sqrt(6.2831853)*P[g,i].sigma); writeln(E_logit[g,i]:6:4); } end else begin P[g,i].mu := ln(m); P[g,i].sigma := (ln(s) - P[g,i].mu)/NORM_95; { E_log[g,i] := exp(P[g,i].mu + 0.5*sqr(P[g,i].sigma)); writeln(E_log[g,i]:7:4); } end; end; end; { readln(ch); } end; for r:=1 to 5 do begin if r=1 then begin C := 0.000001; assign(out[r],'g:\data\chessa\res1.out'); end else if r=2 then begin C := 0.25; age := lft; assign(out[r],'g:\data\chessa\res2.out'); end else if r=3 then begin C := 0.5; age := lft; assign(out[r],'g:\data\chessa\res3.out'); end else if r=4 then begin C := 0.75; age := lft; assign(out[r],'g:\data\chessa\res4.out'); end else if r=5 then begin C := 0.999; age := lft; assign(out[r],'g:\data\chessa\res5.out'); end; rewrite(out[r]); {determine lower triangular matrix aa for correlated variates} CHOLESKY_MATRIX(C,rho,aa); writeln(out[r],'Size of simulation = ', simulatie:3); writeln(out[r],'Correlation among different valve types = ', C:5:3); writeln(out[r],'Positive correlation between var`s 1 and 5 for same valve type = ', plus:5:3); writeln(out[r],'Negative correlation between var`s 12, 9 and 6 for same valve type = ', NC:5:3); writeln(out[r]); repeat {for all ages} begin for g:=1 to aantkleppen do begin rangklep[g].eerste := 0; rangklep[g].tweede := 0; rangklep[g].derde := 0; rangklep[g].vierde := 0; for i := 1 to aantvar do begin rangtussenkleppen[g,i] := 0; rang_plusdelta[g,i] := 0; rang_mindelta[g,i] := 0; end; end; z:=1; while z <= simulatie do begin writeln(z:2,' age = ',age:2); AcceptSim := true; {draw realizations for all variables and all valve types} REALISATIES(aa); tiep:= new(pkans); allo:= new(pkans); for g:= 1 to aantkleppen do begin b:=0; if g = 4 then b:=1; with tiep^ do begin Pe:= a[g,1]; ve:= a[g,2]; riscper:= a[g,3]; Pes:= a[g,4]; Pei:= a[g,5]; Pk:= a[g,6]; Pks:= a[g,7]; Pki:= a[g,8]; Pt:= a[g,9]; Pts:= a[g,10]; Pti:= a[g,11]; Pb:= a[g,12]; Pbs:= a[g,13]; Pbi:= a[g,14]; Pf:=a[g,15]; Pok:= a[g,16]; beta:= a[g,17]; sigmac:= a[g,18]; sigmal:= a[g,19]; eps := a[g,20]; Po[g,age]:= Oper(age,0,eps)*(1+2.0*b)/(1+2.0*b*Oper(age,0,eps)); Pi[g,age]:= (1-Po[g,age])*Po[g,age]; Ph[g,age]:= 1 - Po[g,age] - Pi[g,age]; Phh[g,age]:= 0; Pih[g,age]:= 0; for j:= age+1 to 109 do begin Po[g,j]:= 0; Pi[g,j]:= 0; Ph[g,j]:= 0; Pih[g,j]:= 0; Phh[g,j]:= 0; end; Pk:= Pk/(Pk+Pt+Pb+Pf)*(1-exp(-Pk-Pt-Pb-Pf)); Pt:= Pt/(Pk+Pt+Pb+Pf)*(1-exp(-Pk-Pt-Pb-Pf)); Pb:= Pb/(Pk+Pt+Pb+Pf)*(1-exp(-Pk-Pt-Pb-Pf)); Pf:= Pf/(Pk+Pt+Pb+Pf)*(1-exp(-Pk-Pt-Pb-Pf)); Pns:= Pok*Pf; {kans op niet-structureel klepfalen die leidt tot operatie in een periode} {calculation probability to die from thromboembolism or haemorrhage} Phgs:= Pt*Pts + Pb*Pbs; for i:= 0 to maxvry do for k:= 1 to maxop do begin qj[i,k]:= 0; hqj[i,k]:=0; end; qj[0,1]:= 1; Pos[age]:= Oper(age,1,eps); Poi[age]:= (1-Pos[age])*Pos[age]; allo^.Pos[age]:= 0; allo^.Poi[age]:= 0; for j:= age + 1 to 109 do begin for i:= 0 to maxvry do for k:= 1 to maxop do begin qj_1[i,k]:= qj[i,k]; qj[i,k]:= 0; hqj_1[i,k]:= hqj[i,k]; hqj[i,k]:= 0; end; {calculation probabilities structural deterioration and endocarditis in} {corresponding year based on distribution of time since last operation at} {beginning of year} Ps:= 0; Pen:= 0; allo^.Ps:= 0; allo^.Pen:= 0; for i:= 0 to maxvry do begin Klep:= Klepfalen(j-1,i,beta,sigmac,sigmal); En:= Endo(i,Pe,ve,riscper); if Pe*ve > 1 then AcceptSim := false; if g > 1 then begin Kleph:= Klepfalen(j-1,i,allo^.beta,allo^.sigmac,allo^.sigmal); Enh:= Endo(i,allo^.Pe,allo^.ve,allo^.riscper); if allo^.Pe*allo^.ve > 1 then AcceptSim := false; end; for k:= 1 to maxop do begin if qj_1[i,k] > 0 then begin Ps:= Ps + qj_1[i,k]*Klep; Pen:= Pen + qj_1[i,k]*En; end; if g > 1 then if hqj_1[i,k] > 0 then begin allo^.Ps:= allo^.Ps + hqj_1[i,k]*Kleph; allo^.Pen:= allo^.Pen + hqj_1[i,k]*Enh; end; end; end; if g = 1 then allo^:= tiep^; Pss:= Ps*Pos[j-1]; Pssh:= allo^.Ps*allo^.Pos[j-1]; Psi:= Ps*(1-Pos[j-1]); Psih:= allo^.Ps*(1-allo^.Pos[j-1]); {calculation of probabilities of being in some state in a certain year based} {on probabilities at the end of the preceding year and transition probabilities} {for the present year} Po[g,j]:= Po[g,j-1] + (Ph[g,j-1]+Pi[g,j-1]) *(CBS[j-1] + (1-CBS[j-1])*(Pss+ (1-Pss)*(Pen*Pes+Pk*Pks+Pt*Pts+Pb*Pbs+Pns*Pos[j-1]))) + (Phh[g,j-1]+Pih[g,j-1]) *(CBS[j-1] + (1-CBS[j-1])*(Pssh+(1-Pssh)*(allo^.Pen*allo^.Pes+allo^.Pt*allo^.Pts+allo^.Pb*allo^.Pbs +allo^.Pns*allo^.Pos[j-1]))); Pi[g,j] := Pi[g,j-1]*(1-Pss)*(1-CBS[j-1]); Pi[g,j] := Pi[g,j]*((1-b)*(1-Pen-Pk-Pt-Pb-Pns + Pk*(1-Pks)+Pt*(1-Pts)+Pb*(1-Pbs)+Pns*(1-Pos[j-1])) + b*(1-Pen-Pk-Pt-Pb-Pns-Psi + Pk*(1-Pks)+Pt*(1-Pts)+Pb*(1-Pbs))) + Ph[g,j-1]*(1-Pss)*(1-CBS[j-1])*((1-b)*(Pk*Pki+Pt*Pti+Pb*Pbi+Pns*Poi[j-1]+Psi*(Poi[j-1]/(1-Pos[j-1]))) + b*(Pk*Pki+Pt*Pti+Pb*Pbi)); Pih[g,j]:= Phh[g,j-1]*(1-Pssh)*(1-CBS[j-1])*(allo^.Pen*allo^.Pei+allo^.Pt*allo^.Pti+allo^.Pb*allo^.Pbi +allo^.Pns*allo^.Poi[j-1]+Psih*allo^.Poi[j-1]/(1-allo^.Pos[j-1])) + Pih[g,j-1]*(1-Pssh)*(1-CBS[j-1])*(1-allo^.Pen-allo^.Pt-allo^.Pb-allo^.Pns + allo^.Pen*(1-allo^.Pes) +allo^.Pt*(1-allo^.Pts)+allo^.Pb*(1-allo^.Pbs) +allo^.Pns*(1-allo^.Pos[j-1])) + Ph[g,j-1]*(1-Pss)*(1-CBS[j-1])*((1-b)*(Pen*Pei) +b*(Pen*Pei+Pns*Poi[j-1]+Psi*Poi[j-1]/(1-Pos[j-1]))) + Pi[g,j-1]*(1-Pss)*(1-CBS[j-1])*((1-b)*(Pen*(1-Pes)) +b*(Pen*(1-Pes)+Pns*(1-Pos[j-1])+Psi)); Ph[g,j]:= Ph[g,j-1]*(1-Pss)*(1-CBS[j-1])*((1-b)*(1-Pen-Pk-Pt-Pb-Pns-Psi + Pk*(1-Pks-Pki)+Pt*(1-Pts-Pti)+Pb*(1-Pbs-Pbi) +Pns*(1-Pos[j-1]-Poi[j-1])+Psi*(1-Pos[j-1]-Poi[j-1])/(1-Pos[j-1])) + b*(1-Pen-Pk-Pt-Pb-Pns-Psi + Pk*(1-Pks-Pki)+Pt*(1-Pts-Pti)+Pb*(1-Pbs-Pbi))); Phh[g,j]:= Phh[g,j-1]*(1-Pssh)*(1-CBS[j-1])*(1-allo^.Pen-allo^.Pk-allo^.Pt-allo^.Pb-allo^.Pns-Psih + allo^.Pen*(1-allo^.Pes-allo^.Pei)+allo^.Pt*(1-allo^.Pts-allo^.Pti) +allo^.Pb*(1-allo^.Pbs-allo^.Pbi)+allo^.Pns*(1-allo^.Pos[j-1]-allo^.Poi[j-1]) +Psih*(1-allo^.Pos[j-1]-allo^.Poi[j-1])/(1-allo^.Pos[j-1])) + Ph[g,j-1]*(1-Pss)*(1-CBS[j-1])*((1-b)*(Pen*(1-Pes-Pei)) + b*(Pen*(1-Pes-Pei)+Pns*(1-Pos[j-1]-Poi[j-1]) +Psi*(1-Pos[j-1]-Poi[j-1])/(1-Pos[j-1]))); {calculation of the probability of having a homograft given that the patient} {is still alive} if Pih[g,j-1]+Phh[g,j-1]+Pi[g,j-1]+Ph[g,j-1] <= 0.0 then begin perhh := 0.5; AcceptSim := false end else perhh:= (Pih[g,j-1]+Phh[g,j-1])/(Pih[g,j-1]+Phh[g,j-1]+Pi[g,j-1]+Ph[g,j-1]); {determine the matrices qj en hqj with the aid of a second MARKOV-MODEL} for i:=1 to maxvry do if qj_1[i-1,1] > 0 then qj[i,1]:= qj_1[i-1,1] * (1-Klepfalen(j-1,i-1,beta,sigmac,sigmal)) * (1-Endo(i-1,Pe,ve,riscper)) * (1-Pns) * (1-Pk) * (1-Phgs) * (1-CBS[j-1]) else qj[i,1]:= 0; for k:=2 to maxop do begin Opovl:= 1 - Oper(j-1,k-1,eps); for i:= 0 to maxvry do begin Klep:= Klepfalen(j-1,i,beta,sigmac,sigmal); En:= Endo(i,Pe,ve,riscper); if Pe*ve > 1 then AcceptSim := false; if g > 1 then begin Kleph:= Klepfalen(j-1,i,allo^.beta,allo^.sigmac,allo^.sigmal); Enh:= Endo(i,allo^.Pe,allo^.ve,allo^.riscper); if allo^.Pe*allo^.ve > 1 then AcceptSim := false; end; if qj_1[i,k-1] > 0 then qj[0,k]:= qj[0,k] + qj_1[i,k-1] * (1-b) * (1-En) * (1-Phgs) * (Klep*Opovl + Pns*Opovl + Pk*(1-Pks))* (1-CBS[j-1]); if g > 1 then if (hqj_1[i,k-1] > 0) or (qj_1[i,k-1] > 0) then hqj[0,k]:= hqj[0,k] + perhh * hqj_1[i,k-1] * (1-allo^.Phgs) * (Kleph*Opovl + allo^.Pns*Opovl + Enh*(1-allo^.Pes))* (1-CBS[j-1]) + (1-perhh) * qj_1[i,k-1] * (1-Phgs) * ((1-b) * (1-Klep) * En*(1-Pes) * (1-Pns) + b * (Klep*Opovl + Pns*Opovl + En*(1-Pes))) * (1-CBS[j-1]); end; for i:= 1 to maxvry do begin if qj_1[i-1,k] > 0 then qj[i,k]:= qj_1[i-1,k] * (1-Klepfalen(j-1,i-1,beta,sigmac,sigmal)) * (1-Endo(i-1,Pe,ve,riscper)) * (1-Pns) * (1-Pk) * (1-Phgs)* (1-CBS[j-1]) else qj[i,k]:= 0; if g > 1 then if hqj_1[i-1,k] > 0 then hqj[i,k]:= hqj_1[i-1,k] * perhh * (1-Klepfalen(j-1,i-1,allo^.beta,allo^.sigmac,allo^.sigmal)) * (1-Endo(i-1,allo^.Pe,allo^.ve,allo^.riscper)) * (1-allo^.Pns) * (1-allo^.Phgs)* (1-CBS[j-1]) else hqj[i,k]:= 0; end; end; for k:= 1 to maxop do begin if qj_1[maxvry,k] > 0 then qj[maxvry,k]:= qj[maxvry,k] + qj_1[maxvry,k] * (1-Klepfalen(j-1,maxvry,beta,sigmac,sigmal)) * (1-Endo(maxvry,Pe,ve,riscper)) * (1-Pns) * (1-Pk) * (1-Phgs)* (1-CBS[j-1]); if g > 1 then if hqj_1[maxvry,k] > 0 then hqj[maxvry,k]:= hqj[maxvry,k] + hqj_1[maxvry,k] * perhh * (1-Klepfalen(j-1,maxvry,allo^.beta,allo^.sigmac,allo^.sigmal))* (1-CBS[j-1]) * (1-Endo(maxvry,allo^.Pe,allo^.ve,allo^.riscper)) * (1-allo^.Pns) * (1-allo^.Phgs); end; Opovl:= 1 - Oper(j-1,maxop,eps); for i:= 0 to maxvry do begin Klep:= Klepfalen(j-1,i,beta,sigmac,sigmal); En:= Endo(i,Pe,ve,riscper); if Pe*ve > 1 then AcceptSim := false; if g > 1 then begin Kleph:= Klepfalen(j-1,i,allo^.beta,allo^.sigmac,allo^.sigmal); Enh:= Endo(i,allo^.Pe,allo^.ve,allo^.riscper); if allo^.Pe*allo^.ve > 1 then AcceptSim := false; end; if qj_1[i,maxop] > 0 then qj[0,maxop]:= qj[0,maxop] + qj_1[i,maxop] * (1-b) * (1-En) * (1-Phgs) * (Klep*Opovl + Pns*Opovl + Pk*(1-Pks))* (1-CBS[j-1]); if g > 1 then if (hqj_1[i,maxop] > 0) or (qj_1[i,maxop] > 0) then hqj[0,maxop]:= hqj[0,maxop] + perhh * hqj_1[i,maxop] * (1-allo^.Phgs) * (Kleph*Opovl + allo^.Pns*Opovl + Enh*(1-allo^.Pes))* (1-CBS[j-1]) + (1-perhh) * qj_1[i,maxop] * (1-Phgs) * ((1-b) * (1-Klep) * En*(1-Pes) * (1-Pns) + b * (Klep*Opovl + Pns*Opovl + En*(1-Pes)) )* (1-CBS[j-1]); end; {patients that have died leave the population, so the probabilities do not} {sum up to 1 after the calculations above. We are interested in the probs} {of the #years since operation and the #operations given that the patient} {is alive and given the valve type. For this reason the probabilities in the} {matrix will be normalised} som:= 0; for i:= 0 to maxvry do for k:= 1 to maxop do som:= som + qj[i,k]; for i:= 0 to maxvry do for k:= 1 to maxop do if qj[i,k] > 0 then qj[i,k]:= qj[i,k]/som; if g > 1 then begin som1:= som; som:= 0; for i:= 0 to maxvry do for k:= 1 to maxop do som:= som + hqj[i,k]; for i:= 0 to maxvry do for k:= 1 to maxop do if hqj[i,k] > 0 then hqj[i,k]:= hqj[i,k]/som; som:= som+(1-perhh)*som1; end; {determine operation mortality for next year given the distribution for} {#previous operations} Pos[j]:= 0; allo^.Pos[j]:= 0; for i:= 0 to maxvry do for k:= 1 to maxop do begin if qj[i,k] > 0 then Pos[j]:= Pos[j] + qj[i,k] * Oper(j,k,eps); if g > 1 then if hqj[i,k] > 0 then allo^.Pos[j]:= allo^.Pos[j] + hqj[i,k] * Oper(j,k,eps); end; Poi[j]:= (1-Pos[j])*Pos[j]; if g > 1 then allo^.Poi[j]:= (1-allo^.Pos[j])*allo^.Pos[j]; end; {:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::} {determine life expectancy with and without morbidity} cum:= 0; cumi:= 1; Ph[g,age-1]:= 1; Phh[g,age-1]:= 0; if AcceptSim = true then begin gemov[z,g] := 0; gemgez[z,g]:= 0; for i:= age to 109 do begin p1:= (i-age-0.5)*(Po[g,i]-cum); cum:= Po[g,i]; if p1 < 0 then p1:= 0; gemov[z,g]:= gemov[z,g] + p1; p2:= Ph[g,i-1]+Phh[g,i-1]-Ph[g,i]-Phh[g,i]; {kans op gezond t/m i-e jaar} cumi:= cumi - p2; gemgez[z,g]:= gemgez[z,g] + (i-age-0.5)*p2; if gemgez[z,g] < 0 then gemgez[z,g]:= 0; end; writeln(g:1,' lf exp ',gemov[z,g]:5:3 {,', gez ',gemgez[z,g]:6:2}); end; end; if g = 1 then allo^:= tiep^; end; dispose(tiep); dispose(allo); if AcceptSim = true then begin for g := 1 to aantkleppen do begin k := 1; for j := 1 to aantkleppen do if ( g <> j ) then begin if ( z = 1 ) then begin gem_verschil[g,j] := 0; sigma_verschil[g,j] := 0; end; verschil[g,j] := gemov[z,g] - gemov[z,j]; if ( verschil[g,j] > 0 ) then rangtussenkleppen[g,j] := rangtussenkleppen[g,j] + 1; {How often is valve g preferred to valve j} if ( verschil[g,j] > vrschl_lft ) then rang_plusdelta[g,j] := rang_plusdelta[g,j] + 1; if ( verschil[g,j] < -vrschl_lft ) then rang_mindelta[g,j] := rang_mindelta[g,j] + 1; gem_verschil[g,j] := gem_verschil[g,j] + verschil[g,j]; sigma_verschil[g,j] := sigma_verschil[g,j] + SQR( verschil[g,j] ); if ( verschil[g,j] < 0 ) then k := k + 1; {determine rank order of four valve types} end; if ( k = 1 ) then rangklep[g].eerste := rangklep[g].eerste + 1 else if ( k = 2 ) then rangklep[g].tweede := rangklep[g].tweede + 1 else if ( k = 3 ) then rangklep[g].derde := rangklep[g].derde + 1 else rangklep[g].vierde := rangklep[g].vierde + 1; end; i:=1; for g:=1 to aantkleppen do for j:=g+1 to aantkleppen do begin dif_lifexp[z,i] := gemov[z,g] - gemov[z,j]; i:=i+1; end; z:=z+1; end; end; {end loop simulation} {Determination of medians of the difference in life expectancy between} {two valve types} for g:=1 to 6 do begin for z:=1 to simulatie do for j:=z+1 to simulatie do if dif_lifexp[z,g] > dif_lifexp[j,g] then begin h := dif_lifexp[z,g]; dif_lifexp[z,g] := dif_lifexp[j,g]; dif_lifexp[j,g] := h; end; j := trunc(simulatie/2); median[g] := dif_lifexp[j,g]; end; writeln(out[r],'patient is a ',sex,' of age ',age:2); writeln(out[r]); writeln(out[r],'valve type % first % second % third % fourth'); for g := 1 to aantkleppen do begin write(out[r],g:2,' ',(rangklep[g].eerste/simulatie)*100:4:1,' '); write(out[r],(rangklep[g].tweede/simulatie)*100:4:1,' '); write(out[r],(rangklep[g].derde/simulatie)*100:4:1,' '); writeln(out[r],(rangklep[g].vierde/simulatie)*100:4:1); end; writeln(out[r]); for g := 1 to aantkleppen do {determine means and standard deviations} begin for j := g+1 to aantkleppen do begin if ( g = j ) then gem_verschil[g,j] := 0 else begin gem_verschil[g,j] := gem_verschil[g,j] / simulatie; sigma_verschil[g,j] := SQRT( ( sigma_verschil[g,j] - simulatie * SQR( gem_verschil[g,j] ) ) / ( simulatie - 1 ) ); end; end; end; writeln(out[r]); writeln(out[r],'delta = ',vrschl_lft:5:3); writeln(out[r],' mu , sigma , > 0 , <-delta , -delta< delta '); for g := 1 to aantkleppen do { printen verschil met drempelwaarde } for j := g+1 to aantkleppen do begin writeln(out[r],g:2,j:2,' ',gem_verschil[g,j]:6:3,' , ',sigma_verschil[g,j]:5:3,' , ', (rangtussenkleppen[g,j]/simulatie)*100:5:1,' , ', (rang_mindelta[g,j]/simulatie)*100:5:1,' , ', ( 100 - (rang_mindelta[g,j]/simulatie)*100 - (rang_plusdelta[g,j]/simulatie)*100 ):6:1,' , ', (rang_plusdelta[g,j]/simulatie)*100:5:1); end; writeln(out[r]); write(out[r],'Medians of the differences between the life expectancies'); writeln(out[r],' of the valve pairs:'); for g:=1 to 6 do writeln(out[r],median[g]:5:3); writeln(out[r]); if age = 20 then begin age := 40; vrschl_lft := 0.91; end else if age = 40 then begin age := 55; vrschl_lft := 0.525; end else if age = 55 then begin age := 70; vrschl_lft := 0.14; end else if age = 70 then age := 109 else begin writeln('choose from ages: 20, 40, 55 of 70.'); readln(age); end; end until ( age = 109 ); close(out[r]); end; end; (* Bereken *) (* ........................ MAIN PROGRAM ........................... *) begin bereken; end.