Chisquared on categorical data or a Wilcoxon test on the counts hereof?

Question

This question is related to a question I had in R on SO here

The background of my question is fairly simple. I was given two "databases" in the form of data.frames that each contain about 100 patients each from two different hospitals. Each patient received multiple different antibiotics during the course of their treatment. I now wanted to calculcate all of the different classes of antibiotics given and then compare whether antibiotic "xyz" was given more frequently in hospital 1 versus hospital 2 (or vice a versa; two-tailed).

Now a bit of a catch is that the documentation in hospital 1 is better than in hospital 2, so ideally, you'd want to compare proportions of the total number of antibiotics (given/documented) and not "raw counts" in an attempt to "correct" for this difference in accuracy of the two databases.

So I thought I'd have to summarize my data in contingency tables and then run a chisq.test on it, but the way the data is summarized below, makes me think I need to run a wilcox.test.

May I ask of you experts here if

• any of you have had to deal with this kind of issue before and
• what you suggest is the best way of dealing with it?

Below you can find selected columns of the interim-summary of the respective data.frames. Hospital 1 is hosp1, hospital 2 is hosp2. The data can be pulled from the link provided above

>hosp1 ### this is the summary of hospital 1

                      id total perc
1  1st gen Cephalosporin     6  1.9
2  3rd gen Cephalosporin    65 20.5
3  4th gen Cephalosporin    10  3.2
4         Aminoglycoside    31  9.8
5           Glycopeptide    55 17.4
6            Lincosamide     2  0.6
7              Macrolide     3  0.9
8             Penicillin    36 11.4
9           Tetracycline     2  0.6
10          Trimethoprim     2  0.6
11      Ureidopenicillin    46 14.5
12            Carbapenem    19  6.0
13       Fluoroquinolone    17  5.4
14        Nitroimidazole    12  3.8
15            Antifungal     6  1.9
16         Oxazolidinone     2  0.6
17             Rifamycin     1  0.3
18           Polypeptide     1  0.3
19          Lipopeptide      1  0.3

> hosp2 ### this is the summary of hosp2

                      id total perc
1  3rd gen Cephalosporin    19  9.4
2             Carbapenem    37 18.2
3        Fluoroquinolone    24 11.8
4           Glycopeptide    32 15.8
5             Penicillin    29 14.3
6       Ureidopenicillin    36 17.7
7            Lipopeptide     4  2.0
8              Macrolide     2  1.0
9         Aminoglycoside     9  4.4
10           Polypeptide     1  0.5
11             Rifamycin     1  0.5
12          Tetracycline     1  0.5
13           Lincosamide     1  0.5
14             Quinolone     2  1.0
15           Sulfonamide     2  1.0
16        Nitroimidazole     1  0.5
17            Polymyxine     1  0.5
18              Colistin     1  0.5

Perhaps the merged data makes more sense as to what I'm aiming to compare:

new_df2 <- merge(hosp1, hosp2, by=id, all=TRUE)
                      id total.x perc.x total.y perc.y
1  1st gen Cephalosporin       6    1.9      NA     NA
2  3rd gen Cephalosporin      65   20.5      19    9.4
3  4th gen Cephalosporin      10    3.2      NA     NA
4         Aminoglycoside      31    9.8       9    4.4
5           Glycopeptide      55   17.4      32   15.8
6            Lincosamide       2    0.6       1    0.5
7              Macrolide       3    0.9       2    1.0
8             Penicillin      36   11.4      29   14.3
9           Tetracycline       2    0.6       1    0.5
10          Trimethoprim       2    0.6      NA     NA
11      Ureidopenicillin      46   14.5      36   17.7
12            Carbapenem      19    6.0      37   18.2
13       Fluoroquinolone      17    5.4      24   11.8
14        Nitroimidazole      12    3.8       1    0.5
15            Antifungal       6    1.9      NA     NA
16         Oxazolidinone       2    0.6      NA     NA
17             Rifamycin       1    0.3       1    0.5
18           Polypeptide       1    0.3       1    0.5
19          Lipopeptide        1    0.3      NA     NA
20           Lipopeptide      NA     NA       4    2.0
21             Quinolone      NA     NA       2    1.0
22           Sulfonamide      NA     NA       2    1.0
23            Polymyxine      NA     NA       1    0.5
24              Colistin      NA     NA       1    0.5

And then basically run something like:

with(new_df2(chisq.test(total.x[id=="1st gen Cephaolosporin], total.y[id=="1st gen Cephaolosporin])

Would the wilcox.test() only apply if I were comparing total counts along the whole column? I'm just getting a bit confused here, because effectively the number in the columns total.x and total.y represent counts and for some reason I'm now thinking of a wilcoxon test...but it's all based on categorical data.

Thus, making a contingency table out of this would actually require you to reformat the table into something along the lines of this (to perform for example a comparison of "3rd generation Cephalosporins")

thirdgenhosp1 <- rep(c("Yes", "No"), times=c(65, (sum(new.df2$total.x)-65)))
thirdgenhosp2 <- rep(c("Yes", "No"), times=c(19, (sum(new.df2$total.y)-19)))

### combine the two and try to "correct" for the difference in accuracy of documentation

thirdgen_all <- cbind(thirdgenhosp1, 
c(thirdgenhosp2, rep(NA, length(thirdgenhosp1)-length(thirdgenhosp2)))) 

### then make a data.frame out of this to be able to analyse it

thirdgen_all_df <- data.frame(thirdgen_all)
names(thirdgen_all_df)[2] <- "thirdgenhosp2"

# then perform the comparison
with(thirdgen_all_df, chisq.test(thirdgenhosp1, thirdgenhosp2, correct=F))

Would there be a more efficient way of doing this? And am I actually doing the right thing?

Thanks for any help whatsoever. This has now been bountied :).

Answer 1

It seems that you want to analyze each antibiotic separately. So, above table as a whole is not important. What you are asking is how to analyze each row separately. For this you may use prop.test in R which compares 2 proportions. You can replace NAs by 0. For example, for 1st gen cephalosporin:

> sum.x= sum(mydf$total.x)
> sum.y= sum(mydf$total.y)
> 
> sum.x
[1] 317
> sum.y
[1] 203
> 
> prop.test(c(6,0), c(sum.x, sum.y))

        2-sample test for equality of proportions with continuity correction

data:  c(6, 0) out of c(sum.x, sum.y)
X-squared = 2.4047, df = 1, p-value = 0.121
alternative hypothesis: two.sided
95 percent confidence interval:
 -0.0001137171  0.0379686067
sample estimates:
    prop 1     prop 2 
0.01892744 0.00000000 

Warning message:
In prop.test(c(6, 0), c(sum.x, sum.y)) :
  Chi-squared approximation may be incorrect

You can get all P values in one go using following code:

for(i in 1:nrow(mydf)){
    cat("",as.character(mydf[i,'id']),
       ': P=',prop.test(c(mydf[i,'total.x'],mydf[i,'total.y']), 
       c(sum.x, sum.y))$p.value,'\n')
}

 1g_Cephalosporin : P= 0.1209697 
 3g_Cephalosporin : P= 0.001167431 
 4g_Cephalosporin : P= 0.02588154 
 Aminoglycoside : P= 0.03911204 
 Glycopeptide : P= 0.7244962 
 Lincosamide : P= 1 
 Macrolide : P= 1 
 Penicillin : P= 0.3956558 
 Tetracycline : P= 1 
 Trimethoprim : P= 0.6834551 
 Ureidopenicillin : P= 0.3895501 
 Carbapenem : P= 2.186687e-05 
 Fluoroquinolone : P= 0.01242847 
 Nitroimidazole : P= 0.03955267 
 Antifungal : P= 0.1209697 
 Oxazolidinone : P= 0.6834551 
 Rifamycin : P= 1 
 Polypeptide : P= 1 
 Lipopeptide : P= 1 
 Lipopeptide : P= 0.04609961 
 Quinolone : P= 0.296246 
 Sulfonamide : P= 0.296246 
 Polymyxine : P= 0.8220467 
 Colistin : P= 0.8220467 

Lack of documentation should not affect the proportions obtained.

Answer 2

I am unclear if the units of interest are patients within hospitals or purely aggregated hospital data. The former seems a more natural structure for the data. The recording bias is probably not resolved by a percentage adjustment. If a person was definitely admitted to a hospital with an infection, some antibiotic would be given. If an admission occurred and no antibiotic was given, one hospital is certainly better than the other. A statistical test may not be required.

Overall, measuring the association between hospital and use of any given antibiotic or aggregate classes of antibiotics requires adjustment for types of infections and identified organisms.

I do not use R, so my comments below don't have the benefit of working with your provided data.

If you would like to compare the number of antibiotics classes administered to a single patient, you should use poisson regression, as the use of an antibiotic is neither binary, nor rank-ordered. You could then regress the count of antibiotic classes used per patient using hospital as the variable of interest. Confounder adjustment could then be handled straightforwardly. Overall, I believe this would be the best approach and would make the best use of all data available.

If you would like to assess the likelihood of a single class of antibiotic being used for any given patient by hospital, this would also be poisson data, as the denominator is not clearly defined (many different types of infections, allergies to antibiotics, etc) and the use of a particular antibiotic for any given case is not a Bernoulli trial.

Todd