Quality Adjusted Life Years calculation method - Area under the curve, change from baseline, regression methods

[Tabitha Green]

I need to calculate quality adjusted life years. I understand that the most common way to do that is to calculate the area under the curve. I have written the code and it works beautifully.

* Calculate QALYs for each participant
gen qaly0 = eq5d5lscore0*(0.5)  // baseline QALY score
gen qaly6 = eq5d5lscore6*(0.5) // QALY score at 6 months
gen qaly12 = eq5d5lscore12*(0.5)  // QALY score at 12 months

* Account for lost to follow-up
replace qaly6 = . if ltfu6mth == 1  // assume QALY score at 6 months is same as baseline if lost to follow-up
replace qaly12 = . if ltfuv12mth == 1  // assume QALY score at 12 months is same as baseline if lost to follow-up

* Account for withdrawals
replace qaly6 = . if wd6mth == 1  // set QALY score at 6 months to missing if participant withdrew
replace qaly12 = . if wd12mth == 1  // set QALY score at 12 months to missing if participant withdrew

* Account for deaths
replace qaly6 = . if died6mth == 1  // set QALY score at 6 months to 0 if participant died
replace qaly12 = . if died12mth == 1  // set QALY score at 12 months to 0 if participant died

* Calculate area under the curve to get QALY
gen aucqaly = (qaly0 + qaly6)/2 + (qaly6 + qaly12)/2 if !missing(qaly6, qaly12)

* Calculate total QALYs for the study
sum aucqaly, meanonly

However, I need to also calculate it using two other methods, the change from baseline and the regression method. I have been trying to work out from the equations provided in a paper what the code would be, but I can't work it out. These two papers gives the three techniques: Manca et al. Estimating mean QALYs in trial-based cost-effectiveness analysis: The importance of controlling for baseline utility. Health Economics. 2005;14(5):487-96 and McClure et al. Modifying the quality-adjusted life year calculation to account for meaningful change in health-related quality of life: insights from a pragmatic clinical trial. Eur J Health Econ. 2021;22(9):1441-51


The equation for calculating the QALY-AUC (which is my code above, which I've got sorted out) is:

But I can't work out how I would write the code to calculate the QALY-CFB or QALY-R versions. I just can't tend to "see" how the equations relate to the code. Is anyone here able to help me with what the code would be to calculate the QALY-CFB and the QALY-R?

Variables:
eq5d5lscore0 eq5d5lscore6 eq5d5lscore12 utilities at baseline, 6 months and 12 months.
ltfu6mth ltfu12mth: Lost to follow up by 6 months and 12 months respectively
wd6mth wd12mth: Withdrawn by 6 month or 12 month follow up respectively respectively
died6mth died12mth: Died by the 6 month or 12 month follow up respectively

* Example generated by -dataex-. For more info, type help dataex
clear
input float pat_id long group float(eq5d5lscore0 eq5d5lscore6 eq5d5lscore12) byte(wd6mth ltfu6mth died6mth wd12mth ltfuv12mth died12mth)
 1 1       .658      .917 .797 . . . . . .
 2 2 -.04599999      .861 .518 . . . . . .
 3 2       .385      .852  .88 . . . . . .
 4 2      -.124         1    1 . . . . . .
 5 1       .917      .852 .852 . . . . . .
 6 1       .917      .859 .857 . . . . . .
 7 2 -.04599999         .    . . 1 . . . .
 8 2       .832      .934 .956 . . . . . .
 9 1       .588         .    . . . . . . 1
10 1        .85      .917 .961 . . . . . .
11 2       .852      .956    . . . . . . .
12 2       .457       .66 .739 . . . . . .
13 1       .385         .    . . . 1 . . .
14 1       .754         .    . 1 . . . . .
15 1          1         1 .956 . . . . . .
16 2       .163         .    . . . 1 . . .
17 1        .88         .    . . . . 1 . .
18 2       .917      .917 .889 . . . . . .
19 1        .88         . .887 . . . . . .
20 1       .719      .956 .889 . . . . . .
21 2       .848      .396    . . . . . . .
22 2        .97         .    . 1 . . . . .
23 1       .961         . .731 . . . . . .
24 2       .039         . .794 . . . . . .
25 1       .146         .    . . 1 . . . .
26 1       .623         .    . . 1 . . . .
27 2       .822         . .933 . . . . . .
28 2       .871      .961 .933 . . . . . .
29 2       .222      .208 .887 . . . . . .
30 1       .623         .    . . . 1 . . .
31 2       .871      .536  .88 . . . . . .
32 1       .776         .    . 1 . . . . .
33 1       .428      .735 .178 . . . . . .
34 2       .057      .851 .784 . . . . . .
35 2       .862         . .685 . . . . . .
36 2       .724      .934 .968 . . . . . .
37 2       .296      .889 .901 . . . . . .
38 1       .499      .831 .831 . . . . . .
39 1       .484 .04200001    . . . . . . .
40 2       .432         . .673 . . . . . .
41 1          .      .942 .667 . . . . . .
42 1       .889         1    1 . . . . . .
43 1  .26700002         .    . . . 1 . . .
44 2       .377         1 .933 . . . . . .
45 2       .766      .875 .804 . . . . . .
46 1       .697      .956 .956 . . . . . .
47 1       .961      .624 .831 . . . . . .
48 2       .335      .739 .887 . . . . . .
49 1       .802      .852 .917 . . . . . .
50 1       .889         .    . . 1 . . . .
51 2       .496         .    . . . . . . 1
52 1       .933         .    . 1 . . . . .
53 2       .374         .    . . . . 1 . .
54 1       .735      .831 .901 . . . . . .
55 2       .697      .216 .466 . . . . . .
56 2          .         . .682 . . . . . .
57 2       .903      .804    . . . . . . .
58 2          .      .885 .919 . . . . . .
59 1       .685      .852 .733 . . . . . .
60 2          1      .956    1 . . . . . .
61 1          1         .    . . . 1 . . .
62 2       .852      .887    1 . . . . . .
63 1       .766         . .701 . . . . . .
64 2       .325      .396 .822 . . . . . .
65 1       .596      .625 .616 . . . . . .
66 1       .555      .308    . . . . . . .
67 1       .309         . .956 . . . . . .
68 2       .495      .673 .117 . . . . . .
69 2       .396         1 .834 . . . . . .
70 2        .83         .    . . . 1 . . .
71 1       .765       .97 .956 . . . . . .
72 1       .618      .956 .933 . . . . . .
73 2       .871         .    . . 1 . . . .
74 2       .253         .    . 1 . . . . .
75 2       .552      .673    . . . . . . 1
76 2       .744      .885 .756 . . . . . .
77 1       .318      .895 .846 . . . . . .
78 2          1         .    . . 1 . . . .
79 1       .956      .917 .917 . . . . . .
80 2       .851      .823 .424 . . . . . .
81 2       .379      .113 .889 . . . . . .
82 1       .159      .809 .178 . . . . . .
83 1       .366      .848 .255 . . . . . .
84 2       .961      .933 .961 . . . . . .
85 1       .887      .543 .704 . . . . . .
86 1       .572         1 .794 . . . . . .
87 1  .10900001      .961 .917 . . . . . .
88 1       .848      .906    . . . . . . .
89 1        .52      .731 .219 . . . . . .
90 1       .837         .    . . . 1 . . .
91 1       .306         .    . 1 . . . . .
92 2       .681         .    . 1 . . . . .
93 1       .844      .762 .887 . . . . . .
94 2       .739         .    . . 1 . . . .
95 2       .795         . .472 . . . . . .
96 1       .784         1 .919 . . . . . .
97 2       .572         . .219 . . . . . .
98 2       .665      .487 .495 . . . . . .
end
label values group group_numeric
label def group_numeric 1 "GROUP A", modify
label def group_numeric 2 "GROUP B", modify

[Joseph Coveney]

I need to calculate quality adjusted life years. . . . I have written the code and it works beautifully.

Are you sure? The code you show for imputing values for six- and twelve-month visits for LTFU doesn't seem to do what its comment says it's supposed to. Ditto for deaths.

But I can't work out how I would write the code to calculate the QALY-CFB or QALY-R versions. I just can't tend to "see" how the equations relate to the code. Is anyone here able to help me with what the code would be to calculate the QALY-CFB and the QALY-R?

Here's my crack at it below.

I imputed missing values as last observation carried forward instead of the baseline score for LTFU (and simple missed visits). And I didn't divide all of the QOL scores by two as you did, but otherwise the code should be okay.

I also shortened the variable names for laziness. Both do-file and log file are attached if you're interested further.

version 18.0

clear *

quietly input byte(pat_id group) ///
    double(eq5d5lscore0 eq5d5lscore6 eq5d5lscore12) ///
    byte(wd6mth ltfu6mth died6mth wd12mth ltfuv12mth died12mth)
<redacted for brevity>
end

rename pat_id pid

assert inlist(group, 1, 2)
generate byte trt = group == 2
drop group
label define Groups 0 A 1 B
label values trt Groups

rename eq5d5lscore* sco*
rename wd*mth wid*
rename ltfu*mth ltf*
rename ltfv12 ltf12 // n.b.
rename died*mth dea*

foreach var of varlist w* l* d* {
    assert inlist(`var', 1, .)
    quietly replace `var' = 0 if missing(`var')
}

*
* Begin here
*
// Impute last observation carried forward for LTFU or missed-visit
quietly generate double isc0 = sco0
quietly generate double isc6 = cond(missing(sco6) & !wid6 & !dea6, sco0, sco6)
quietly generate double isc12 = cond(missing(sco12) & !wid12 & !dea12, isc6, sco12)

// Impute zero for death
quietly replace isc6 = 0 if dea6
quietly replace isc12 = 0 if dea6 | dea12

// QALY-CFB
quietly generate double cfb = (isc6 - sco0) / 2 + (isc12 - sco0) / 4
quietly ttest cfb, by(trt)
display in smcl as text "Incremental QALY-CFB: " as result %05.3f ///
    r(mu_1) - r(mu_2) " ± " %05.3f r(se) ///
    as text ", n = " as result r(N_1) ", " r(N_2)

// QALY-AUC
quietly reshape long sco isc dea ltf wid, i(pid) j(tim)
bysort pid (tim): integ isc tim, double generate(auc) trapezoid
quietly by pid: replace auc = . if mi(sco[1])
quietly replace auc = auc / 12 // QALMonth → QALY
quietly ttest auc if tim == 12, by(trt)
display in smcl as text "Incremental QALY-AUC: " as result %05.3f ///
    r(mu_1) - r(mu_2) " ± " %05.3f r(se) ///
    as text ", n = " as result r(N_1) ", " r(N_2)

// QALY-R
quietly by pid: generate double sco0 = sco[1]
regress auc i.trt c.sco0 if tim == 12
display in smcl as text "Incremental QALY-R: " as result %05.3f ///
    _b[1.trt] " ± " %05.3f _se[1.trt] ///
    as text ", df = " as result e(df_r)

exit

[Tabitha Green]

Thank you so much, I really appreciate it (and identifying that I hadn't imputed the values correctly in my code. I went back and checked and your were correct).

Can I ask you a quick follow up? With the QALY-CFB, is the ttest line supposed to be at time 12 only too?

ttest cfb, by(trt)
ttest cfb if tim == 12, by(trt)

[Tabitha Green]

Oh sorry, that's because I re-ran the code a second time and missed the wide to long conversion was in the code after that line.