Analysis of HCQ and azithromycin trial
T Nguyen
24/03/2020
Note
Gautret and colleagues have just published a paper on the efficacy of hydroxychloroquine (HCQ) and azithromycin in the treatment of patients infected with SARS-cov-2 (eg Wuhan virus). In the paper, the authors claim that “Despite its small sample size our survey shows that hydroxychloroquine treatment is significantly associated with viral load reduction/disappearance in COVID-19 patients and its effect is reinforced by azithromycin”. The paper can be accessed here:
https://www.sciencedirect.com/science/article/pii/S0924857920300996
I was curious of the conclusion and the idea underlying the study. The two drugs were not designed to treat SARS-cov-2 infection, so the hypothesis behind the study is quite novel. However, the more important question is whether the data are consistent with the authors’ conclusion. I noted that the authors have used a rather simple and descriptive analysis that did not account for the within-patient correlation. That is, to me, a significant weakness of the paper.
Therefore, I have downloaded and re-analyzed the data (which are provided in the paper’s appendix). My analysis is shown in this document. In contrast to the authors’ conclusion, I found no statistically significant effect of HCQ on the rate of negative tests.
Loading packages for analysis
library(lme4)## Loading required package: Matrixlibrary(ggplot2)
library(gridExtra)
library(dplyr)##
## Attaching package: 'dplyr'## The following object is masked from 'package:gridExtra':
##
## combine## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(table1)##
## Attaching package: 'table1'## The following objects are masked from 'package:base':
##
## units, units<-library(epiDisplay)## Loading required package: foreign## Loading required package: survival## Loading required package: MASS##
## Attaching package: 'MASS'## The following object is masked from 'package:dplyr':
##
## select## Loading required package: nnet##
## Attaching package: 'epiDisplay'## The following object is masked from 'package:ggplot2':
##
## alphalibrary(ICCbin) # for intraclass corr
library(sjPlot) # for model plotting## #refugeeswelcomeReading data
df = read.csv("~/Dropbox/Bao chi/Coronavirus articles/HCQ study LONG.csv")
# recoding character variables
df$outc = ifelse(df$outcome=="NEG", 1, 0)
baseline = subset(df, time==0)
final = subset(df, time==6)Some descriptive analyses
# Baseline data
table1(~age + sex + factor(outc) | hcq, data=baseline)| No (n=16) |
Yes (n=20) |
Overall (n=36) |
|
|---|---|---|---|
| age | |||
| Mean (SD) | 37.4 (24.0) | 51.2 (18.7) | 45.1 (22.0) |
| Median [Min, Max] | 32.5 [10.0, 75.0] | 51.5 [20.0, 87.0] | 47.0 [10.0, 87.0] |
| sex | |||
| F | 10 (62.5%) | 11 (55.0%) | 21 (58.3%) |
| M | 6 (37.5%) | 9 (45.0%) | 15 (41.7%) |
| factor(outc) | |||
| 0 | 16 (100%) | 20 (100%) | 36 (100%) |
table1(~age + sex + factor(outc) | azi, data=baseline)| No (n=30) |
Yes (n=6) |
Overall (n=36) |
|
|---|---|---|---|
| age | |||
| Mean (SD) | 44.6 (23.4) | 47.5 (14.1) | 45.1 (22.0) |
| Median [Min, Max] | 44.5 [10.0, 87.0] | 51.5 [20.0, 60.0] | 47.0 [10.0, 87.0] |
| sex | |||
| F | 19 (63.3%) | 2 (33.3%) | 21 (58.3%) |
| M | 11 (36.7%) | 4 (66.7%) | 15 (41.7%) |
| factor(outc) | |||
| 0 | 30 (100%) | 6 (100%) | 36 (100%) |
# Day 6 data
table1(~age + sex + factor(outc) | hcq, data=final)| No (n=16) |
Yes (n=20) |
Overall (n=36) |
|
|---|---|---|---|
| age | |||
| Mean (SD) | 37.4 (24.0) | 51.2 (18.7) | 45.1 (22.0) |
| Median [Min, Max] | 32.5 [10.0, 75.0] | 51.5 [20.0, 87.0] | 47.0 [10.0, 87.0] |
| sex | |||
| F | 10 (62.5%) | 11 (55.0%) | 21 (58.3%) |
| M | 6 (37.5%) | 9 (45.0%) | 15 (41.7%) |
| factor(outc) | |||
| 0 | 14 (87.5%) | 7 (35.0%) | 21 (58.3%) |
| 1 | 2 (12.5%) | 13 (65.0%) | 15 (41.7%) |
table1(~age + sex + factor(outc) | azi, data=final)| No (n=30) |
Yes (n=6) |
Overall (n=36) |
|
|---|---|---|---|
| age | |||
| Mean (SD) | 44.6 (23.4) | 47.5 (14.1) | 45.1 (22.0) |
| Median [Min, Max] | 44.5 [10.0, 87.0] | 51.5 [20.0, 60.0] | 47.0 [10.0, 87.0] |
| sex | |||
| F | 19 (63.3%) | 2 (33.3%) | 21 (58.3%) |
| M | 11 (36.7%) | 4 (66.7%) | 15 (41.7%) |
| factor(outc) | |||
| 0 | 21 (70.0%) | 0 (0%) | 21 (58.3%) |
| 1 | 9 (30.0%) | 6 (100%) | 15 (41.7%) |
# Logistic regression analysis
m6 = glm(outc ~ age + sex + hcq, family=binomial, data=final)
logistic.display(m6)##
## Logistic regression predicting outc
##
## crude OR(95%CI) adj. OR(95%CI) P(Wald's test)
## age (cont. var.) 1 (0.97,1.03) 0.97 (0.93,1.02) 0.252
##
## sex: M vs F 2.29 (0.59,8.91) 1.72 (0.29,10.07) 0.547
##
## hcq: Yes vs No 13 (2.27,74.32) 24.17 (2.66,219.27) 0.005
##
## P(LR-test)
## age (cont. var.) 0.228
##
## sex: M vs F 0.545
##
## hcq: Yes vs No < 0.001
##
## Log-likelihood = -17.6083
## No. of observations = 36
## AIC value = 43.2166m7 = glm(outc ~ age + sex + hcq + azi, family=binomial, data=final)
logistic.display(m7)##
## Logistic regression predicting outc
##
## crude OR(95%CI) adj. OR(95%CI) P(Wald's test)
## age (cont. var.) 1 (0.97,1.03) 0.97 (0.92,1.02) 0.264
##
## sex: M vs F 2.29 (0.59,8.91) 1.18 (0.17,8.37) 0.868
##
## hcq: Yes vs No 13 (2.27,74.32) 12.8 (1.38,118.52) 0.025
##
## azi: Yes vs No 269847183.02 (0,Inf) 98615039.9 (0,Inf) 0.994
##
## P(LR-test)
## age (cont. var.) 0.239
##
## sex: M vs F 0.868
##
## hcq: Yes vs No 0.01
##
## azi: Yes vs No 0.018
##
## Log-likelihood = -14.8101
## No. of observations = 36
## AIC value = 39.6202# Estimate intraclass correlation
iccbin(cid = id, y = outc, data = df)## Warning in iccbin(cid = id, y = outc, data = df): 'cid' has been coerced to
## a factor## Warning in iccbin(cid = id, y = outc, data = df): One or Both of 'Zou and
## Donner's Modified Wald' Confidence Limits Fell Outside of [0, 1]## Warning in iccbin(cid = id, y = outc, data = df): ICC Not Estimable by
## 'Moment with Equal Weights' Method## Warning in iccbin(cid = id, y = outc, data = df): ICC Not Estimable by
## 'Moment with Weights Proportional to Cluster Size' Method## Warning in iccbin(cid = id, y = outc, data = df): ICC Not Estimable by
## 'Modified Moment with Equal Weights' Method## Warning in iccbin(cid = id, y = outc, data = df): ICC Not Estimable by
## 'Modified Moment with Weights Proportional to Cluster Size' Method## $estimates
## Methods
## 1 ANOVA Estimate
## 2 Modified ANOVA Estimate
## 3 Moment Estimate with Equal Weights
## 4 Moment Estimate with Weights Proportional to Cluster Size
## 5 Modified Moment Estimate with Equal Weights
## 6 Modified Moment Estimate with Weights Proportional to Cluster Size
## 7 Stabilized Moment Estimate
## 8 Moment Estimate from Unbiased Estimating Equation
## 9 Fleiss-Cuzick Kappa Type Estimate
## 10 Mak's Unweighted Average Estimate
## 11 Correlation Estimate with Equal Weight to Every Pair of Observations
## 12 Correlation Estimate with Equal Weight to Each Cluster Irrespective of Size
## 13 Correlation Estimate with Weighting Each Pair According to Number of Pairs individuals Appear
## 14 Resampling Estimate
## 15 First-order Model Linearized Estimate
## 16 Monte Carlo Simulation Estimate
## ICC
## 1 0.370367364605313
## 2 0.362212276214834
## 3 -
## 4 -
## 5 -
## 6 -
## 7 0.452323103154305
## 8 0.370367364605313
## 9 0.362212276214834
## 10 0.370367364605313
## 11 0.362212276214834
## 12 0.362212276214834
## 13 0.362212276214834
## 14 0.375164413591523
## 15 0.354513821810483
## 16 0.398377088822742
##
## $ci
## Type LowerCI UpperCI
## 1 Smith's Large Sample Confidence Interval 0.2239875 0.5167472
## 2 Zou and Donner's Modified Wald Confidence Interval 0.0000000 1.0000000
## 3 Fleiss-Cuzick Confidence Interval 0.1505523 0.5738722
## 4 Pearson Correlation Type Confidence Interval 0.1473362 0.5770884
## 5 Resampling Based Confidence Interval 0.1567670 0.5935618Mixed-effects logistic regression analysis
m1 = glmer(outc ~ hcq + azi + time + hcq:time + azi:time + (1 | id), data=df, family = binomial(link="logit"))
summary(m1)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: outc ~ hcq + azi + time + hcq:time + azi:time + (1 | id)
## Data: df
##
## AIC BIC logLik deviance df.resid
## 174.4 199.1 -80.2 160.4 245
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.6858 -0.1877 -0.0817 0.0055 2.6802
##
## Random effects:
## Groups Name Variance Std.Dev.
## id (Intercept) 8.131 2.851
## Number of obs: 252, groups: id, 36
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -5.8256 1.5973 -3.647 0.000265 ***
## hcqYes 1.7120 1.7309 0.989 0.322641
## aziYes -5.0888 3.4290 -1.484 0.137800
## time 0.4545 0.2288 1.986 0.047007 *
## hcqYes:time 0.2753 0.2980 0.924 0.355525
## aziYes:time 2.7199 1.1586 2.348 0.018897 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) hcqYes aziYes time hcqYs:
## hcqYes -0.713
## aziYes 0.098 -0.240
## time -0.668 0.550 -0.031
## hcqYes:time 0.366 -0.625 0.110 -0.721
## aziYes:time -0.188 0.126 -0.871 0.060 -0.076# Plotting data
summary.data = df %>% group_by(hcq, time) %>% summarise(outcome=mean(outc))
p1 = plot_model(m1, type="pred", terms = c("time", "hcq"))
p1 = p1 + geom_point(data=summary.data, aes(x=time, y=outcome, color=hcq), inherit.aes = F) + scale_color_manual(values=c("blue", "red")) + scale_fill_manual(values=c("blue", "red")) + labs(title="Effect of hydroxychloroquine", x="Day", y="Probability of Negative Tests") + theme(legend.position="top")## Scale for 'colour' is already present. Adding another scale for
## 'colour', which will replace the existing scale.## Scale for 'fill' is already present. Adding another scale for 'fill',
## which will replace the existing scale.summary.data = df %>% group_by(azi, time) %>% summarise(outcome=mean(outc))
p2 = plot_model(m1, type="pred", terms = c("time", "azi"))
p2 = p2 + geom_point(data=summary.data, aes(x=time, y=outcome, color=azi), inherit.aes = F) + scale_color_manual(values=c('blue','red')) + scale_fill_manual(values=c('blue','red')) + labs(title="Effect of azithromycin", x="Day", y="Probability of Negative Tests") + theme(legend.position="top")## Scale for 'colour' is already present. Adding another scale for
## 'colour', which will replace the existing scale.
## Scale for 'fill' is already present. Adding another scale for 'fill',
## which will replace the existing scale.grid.arrange(p1, p2, ncol=2)