# Chi-squared for repeated measures for 2 by n

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))

#   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?

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:

• 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) Commented Sep 8, 2022 at 20:43
• 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. Commented Sep 8, 2022 at 21:11