# Pseudomonas aeruginosa VAP 8 vs 15 days

This is an attempt to make sense of the statistics of the paper “Comparison of 8 versus 15 days of antibiotic therapy for Pseudomonas aeruginosa ventilator-associated pneumonia in adults: a randomized, controlled, open-label trial” published in Intensive Care Med

The authors claim to show no difference between short and long-duration treatments, but their evidence contradicts their claims.

```
library(magrittr)
library(scales)
set.seed(2022-05-13)
```

I will focus on the authors bayesian analysis, as it provides clearer insight, in my opinion.

The authors perform a betabinomial analysis with equal beta priors for the event rate in both groups centered on previously described event rates, mean = 35,7%. Then vary their effective samples size to generate variance.

Lets apply some sanity checks to their priors.

Effective sample size 5

```
samples <- 10000
a <- 1.785
b <- 3.215
d <- rbeta(samples, a, b)
plot(density(d), xlim = c(0, 1), xlab = "Event rate", main = "Prior event rate for ESS = 5")
```

`qbeta(c(0.025, 0.975), a, b)`

`## [1] 0.04875906 0.76979588`

`qbeta(c(0.1, 0.9), a, b)`

`## [1] 0.1116929 0.6338099`

Seems overly broad for a prior, we have a priori knowledge from other studies that event rates are around 20-50%. This prior places 80% of the probability between 11% and 63%.

Also the prior probability that the event is rare (<5% event rate) is 2.59% and the prior probability that the event is near inevitable (>95% event rate) is 0.09%

So ESS = 5 is out.

Lets see if ESS=10 is better

```
a <- 3.57
b <- 6.43
d <- rbeta(samples, a, b)
plot(density(d), xlim = c(0, 1), xlab = "Event rate", main = "Prior event rate for ESS = 10")
```

`mean(d < 0.05) %>% percent(0.01)`

`## [1] "0.26%"`

`mean(d > 0.95) %>% percent(0.01)`

`## [1] "0.00%"`

`qbeta(c(0.025, 0.975), a, b)`

`## [1] 0.1086034 0.6592063`

`qbeta(c(0.1, 0.9), a, b)`

`## [1] 0.1743354 0.5536904`

Seems quite a bit better, still a little broader than i would prefer but fair enough.

The prior probability that the event is rare (<5% event rate) is now 0.26% and the prior probability that the event is near inevitable (>95% event rate) is now 0.00%.

I’ll accept this as a *not not sane* prior.

Not lets get into the meat of the ITT analysis, but broaden our scope to look at the probability of one treatment group being worse off by any amount. We’ll keep the priors from ESS=10.

```
c <- a + 25 # 15 day with event + prior
d <- b + 98 - 25 # 15 day with no event + prior
d <- rbeta(samples, c, d)
plot(density(d), xlim = c(0, 1), xlab = "Event rate", main = "Comparing posteriors\n 15 days in black, 8 days in red")
e <- a + 31 # 8 day with event + prior
f <- b + 88 - 31 # 8 day with no event + prior
d2 <- rbeta(samples, e, f)
lines(density(d2), col = "red")
```

Visually there is clear distinction, but what is statistics if not a way to stop us from eyeballing important things.

`mean(d2-d > 0) %>% percent(0.01)`

`## [1] "91.01%"`

There is a **91% probability that 8 days is worse by any amount**.

Finally let’s reproduce their estimate from supplementary 2

`mean(d2 - d < 0.1) %>% percent(0.01)`

`## [1] "57.22%"`

Close enough. I’m not using exact methods here so a perfect concurrence isn’t expected on the trailing digits.

Citation: Bouglé A, Tuffet S, Federici L, et al (2022) Comparison of 8 versus 15 days of antibiotic therapy for Pseudomonas aeruginosa ventilator-associated pneumonia in adults: a randomized, controlled, open-label trial. Intensive Care Med. https://doi.org/10.1007/s00134-022-06690-5