Ước tính Estrogen | juliablogger.com – Beragampengetahuan

Bởi: oxinabox.github.io

Đăng lại từ:

Là một người chuyển giới nữ tham gia HRT, tôi muốn biết nồng độ estradiol trong máu của mình vào tất cả các giờ trong ngày.
Điều này rất hữu ích vì đỉnh, đáy và mức trung bình đều có tác dụng.
Tuy nhiên, tôi chỉ được xét nghiệm máu một số lần hữu hạn mỗi ngày – thường là một lần.
Tôi không phải là bác sĩ y khoa, nhưng tôi là loại bác sĩ có thể áp dụng mô hình khoa học vào nhiệm vụ ước tính các đường cong dựa trên những quan sát hạn chế.
Tôi thành thật ngạc nhiên không ai đã làm điều này.
Sự giao thoa giữa máy tính khoa học và dân gian xuyên không là không tầm thường.
Rốt cuộc, vấn đề khó khăn nhất trong khoa học máy tính là phiền muộn giới tính.

Để làm điều này, chúng ta sẽ sử dụng lập trình xác suất, để phân phối các đường cong mức có thể.
Đây là một trường hợp sử dụng tuyệt vời cho lập trình xác suất.
Chúng tôi có rất nhiều kiến ​​thức về miền, nhưng kiến ​​thức về miền đó có một vài tham số mà chúng tôi không biết và chúng tôi chỉ có một ít dữ liệu.
Và điều quan trọng là chúng tôi hoàn toàn ổn với việc phân phối các câu trả lời, thay vì một câu trả lời duy nhất.

Đây là lần đầu tiên tôi sử dụng lập trình xác suất, mặc dù tôi đã biết về nó trong nhiều năm.
Tôi chỉ hiếm khi gặp một vấn đề mà nó sẽ trở nên tuyệt vời.
Nhưng đây là một.
Bài đăng trên blog này là cuối cùng tôi đã học và sử dụng lập trình xác suất.

Chúng tôi bắt đầu bằng cách xem xét hình dạng của đường cong.
Cách các cấp độ thay đổi theo thời gian là phần kiến ​​thức miền lớn nhất mà chúng ta có.
Chúng tôi muốn xây dựng một đường cong được tham số hóa trong đó hình dạng có thể thay đổi bằng cách điều chỉnh các tham số.

Trong bài đăng trên blog này về Estrogen ngậm dưới lưỡi so với uống, họ đã tính gần đúng hàm estradiol với một tuyến tính đến cực đại sau đó là sự phân rã theo cấp số nhân.

Tôi đã vẽ các điểm từ đường cong estradiol dưới lưỡi và đưa ra ước tính về hàm estradiol từ 1 giờ đến 24 giờ như 350.54*(HOURS^-0.907). Từ 0 đến 1 giờ, tôi ước tính tuyến tính mức estradiol, như 451*HOURS.

Hôm nay tôi quan tâm đến chức năng estradiol cho gel estradiol.
Nó có hình dạng giống nhau, vì vậy nên áp dụng một chiến lược tương tự.

Järvinen A, Granander M, Nykänen S, Laine T, Geurts P, Viitanen A (tháng 11 năm 1997). “Dược động học ở trạng thái ổn định của gel oestradiol ở phụ nữ sau mãn kinh: tác dụng của vùng bôi và rửa”. Br J Sản khoa Gynaecol. 104 Bổ sung 16:14–8. doi:10.111/j.1471-0528.1997.tb11562.x. PMID 9389778. S2CID 36677042.
đưa ra 3 đường cong cho liều duy nhất.
Cũng có sẵn từ Wikimedia Commons. Lô đất dưới đây:

3 đường cong này trong nghiên cứu đó là từ việc bôi gel lên các khu vực lớn hơn hoặc nhỏ hơn.
Đối với mục đích của chúng tôi mà không quan trọng.
Thay vào đó, chúng ta có thể coi chúng là 3 đường cong thực tế khác nhau, có một số khác biệt dựa trên các nguyên nhân khác nhau (có thể là khu vực ứng dụng, vị trí, sinh lý học, v.v.).
Sau đó chúng ta sẽ ước lượng các tham số của các đường cong đó.

Tuy nhiên, trước khi chúng tôi làm bất cứ điều gì khác, chúng tôi bắt đầu bằng cách thiết lập môi trường Julia của mình.

using Pkg: @pkg_str
pkg"activate --temp"
pkg"add Tables@1 Plots@1 Turing@0.21.12 MCMCChains@5.5 LsqFit@0.13.0 StatsPlots@0.15.4"

Đường cong chúng tôi muốn xác định nồng độ estradiol trong máu hiện tại c tại thời điểm t giờ sau khi bôi gel.
Nó được mô tả bởi 3 thông số:

• c_max: nồng độ đỉnh.
• t_max: thời gian cần thiết để đạt được nồng độ cực đại.
• halflife: thời gian để nồng độ giảm còn một nửa sau khi đạt cực đại.

(Cùng một đường cong, với các giá trị tham số khác nhau có thể hoạt động đối với các phương pháp quản trị khác)
Chúng tôi biểu thị điều này dưới dạng hàm bậc cao hơn của 3 tham số, trả về hàm thời gian t.

function single_dose(c_max, halflife, t_max)
    function(t)
        if t < t_max
            c_max/t_max * t
        else
            c_max * 2^(-(t-t_max)/halflife)
        end
    end
end

Hãy xem chúng tôi đã làm như thế nào.
Tôi sẽ vẽ biểu đồ dữ liệu từ Järvinen et al đối với các đường cong bằng cách sử dụng công thức của tôi, phù hợp nhất với sự kiểm tra của riêng tôi.
Tôi đang giảm tất cả dữ liệu từ Järvinen và cộng sự xuống 25 pg/mL, vì dữ liệu đó là từ những phụ nữ cis sau mãn kinh, những người tự sản xuất khoảng 25 pg/mL estradiol trước khi bạn tính đến HRT.
Chúng tôi chỉ muốn lập mô hình thành phần HRT.

using Plots
plot(layout=(1,3))

plot!(
    0:0.1:24, single_dose(132, 3.5, 3),
    label="A200 predicted"; linecolor=:red, subplot=1, yrange=(0,130),
)
scatter!(
    [0,1,2,3,4,6,8,10,12,16,24], [0,25,100,132,90,82,60,55,32,15,4],
    label="A200 actual", markercolor=:red, subplot=1, yrange=(0,130),
)

plot!(
    single_dose(100, 3.5, 2.5),
    label="A400 predicted", linecolor=:magenta, subplot=2, yrange=(0,130),
)
scatter!(
    [0,1,2,3,4,6,8,10,12,16,24], [25,35,70,75,55,45,35,32,22,15,4];
    label="A400 actual", markercolor=:magenta, subplot=2, yrange=(0,130),
)

plot!(
    single_dose(20, 3.5, 2.7),
    label="Amax predicted", linecolor=:blue, subplot=3, yrange=(0,130),
)
scatter!(
    [0,1,2,3,4,6,8,10,12,16,24], [5,14,17,20,12,10,5,2,5,5,4];
    label="Amax actual", markercolor=:blue, subplot=3, yrange=(0,130),
)

Bằng cách nhìn vào những ô này, nó có vẻ là một mô hình khá tốt.
Tất nhiên với đủ bậc tự do, bạn có thể nhét vừa một con voi.
Tuy nhiên, chúng tôi có 10 điểm và chỉ có 3 bậc tự do, trong đó chúng tôi chỉ thay đổi 2 trong số chúng trên 3 bộ dữ liệu.
Vì vậy, có vẻ như chúng tôi là tốt.

Nói chung là hợp lý về mặt sinh học.
Chúng tôi mong đợi sự hấp thụ ban đầu nhanh chóng, điều đó sẽ kết thúc vào một thời điểm nào đó sau vài giờ nữa.
Vì nó nhanh và ngắn, nên việc chúng ta mô hình hóa nó bằng cái gì không thực sự quan trọng, vì vậy tuyến tính là tốt.
Sau đó, chúng tôi mong đợi một cái đuôi tắt khi nó được tiêu thụ.
Điều hợp lý là tốc độ hấp thụ có liên quan đến số lượng còn lại – điều này cho thấy một số cấp số nhân.
Chúng tôi mong đợi điều này bởi vì càng nhiều estradiol xung quanh thì khả năng một trong số chúng chạm vào thụ thể và bị tiêu thụ càng cao. Sau đó, khi estradiol và các thụ thể được tiêu thụ, cơ hội va chạm sẽ giảm đi.
Điều này có thể là vô nghĩa, tôi không phải là nhà sinh vật học hệ thống.

Bây giờ tôi chỉ phù hợp với những đường cong bằng mắt.
Chúng ta có thể tìm thấy các tham số có khả năng nhất thông qua hồi quy bình phương nhỏ nhất.
Đối với điều này, chúng ta có thể sử dụng LsqFit.jl, thông thường tôi sẽ sử dụng Optim.jl và tự mình viết bài toán bình phương nhỏ nhất và sử dụng L-BFGS để giải quyết vấn đề đó.
Nhưng thuật toán Levenberg-Marquardt của LsqFit.jl phù hợp hơn.

using LsqFit

fit = curve_fit(
    (t,p)->single_dose(p...).(t),
    [0,1,2,3,4,6,8,10,12,16,24], # t values
    [0,25,100,132,90,82,60,55,32,15,4], # c
    [132, 3.5, 3]  # initial guess at parameters
)

c_max, halflife, t_max = fit.param
@show c_max halflife t_max

plot(
    0:0.1:24, single_dose(c_max, halflife, t_max),
    label="A200 predicted"; linecolor=:red,
)
scatter!(
    [0,1,2,3,4,6,8,10,12,16,24], [0,25,100,132,90,82,60,55,32,15,4],
    label="A200 actual", markercolor=:red,
)

c_max = 125.66238991436627
halflife = 4.883455201973709
t_max = 2.7924975536499423

Đó thực sự là một sự phù hợp tốt đẹp.

Độ khớp bình phương nhỏ nhất này là ước tính khả năng tối đa (MLE), nó là đường cong tối đa hóa khả năng xảy ra của các tham số.
Nhưng thực sự chúng ta không theo đuổi một đường cong nào cả.
Chúng tôi quan tâm đến phân phối trên các đường cong có thể, dựa trên các quan sát.
Những điều này cho biết sử dụng những thực tế có thể xảy ra để giải thích những gì chúng ta đang thấy.

Để làm điều này, chúng tôi sẽ sử dụng lập trình xác suất.
Đặc biệt là thư viện Turing.jl.
Đây là một trường hợp sử dụng lý tưởng cho lập trình xác suất.
Chúng tôi có rất nhiều kiến ​​​​thức về miền để cung cấp cho các linh mục.
Chúng tôi chỉ có một số lượng rất nhỏ các quan sát.
và chúng tôi muốn thực hiện suy luận để xác định phân phối trên một số lượng nhỏ tham số.

Để bắt đầu, hãy nghĩ về những điều trước đây của chúng tôi. Đây là niềm tin của chúng tôi về các giá trị mà các tham số có thể nhận trước khi chúng tôi xem xét dữ liệu.

• c_max nằm trong khoảng từ 0 đến 500 pg/mL (tức là 0-1835 pmol/L).
  Nếu E2 của bạn cao hơn thì có gì đó không ổn. Bây giờ, chúng ta đừng giả sử bất cứ điều gì nữa và chỉ cần đi với một Uniform phân bổ. Mặc dù có lẽ chúng ta có thể làm điều gì đó thông minh hơn bằng cách chọn một thứ gì đó kết thúc tốt đẹp về phía cuối. ví dụ 500Beta(1.01,1.01).
• t_max là khoảng từ 1 đến 4 giờ, chúng tôi biết điều này vì hướng dẫn nói rằng đừng để bất kỳ ai chạm vào bạn trong giờ đầu tiên (vì vậy sau đó nó chắc chắn vẫn còn thấm nước) và thông thường là không rửa khu vực đó trong ít nhất 4 giờ – thế là phải làm nhưng rồi. Nếu chúng ta sử dụng một 4Beta(2,2)+1 phân phối nó có một số lực đẩy về phía trung tâm.
• halflife, chúng tôi biết điều này phải tích cực, vì nếu không thì nó sẽ không phân rã. Nhật ký bình thường có ý nghĩa vì nó xuất hiện theo cấp số nhân. Chúng tôi muốn nó có chế độ 3.5 vì đó là những gì bằng mắt thường mà chúng tôi thấy phù hợp với các đường cong một cách độc đáo (có lẽ là gian lận Bayesian tồi tệ ở đây) và bởi vì điều đó có nghĩa là nó hầu như bị phân hủy hoàn toàn sau 24 giờ – thông thường nó không thể cao hơn nhiều vì nếu không thì sẽ không cần liều lượng hàng ngày , cũng không thấp hơn nhiều vì trong trường hợp đó sẽ cần nhiều liều mỗi ngày. Để đặt chế độ thành 3.5, chúng tôi sử dụng LogNormal(log(3.5)+1, 1)

using Turing, StatsPlots

plot(layout=(1,3), legend=false, )
plot!(Uniform(0, 500), title="c_max", subplot=1, linewidth=2)
plot!(3Beta(2,2) + 1, title="t_max", subplot=2, linewidth=2)
plot!(LogNormal(log(3.5)+1, 1), title="halflife", xrange=(0,50), subplot=3, linewidth=2)

Thành phần khác mà chúng tôi sẽ muốn là một thuật ngữ lỗi.
Chúng tôi muốn thể hiện các quan sát của mình về nồng độ dưới dạng các mẫu nhiễu từ phân phối chuẩn có tâm trên đường cong thực tế mà chúng tôi đang ước tính.
Vì vậy, chúng tôi cần một thuật ngữ sai số sẽ cho phép một số chỗ ngọ nguậy về đường cong đó mà không loại bỏ suy luận về các tham số thực.
Chúng ta sẽ định nghĩa một biến gọi là err đó là độ lệch chuẩn của thuật ngữ lỗi này.
Tiên đoán của chúng ta về thuật ngữ lỗi này phải dương với đỉnh bằng 0 và nhanh chóng tắt dần.
Gamma(1, 1) đáp ứng yêu cầu của chúng tôi.

plot(Gamma(1,1), title="err", legend=false)

@model function single_dose_model()
    c_max ~ Uniform(0, 500)
    t_max ~ 3Beta(2,2) + 1
    halflife ~ LogNormal(log(3.5)+1, 1)
    
    dose_f = single_dose(c_max, halflife, t_max)  

    err ~ Gamma(1, 1)
    
    # There is probably a smarter way to do this, but we want to allow observations at different hours
    c1 ~ Normal(dose_f(1), err)
    c2 ~ Normal(dose_f(2), err)
    c3 ~ Normal(dose_f(3), err)
    c4 ~ Normal(dose_f(4), err)
    c6 ~ Normal(dose_f(6), err)
    c8 ~ Normal(dose_f(8), err)
    c10 ~ Normal(dose_f(10), err)
    c12 ~ Normal(dose_f(12), err)
    c16 ~ Normal(dose_f(16), err)
    c24 ~ Normal(dose_f(24), err)
    
    return (c1, c2, c3, c4, c6, c8, c10, c12, c16, c24)
end

model = single_dose_model() | (;c1=25,c2=100,c3=132,c4=90,c6=82,c8=60,c10=55,c12=32,c16=15,c24=4)
chain=sample(model, NUTS(), 4_000)

using StatsPlots
plot(chain)

Vì vậy, hãy xem phân phối trên các đường cong (như được biểu thị bằng các mẫu).

using Tables

gr(fmt=:png)  # make sure not to plot svg or will crash browser
scatter(
    [0,1,2,3,4,6,8,10,12,16,24], [0,25,100,132,90,82,60,55,32,15,4],
    label="A200 actual", markercolor=:red,
    legend=false
)

for samp in rowtable(chain)
    f = single_dose(samp.c_max, samp.halflife, samp.t_max)
    plot!(
        0:0.1:24, f,
        linewidth=0.5,
        linealpha=0.005, linecolor=:red, yrange=(0,200),
    )
end
plot!()

Chúng tôi thấy phạm vi giá trị khá nhỏ và rõ ràng này cho các tham số: c_max, t_max, halflife.
Thuật ngữ lỗi, errkhá lớn

Bây giờ chúng tôi đã chỉ ra rằng chúng tôi có thể suy luận để tìm các phân phối trên các tham số phù hợp với dữ liệu, hãy bắt tay vào một nhiệm vụ thực tế hơn.
Không ai được xét nghiệm máu vài giờ một lần ngoài bài tập thu thập dữ liệu thử nghiệm.
Các xét nghiệm máu thường xuyên nhất mà tôi được biết là 2 tuần một lần và hầu hết là 3-6 tháng một lần.
Vì vậy, những gì chúng tôi thực sự quan tâm là suy luận điều gì có thể xảy ra với nồng độ trong máu từ một lần quan sát.

Trường hợp thực tế là bôi gel vào buổi sáng trước khi đi làm và xét nghiệm máu vào cuối ngày làm việc ngay trước khi máy lấy máu đóng cửa.
Vì vậy, chúng tôi sẽ xem xét mẫu cho c8nồng độ trong máu 8 giờ sau khi dùng thuốc.
Và chúng tôi muốn thấy sự phân phối trên tất cả các chức năng mức độ estrogen có thể dẫn đến quan sát đó.
Đây là những gì nó thực sự là tất cả về.

model = single_dose_model() | (;c8=60)
chain=sample(model, NUTS(), 4_000)

scatter(
    [8], [60],
    label="A200 actual", markercolor=:red,
    legend=false
)

for samp in rowtable(chain)
    f = single_dose(samp.c_max, samp.halflife, samp.t_max)
    plot!(
        0:0.1:24, f,
        linewidth=1,
        linealpha=0.005, linecolor=:red, yrange=(0,200),
    )
end
plot!()

Vì vậy, đó là thực sự thực sự thông tin.
Có một loạt các giải thích có thể.
Từ rất nhỏ t_max và lớn c_max có nghĩa là nó đạt đỉnh sớm và giảm đi rất nhiều,
đến những cái có nhiều khả năng trông giống với loại đường cong mà chúng ta đã thấy dựa trên dữ liệu thực nghiệm với các phép đo thường xuyên hơn.
Điều này thực sự tuyệt vời.

Chúng ta có thể thêm nhiều điểm quan sát hơn và cắt giảm số lượng thực tế mà chúng ta có thể tham gia.
Điều này thực tế thực sự là một điều thiết thực để làm.
Mỗi điểm là một lần xét nghiệm máu – ở Luân Đôn, người chuyển giới có thể lấy máu miễn phí vào buổi tối vài ngày một tuần tại 56T Dean St hoặc CliniQ.
Về lý thuyết, một số bác sĩ đa khoa của NHS cũng sẽ làm chúng cho người chuyển giới, nhưng nhiều người (bao gồm cả tôi) sẽ không làm vì NHS kỵ khí một cách có hệ thống và các bác sĩ cho rằng người chuyển giới quá phức tạp đối với họ. Vì vậy, để có được một ít trong một ngày, bạn sẽ phải trả riêng khoảng 50 bảng Anh. Thật tuyệt khi chúng ta biết chính xác mỗi lần quan sát tốn bao nhiêu tiền.

Vì vậy, chúng ta có thể mô phỏng thêm một quan sát khác.
Bạn có thể thấy trong các biểu đồ sau rằng nếu chúng ta thêm số đọc là 60 vào 8 giờ sau khi ứng dụng, chúng ta sẽ chia các vũ trụ khả dĩ thành hai nhóm giải thích khả dĩ. Một bộ có số đọc 3 giờ khi nó vẫn đang tăng và một bộ khi nó đang giảm.

model = single_dose_model() | (;c3=100, c8=60)
chain=sample(model, NUTS(), 4_000)

scatter(
    [3, 8], [100, 60],
    label="A200 actual", markercolor=:red,
    legend=false
)

for samp in rowtable(chain)
    f = single_dose(samp.c_max, samp.halflife, samp.t_max)
    plot!(
        0:0.1:24, f,
        linewidth=1,
        linealpha=0.005, linecolor=:red, yrange=(0,200),
    )
end
plot!()

model = single_dose_model() | (;c3=100, c8=60, c1=50)
chain=sample(model, NUTS(), 4_000)

scatter(
    [3, 8, 1], [100, 60, 50],
    label="A200 actual", markercolor=:red,
    legend=false
)

for samp in rowtable(chain)
    f = single_dose(samp.c_max, samp.halflife, samp.t_max)
    plot!(
        0:0.1:24, f,
        linewidth=1,
        linealpha=0.005, linecolor=:red, yrange=(0,200),
    )
end
plot!()

Điều đáng chú ý là mặc dù 3 điểm này dẫn đến một tập hợp tương đối nhỏ các thế giới có thể xảy ra, nhưng điều này không phải lúc nào cũng đúng.
Không giống như để khớp một tuyến tính đa biến (hoặc một số đường cong khác), một điểm trên mỗi bậc tự do là không đủ.
Ví dụ: nếu tất cả các điểm có thể nằm trong phân đoạn theo từng phần đầu tiên, thì chúng ta thậm chí còn biết ít hơn về tình trạng của thế giới so với việc chúng ta có 1 điểm trong phân đoạn nối đuôi (đặc biệt nếu những điểm đó gợi ý một lỗi lớn trong kết quả đọc err)

model = single_dose_model() | (;c1=30, c2=90, c3=90)
chain=sample(model, NUTS(), 4_000)

scatter(
    [1,2,3], [30, 90, 90],
    label="A200 actual", markercolor=:blue,
    legend=false
)

for samp in rowtable(chain)
    f = single_dose(samp.c_max, samp.halflife, samp.t_max)
    plot!(
        0:0.1:24, f,
        linewidth=1,
        linealpha=0.005, linecolor=:blue, yrange=(0,200),
    )
end
plot!()

Đây chỉ là cái nhìn đầu tiên về chủ đề này.
Tôi tưởng tượng tôi có thể trở lại với nó một lần nữa trong tương lai.
Dưới đây là một số điều bổ sung mà chúng tôi có thể muốn xem xét:

• Xác định thời điểm tối ưu để kiểm tra: như đã thảo luận, 3 lần đọc sẽ không phải lúc nào cũng nắm bắt được đường cong, những thời điểm khác nhau có thể mang lại nhiều thông tin hơn những thời điểm khác, đặc biệt là khi chúng tôi xem xét mức lỗi so với mức tín hiệu.
• Mức trung bình: phân bố của mức trung bình có khả năng bị sụp đổ khá nhiều – nhiều bộ giá trị tham số khác nhau có thể dẫn đến cùng một mức trung bình.
• Nhiều ngày: Vì estrogen không đạt đến mức 0 sau 24 giờ nên có thể lập mô hình theo ngày. Cũng có thể bao gồm một thuật ngữ cho sự thay đổi khi nó được áp dụng trong ngày vì mọi người không nhất quán. Thứ sáu là rất quan trọng để làm cho mô hình thực tế
• Thay đổi liều lượng: kéo dài đến thứ sáu, mọi người thay đổi ở đó và chúng tôi biết liều lượng cao hơn dẫn đến mức cao hơn nên chúng tôi có thể thêm kiến ​​thức trước đó.

Lập trình xác suất là một kỹ thuật tuyệt vời để nghiên cứu về dược lực học.
Nó cho phép chúng ta xử lý thực tế là chúng ta có nhiều điều chưa biết về sinh học cá nhân của mọi người, trong khi vẫn thu hẹp phạm vi thế giới khả dĩ mà họ có thể sống.