Search

Simulation of hemophilia bleeds

A SAS simulation study describing bleeds in hemophilia type A and B patients. Based on the “proof-of-concept study” Randomized, prospective clinical trial of recombinant factor VIIa for secondary prophylaxis in hemophilia patients with inhibitors by Konkle et. al. published in Journal of Thrombosis and Haemostasis, 5: 1904-1913 in 2007.

The epidemiology and design of the study are not in accordance with CONSORT and SPIRIT as of December 2019. We observe 22 included patients (11 in each group) allocated from 20 sites in 11 countries. Neither the randomization nor the analysis is stratified by center, and it is possible to question the logistics of the experiment. People have asked if a condition on the number of bleeds (above 4 for each month in the pre-prophylaxis period) introduces bias: The study is a comparison of the effect of dose-difference not the effect of the drug, although the mean of bleeds in the first period is obviously higher than in the second period. Suspicion should be aimed at analysis methods and exposition of the results. The primary end-point is not even well defined or reported, which is obvious since the number of bleeds per period differs in tables and supplementary graphics of bleeds per month per patient in each group. Serum level FVII is discussed briefly and not compared to the administration of rFVIIa. Outliers are discussed based on clinical knowledge although there is no identification of outliers based on pattern analysis.

NovoSeven is a well-discussed medication. The active ingredient has been on the market many years before the introduction of the Novo Nordisk product, which was abandoned in 2012 for use outside clinical trials. This may be the reason behind poor design and loose exposition of results in the above-mentioned paper and maybe other studies and consecutive guidelines related to the application of the product. Updated TD or LD studies and strict adherence to modern CONSORT and SPIRIT guidelines would leave other more precise conclusions and recommendations based on an updated design for relevant dose level effect studies, more appropriate selection, and allocation algorithms and statistical analysis methods. The latter includes the use of generalized linear mixed models and statistical modeling of conditional probabilities for the exact unbiased estimation of difference in effect between pre-prophylaxis period and prophylaxis period.

The simulation and estimation routine:

Time series are simulated in a balanced RCT based on a Markov chain assumption. Results from likelihood ratio tests for group effect are analyzed with GLIMMIX and output to a file sim.xlsx

A professional simulation would include graphics with power calculations using different sets of input values, and a qualitative assessment of the simulated data compared to clinical settings.

The sketch do not contain tests on input values and output tables. With a 0.1 (10%) difference between groups vi observe a power close to 0.85 (85%) at alpha level 0.01 (1%). This is not a replica of the conditions prior to or results from the study mentioned above for which I have no raw data. The sketch provides a framework to simulate and analyze the experiment.


Discussion of mean-values lead to a rather straight forward calculation of conditional mean values. Calculation of unconditional mean values are not of interest themselves. With the given design we expect a sufficient number of bleeds in the evaluation of the difference in effect related to drug dose. We fit a regression line to the data, and observe that the number of bleeds is decreasing due to other factors than the effect of the drug in question. Thus the effect of the drug may be introduction of increased variation seen as unstable paths in the depiction of monthly number of bleeds. Making it difficult to discuss an effect of the drug itself compared to placebo. A much more positive effect on patients with a lower number of bleeds may be present, and we also observed a few cases of a complete stop in the prophylaxis period.


CODEBLOCK 1 (first graph):

%macro graph;
data graph; set _NULL_; run;
%do day=1 %to 270;
%do group=1 %to 2;
data pred;
set Ests;
retain eff day group;
day=&day;
group=&group;
if Effect EQ "intervention" AND intervention EQ 0 AND &day LE 90 then eff=Estimate;
if Effect EQ "intervention" AND intervention EQ 1 AND &day GT 90 then eff=Estimate;
if &group EQ 2 AND grp EQ 2 AND Effect EQ "grp" then eff=eff+Estimate;
if &day LE 90 AND Effect EQ "pre" then eff=eff+&day*Estimate;
if &day GT 90 AND Effect EQ "pro" then eff=eff+(&day-90)*Estimate;
if &group EQ 2 AND grp EQ 2 AND &day GT 90 AND Effect EQ "pro*grp" then eff=eff+(&day-90)*Estimate;
if &group EQ 1 AND grp EQ 1 AND &day GT 180 AND Effect EQ "post*grp" then eff=eff+(&day-180)*Estimate;
if &group EQ 2 AND grp EQ 2 AND &day GT 180 AND Effect EQ "post*grp" then eff=eff+(&day-180)*Estimate;
eff=eff;
if Effect EQ "post*grp" AND grp EQ 1 then output;
run;
data graph;
set graph pred;
run;
%end;
%end;
%mend;

%macro bleedsim(pbleed_1,pbleed_2,inc,cut,rate,retention,nmb_patients,grp);
*Simulation of bleed pr day for each period;
*Assuming a reduction in first month of period 2;
*i.e. 'pbleed_1' is less than 'pbleed_2';
*and both are numbers between 0 and 1 given proportion of non-bleeds;
*Reduction is gradually cancelled (using a linear trend);
*from day 'cut' in period two;
*'rate' is a factor describing speed of cancellation;
*We condition number of bleeds pr month in period 1 >= 4;
*'retention' is retention probability (drug was administrated on the given day);
*'inc' is increase in bleeds given as decimal number;
*if bleed was observed the previous day;

data pt;
retain bin choice bleed month totpat ok b1p ptnr;
grp=1;
totpat=0; ptnr=0;
do while (1);
if totpat GT &nmb_patients/2 then grp=2;
ptnr=ptnr+1;
ok=1; b1p=0; month=0;
do sample=1 to 90;
period=1;
day=sample;
month=floor((day-1)/30)+1;
call streaminit(0);
seed = ceil( (2**31 - 1)*rand("uniform") );
retention=(ranuni(seed) LE &retention);
if bleed EQ 0 then do; 
call streaminit(0);
seed = ceil( (2**31 - 1)*rand("uniform") ); 
choice=0; bin=ranuni(seed); bleed=(bin GE &pbleed_1); end;
else do; call streaminit(0);
seed = ceil( (2**31 - 1)*rand("uniform") );
choice=1; bin=ranuni(seed); bleed=(bin GE (&pbleed_1-&inc)); end;
if retention then do;
b1p=b1p+bleed;
end;
if day in (30 60 90) then do; 
if b1p LT 4 then ok=0;
b1p=0;
end;
if retention then output;
ok=ok;
end;
if ok GT 0 then totpat=totpat+1;
ok=ok;
if totpat LE &nmb_patients then do;
do sample=1 to &cut;
day=sample;
call streaminit(0);
seed = ceil( (2**31 - 1)*rand("uniform") );
retention=(ranuni(seed) LE &retention);
if bleed EQ 0 then do; 
call streaminit(0);
seed = ceil( (2**31 - 1)*rand("uniform") );
choice=0; bin=ranuni(seed); bleed=(bin GE (&pbleed_2-&grp*(grp EQ 1))); end;
else do; 
call streaminit(0);
seed = ceil( (2**31 - 1)*rand("uniform") );
choice=1; bin=ranuni(seed); bleed=(bin GE &pbleed_2-&grp*(grp EQ 1)-&inc); end;
period=2;
if sample GT 90 then period=3;
month=floor((day-1)/30)+1;
if retention then output;
end;
do sample=(&cut+1) to 180;
day=sample;
call streaminit(0);
seed = ceil( (2**31 - 1)*rand("uniform") );
retention=(ranuni(seed) LE &retention);
if bleed EQ 0 then do; 
call streaminit(0);
seed = ceil( (2**31 - 1)*rand("uniform") );
choice=0; bin=ranuni(seed); bleed=(bin GE (&pbleed_1+(1-&rate)*((&pbleed_2-&grp*(grp EQ 1))-&pbleed_1)*((190-&cut)-day+&cut)/(190-&cut))); end;
else do; 
call streaminit(0);
seed = ceil( (2**31 - 1)*rand("uniform") );
choice=1; bin=ranuni(seed); bleed=(bin GE ((&pbleed_1-&inc)+(1-&rate)*((&pbleed_2-&grp*(grp EQ 1))-&pbleed_1)*((190-&cut)-day+&cut)/(190-&cut))); end;
period=2;
if sample GT 90 then do; period=3; day=day-90; end;
month=floor((day-1)/30)+1;
if retention then output;
end;
end;
if totpat EQ &nmb_patients then leave;
end;
run;
proc sort data=pt out=patlst; by ptnr ok; run;
proc sort data=patlst nodupkey; by ptnr; run;
data patlst; set patlst; if ok EQ 0 then delete; run;
proc sort data=pt; by ptnr period month; run;
data patienter;
merge patlst(in=l) pt(in=p);
by ptnr;
if l EQ 1 AND p EQ 1 then output;
run;
proc means data=patienter sum;
class period month;
var bleed;
run;
proc means data=patienter sum noprint;
by ptnr;
class period month;
var bleed;
output out=table;
id grp;
run;
proc format library=work;
value timevar
11=1
12=2
13=3
21=4
22=5
23=6
31=7
32=8
33=9
;
run;
data table;
set table;
trials=_FREQ_;
bleed=bleed*_FREQ_;
time=period*10+month;
timevar=input(put(time,timevar.),8.);
if missing(period) OR missing(month) OR _STAT_ NE "MEAN" then delete;
run;
ods graphics on;
proc sgplot data=table;
series x=timevar y=bleed / group=ptnr;
run;
ods graphics off;
proc sort data=patienter; by ptnr period month day; run;

data analysedata;
set patienter;
pre=0;
if period EQ 1 then pre=day;
pro=0;
if period EQ 2 then pro=day;
post=0;
if period EQ 3 then do;
pro=day+90;
post=day;
end;
intervention=0;
if period GT 1 then intervention=1;
run;
ods exclude all;
proc glimmix data=analysedata method=mspl;
class ptnr grp(ref="1") intervention;
model bleed=intervention grp pre pro pro*grp post*grp /noint SOLUTION CHISQ CORRB;
random _RESIDUAL_ / subject = ptnr type = ARMA(1,1); 
ods output ParameterEstimates=Ests FitStatistics=Fit1;
run;
%graph;
/*proc export data=graph DBMS=xlsx
   outfile="/folders/myfolders/graph.xlsx" replace;
run;*/
proc sgplot data=graph;
series x=day y=eff / group=group;
run;
/*proc export data=Ests DBMS=xlsx
   outfile="/folders/myfolders/ests.xlsx" replace;
run;*/
proc glimmix data=analysedata method=mspl;
class ptnr intervention;
model bleed=intervention pre pro /noint SOLUTION CHISQ CORRB;
random _RESIDUAL_ / subject = ptnr type = ARMA(1,1); 
ods output FitStatistics=Fit2;
run;
ods exclude none;
data Fit;
set Fit1(in=f) Fit2(in=r);
if f EQ 1 then model=1;
if r EQ 1 then model=2;
if Descr EQ "-2 Log Likelihood" then output;
run;
data Fit(keep=ll1 ll2 chisq p);
set Fit;
retain ll1;
if model EQ 1 then ll1=Value;
if model EQ 2 then do;
ll2=Value;
chisq=ll2-ll1;
format chisq 13.2;
p=1-CDF('chisquare',chisq,3);
format p 13.2;
output;
end;
run;
data sim;
set sim Fit;
run;
%mend;

%macro powersim(nmb);
data sim; set _NULL_; run;
%do i=1 %to &nmb;
%bleedsim(0.75,0.95,0.20,60,0.1,0.85,22,0.15);
%end;
%mend;
%powersim(100);
proc export data=sim DBMS=xlsx
   outfile="/folders/myfolders/sim.xlsx" replace;
run;


CODEBLOCK 2 (second graph):


/* Generated Code (IMPORT) */
/* Source File: novoint.csv */
/* Source Path: /folders/myfolders/Novo */
/* Code generated on: 12/21/19, 4:27 PM */
%web_drop_table(WORK.IMPORT);
FILENAME REFFILE '/folders/myfolders/Novo/novoint.csv';
PROC IMPORT DATAFILE=REFFILE
 DBMS=CSV
 OUT=WORK.IMPORT;
 GETNAMES=YES;
 delimiter=";";
RUN;
PROC CONTENTS DATA=WORK.IMPORT; RUN;
%web_open_table(WORK.IMPORT);
data import2;
set import;
grp=1;
if substr(ptnr,1,3) EQ "270" then grp=2;
cond_1=0;cond_2=0;cond_3=0;cond_4=0;
cond_5=0;cond_6=0;cond_7=0;cond_8=0;cond_9=0;
o_1=0;o_2=0;o_3=0;o_4=0;o_5=0;
o_6=0;o_7=0;o_8=0;o_9=0;
if _1 GE 4 then cond_1=1;
if _2 GE 4 then cond_2=1;
if _3 GE 4 then cond_3=1;
if _4 GE 4 then cond_4=1;
if _5 GE 4 then cond_5=1;
if _6 GE 4 then cond_6=1;
if _7 GE 4 then cond_7=1;
if _8 GE 4 then cond_8=1;
if _9 GE 4 then cond_9=1;
if _1 GE 4 then o_1=_1;
if _2 GE 4 then o_2=_2;
if _3 GE 4 then o_3=_3;
if _4 GE 4 then o_4=_4;
if _5 GE 4 then o_5=_5;
if _6 GE 4 then o_6=_6;
if _7 GE 4 then o_7=_7;
if _8 GE 4 then o_8=_8;
if _9 GE 4 then o_9=_9;
run;
proc means data=import2 noprint;
class grp;
var _1 _2 _3 _4 _5 _6 _7 _8 _9 
cond_1 cond_2 cond_3 cond_4 cond_5 cond_6 cond_7 cond_8 cond_9
o_1 o_2 o_3 o_4 o_5 o_6 o_7 o_8 o_9;
output out=table(where=(grp in (1 2) AND _STAT_ in ("MEAN")));
run;
data estimates(keep=grp cm_1-cm_9 month bleed);
set table;
_1=_FREQ_*_1;_2=_FREQ_*_2;_3=_FREQ_*_3;_4=_FREQ_*_4;_5=_FREQ_*_5;_6=_FREQ_*_6;_7=_FREQ_*_7;
_8=_FREQ_*_8;_9=_FREQ_*_9;
cm_1=o_1/cond_1;cm_2=o_2/cond_2;cm_3=o_3/cond_3;cm_4=o_4/cond_4;
cm_5=o_5/cond_5;cm_6=o_6/cond_6;cm_7=o_7/cond_7;cm_8=o_8/cond_8;cm_9=o_9/cond_9;
if missing(cm_1) then cm_1=0; 
if missing(cm_2) then cm_2=0;
if missing(cm_3) then cm_3=0;
if missing(cm_4) then cm_4=0;
if missing(cm_5) then cm_5=0;
if missing(cm_6) then cm_6=0;
if missing(cm_7) then cm_7=0;
if missing(cm_8) then cm_8=0;
if missing(cm_9) then cm_9=0;
month=1; bleed=cm_1;output;
month=2; bleed=cm_2;output;
month=3; bleed=cm_3;output;
month=4; bleed=cm_4;output;
month=5; bleed=cm_5;output;
month=6; bleed=cm_6;output;
month=7; bleed=cm_7;output;
month=8; bleed=cm_8;output;
month=9; bleed=cm_9;output;
run;
proc sgplot data=estimates;
series x=month y=bleed /group=grp;
reg x=month y=bleed / ;
run;



1 view

Recent Posts

See All

©2020 by Danish Institute for Data Science. Proudly created with Wix.com