1
$\begingroup$

I have data on patients who have two gastric samples drawn to culture bacteria. I would like to see if gastric pH (categorized as acidic, neutral, or basic) has an influence on the dichotomous outcome of culture (growth or no growth). The two samples from each patient will have the same pH, but may have discordant culture outcomes.

Sample data look like this:

sample_data <- data.frame(studyid = rep(1:20, each = 2),
                          sample_number = rep(1:2, times = 20),
                          result = rep(c(0,0,1,0,1,0,1,1,0,0,0,1,1,1,0,1,1,0,1,1),2),
                          pH_category = rep(c(3,2,2,1,1,2,1,2,3,3,1,2,3,3,2,3,1,1,1,2), each = 2))

head(sample_data)
#   studyid sample_number result pH_category
# 1       1             1      0           3
# 2       1             2      0           3
# 3       2             1      1           2
# 4       2             2      0           2
# 5       3             1      1           2
# 6       3             2      0           2

The sample table looks like this:

# Table
sample_table <- table(culture = sample_data$result,
                      ph_category = sample_data$pH_category)
sample_table

#       ph_category
#culture 1 2 3
#      0 6 7 5
#      1 8 7 7

I assume a Chi-squared test is not appropriate since they are paired data, however the McNemar Test and McNemar-Bowker test are also not appropriate since I am not comparing the first sample to the second.

Is there another test to examine this table's distribution? Despite the frequency of this question, I can't find any other answers on CV (here, here, or here).

Would running a GEE logistic model to account for within-subject correlation be the best alternative?

$\endgroup$
1

1 Answer 1

2
$\begingroup$

One possibility is to use a random intercepts logistic regression model, which can be fitted in R with lme4. Something like

library(lme4)
sample_data <- data.frame(studyid = rep(1:20, each = 2),
                          sample_number = rep(1:2, times = 20),
                          result = rep(c(0,0,1,0,1,0,1,1,0,0,0,1,1,1,0,1,1,0,1,1),2),
                          pH_category = rep(c(3,2,2,1,1,2,1,2,3,3,1,2,3,3,2,3,1,1,1,2), each = 2))

sample_data$pH_category <- factor(sample_data$pH_category)   

mod0 <- glmer(result  ~  0 + pH_category  + (1 | studyid), data=sample_data,
              family=binomial, nAGQ=9)

but this gives a singular fit, which might be OK, but let us try with some regularization:

library(blme)

mod0.b <- bglmer(result  ~  0 + pH_category  + (1 | studyid), 
    data=sample_data, family=binomial, nAGQ=9)

 summary(mod0.b)
Cov prior  : studyid ~ wishart(df = 3.5, scale = Inf, posterior.scale = cov, common.scale = TRUE)
Prior dev  : -0.6106

Generalized linear mixed model fit by maximum likelihood (Adaptive
  Gauss-Hermite Quadrature, nAGQ = 9) [bglmerMod]
 Family: binomial  ( logit )
Formula: result ~ 0 + pH_category + (1 | studyid)
   Data: sample_data

     AIC      BIC   logLik deviance df.resid 
    64.0     70.8    -28.0     56.0       36 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-1.1312 -0.8000  0.5621  0.7039  1.0000 

Random effects:
 Groups  Name        Variance Std.Dev.
 studyid (Intercept) 1.502    1.226   
Number of obs: 40, groups:  studyid, 20

Fixed effects:
              Estimate Std. Error z value Pr(>|z|)
pH_category1 3.763e-01  7.818e-01   0.481    0.630
pH_category2 3.138e-05  7.569e-01   0.000    1.000
pH_category3 4.307e-01  8.370e-01   0.515    0.607

Correlation of Fixed Effects:
            pH_ct1 pH_ct2
pH_categry2 0.000        
pH_categry3 0.011  0.000 

We can also show the (fixed effects part of) the fit with package visreg

library(visreg)
visreg(mod0.b, scale="response", alpha=0.01)

resulting in:

visreg plot of fixed effects

$\endgroup$
2
  • $\begingroup$ Thank you - can you clarify why you reduced the intercept to 0? If you had not graciously provided code, I would have specified the formula as result ~ pH_category + (1 | studyid) $\endgroup$
    – jpsmith
    Commented Sep 8, 2022 at 20:43
  • 1
    $\begingroup$ When there is only one factor variable, removing the intercept does not change the model, it only makes for more clear interpretation. If you compare the fit with the results from your model in the comment (do it!), you will see that fitted values or predicted values) are the same, AIC, BIC, deviance etc are all the same, and the differences between the fitted coefs for the three pH levels are the same. $\endgroup$ Commented Sep 8, 2022 at 21:11

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.