Selection bias in clinical sciences: Difference between revisions
Selection bias in clinical sciences (view source)
Revision as of 08:52, 28 September 2022
, 1 year agojulia example
m (wording) |
(julia example) |
||
Line 54:
<br /><br /><br />
=={{header|Julia}}==
<syntaxhighlight lang="julia">using HypothesisTests
@enum TreatmentClass Untreated Irregular Regular
mutable struct Subject
cum_dose::Float64
treatment_class::TreatmentClass
had_covid::Bool
end
function update!(subjects::Vector{Subject}, pcovid = 0.001, pstart = 0.005, pdosing = 0.25, drange = 3:3:9)
for subj in subjects
if subj.had_covid
continue
elseif rand() < pcovid
subj.had_covid = true
elseif subj.cum_dose > 0 && rand() <= pdosing * 10 || subj.cum_dose == 0 && rand() <= pstart
subj.cum_dose += rand(drange)
subj.treatment_class =
subj.cum_dose == 0 ? Untreated : subj.cum_dose >= 100 ? Regular : Irregular
end
end
end
function run_study(N = 10_000, duration = 180)
population = [Subject(0.0, Untreated, false) for _ in 1:N]
unt, unt_covid, irr, irr_covid, reg, reg_covid = 0, 0, 0, 0, 0, 0
println("Population size $N, daily infection risk 0.1%")
for day in 1:duration
update!(population)
if day % 30 == 0
println("\nDay $day:")
unt = count(s -> s.treatment_class == Untreated, population)
unt_covid = count(s -> (s.treatment_class == Untreated) && s.had_covid, population)
println("Untreated: N = $unt, with infection = $unt_covid")
irr = count(s -> s.treatment_class == Irregular, population)
irr_covid = count(s -> (s.treatment_class == Irregular) && s.had_covid, population)
println("Irregular Use: N = $irr, with infection = $irr_covid")
reg = count(s -> s.treatment_class == Regular, population)
reg_covid = count(s -> (s.treatment_class == Regular) && s.had_covid, population)
println("Regular Use: N = $reg, with infection = $reg_covid")
end
if day == 90
println("\nAt midpoint, Infection case percentages are:")
println(" Untreated : ", Float16(100 * unt_covid / unt))
println(" Irregulars: ", Float16(100 * irr_covid / irr))
println(" Regulars : ", Float16(100 * reg_covid / reg))
end
end
println("\nAt study end, Infection case percentages are:")
println(" Untreated : ", Float16(100 * unt_covid / unt), " of group size of $unt")
println(" Irregulars: ", Float16(100 * irr_covid / irr), " of group size of $irr")
println(" Regulars : ", Float16(100 * reg_covid / reg), " of group size of $reg")
untreated = [s.had_covid for s in population if s.treatment_class == Untreated]
irregular = [s.had_covid for s in population if s.treatment_class == Irregular]
regular = [s.had_covid for s in population if s.treatment_class == Regular]
println("\n\n Final statistics:\n")
@show KruskalWallisTest(untreated, irregular, regular)
end
run_study()
</syntaxhighlight>{{out}}
<pre>
Population size 10000, daily infection risk 0.1%
Day 30:
Untreated: N = 8671, with infection = 263
Irregular Use: N = 694, with infection = 16
Regular Use: N = 635, with infection = 2
Day 60:
Untreated: N = 7529, with infection = 505
Irregular Use: N = 620, with infection = 43
Regular Use: N = 1851, with infection = 40
Day 90:
Untreated: N = 6583, with infection = 702
Irregular Use: N = 578, with infection = 60
Regular Use: N = 2839, with infection = 113
At midpoint, Infection case percentages are:
Untreated : 10.664
Irregulars: 10.38
Regulars : 3.98
Day 120:
Untreated: N = 5763, with infection = 862
Irregular Use: N = 526, with infection = 78
Regular Use: N = 3711, with infection = 208
Day 150:
Untreated: N = 5106, with infection = 986
Irregular Use: N = 420, with infection = 91
Regular Use: N = 4474, with infection = 327
Day 180:
Untreated: N = 4552, with infection = 1099
Irregular Use: N = 369, with infection = 97
Regular Use: N = 5079, with infection = 472
At study end, Infection case percentages are:
Untreated : 24.14 of group size of 4552
Irregulars: 26.28 of group size of 369
Regulars : 9.3 of group size of 5079
Final statistics:
KruskalWallisTest(untreated, irregular, regular) = Kruskal-Wallis rank sum test (chi-square approximation)
-------------------------------------------------------
Population details:
parameter of interest: Location parameters
value under h_0: "all equal"
point estimate: NaN
Test summary:
outcome with 95% confidence: reject h_0
one-sided p-value: <1e-88
Details:
number of observation in each group: [4552, 369, 5079]
χ²-statistic: 406.313
rank sums: [2.44609e7, 2.02244e6, 2.35217e7]
degrees of freedom: 2
adjustment for ties: 0.416933
Kruskal-Wallis rank sum test (chi-square approximation)
-------------------------------------------------------
Population details:
parameter of interest: Location parameters
value under h_0: "all equal"
point estimate: NaN
Test summary:
outcome with 95% confidence: reject h_0
one-sided p-value: <1e-88
Details:
number of observation in each group: [4552, 369, 5079]
χ²-statistic: 406.313
rank sums: [2.44609e7, 2.02244e6, 2.35217e7]
degrees of freedom: 2
adjustment for ties: 0.416933
</pre>
=={{header|Python}}==
| |||