Selection bias in clinical sciences: Difference between revisions
m
→{{header|jq}}: typo
m (→{{header|Wren}}: Changed to Wren S/H) |
m (→{{header|jq}}: typo) |
||
| (3 intermediate revisions by 2 users not shown) | |||
Line 53:
<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}}==
| |||