Load data (not shown here)

Descriptive plots

Look at the distribution of age and gender across treatment group

Reproduce figure 1

Which compares patients receiving chloroquine to the control group. Note that this differs slightly from the original plot, as we did not consider missing data. The general trend seems unchanged.

Fisher test comparing patient counts at day 6

We still get a significant effect of treatment, even if of lower magnitude

d = data[data$time_c==6,]
t2 = table(d$chloroquine,d$outcome)
fisher.test(t2)
## 
##  Fisher's Exact Test for Count Data
## 
## data:  t2
## p-value = 0.02094
## alternative hypothesis: true odds ratio is not equal to 1
## 95 percent confidence interval:
##  0.009066969 0.771905850
## sample estimates:
## odds ratio 
##  0.1118726

Fit figure 1

Using a mixed-effects logistic regression of outcome, with chloroquine and time as fixed effects, and a random intercept per patient. This shows a main effect of time and an interaction between time and chloroquine, indicating an increased effect of treatment over time. Note that the goodness of fit is rather poor..

m1 = glmer(outcome ~ chloroquine * time_c  + (1| patient), 
           data = data, family = binomial) # more complex models fail to converge

summary(m1)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: outcome ~ chloroquine * time_c + (1 | patient)
##    Data: data
## 
##      AIC      BIC   logLik deviance df.resid 
##    169.4    186.0    -79.7    159.4      198 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.1286 -0.1846  0.0687  0.2670  1.5357 
## 
## Random effects:
##  Groups  Name        Variance Std.Dev.
##  patient (Intercept) 12.2     3.494   
## Number of obs: 203, groups:  patient, 36
## 
## Fixed effects:
##                       Estimate Std. Error z value Pr(>|z|)    
## (Intercept)             6.0140     1.7999   3.341 0.000834 ***
## chloroquineYes         -0.4786     1.8621  -0.257 0.797148    
## time_c                 -0.4814     0.2454  -1.962 0.049808 *  
## chloroquineYes:time_c  -0.9028     0.3710  -2.433 0.014963 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) chlrqY time_c
## chloroqunYs -0.685              
## time_c      -0.631  0.538       
## chlrqnYs:t_  0.111 -0.565 -0.583
avg_data = data %>% group_by(chloroquine,time_c,time) %>% summarise(outcome=mean(outcome))

plot_m1 = plot_model(m1,type='pred',terms = c('time_c','chloroquine'))
plot_m1 + geom_point(data=avg_data,aes(time_c,outcome,color=chloroquine),inherit.aes = F)+ 
  scale_color_manual(values=c('black','red'))+
  scale_fill_manual(values=c('black','red'))

Same with age as covariate of no interest

Results remain qualitatively similar

m1_a = glmer(outcome ~ chloroquine * time_c + age + (1| patient), 
           data = data, family = binomial) # more complex models fail to converge

summary(m1_a)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: outcome ~ chloroquine * time_c + age + (1 | patient)
##    Data: data
## 
##      AIC      BIC   logLik deviance df.resid 
##    169.4    189.3    -78.7    157.4      197 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.96578 -0.18816  0.05496  0.25056  1.63779 
## 
## Random effects:
##  Groups  Name        Variance Std.Dev.
##  patient (Intercept) 11.7     3.421   
## Number of obs: 203, groups:  patient, 36
## 
## Fixed effects:
##                       Estimate Std. Error z value Pr(>|z|)  
## (Intercept)            4.49112    1.90162   2.362   0.0182 *
## chloroquineYes        -1.44067    2.02867  -0.710   0.4776  
## time_c                -0.50861    0.25473  -1.997   0.0459 *
## age                    0.04663    0.03445   1.354   0.1758  
## chloroquineYes:time_c -0.85218    0.37165  -2.293   0.0218 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) chlrqY time_c age   
## chloroqunYs -0.420                     
## time_c      -0.540  0.557              
## age         -0.422 -0.387 -0.130       
## chlrqnYs:t_  0.187 -0.564 -0.612 -0.008

Same model as m1 applied on symptomatic patients only

The interaction between time and chloroquine becomes marginal, although we lose quite some data

m1_s = glmer(outcome ~ chloroquine * time_c  + (1| patient) , 
           data = data[!is.na(data$outcome) &data$status!="Asymptomatic",], family = binomial) # more complex models fail to converge
summary(m1_s)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: outcome ~ chloroquine * time_c + (1 | patient)
##    Data: data[!is.na(data$outcome) & data$status != "Asymptomatic", ]
## 
##      AIC      BIC   logLik deviance df.resid 
##    114.1    129.5    -52.0    104.1      157 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.2385 -0.1048  0.0427  0.1837  1.3816 
## 
## Random effects:
##  Groups  Name        Variance Std.Dev.
##  patient (Intercept) 10.77    3.282   
## Number of obs: 162, groups:  patient, 30
## 
## Fixed effects:
##                       Estimate Std. Error z value Pr(>|z|)  
## (Intercept)             7.4250     3.0582   2.428   0.0152 *
## chloroquineYes         -1.3640     3.1565  -0.432   0.6657  
## time_c                 -0.2981     0.5783  -0.515   0.6063  
## chloroquineYes:time_c  -1.0743     0.6538  -1.643   0.1004  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) chlrqY time_c
## chloroqunYs -0.867              
## time_c      -0.781  0.747       
## chlrqnYs:t_  0.597 -0.772 -0.876
avg_data = data %>% filter(status!='Asymptomatic') %>% 
  group_by(chloroquine,time_c,time) %>% summarise(outcome=mean(outcome))

plot_m1_s = plot_model(m1_s,type='pred',terms = c('time_c','chloroquine'))
plot_m1_s + geom_point(data=avg_data,aes(time_c,outcome,color=chloroquine),inherit.aes = F)+ 
  scale_color_manual(values=c('black','red'))+
  scale_fill_manual(values=c('black','red'))

Reproduce figure 2

This time comparing three groups: no treatment, chloroquine only, and chloroquine + azithromycin

ggplot(data,aes(time,outcome,group=treatment,color=treatment,fill=treatment)) +
  stat_summary(fun.data = mean_cl_boot,geom='ribbon',alpha=.3,color=NA) + 
  stat_summary(fun.y = mean,geom='line',size=1.2)+
  stat_summary(fun.y = mean,geom='point',size=4)+
  theme(legend.position='bottom') + 
  labs(y = 'Percentage of patients with PCR positive samples', x = 'Day')

Fit figure 2

Using a mixed-effects logistic regression of outcome, with treatment and time as fixed effects, and a random intercept per patient. This shows a main effect of time, and an interaction between time and treatment. Contrasts show a difference over time between [chloroquine + azithromycin] and [chloroquine].

m2 = glmer(outcome ~ treatment * time_c  + (1| patient) , data = data[!is.na(data$outcome),], family = binomial) # more complex models fail to converge
Anova(m2)
## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: outcome
##                    Chisq Df Pr(>Chisq)    
## treatment         3.6512  2  0.1611245    
## time_c           12.2128  1  0.0004746 ***
## treatment:time_c  6.7888  2  0.0335600 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(m2)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: outcome ~ treatment * time_c + (1 | patient)
##    Data: data[!is.na(data$outcome), ]
## 
##      AIC      BIC   logLik deviance df.resid 
##    161.5    184.6    -73.7    147.5      196 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.5631 -0.0775  0.0654  0.2457  3.6834 
## 
## Random effects:
##  Groups  Name        Variance Std.Dev.
##  patient (Intercept) 9.255    3.042   
## Number of obs: 203, groups:  patient, 36
## 
## Fixed effects:
##                                   Estimate Std. Error z value Pr(>|z|)   
## (Intercept)                          9.385      3.631   2.584  0.00976 **
## treatmentcholoroquine only          -4.516      3.604  -1.253  0.21016   
## treatmentnone                       -3.881      3.630  -1.069  0.28497   
## time_c                              -3.523      1.243  -2.835  0.00459 **
## treatmentcholoroquine only:time_c    2.552      1.211   2.108  0.03504 * 
## treatmentnone:time_c                 3.061      1.249   2.451  0.01423 * 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) trtmno trtmnt time_c tonl:_
## trtmntchlro -0.911                            
## treatmentnn -0.891  0.872                     
## time_c      -0.918  0.820  0.800              
## trtonly:tm_  0.895 -0.867 -0.808 -0.975       
## trtmntnn:t_  0.901 -0.812 -0.842 -0.982  0.960
avg_data = data %>% filter(!is.na(outcome)) %>% 
  group_by(treatment,time_c,time) %>% summarise(outcome=mean(outcome))

plot_m2= plot_model(m2,type='pred',terms = c('time_c','treatment'))
plot_m2 + geom_point(data=avg_data,aes(time_c,outcome,color=treatment),inherit.aes = F)