Plotting our response variables

Figure 1: Antibiotic duration vs length of stay

The marginal histograms can be further broken down by their PCR era below:

Figure 2: Histogram of duration on antibiotics and length of stay

Distribution

The LOS looks pretty good, but we have an large cluster of patients who had zero hours of antibiotics (ABX) regardless of how long their LOS was. On one hand this is actually good news, since we have an unexpectedly large number of subjects in the Post-PCR era who had zero hours of antibiotics.

But while this might be good for stewardship, this presents a bit of a problem for me in modeling the data. As shown in figure 2, the distribution of ABX is bimodal, with a peak at zero and a peak farther out around ~40 hours. This means our model needs to account for why some folks didn’t get any antibiotics (regardless of PCR era) and other folks did.

To make matters even more complicated, the distributions of ABX between the Pre-PCR and Post-PCR eras aren’t the same (figure 3, below). Again, this is good because it indicates the PCR is having an effect on stewardship, but could mean that we’ll have unequal variance in our models.

Figure 3: Violin plot of antibiotic duration, split by era

For what it’s worth, accounting for the PCR era alone doesn’t fix the distribution. This can be illustrated by the QQ plot of the model: \(Y_{i} = \beta_{0} + \beta_{1}\ {\sf PCR\ era}_{i} + \epsilon_{i}\)

Figure 4: QQ plot for linear model of PCR era alone

The QQ plot shows that it’s not normally distributed, and is strongly light-tailed¹.

Balanced data

Additionally, our data is not balanced, which means that the factorial ANOVA is less resistant to violations of the assumptions of the normality². This means for our models we’ll have to use a type II or III sums of squares in our model. This can be seen below by tallying the data (see column n).

And this problem only worsens as you stratify by additional predictors in the model. For my referance, this Cross Validated post explains the difference among these methods well, and this paper outlines the pros/cons on each type³.

Additionally, here’s a blog post showing how to easily change them in R. This post shows how to use the contrast car::contr.Sum().

Factors associated with no ABX

In figure 5 and 6, we can see that Tmax and Age are both associated with not receiving any antibiotics (illustrated by the density plots), but not strongly associated with the duration of antibiotics (illustrated by the scatterplots). Specifically, a low Tmax was associated with not getting antibiotics (figure 5) while a young age was associated with getting some antibiotics (figure 6).

Figure 5: Duration of antibiotics by Tmax

Figure 6: Duration of antibiotics by age

It’s also worth mentioning that in the previous models we had already (partially) accounted for the folks who didn’t get antibiotics when we considered their LP status. Figure 7 illustrates that not getting an LP was associated with not getting any antibiotics.

Figure 7: Duration of antibiotics by LP status

Modeling these factors

This model

ABX ~ Era + Infxn + LP_Done + age

Model diagnostics: This model is great (as far as the Adjusted \(R^2\) is concerned), but is problematic because it has a singularity on LP_DoneTRUE:InfxnViral:age