user
In epidemiology, retrospective analyses have well-known limitations compared to prospective studies.

One such limitation is the occurrence of *selection bias* in the choice of subjects between treated
and untreated groups about whom the data is collected. For example, a treatment may have only been
given to persons who were less severely ill, which would bias the results in favor of such subjects
appearing to have done better because of the treatment when the biased group is then compared to those
who who did not receive the study treatment. Or, in a retrospective study, there may a choice to place
subjects in a particular study group using a method which is inadvertently biased by the outcome being
measured. Creating a programming example of a simulation of such selection bias in the design of a
retrospective study is the topic of this task.

The genuine, historical example (only partially approximated in this task) is of a study done of persons who, over a course
of 180 days, may or may not have become infected with Covid-19. Prior to becoming ill, these subjects may or may not have taken an available medication, which was actually taken on a particular schedule not used here, but is approximated by stating the medication was taken in doses of 3, 6, or 9 mg daily. The historical study then divided its subjects into three groups based on their cumulative dosage of the study medication:

* Group UNTREATED were those who did not take the study medication at all before they got Covid-19, including those who exited the study period without Covid-19 and having never taken the study medication.

* Group IRREGULAR is those who took the study medication but whose cumulative dose was less than a certain amount (approximated for our purposes as 100 mg) before they either came down with Covid-19 during the study or the study period ended.

* Group REGULAR is those who took (our approximation is >= 100 mg) of the study medication either before they came down with Covid-19 or took >= 100 mg by the end of the study and never became infected during the study.

Assumptions for the study simulation programming task

* Daily risk of getting Covid-19 infection for each subject was 0.1% per day, or 18% over the 180 cumulative days of the study.

* The probability of starting treatment medication for anyone not already taking it was 0.5% per day. For those who started medication, the chance of continuing the treatment was increased 50-fold to 25% each day, since most who started the medication continued to take it to some extent.

* Study dose per day is random between the approximation for the simulation of 3, 6 and 9 mg. The daily cumulative dosage is used to determine the group the subject is in, unless a subject develops Covid-19. If a subject was diagnosed with Covid-19, their group at the time of that diagnosis is used in the statistical analysis of that group.

Task

* Create a simulation of the subjects, keeping track of their medication dosages, group membership, and Covid-19 status during the study.

* 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 Kruskal 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 Kruskal-Wallis test.

* You should get a statistical result highly favoring the REGULAR group.

Stretch task

* Show monthly outcomes.

A note regarding outcome: Note that by simulation design all subjects must have an IDENTICAL risk, that is 0.1 per cent or p = 0.001 per day, of developing Covid-19. Because of the design, any statistical differences between the groups CANNOT come from an influence of the treatment on that risk, but must come from some other feature of the study design.

See also
:   * Wikipedia - Selection Bias * Wikipedia - Retrospective studies * UCDavis (pdf) Challenges in Retrospective studies
assistant
**Translation of**: Phix

```
# 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 )

```

Output:

```
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
```