Hi, I am trying to create a Bland Altman plot for dataset. I am comparing continuous values from two devices and my data contains maximum of 2 observations per patient. How do I create a Bland Altman plot while accounting for correlation at patient level?
I have written the following code and tried to use the bootstrap technique but not sure if it is correct. Can someone please help me through this? Thanks!
* Calculate the average and difference (bias)
gen avg = (ctrl_thickness_cor + ctrl_thickness_pent) / 2
gen bias = ctrl_thickness_cor - ctrl_thickness_pent
* Summarize the bias and calculate LoA
summarize bias, detail
scalar mean_bias = r(mean)
scalar sd_bias = r(sd)
scalar ul_bias = r(mean) + 1.96 * r(sd)
scalar ll_bias = r(mean) - 1.96 * r(sd)
* Define a program for bootstrapping
capture program drop blandaltman
program blandaltman, rclass
summarize bias, detail
return scalar mean_bias = r(mean)
return scalar ul_bias = r(mean) + 1.96 * r(sd)
return scalar ll_bias = r(mean) - 1.96 * r(sd)
end
* Run the bootstrap with clustering by record_id, using BCa CIs
bootstrap mean_bias=r(mean_bias) ul_bias=r(ul_bias) ll_bias=r(ll_bias), ///
cluster(record_id) reps(1000) saving(C:\CCT, replace) seed(1234): blandaltman
* Store bootstrap results
estat bootstrap, all
ereturn list
matrix list = e(b)
local actual_mean_bias: di %5.2f e(b)[1,1]
local actual_ll: di %5.2f e(b)[1,2]
local actual_ul: di %5.2f e(b)[1,3]
matrix list = e(ci_percentile)
local bs_meanll: di %5.2f e(ci_percentile)[1,1]
local bs_meanul: di %5.2f e(ci_percentile)[2,1]
local bs_lower_ll: di %5.2f e(ci_percentile)[1,2]
local bs_lower_ul: di %5.2f e(ci_percentile)[2,2]
local bs_upper_ll: di %5.2f e(ci_percentile)[1,3]
local bs_upper_ul: di %5.2f e(ci_percentile)[2,3]
* Debug: Display the CI values to check their width
di "Mean Bias: " `actual_mean_bias' " (" `bs_meanll' " to " `bs_meanul' ")"
di "LoA Upper: " `actual_ll' " (" `bs_lower_ll' " to " `bs_lower_ul' ")"
di "LoA Lower: " `actual_ul' " (" `bs_upper_ll' " to " `bs_upper_ul' ")"
twoway (scatter bias avg, mcolor("10 90 255") msize(small)) ///
(function y = `actual_mean_bias', range(400 700) lcolor(red) lwidth(medium) lpattern(shortdash)) /// * Mean line in bold red
(function y = `bs_meanll', range(400 700) lcolor(red) lwidth(medium) lpattern(dot)) /// * Mean ll line in bold red
(function y = `bs_meanul', range(400 700) lcolor(red) lwidth(medium) lpattern(dot)) /// * Mean ul line in bold red
(function y = `actual_ll', range(400 700) lcolor(blue) lwidth(medium) lpattern(shortdash)) /// * actual lowelimiy
(function y = `bs_lower_ll', range(400 700) lcolor(blue) lwidth(medium) lpattern(dot)) /// * Upper limit in bold blue
(function y = `bs_lower_ul', range(400 700) lcolor(blue) lwidth(medium) lpattern(dot)) /// * Lower limit in bold blue
(function y = `actual_ul', range(400 700) lcolor(blue) lwidth(medium) lpattern(shortdash)) /// * actual upper limit
(function y = `bs_upper_ll', range(400 700) lcolor(blue) lwidth(medium) lpattern(dot)) /// * Upper limit in bold blue
(function y = `bs_upper_ul', range(400 700) lcolor(blue) lwidth(medium) lpattern(dot)) ///
(function y = 0, range(400 700) lpattern(solid) lcolor(darkgray)) /// * Dotted line at y = 0
, yscale(range(-80 80)) /// * Set the y-axis range
ylabel(-80(20)80) /// * Adjust y-axis labels
xlabel(400(50)700)
I have written the following code and tried to use the bootstrap technique but not sure if it is correct. Can someone please help me through this? Thanks!
* Calculate the average and difference (bias)
gen avg = (ctrl_thickness_cor + ctrl_thickness_pent) / 2
gen bias = ctrl_thickness_cor - ctrl_thickness_pent
* Summarize the bias and calculate LoA
summarize bias, detail
scalar mean_bias = r(mean)
scalar sd_bias = r(sd)
scalar ul_bias = r(mean) + 1.96 * r(sd)
scalar ll_bias = r(mean) - 1.96 * r(sd)
* Define a program for bootstrapping
capture program drop blandaltman
program blandaltman, rclass
summarize bias, detail
return scalar mean_bias = r(mean)
return scalar ul_bias = r(mean) + 1.96 * r(sd)
return scalar ll_bias = r(mean) - 1.96 * r(sd)
end
* Run the bootstrap with clustering by record_id, using BCa CIs
bootstrap mean_bias=r(mean_bias) ul_bias=r(ul_bias) ll_bias=r(ll_bias), ///
cluster(record_id) reps(1000) saving(C:\CCT, replace) seed(1234): blandaltman
* Store bootstrap results
estat bootstrap, all
ereturn list
matrix list = e(b)
local actual_mean_bias: di %5.2f e(b)[1,1]
local actual_ll: di %5.2f e(b)[1,2]
local actual_ul: di %5.2f e(b)[1,3]
matrix list = e(ci_percentile)
local bs_meanll: di %5.2f e(ci_percentile)[1,1]
local bs_meanul: di %5.2f e(ci_percentile)[2,1]
local bs_lower_ll: di %5.2f e(ci_percentile)[1,2]
local bs_lower_ul: di %5.2f e(ci_percentile)[2,2]
local bs_upper_ll: di %5.2f e(ci_percentile)[1,3]
local bs_upper_ul: di %5.2f e(ci_percentile)[2,3]
* Debug: Display the CI values to check their width
di "Mean Bias: " `actual_mean_bias' " (" `bs_meanll' " to " `bs_meanul' ")"
di "LoA Upper: " `actual_ll' " (" `bs_lower_ll' " to " `bs_lower_ul' ")"
di "LoA Lower: " `actual_ul' " (" `bs_upper_ll' " to " `bs_upper_ul' ")"
twoway (scatter bias avg, mcolor("10 90 255") msize(small)) ///
(function y = `actual_mean_bias', range(400 700) lcolor(red) lwidth(medium) lpattern(shortdash)) /// * Mean line in bold red
(function y = `bs_meanll', range(400 700) lcolor(red) lwidth(medium) lpattern(dot)) /// * Mean ll line in bold red
(function y = `bs_meanul', range(400 700) lcolor(red) lwidth(medium) lpattern(dot)) /// * Mean ul line in bold red
(function y = `actual_ll', range(400 700) lcolor(blue) lwidth(medium) lpattern(shortdash)) /// * actual lowelimiy
(function y = `bs_lower_ll', range(400 700) lcolor(blue) lwidth(medium) lpattern(dot)) /// * Upper limit in bold blue
(function y = `bs_lower_ul', range(400 700) lcolor(blue) lwidth(medium) lpattern(dot)) /// * Lower limit in bold blue
(function y = `actual_ul', range(400 700) lcolor(blue) lwidth(medium) lpattern(shortdash)) /// * actual upper limit
(function y = `bs_upper_ll', range(400 700) lcolor(blue) lwidth(medium) lpattern(dot)) /// * Upper limit in bold blue
(function y = `bs_upper_ul', range(400 700) lcolor(blue) lwidth(medium) lpattern(dot)) ///
(function y = 0, range(400 700) lpattern(solid) lcolor(darkgray)) /// * Dotted line at y = 0
, yscale(range(-80 80)) /// * Set the y-axis range
ylabel(-80(20)80) /// * Adjust y-axis labels
xlabel(400(50)700)
