Rosetta Code

Selection bias in clinical sciences: Difference between revisions

Page Discussion

Selection bias in clinical sciences (view source)

Revision as of 04:22, 7 May 2024

57,875 bytes added ,  28 days ago
m
→‎{{header|jq}}: typo
m (sp)
 
(17 intermediate revisions by 9 users not shown)
Line 1:
{{draft task|Probability and statistics}}
 
In epidemiology, retrospective analyses have well-known limitations compared to prospective studies.
Line 36:
* Use at least 1000 subjects in the simulation over the 180 days (historically, the study size was 80,000).
 
* Statistics used are to be the KruscalKruskal statistic for the analysis of multiple groups, with the boolean study outcome variable whether the subject got Covid-19 during the study period, analyzed versus category. If the programming language does not have a statistical package for such a test, calculation of the Kruskal H statistic may be used, since any values of H over about 50 are considered very highly indicative of significant group difference. A description of that test is found at [[wp:Kruskal%E2%80%93Wallis_one-way_analysis_of_variance|Kruskal-Wallis test]].
 
* You should get a statistical result highly favoring the REGULAR group.
Line 48:
 
;See also:
:;*[[httpswp://en.wikipedia.org/wiki/Selection_bias|Wikipedia - Selection Bias]]
:;*[[httpswp://en.wikipedia.org/wiki/Retrospective_cohort_study |Wikipedia entry on- retrospectiveRetrospective studies]]
:;*[[https://health.ucdavis.edu/ctsc/area/Resource_Library/documents/Challenges-of-Retrospective-Observational-Studies_8March2017_Kim.pdf UCDavis (pdf) Challenges in Retrospective studies]]
<br /><br /><br />
 
 
=={{header|FreeBASIC}}==
{{trans|Wren}}
<syntaxhighlight lang="vbnet">Const UNTREATED = 0
Const IRREGULAR = 1
Const REGULAR = 2
Const DOSE_FOR_REGULAR = 100
 
' A subject for the study
Type Subject
cumDose As Double
category As Integer
hadCovid As Boolean
updateCount As Integer
End Type
 
'Set treatment category based on cumulative treatment taken.
Sub categorize (Byref s As Subject)
If s.cumDose = 0 Then
s.category = UNTREATED
Elseif s.cumDose >= DOSE_FOR_REGULAR Then
s.category = REGULAR
Else
s.category = IRREGULAR
End If
End Sub
 
'Daily update on the subject to check for infection and randomly dose.
Sub update (Byref s As Subject, pCovid As Double, pStartingTreatment As Double, pRedose As Double, dRange As Double)
If Not s.hadCovid Then
If Rnd < pCovid Then
s.hadCovid = True
Elseif (s.cumDose = 0 And Rnd < pStartingTreatment) Or (s.cumDose > 0 And Rnd < pRedose) Then
s.cumDose += Int(Rnd * dRange)
categorize(s)
End If
End If
s.updateCount += 1
End Sub
 
Function sumRanks(a() As Integer) As Integer
Dim As Integer sum = 0
For i As Integer = Lbound(a) To Ubound(a)
sum += a(i)
Next i
Return sum
End Function
 
Function kruskal(a() As Integer, b() As Integer, c() As Integer) As Double
Dim As Integer n = Ubound(a) + Ubound(b) + Ubound(c) + 3
Dim As Double sra = sumRanks(a())
Dim As Double srb = sumRanks(b())
Dim As Double src = sumRanks(c())
Dim As Double H = 12/(n*(n+1)) * (sra*sra/Ubound(a) + srb*srb/Ubound(b) + src*src/Ubound(c)) - 3*(n+1)
Return H
End Function
 
'Run the study using the population of size 'N' for 'duration' days
Sub runStudy(numSubjects As Integer, duration As Integer, interval As Integer)
Dim As Integer i
Dim As Subject population(numSubjects - 1)
For i = 0 To numSubjects - 1
population(i).cumDose = 0
population(i).category = UNTREATED
population(i).hadCovid = False
population(i).updateCount = 0
Next i
Dim As Integer unt, untCovid, irr, irrCovid, reg, regCovid, dia
Print "Total subjects:"; numSubjects
For dia = 0 To duration - 1
For i = 0 To numSubjects - 1
update(population(i), 0.001, 0.005, 0.25, Int(Rnd * 9) + 3)
Next i
If (dia + 1) Mod interval = 0 Then
Print !"\nDay"; dia + 1
unt = 0
untCovid = 0
irr = 0
irrCovid = 0
reg = 0
regCovid = 0
For i = 0 To numSubjects - 1
If population(i).category = UNTREATED Then
unt += 1
If population(i).hadCovid Then untCovid += 1
Elseif population(i).category = IRREGULAR Then
irr += 1
If population(i).hadCovid Then irrCovid += 1
Elseif population(i).category = REGULAR Then
reg += 1
If population(i).hadCovid Then regCovid += 1
End If
Next i
Print " Untreated: N ="; unt; ", with infection ="; untCovid
Print "Irregular Use: N ="; irr; ", with infection ="; irrCovid
Print " Regular Use: N ="; reg; ", with infection ="; regCovid
End If
 
If dia = duration \ 2 - 1 Then
Print !"\nAt midpoint, Infection case percentages are:"
Print " Untreated :"; 100 * untCovid / unt
Print " Irregulars:"; 100 * irrCovid / irr
Print " Regulars :"; 100 * regCovid / reg
End If
Next dia
Print !"\nAt study end, Infection case percentages are:"
Print " Untreated :"; 100 * untCovid / unt; " of group size of"; unt
Print " Irregulars:"; 100 * irrCovid / irr; " of group size of"; irr
Print " Regulars :"; 100 * regCovid / reg; " of group size of"; reg
Dim As Integer sinTratar(numSubjects - 1)
Dim As Integer esIrregular(numSubjects - 1)
Dim As Integer esRegular(numSubjects - 1)
For i = 0 To numSubjects - 1
If population(i).category = UNTREATED Then sinTratar(i) = population(i).hadCovid
If population(i).category = IRREGULAR Then esIrregular(i) = population(i).hadCovid
If population(i).category = REGULAR Then esRegular(i) = population(i).hadCovid
Next i
Print !"\nFinal statistics: H = "; kruskal(sinTratar(), esIrregular(), esRegular())
End Sub
 
Randomize Timer
runStudy(1000, 180, 30)
 
Sleep</syntaxhighlight>
{{out}}
<pre>Total subjects: 1000
 
Day 30
Untreated: N = 908, with infection = 41
Irregular Use: N = 92, with infection = 1
Regular Use: N = 0, with infection = 0
 
Day 60
Untreated: N = 794, with infection = 73
Irregular Use: N = 206, with infection = 2
Regular Use: N = 0, with infection = 0
 
Day 90
Untreated: N = 718, with infection = 92
Irregular Use: N = 281, with infection = 6
Regular Use: N = 1, with infection = 0
 
At midpoint, Infection case percentages are:
Untreated : 12.81337047353761
Irregulars: 2.135231316725978
Regulars : 0
 
Day 120
Untreated: N = 640, with infection = 101
Irregular Use: N = 347, with infection = 11
Regular Use: N = 13, with infection = 0
 
Day 150
Untreated: N = 589, with infection = 116
Irregular Use: N = 349, with infection = 27
Regular Use: N = 62, with infection = 0
 
Day 180
Untreated: N = 528, with infection = 131
Irregular Use: N = 331, with infection = 36
Regular Use: N = 141, with infection = 4
 
At study end, Infection case percentages are:
Untreated : 24.81060606060606 of group size of 528
Irregulars: 10.8761329305136 of group size of 331
Regulars : 2.836879432624114 of group size of 141
 
Final statistics: H = -9002.999975352894</pre>
 
=={{header|jq}}==
'''Adapted from [[#Wren|Wren]]'''
 
'''Works with jq, the C implementation of jq'''
 
'''Works with gojq, the Go implementation of jq'''
 
Since jq does not include a PRN generator, an external source of entropy such as /dev/urandom is assumed.
A suitable invocation of jq would be along the lines of:
<pre>
< /dev/urandom tr -cd '0-9' | fold -w 1 | $jq -Rcnr --rawfile text alice_oz.txt -f mc.jq
</pre>
<syntaxhighlight lang="jq">
### Generic utilities
def lpad($len): tostring | ($len - length) as $l | (" " * $l) + .;
 
def sum(s): reduce s as $x (0; . + $x);
 
# For readability:
def count(s): reduce s as $x (0; . + 1);
 
### Kruskal-Wallis statistic
 
# The arguments are assumed to be arrays of boolean values.
def kruskal($a; $b; $c):
def tobits: map(if . then 1 else 0 end);
# map the boolean values to 1 (true) or 0 (false).
($a|tobits) as $aa
| ($b|tobits) as $bb
| ($c|tobits) as $cc
# aggregate and sort them
| (($aa + $bb + $cc)|sort) as $ss
# find rank of first occurrence of 1
| (($ss|index(1)) + 1) as $ix
# calculate average ranks for 0 and 1
| ($ix / 2) as $arf
| ($ss|length) as $n
| (($ix + $n) / 2) as $art
# calculate sum of ranks for each list
| if any($a,$b,$c; length == 0) then null
else
sum($a[] | if . then $art else $arf end) as $sra
| sum($b[] | if . then $art else $arf end) as $srb
| sum($c[] | if . then $art else $arf end) as $src
# calculate H
| (12/($n*($n+1))) * ($sra*$sra/($a|length) + $srb*$srb/($b|length) + $src*$src/($c|length)) - 3 * ($n + 1)
end;
 
### Pseuo-random numbers
 
# Output: a prn in range(0;$n) where $n is `.`
def prn:
if . == 1 then 0
else . as $n
| ([1, (($n-1)|tostring|length)]|max) as $w
| [limit($w; inputs)] | join("") | tonumber
| if . < $n then . else ($n | prn) end
end;
 
def sample:
if length == 0 # e.g. null or []
then null
else .[length|prn]
end;
 
def randFloat:
(1000|prn) / 1000;
 
### Simulation
 
def UNTREATED: 0;
def IRREGULAR: 1;
def REGULAR: 2;
def DOSE_FOR_REGULAR: 100;
 
def Subject:
{cumDose: 0,
category: UNTREATED,
hadCovid: false,
updateCount: 0 };
 
def Subject($cumDose; $category; $hadCovid; $updateCount):
{$cumDose, $category, $hadCovid, $updateCount};
 
# Set treatment category based on cumulative treatment taken.
def categorize:
if .cumDose == 0 then UNTREATED
elif .cumDose >= DOSE_FOR_REGULAR then REGULAR
else IRREGULAR
end;
 
# Daily update on the input Subject to check for infection and randomly dose.
def update($pCovid; $pStartingTreatment; $pRedose; $dRange):
if (.hadCovid|not)
then if (randFloat < $pCovid)
then .hadCovid = true
else if (.cumDose == 0 and randFloat < $pStartingTreatment) or
(.cumDose > 0 and randFloat < $pRedose)
then .cumDose += ($dRange|sample)
| .category = categorize
end
end
| .updateCount += 1
end;
 
# Update using default parameters.
def update: update(0.001; 0.005; 0.25; [3, 6, 9]);
 
# Run the study using the population of size $numSubjects for $duration days.
def runStudy($numSubjects; $duration; $interval):
def r: . * 1000 | round / 1000;
def div($i; $j): if $j == 0 then 0 else $i / $j | r | lpad(6) end;
{ population: [range(0; $numSubjects) | Subject],
unt: 0,
untCovid: 0,
irr: 0,
irrCovid: 0,
reg: 0,
regCovid: 0
}
| "Total subjects: \($numSubjects)",
(foreach range(0; $duration) as $day (.;
reduce range(0; $numSubjects) as $i (.; .population[$i] |= update) ;
if ( ($day + 1) % $interval == 0 ) or $day == (($duration/2)|floor) - 1 or $day == $duration-1
then
.unt = count(.population[] | select(.category == UNTREATED))
| .untCovid = count(.population[] | select( .category == UNTREATED and .hadCovid ))
| .irr = count(.population[] | select(.category == IRREGULAR ))
| .irrCovid = count(.population[] | select(.category == IRREGULAR and .hadCovid))
| .reg = count(.population[] | select (.category == REGULAR ))
| .regCovid = count(.population[] | select(.category == REGULAR and .hadCovid))
| (select( ($day + 1) % $interval == 0 )
| "Day \($day + 1):",
"Untreated: N = \(.unt), with infection = \(.untCovid)",
"Regular Use: N = \(.reg), with infection = \(.regCovid)" ),
 
(select($day == (($duration/2)|floor) - 1)
| "\nAt midpoint, Infection case percentages are:",
" Untreated : \(div(100 * .untCovid ; .unt))",
" Irregulars: \(div(100 * .irrCovid ; .irr))",
" Regulars : \(div(100 * .regCovid ; .reg))" ),
(select($day == $duration-1)
| "\nAt study end, infection case percentages are:",
" Untreated : \(div(100 * .untCovid ; .unt)) of group size \(.unt)",
" Irregulars: \(div(100 * .irrCovid ; .irr)) of group size \(.irr)",
" Regulars : \(div(100 * .regCovid ; .reg)) of group size \(.reg)",
( (.population | map(select(.category == UNTREATED)) | map(.hadCovid)) as $untreated
| (.population | map(select(.category == IRREGULAR)) | map(.hadCovid)) as $irregular
| (.population | map(select(.category == REGULAR)) | map(.hadCovid)) as $regular
| "\nFinal statistics: H = \(kruskal($untreated; $irregular; $regular)|r)" ) )
else empty
end ) );
 
runStudy(10000; 180; 60)
</syntaxhighlight>
{{output}}
<pre>
Total subjects: 10000
 
Day 60:
Untreated: N = 7471, with infection = 508
Regular Use: N = 185, with infection = 0
 
At midpoint, Infection case percentages are:
Untreated : 10.74
Irregulars: 6.2
Regulars : 1.43
 
Day 120:
Untreated: N = 5780, with infection = 880
Regular Use: N = 2139, with infection = 75
 
Day 180:
Untreated: N = 4583, with infection = 1127
Regular Use: N = 3823, with infection = 249
 
At study end, infection case percentages are:
Untreated : 24.591 of group size 4583
Irregulars: 18.068 of group size 1594
Regulars : 6.513 of group size 3823
 
Final statistics: H = 205.487
</pre>
 
=={{header|Julia}}==
Line 198 ⟶ 557:
degrees of freedom: 2
adjustment for ties: 0.417533
</pre>
 
=={{header|Nim}}==
{{trans|Python}}
{{trans|Wren}}
<syntaxhighlight lang="Nim">import std/[algorithm, math, random, sequtils, strformat, sugar]
 
type
 
TreatmentClass {.pure.} = enum Untreated, Irregular, Regular
 
Subject = object
cumDose: float
category: TreatmentClass
hadCovid: bool
 
const DoseForRegular = 100
 
func categorize(subject: var Subject) =
## Set treatment category based on cumulative treatment taken.
subject.category = if subject.cumDose == 0: Untreated
elif subject.cumDose >= DoseForRegular: Regular
else: Irregular
 
proc update(subject: var Subject; pCovid = 0.001;
pStart = 0.005; pDosing = 0.25; doses = @[3, 6, 9]) =
## Daily update on the subject to check for infection and randomly dose.
if not subject.hadCovid:
if rand(1.0) < pCovid:
subject.hadCovid = true
elif (subject.cumDose == 0 and rand(1.0) < pStart) or
(subject.cumDose > 0 and rand(1.0) < pDosing):
subject.cumDose += sample(doses).toFloat
subject.categorize()
 
func kruskalWallis(a, b, c: seq[bool]): float =
 
# Aggregate and sort.
let s = sorted(a & b & c)
# Find rank of first occurrence of "true".
let ix = s.find(true) + 1
# Calculate average ranks for "false" and "true".
let n = s.len
let arf = ix / 2
let art = (ix + n) / 2
# Calculate sum of ranks for each list.
let sra = sum(a.mapIt(if it: art else: arf))
let srb = sum(b.mapIt(if it: art else: arf))
let src = sum(c.mapIt(if it: art else: arf))
# Calculate H.
result = 12 / (n * (n + 1)) * (sra * sra / a.len.toFloat + srb * srb / b.len.toFloat +
src * src / c.len.toFloat) - 3 * float(n + 1)
 
proc runStudy(numSubjects = 1000; duration = 180; interval= 30) =
## Run the study using the population of size "numSubjects" for "duration" days.
var population = newSeqWith(numSubjects, Subject())
var unt, untCovid, irr, irrCovid, reg, regCovid = 0
echo &"Total subjects: {num_subjects}"
 
for day in 1..duration:
for subj in population.mitems:
subj.update()
 
if day mod interval == 0:
echo &"\nDay {day}:"
unt = population.countIt(it.category == Untreated)
untCovid = population.countIt(it.category == Untreated and it.hadCovid)
echo &"Untreated: N = {unt}, with infection = {untCovid}"
irr = population.countIt(it.category == IRREGULAR)
irrCovid = population.countIt(it.category == IRREGULAR and it.hadCovid)
echo &"Irregular Use: N = {irr}, with infection = {irrCovid}"
reg = population.countIt(it.category == REGULAR)
regCovid = population.countIt(it.category == REGULAR and it.hadCovid)
echo &"Regular Use: N = {reg}, with infection = {reg_covid}"
 
if day == duration div 2:
echo "\nAt midpoint, infection case percentages are:"
echo " Untreated : ", 100 * untCovid / unt
echo " Irregulars: ", 100 * irrCovid / irr
echo " Regulars : ", 100 * regCovid / reg
 
echo "\nAt study end, infection case percentages are:"
echo &" Untreated : {100 * untCovid / unt} of group size of {unt}"
echo &" Irregulars: {100 * irrCovid / irr} of group size of {irr}"
echo &" Regulars : {100 * regCovid / reg} of group size of {reg}"
let untreated = collect:
for s in population:
if s.category == Untreated:
s.hadCovid
let irregular = collect:
for s in population:
if s.category == Irregular:
s.hadCovid
let regular = collect:
for s in population:
if s.category == Regular:
s.hadCovid
echo "\nFinal statistics: ", kruskalWallis(untreated, irregular, regular)
 
randomize()
runStudy(10_000)
</syntaxhighlight>
 
{{out}}
<pre>Total subjects: 10000
 
Day 30:
Untreated: N = 8619, with infection = 283
Irregular Use: N = 1379, with infection = 27
Regular Use: N = 2, with infection = 0
 
Day 60:
Untreated: N = 7480, with infection = 498
Irregular Use: N = 2362, with infection = 77
Regular Use: N = 158, with infection = 3
 
Day 90:
Untreated: N = 6479, with infection = 712
Irregular Use: N = 2440, with infection = 151
Regular Use: N = 1081, with infection = 24
 
At midpoint, infection case percentages are:
Untreated : 10.98935020836549
Irregulars: 6.188524590163935
Regulars : 2.220166512488436
 
Day 120:
Untreated: N = 5655, with infection = 873
Irregular Use: N = 2199, with infection = 210
Regular Use: N = 2146, with infection = 70
 
Day 150:
Untreated: N = 4998, with infection = 1007
Irregular Use: N = 1910, with infection = 257
Regular Use: N = 3092, with infection = 135
 
Day 180:
Untreated: N = 4479, with infection = 1122
Irregular Use: N = 1575, with infection = 302
Regular Use: N = 3946, with infection = 232
 
At study end, infection case percentages are:
Untreated : 25.05023442732753 of group size of 4479
Irregulars: 19.17460317460317 of group size of 1575
Regulars : 5.879371515458693 of group size of 3946
 
Final statistics: 235.1089104422972
</pre>
 
=={{header|Perl}}==
{{trans|Raku}}
<syntaxhighlight lang="perl" line>use v5.036;
use List::AllUtils <sum0 indexes firstidx>;
use enum qw<False True UNTREATED REGULAR IRREGULAR>;
use constant DOSE_FOR_REGULAR => 100;
 
my ($nSubjects,$duration,$interval) = (10000, 180, 30);
my (@dosage) = (0) x $nSubjects;
my (@category) = (UNTREATED) x $nSubjects;
my (@hadcovid) = (False) x $nSubjects;
 
sub update {
my $pCovid = 1e-3;
my $pStartTreatment = 5e-3;
my $pRedose = 1/4;
my @dRange = <3 6 9>;
for my $i (0 .. $nSubjects-1) {
unless ($hadcovid[$i]) {
if (rand() < $pCovid) {
$hadcovid[$i] = True
} else {
my $dose = $dosage[$i];
if ($dose==0 && rand() < $pStartTreatment or $dose > 0 && rand() < $pRedose) {
$dosage[$i] = $dose += $dRange[3*rand()];
$category[$i] = ($dose > DOSE_FOR_REGULAR) ? REGULAR : IRREGULAR
}
}
}
}
}
 
sub kruskal (@sets) {
my @ranked = sort @{$sets[0]}, @{$sets[1]}, @{$sets[2]};
my $n = @ranked;
my @sr = (0) x @sets;
my $ix = firstidx { $_ == 1 } @ranked;
my ($arf,$art) = ($ix/2, ($ix+$n)/2);
 
for my $i (0..2) {
$sr[$i] += $_ ? $art : $arf for @{$sets[$i]}
}
 
my $H = sum0 map { my $s = $sr[$_]; $s**2 / @{$sets[$_]} } 0..$#sr;
12/($n*($n+1)) * $H - 3 * ($n + 1)
}
 
my($unt,$irr,$reg,$hunt,$hirr,$hreg,@sunt,@sirr,@sreg);
my $midpoint = int $duration / 2;
 
say "Total subjects: $nSubjects\n";
 
for my $day (1 .. $duration) {
update();
if (0 == $day % $interval or $day == $duration or $day == $midpoint) {
@sunt = @hadcovid[ indexes { $_ == UNTREATED } @category];
@sirr = @hadcovid[ indexes { $_ == IRREGULAR } @category];
@sreg = @hadcovid[ indexes { $_ == REGULAR } @category];
( $unt, $irr, $reg) = (scalar(@sunt), scalar(@sirr), scalar(@sreg));
($hunt,$hirr,$hreg) = ( sum0(@sunt), sum0(@sirr), sum0(@sreg));
}
 
if (0 == $day % $interval) {
printf "Day %d:\n", $day;
printf "Untreated: N = %4d, with infection = %4d\n", $unt,$hunt;
printf "Irregular Use: N = %4d, with infection = %4d\n", $irr,$hirr;
printf "Regular Use: N = %4d, with infection = %4d\n\n",$reg,$hreg;
}
 
if ($day == $midpoint or $day == $duration) {
my $stage = $day == $midpoint ? 'midpoint' : 'study end';
printf "At $stage, Infection case percentages are:\n";
printf " Untreated : %6.2f\n", 100*$hunt/$unt;
printf " Irregulars: %6.2f\n", 100*$hirr/$irr;
printf " Regulars : %6.2f\n\n", 100*$hreg/$reg;
}
}
 
printf "Final statistics: H = %.2f\n", kruskal ( \@sunt, \@sirr, \@sreg );
 
</syntaxhighlight>
{{out}}
<pre>Total subjects: 10000
 
Day 30:
Untreated: N = 8640, with infection = 283
Irregular Use: N = 1360, with infection = 18
Regular Use: N = 0, with infection = 0
 
Day 60:
Untreated: N = 7499, with infection = 532
Irregular Use: N = 2338, with infection = 75
Regular Use: N = 163, with infection = 1
 
Day 90:
Untreated: N = 6537, with infection = 709
Irregular Use: N = 2382, with infection = 144
Regular Use: N = 1081, with infection = 12
 
At midpoint, Infection case percentages are:
Untreated : 10.85
Irregulars: 6.05
Regulars : 1.11
 
Day 120:
Untreated: N = 5739, with infection = 855
Irregular Use: N = 2102, with infection = 216
Regular Use: N = 2159, with infection = 52
 
Day 150:
Untreated: N = 5091, with infection = 1012
Irregular Use: N = 1792, with infection = 270
Regular Use: N = 3117, with infection = 134
 
Day 180:
Untreated: N = 4468, with infection = 1126
Irregular Use: N = 1669, with infection = 315
Regular Use: N = 3863, with infection = 243
 
At study end, Infection case percentages are:
Untreated : 25.20
Irregulars: 18.87
Regulars : 6.29
 
Final statistics: H = 218.74</pre>
 
=={{header|Phix}}==
{{trans|Wren}}
<!--<syntaxhighlight lang="phix">(phixonline)-->
<span style="color: #008080;">with</span> <span style="color: #008080;">javascript_semantics</span>
<span style="color: #004080;">sequence</span> <span style="color: #000000;">dosage</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">category</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">hadcovid</span>
<span style="color: #008080;">enum</span> <span style="color: #000000;">UNTREATED</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">REGULAR</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">IRREGULAR</span>
<span style="color: #008080;">constant</span> <span style="color: #000000;">DOSE_FOR_REGULAR</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">100</span>
<span style="color: #008080;">procedure</span> <span style="color: #000000;">update</span><span style="color: #0000FF;">(</span><span style="color: #004080;">atom</span> <span style="color: #000000;">pCovid</span><span style="color: #0000FF;">=</span><span style="color: #000000;">0.001</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">pStartTreatment</span><span style="color: #0000FF;">=</span><span style="color: #000000;">0.005</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">pRedose</span><span style="color: #0000FF;">=</span><span style="color: #000000;">0.25</span><span style="color: #0000FF;">,</span> <span style="color: #004080;">sequence</span> <span style="color: #000000;">dRange</span><span style="color: #0000FF;">={</span><span style="color: #000000;">3</span><span style="color: #0000FF;">,</span><span style="color: #000000;">6</span><span style="color: #0000FF;">,</span><span style="color: #000000;">9</span><span style="color: #0000FF;">})</span>
<span style="color: #008080;">for</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">dosage</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">do</span>
<span style="color: #008080;">if</span> <span style="color: #008080;">not</span> <span style="color: #000000;">hadcovid</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">]</span> <span style="color: #008080;">then</span>
<span style="color: #008080;">if</span> <span style="color: #7060A8;">rnd</span><span style="color: #0000FF;">()<</span><span style="color: #000000;">pCovid</span> <span style="color: #008080;">then</span>
<span style="color: #000000;">hadcovid</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">]</span> <span style="color: #0000FF;">=</span> <span style="color: #004600;">true</span>
<span style="color: #008080;">else</span>
<span style="color: #004080;">atom</span> <span style="color: #000000;">dose</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">dosage</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">]</span>
<span style="color: #008080;">if</span> <span style="color: #0000FF;">(</span><span style="color: #000000;">dose</span><span style="color: #0000FF;">=</span><span style="color: #000000;">0</span> <span style="color: #008080;">and</span> <span style="color: #7060A8;">rnd</span><span style="color: #0000FF;">()<</span><span style="color: #000000;">pStartTreatment</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">or</span> <span style="color: #0000FF;">(</span><span style="color: #000000;">dose</span><span style="color: #0000FF;">></span><span style="color: #000000;">0</span> <span style="color: #008080;">and</span> <span style="color: #7060A8;">rnd</span><span style="color: #0000FF;">()<</span><span style="color: #000000;">pRedose</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">then</span>
<span style="color: #000000;">dose</span> <span style="color: #0000FF;">+=</span> <span style="color: #000000;">dRange</span><span style="color: #0000FF;">[</span><span style="color: #7060A8;">rand</span><span style="color: #0000FF;">(</span><span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">dRange</span><span style="color: #0000FF;">))]</span>
<span style="color: #000000;">dosage</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">]</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">dose</span>
<span style="color: #000000;">category</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">]</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">iff</span><span style="color: #0000FF;">(</span><span style="color: #000000;">dose</span><span style="color: #0000FF;">></span><span style="color: #000000;">DOSE_FOR_REGULAR</span><span style="color: #0000FF;">?</span><span style="color: #000000;">REGULAR</span><span style="color: #0000FF;">:</span><span style="color: #000000;">IRREGULAR</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">procedure</span>
<span style="color: #008080;">function</span> <span style="color: #000000;">kruskal</span><span style="color: #0000FF;">(</span><span style="color: #004080;">sequence</span> <span style="color: #000000;">sets</span><span style="color: #0000FF;">)</span>
<span style="color: #004080;">sequence</span> <span style="color: #000000;">ranked</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">sort</span><span style="color: #0000FF;">(</span><span style="color: #7060A8;">flatten</span><span style="color: #0000FF;">(</span><span style="color: #000000;">sets</span><span style="color: #0000FF;">)),</span>
<span style="color: #000000;">sr</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">repeat</span><span style="color: #0000FF;">(</span><span style="color: #000000;">0</span><span style="color: #0000FF;">,</span><span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">sets</span><span style="color: #0000FF;">))</span>
<span style="color: #004080;">integer</span> <span style="color: #000000;">ix</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">find</span><span style="color: #0000FF;">(</span><span style="color: #004600;">true</span><span style="color: #0000FF;">,</span><span style="color: #000000;">ranked</span><span style="color: #0000FF;">),</span>
<span style="color: #000000;">n</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">ranked</span><span style="color: #0000FF;">)</span>
<span style="color: #004080;">atom</span> <span style="color: #000000;">arf</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">ix</span><span style="color: #0000FF;">/</span><span style="color: #000000;">2</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">art</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">(</span><span style="color: #000000;">ix</span><span style="color: #0000FF;">+</span><span style="color: #000000;">n</span><span style="color: #0000FF;">)/</span><span style="color: #000000;">2</span>
<span style="color: #008080;">for</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">,</span><span style="color: #000000;">set</span> <span style="color: #008080;">in</span> <span style="color: #000000;">sets</span> <span style="color: #008080;">do</span>
<span style="color: #008080;">for</span> <span style="color: #000000;">b</span> <span style="color: #008080;">in</span> <span style="color: #000000;">set</span> <span style="color: #008080;">do</span>
<span style="color: #000000;">sr</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">]</span> <span style="color: #0000FF;">+=</span> <span style="color: #7060A8;">iff</span><span style="color: #0000FF;">(</span><span style="color: #000000;">b</span><span style="color: #0000FF;">?</span><span style="color: #000000;">art</span><span style="color: #0000FF;">:</span><span style="color: #000000;">arf</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
<span style="color: #004080;">atom</span> <span style="color: #000000;">H</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">0</span>
<span style="color: #008080;">for</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">,</span><span style="color: #000000;">s</span> <span style="color: #008080;">in</span> <span style="color: #000000;">sr</span> <span style="color: #008080;">do</span>
<span style="color: #000000;">H</span> <span style="color: #0000FF;">+=</span> <span style="color: #000000;">s</span><span style="color: #0000FF;">*</span><span style="color: #000000;">s</span><span style="color: #0000FF;">/</span><span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">sets</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">])</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
<span style="color: #000000;">H</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">12</span><span style="color: #0000FF;">/(</span><span style="color: #000000;">n</span><span style="color: #0000FF;">*(</span><span style="color: #000000;">n</span><span style="color: #0000FF;">+</span><span style="color: #000000;">1</span><span style="color: #0000FF;">))</span> <span style="color: #0000FF;">*</span> <span style="color: #000000;">H</span> <span style="color: #0000FF;">-</span> <span style="color: #000000;">3</span> <span style="color: #0000FF;">*</span> <span style="color: #0000FF;">(</span><span style="color: #000000;">n</span> <span style="color: #0000FF;">+</span> <span style="color: #000000;">1</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">return</span> <span style="color: #000000;">H</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
<span style="color: #008080;">procedure</span> <span style="color: #000000;">run_study</span><span style="color: #0000FF;">(</span><span style="color: #004080;">integer</span> <span style="color: #000000;">nSubjects</span><span style="color: #0000FF;">=</span><span style="color: #000000;">10000</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">duration</span><span style="color: #0000FF;">=</span><span style="color: #000000;">180</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">interval</span><span style="color: #0000FF;">=</span><span style="color: #000000;">30</span><span style="color: #0000FF;">)</span>
<span style="color: #000000;">dosage</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">repeat</span><span style="color: #0000FF;">(</span><span style="color: #000000;">0</span><span style="color: #0000FF;">,</span><span style="color: #000000;">nSubjects</span><span style="color: #0000FF;">)</span>
<span style="color: #000000;">category</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">repeat</span><span style="color: #0000FF;">(</span><span style="color: #000000;">UNTREATED</span><span style="color: #0000FF;">,</span><span style="color: #000000;">nSubjects</span><span style="color: #0000FF;">)</span>
<span style="color: #000000;">hadcovid</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">repeat</span><span style="color: #0000FF;">(</span><span style="color: #004600;">false</span><span style="color: #0000FF;">,</span><span style="color: #000000;">nSubjects</span><span style="color: #0000FF;">)</span>
<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"Total subjects: %d\n\n"</span><span style="color: #0000FF;">,</span><span style="color: #000000;">nSubjects</span><span style="color: #0000FF;">)</span>
<span style="color: #004080;">sequence</span> <span style="color: #000000;">sunt</span><span style="color: #0000FF;">,</span><span style="color: #000000;">sirr</span><span style="color: #0000FF;">,</span><span style="color: #000000;">sreg</span>
<span style="color: #004080;">integer</span> <span style="color: #000000;">unt</span><span style="color: #0000FF;">,</span><span style="color: #000000;">irr</span><span style="color: #0000FF;">,</span><span style="color: #000000;">reg</span><span style="color: #0000FF;">,</span><span style="color: #000000;">hunt</span><span style="color: #0000FF;">,</span><span style="color: #000000;">hirr</span><span style="color: #0000FF;">,</span><span style="color: #000000;">hreg</span><span style="color: #0000FF;">,</span><span style="color: #000000;">midpoint</span><span style="color: #0000FF;">=</span><span style="color: #7060A8;">floor</span><span style="color: #0000FF;">(</span><span style="color: #000000;">duration</span><span style="color: #0000FF;">/</span><span style="color: #000000;">2</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">for</span> <span style="color: #000000;">day</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #000000;">duration</span> <span style="color: #008080;">do</span>
<span style="color: #000000;">update</span><span style="color: #0000FF;">()</span>
<span style="color: #008080;">if</span> <span style="color: #7060A8;">rmdr</span><span style="color: #0000FF;">(</span><span style="color: #000000;">day</span><span style="color: #0000FF;">,</span><span style="color: #000000;">interval</span><span style="color: #0000FF;">)=</span><span style="color: #000000;">0</span> <span style="color: #008080;">or</span> <span style="color: #000000;">day</span><span style="color: #0000FF;">=</span><span style="color: #000000;">duration</span> <span style="color: #008080;">or</span> <span style="color: #000000;">day</span><span style="color: #0000FF;">=</span><span style="color: #000000;">midpoint</span> <span style="color: #008080;">then</span>
<span style="color: #000000;">sunt</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">extract</span><span style="color: #0000FF;">(</span><span style="color: #000000;">hadcovid</span><span style="color: #0000FF;">,</span><span style="color: #7060A8;">find_all</span><span style="color: #0000FF;">(</span><span style="color: #000000;">UNTREATED</span><span style="color: #0000FF;">,</span><span style="color: #000000;">category</span><span style="color: #0000FF;">))</span>
<span style="color: #000000;">sirr</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">extract</span><span style="color: #0000FF;">(</span><span style="color: #000000;">hadcovid</span><span style="color: #0000FF;">,</span><span style="color: #7060A8;">find_all</span><span style="color: #0000FF;">(</span><span style="color: #000000;">IRREGULAR</span><span style="color: #0000FF;">,</span><span style="color: #000000;">category</span><span style="color: #0000FF;">))</span>
<span style="color: #000000;">sreg</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">extract</span><span style="color: #0000FF;">(</span><span style="color: #000000;">hadcovid</span><span style="color: #0000FF;">,</span><span style="color: #7060A8;">find_all</span><span style="color: #0000FF;">(</span><span style="color: #000000;">REGULAR</span><span style="color: #0000FF;">,</span><span style="color: #000000;">category</span><span style="color: #0000FF;">))</span>
<span style="color: #000000;">unt</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">sunt</span><span style="color: #0000FF;">);</span> <span style="color: #000000;">hunt</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">sum</span><span style="color: #0000FF;">(</span><span style="color: #000000;">sunt</span><span style="color: #0000FF;">)</span>
<span style="color: #000000;">irr</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">sirr</span><span style="color: #0000FF;">);</span> <span style="color: #000000;">hirr</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">sum</span><span style="color: #0000FF;">(</span><span style="color: #000000;">sirr</span><span style="color: #0000FF;">)</span>
<span style="color: #000000;">reg</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">sreg</span><span style="color: #0000FF;">);</span> <span style="color: #000000;">hreg</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">sum</span><span style="color: #0000FF;">(</span><span style="color: #000000;">sreg</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
<span style="color: #008080;">if</span> <span style="color: #7060A8;">rmdr</span><span style="color: #0000FF;">(</span><span style="color: #000000;">day</span><span style="color: #0000FF;">,</span><span style="color: #000000;">interval</span><span style="color: #0000FF;">)=</span><span style="color: #000000;">0</span> <span style="color: #008080;">then</span>
<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"Day %d:\n"</span><span style="color: #0000FF;">,</span><span style="color: #000000;">day</span><span style="color: #0000FF;">)</span>
<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"Untreated: N = %d, with infection = %d\n"</span><span style="color: #0000FF;">,{</span><span style="color: #000000;">unt</span><span style="color: #0000FF;">,</span><span style="color: #000000;">hunt</span><span style="color: #0000FF;">})</span>
<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"Irregular Use: N = %d, with infection = %d\n"</span><span style="color: #0000FF;">,{</span><span style="color: #000000;">irr</span><span style="color: #0000FF;">,</span><span style="color: #000000;">hirr</span><span style="color: #0000FF;">})</span>
<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"Regular Use: N = %d, with infection = %d\n\n"</span><span style="color: #0000FF;">,{</span><span style="color: #000000;">reg</span><span style="color: #0000FF;">,</span><span style="color: #000000;">hreg</span><span style="color: #0000FF;">})</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
<span style="color: #008080;">if</span> <span style="color: #000000;">day</span><span style="color: #0000FF;">=</span><span style="color: #000000;">midpoint</span> <span style="color: #008080;">then</span>
<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"At midpoint, Infection case percentages are:\n"</span><span style="color: #0000FF;">)</span>
<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">" Untreated : %f\n"</span><span style="color: #0000FF;">,</span><span style="color: #000000;">100</span><span style="color: #0000FF;">*</span><span style="color: #000000;">hunt</span><span style="color: #0000FF;">/</span><span style="color: #000000;">unt</span><span style="color: #0000FF;">)</span>
<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">" Irregulars: %f\n"</span><span style="color: #0000FF;">,</span><span style="color: #000000;">100</span><span style="color: #0000FF;">*</span><span style="color: #000000;">hirr</span><span style="color: #0000FF;">/</span><span style="color: #000000;">irr</span><span style="color: #0000FF;">)</span>
<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">" Regulars : %f\n\n"</span><span style="color: #0000FF;">,</span><span style="color: #000000;">100</span><span style="color: #0000FF;">*</span><span style="color: #000000;">hreg</span><span style="color: #0000FF;">/</span><span style="color: #000000;">reg</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"At study end, Infection case percentages are:\n"</span><span style="color: #0000FF;">)</span>
<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">" Untreated : %f\n"</span><span style="color: #0000FF;">,</span><span style="color: #000000;">100</span><span style="color: #0000FF;">*</span><span style="color: #000000;">hunt</span><span style="color: #0000FF;">/</span><span style="color: #000000;">unt</span><span style="color: #0000FF;">)</span>
<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">" Irregulars: %f\n"</span><span style="color: #0000FF;">,</span><span style="color: #000000;">100</span><span style="color: #0000FF;">*</span><span style="color: #000000;">hirr</span><span style="color: #0000FF;">/</span><span style="color: #000000;">irr</span><span style="color: #0000FF;">)</span>
<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">" Regulars : %f\n\n"</span><span style="color: #0000FF;">,</span><span style="color: #000000;">100</span><span style="color: #0000FF;">*</span><span style="color: #000000;">hreg</span><span style="color: #0000FF;">/</span><span style="color: #000000;">reg</span><span style="color: #0000FF;">)</span>
<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"Final statistics: H = %f\n"</span><span style="color: #0000FF;">,</span><span style="color: #000000;">kruskal</span><span style="color: #0000FF;">({</span><span style="color: #000000;">sunt</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">sirr</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">sreg</span><span style="color: #0000FF;">}))</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">procedure</span>
<span style="color: #000000;">run_study</span><span style="color: #0000FF;">()</span>
<!--</syntaxhighlight>-->
<small>Aside: the best explanation of the Kruskal-Wallis Test I found was http://www.statisticslectures.com/topics/kruskalwallis/ but I didn't implement anything that I found there.</small>
{{out}}
<pre>
Total subjects: 10000
 
Day 30:
Untreated: N = 8628, with infection = 285
Irregular Use: N = 1372, with infection = 23
Regular Use: N = 0, with infection = 0
 
Day 60:
Untreated: N = 7460, with infection = 518
Irregular Use: N = 2366, with infection = 72
Regular Use: N = 174, with infection = 0
 
Day 90:
Untreated: N = 6497, with infection = 730
Irregular Use: N = 2399, with infection = 151
Regular Use: N = 1104, with infection = 17
 
At midpoint, Infection case percentages are:
Untreated : 11.235955
Irregulars: 6.294289
Regulars : 1.539855
 
Day 120:
Untreated: N = 5709, with infection = 879
Irregular Use: N = 2102, with infection = 220
Regular Use: N = 2189, with infection = 66
 
Day 150:
Untreated: N = 5058, with infection = 1005
Irregular Use: N = 1846, with infection = 273
Regular Use: N = 3096, with infection = 152
 
Day 180:
Untreated: N = 4486, with infection = 1109
Irregular Use: N = 1624, with infection = 315
Regular Use: N = 3890, with infection = 253
 
At study end, Infection case percentages are:
Untreated : 24.721355
Irregulars: 19.396552
Regulars : 6.503856
 
Final statistics: H = 211.421331
</pre>
 
Line 329 ⟶ 1,093:
 
Final statistics: KruskalResult(statistic=55.48204323818349, pvalue=8.95833684545873e-13)
</pre>
 
=={{header|Raku}}==
{{trans|Phix}}
<syntaxhighlight lang="raku" line># 20221025 Raku programming solution
 
enum <UNTREATED REGULAR IRREGULAR>;
my \DOSE_FOR_REGULAR = 100;
my ($nSubjects,$duration,$interval) = 10000, 180, 30;
my (@dosage,@category,@hadcovid) := (0,UNTREATED,False)>>.&{($_ xx $nSubjects).Array};
 
sub update($pCovid=1e-3, $pStartTreatment=5e-3, $pRedose=¼, @dRange=<3 6 9>) {
for 0 ..^ @dosage.elems -> \i {
unless @hadcovid[i] {
if rand < $pCovid {
@hadcovid[i] = True
} else {
my $dose = @dosage[i];
if $dose==0 && rand < $pStartTreatment or $dose > 0 && rand < $pRedose {
@dosage[i] = $dose += @dRange.roll;
@category[i] = ($dose > DOSE_FOR_REGULAR) ?? REGULAR !! IRREGULAR
}
}
}
}
}
 
sub kruskal (@sets) {
my $n = ( my @ranked = @sets>>.List.flat.sort ).elems;
my @sr = 0 xx @sets.elems;
my $ix = (@ranked.first: * == True, :k)+1,
my ($arf,$art) = ($ix, $ix+$n) >>/>> 2;
 
for @sets.kv -> \i,@set { for @set -> $b { @sr[i] += $b ?? $art !! $arf } }
 
my $H = [+] @sr.kv.map: -> \i,\s { s*s/@sets[i].elems }
return 12/($n*($n+1)) * $H - 3 * ($n + 1)
}
 
say "Total subjects: $nSubjects\n";
 
my ($midpoint,$unt,$irr,$reg,$hunt,$hirr,$hreg,@sunt,@sirr,@sreg)=$duration div 2;
 
for 1 .. $duration -> \day {
update();
if day %% $interval or day == $duration or day == $midpoint {
@sunt = @hadcovid[ @category.grep: UNTREATED,:k ];
@sirr = @hadcovid[ @category.grep: IRREGULAR,:k ];
@sreg = @hadcovid[ @category.grep: REGULAR, :k ];
($unt,$hunt,$irr,$hirr,$reg,$hreg)=(@sunt,@sirr,@sreg).map:{|(.elems,.sum)}
}
if day %% $interval {
printf "Day %d:\n",day;
printf "Untreated: N = %4d, with infection = %4d\n", $unt,$hunt;
printf "Irregular Use: N = %4d, with infection = %4d\n", $irr,$hirr;
printf "Regular Use: N = %4d, with infection = %4d\n\n",$reg,$hreg
}
if day == $midpoint | $duration {
my $stage = day == $midpoint ?? 'midpoint' !! 'study end';
printf "At $stage, Infection case percentages are:\n";
printf " Untreated : %f\n", 100*$hunt/$unt;
printf " Irregulars: %f\n", 100*$hirr/$irr;
printf " Regulars : %f\n\n",100*$hreg/$reg
}
}
 
printf "Final statistics: H = %f\n", kruskal ( @sunt, @sirr, @sreg )</syntaxhighlight>
{{out}}
<pre>Total subjects: 10000
 
Day 30:
Untreated: N = 8616, with infection = 314
Irregular Use: N = 1383, with infection = 27
Regular Use: N = 1, with infection = 0
 
Day 60:
Untreated: N = 7511, with infection = 561
Irregular Use: N = 2308, with infection = 87
Regular Use: N = 181, with infection = 1
 
Day 90:
Untreated: N = 6571, with infection = 760
Irregular Use: N = 2374, with infection = 141
Regular Use: N = 1055, with infection = 19
 
At midpoint, Infection case percentages are:
Untreated : 11.565972
Irregulars: 5.939343
Regulars : 1.800948
 
Day 120:
Untreated: N = 5757, with infection = 909
Irregular Use: N = 2097, with infection = 189
Regular Use: N = 2146, with infection = 61
 
Day 150:
Untreated: N = 5092, with infection = 1054
Irregular Use: N = 1832, with infection = 238
Regular Use: N = 3076, with infection = 132
 
Day 180:
Untreated: N = 4530, with infection = 1172
Irregular Use: N = 1615, with infection = 282
Regular Use: N = 3855, with infection = 229
 
At study end, Infection case percentages are:
Untreated : 25.871965
Irregulars: 17.461300
Regulars : 5.940337
 
Final statistics: H = 248.419454</pre>
 
=={{header|Wren}}==
{{trans|Python}}
{{libheader|Wren-math}}
<syntaxhighlight lang="wren">import "random" for Random
import "./math" for Nums
 
var Rand = Random.new()
 
var UNTREATED = 0
var IRREGULAR = 1
var REGULAR = 2
var DOSE_FOR_REGULAR = 100
 
/* A subject for the study. */
class Subject {
construct new() {
_cumDose = 0
_category = UNTREATED
_hadCovid = false
_updateCount = 0
}
 
cumDose { _cumDose }
category { _category }
hadCovid { _hadCovid }
updateCount { _updateCount }
 
cumDose(d) { _cumDose = d }
category(c) { _category = c }
hadCovid(h) { _hadCovid = h }
updateCount(u) { _updateCount = u }
 
// Daily update on the subject to check for infection and randomly dose.
update(pCovid, pStartingTreatment, pRedose, dRange) {
if (!_hadCovid) {
if (Rand.float() < pCovid) {
_hadCovid = true
} else if ((_cumDose == 0 && Rand.float() < pStartingTreatment) ||
(_cumDose > 0 && Rand.float() < pRedose)) {
_cumDose = _cumDose + Rand.sample(dRange)
categorize()
}
}
_updateCount = _updateCount + 1
}
 
// Update using default parameters.
update() { update(0.001, 0.005, 0.25, [3, 6, 9]) }
 
// Set treatment category based on cumulative treatment taken.
categorize() {
_category = (_cumDose == 0) ? UNTREATED :
(_cumDose >= DOSE_FOR_REGULAR) ? REGULAR : IRREGULAR
return _category
}
}
 
// a, b and c are assumed to be lists of boolean values.
var kruskal = Fn.new { |a, b, c|
// map the bool values to 1 (true) or 0 (false).
var aa = a.map { |e| e ? 1: 0 }.toList
var bb = b.map { |e| e ? 1: 0 }.toList
var cc = b.map { |e| e ? 1: 0 }.toList
// aggregate and sort them
var ss = (aa + bb + cc).sort()
// find rank of first occurrence of 1
var ix = ss.indexOf(1) + 1
// calculate average ranks for 0 and 1
var arf = (1 + ix - 1) / 2
var n = ss.count
var art = (ix + n) / 2
// calculate sum of ranks for each list
var sra = Nums.sum(a.map { |e| e ? art : arf })
var srb = Nums.sum(b.map { |e| e ? art : arf })
var src = Nums.sum(c.map { |e| e ? art : arf })
// calculate H
var H = 12/(n*(n+1)) * (sra*sra/a.count + srb*srb/b.count + src*src/c.count) - 3 * (n + 1)
return H
}
 
// Run the study using the population of size 'N' for 'duration' days.
var runStudy = Fn.new { |numSubjects, duration, interval|
var population = List.filled(numSubjects, null)
for (i in 0...numSubjects) population[i] = Subject.new()
var unt = 0
var untCovid = 0
var irr = 0
var irrCovid = 0
var reg = 0
var regCovid = 0
System.print("Total subjects: %(numSubjects)")
for (day in 0...duration) {
for (subj in population) subj.update()
if ((day + 1) % interval == 0) {
System.print("\nDay %(day + 1):")
unt = population.count { |s| s.category == UNTREATED }
untCovid = population.count { |s| s.category == UNTREATED && s.hadCovid }
System.print("Untreated: N = %(unt), with infection = %(untCovid)")
irr = population.count { |s| s.category == IRREGULAR }
irrCovid = population.count { |s| s.category == IRREGULAR && s.hadCovid }
reg = population.count { |s| s.category == REGULAR }
regCovid = population.count { |s| s.category == REGULAR && s.hadCovid }
System.print("Regular Use: N = %(reg), with infection = %(regCovid)")
}
if (day == (duration/2).floor - 1) {
System.print("\nAt midpoint, Infection case percentages are:")
System.print(" Untreated : %(100 * untCovid / unt)")
System.print(" Irregulars: %(100 * irrCovid / irr)")
System.print(" Regulars : %(100 * regCovid / reg)")
}
}
System.print("\nAt study end, Infection case percentages are:")
System.print(" Untreated : %(100 * untCovid / unt) of group size of %(unt)")
System.print(" Irregulars: %(100 * irrCovid / irr) of group size of %(irr)")
System.print(" Regulars : %(100 * regCovid / reg) of group size of %(reg)")
var untreated = population.where { |s| s.category == UNTREATED }.map { |s| s.hadCovid }.toList
var irregular = population.where { |s| s.category == IRREGULAR }.map { |s| s.hadCovid }.toList
var regular = population.where { |s| s.category == REGULAR }.map { |s| s.hadCovid }.toList
System.print("\nFinal statistics: H = %(kruskal.call(untreated, irregular, regular))")
}
 
runStudy.call(1000, 180, 30)</syntaxhighlight>
 
{{out}}
Sample run.
<pre>
Total subjects: 1000
 
Day 30:
Untreated: N = 882, with infection = 32
Regular Use: N = 0, with infection = 0
 
Day 60:
Untreated: N = 783, with infection = 58
Regular Use: N = 13, with infection = 0
 
Day 90:
Untreated: N = 679, with infection = 82
Regular Use: N = 91, with infection = 0
 
At midpoint, Infection case percentages are:
Untreated : 12.076583210604
Irregulars: 4.3478260869565
Regulars : 0
 
Day 120:
Untreated: N = 592, with infection = 95
Regular Use: N = 186, with infection = 5
 
Day 150:
Untreated: N = 517, with infection = 106
Regular Use: N = 286, with infection = 19
 
Day 180:
Untreated: N = 459, with infection = 118
Regular Use: N = 375, with infection = 26
 
At study end, Infection case percentages are:
Untreated : 25.708061002179 of group size of 459
Irregulars: 15.060240963855 of group size of 166
Regulars : 6.9333333333333 of group size of 375
 
Final statistics: H = 395.09607330604
</pre>
Peak
2,462

edits

Retrieved from "https://rosettacode.org/wiki/Special:MobileDiff/330770...364295"