elens.R

# I analyse Dutch GP Rob Elens' data using a script
# which I earlier used for Raoult (Marseilles).
#
# I will perform some alternative and simple frequentist and Bayesian analyses.
# 
# Dutch family doctor Elens first had 25 patients whom he gave the
# then standard treatment: 12 of 25 died.
# Then he gave the next 10 chloroquine, they all 10 recovered.
# He was then ordered to stop that treatment. 
# His next patient of course (no chloroquine) ... died. 
# And so his story was reported in the media.
# https://www.elsevierweekblad.nl/domein/americandreamers/kennis/achtergrond/2020/05/trump-prijst-hydroxychloroquine-aan-wat-zijn-de-gevaren-van-het-middel-757234/
# https://www.youtube.com/watch?v=ctALgnCck1g
# Rob Elens said (of course) that his patients were not selected in any way
# He does not tell us anything about sex, age, prior conditions, ...
# If he had not had this experience, he would not have been reported in the media.
# Anyway, it does not surprise me that in the early days of an epidemic 
# the kind of patients who come to the door of their family doctor depends
# very strongly on *time*.
#
# Anyway, in total, Elens has:
# 13 of 26 "control" died, 13 recovered
# 0 of 10 "treatment" died, 10 recovered
#
# I start with using the very conventional "Fisher exact test", see
# https://stat.ethz.ch/R-manual/R-devel/library/stats/html/fisher.test.html.
# The model is: two independent binomials, and we condition on row and column 
# totals, getting a hypergeometric distribution under the null hypothesis.
#
# NB Sometimes he says "10" sometimes he says "8". That does not increase confidence
numbers <- matrix(c(13, 13, 0, 10), 2, 2, byrow = TRUE,
                  dimnames = list(Treatment = c("control", "treatment"),
                                  Outcome = c("bad outcome", "good outcome")))
numbers

##            Outcome
## Treatment   bad outcome good outcome
##   control            13           13
##   treatment           0           10

fisher.test(numbers)

## 
##  Fisher's Exact Test for Count Data
## 
## data:  numbers
## p-value = 0.005848
## alternative hypothesis: true odds ratio is not equal to 1
## 95 percent confidence interval:
##  1.720653      Inf
## sample estimates:
## odds ratio 
##        Inf

# The treatment group does overwhelmingly better than the control group.
#
# Being "exact" has a disadvantage, namely reduction of power. 
# The Barnard test does not have that disadvantage.
# https://www.rdocumentation.org/packages/Barnard/versions/1.6/topics/barnard.test
# The model is simply: two independent binomials.

library(Barnard)
barnard.test(13, 13, 0, 10)

## 
## Barnard's Unconditional Test
## 
##            Treatment I Treatment II
## Outcome I           13           13
## Outcome II           0           10
## 
## Null hypothesis: Treatments have no effect on the outcomes
## Score statistic = -2.79751
## Nuisance parameter = 0.141 (One sided), 0.535 (Two sided)
## P-value = 0.00383363 (One sided), 0.00549331 (Two sided)

# Nothing changed

# Now we go Bayesian

library(abtest)
ab_out <- ab_test(data = list(y1 = 13, n1 = 26, y2 = 10, n2 = 10), 
                  posterior = TRUE)
print(ab_out)

## Bayesian A/B Test Results:
## 
##  Bayes Factors:
## 
##  BF10: 6.997186
##  BF+0: 14.27489
##  BF-0: 0.1953377
## 
##  Prior Probabilities Hypotheses:
## 
##  H+: 0.25
##  H-: 0.25
##  H0: 0.5
## 
##  Posterior Probabilities Hypotheses:
## 
##  H+: 0.8667
##  H-: 0.0119
##  H0: 0.1214

plot(ab_out)

# Note: I see errors in the description of the output of ab_test() !!!!

plot_posterior(ab_out, what = "logor") # This is the posterior, conditional on inequality!

# Posterior probability that odds of a good outcome are more than 20% better:

logOR <- ab_out$post$H1$logor
(1 - ab_out$post_prob["H0"]) * mean(logOR > log(1.2))

##        H0 
## 0.8531781

# I should also compute the prior probability that odds of a good outcome
# are more than 20% better. It's of course

0.50 * pnorm( log(1.2))

## [1] 0.2861674

# So the data certainly changed my opinion. 
# Before I was 50% sure that there was no treatment effect at all, 
# and 25% sure there was a bad effect.
# Now I should be only 13% sure of no effect or a bad effect.
# On the other hand, this was not many patients, no randomisation,
# no attempt at all to control for confounders (not even known ones),
# ... the whole thing is little better than "hearsay" evidence.
# Yes, so do some more careful trials to figure out if this does help
# at all and which patients it helps. There is strong evidence it is
# bad for seriously ill patients. There are reasonable arguments why it
# could well be good *for some people*, 
## *before they actually get sick of Corona infection*.
## The Marseilles group is performing such a trial and are now working with
## good statisticians.