3
$\begingroup$

Would it be possible to replace the McNemar test with regression? What would be the regression formula? What type of regression should be used (logistic?, some survival?)?

Let's have a made-up example with one/matched/paired study population. In the example, the pre-intervention male proportion is 50% and that dropped to 40% after the intervention. Thus, there was a drop of 10 percentage points.

enter image description here

Woult it be correct to use logistic regression with logit-link for analysis? Y variable would be "Sex-Changed" and intercept-only model will be run?

Sex-Changed ~ 1
$\endgroup$

2 Answers 2

4
$\begingroup$

The following example will work to model with logistic regression a test similar to McNemar's test on a 2 x 2 table.

Essentially the dependent variable for a binomial logistic model takes on one value if the observation changed in one direction, and another value if the observation changed in the other direction. Concordant observations are ignored, like in McNemar's test.

At the time of writing, I don't know how to extend this model to tables larger than 2 x 2.

Matrix =as.matrix(read.table(header=TRUE, row.names=1, text="
Before  Yes   No
Yes      9     5
No      17    15
"))

Matrix

mcnemar.test(Matrix)

   ### McNemar's Chi-squared test with continuity correction
   ### McNemar's chi-squared = 5.5, df = 1, p-value = 0.01902

binom.test(17, (17+5))

   ### Exact binomial test
   ### number of successes = 17, number of trials = 22, p-value = 0.0169


### Code to convert matrix to long format adapted from
###  https://rcompanion.org/handbook/H_01.html

Counts = as.data.frame(as.table(Matrix))

colnames(Counts) = c("Before", "After", "Freq")

Long = Counts[rep(row.names(Counts), Counts$Freq), c("Before", "After")]

rownames(Long) = seq(1:nrow(Long))

### Create a new variable with values -1 for change in one direction,
###  0 for no change, and 1 for change in the other direction.

Long$ChangeSign = as.numeric(Long$Before) - as.numeric(Long$After)

Long$ChangeSign = factor(Long$ChangeSign)

Long

### Create a new data frame with only those observations where
###   the result changed.

Long2 = Long[Long$ChangeSign!=0,]

Long2

   ###    Before After ChangeSign
   ### <snip>
   ### 23     No   Yes          1
   ### 24     No   Yes          1
   ### 25     No   Yes          1
   ### 26     No   Yes          1
   ### 27    Yes    No         -1
   ### 28    Yes    No         -1
   ### 29    Yes    No         -1
   ### 30    Yes    No         -1
   ### <snip>

model = glm(ChangeSign ~ 1, data = Long2, family = binomial)

summary(model)

### Coefficients:
###   Estimate Std. Error z value Pr(>|z|)  
### (Intercept)  -1.2238     0.5087  -2.405   0.0162 *
$\endgroup$
2
  • $\begingroup$ you're absolutely right, sorry. I would appreciate your thoughts on my updated answer where I try to figure out what the difference between the approaches is. $\endgroup$
    – Eoin
    Jan 13 at 10:29
  • $\begingroup$ I don't have any deep insights into the connection between McNemar and logistic regression... The original poster had the right idea, that essentially the McNemar is about the model Change ~ 1. Except that Change needs to take on different values if the result changed from Result 1 to Result 2, and another value for the opposite change. Also, McNemar completely ignores observations where there is no change. It doesn't count them as part of the sample size. I have heard some criticisms of this. I could see another test that takes into account observations with no change. $\endgroup$ Jan 13 at 12:29
0
$\begingroup$

Update

This answer is wrong, but the reason for my mistake is illustrative enough that I think it's worth leaving up and explaining. In short, I believe that McNemar's test has the null hypothesis $P(r_1 = y\ \&\ r_2 = n) = P(r_1 = n\ \&\ r_2 = y)$, while the logistic model has the null $P(r_1 = y | r_2 = n) = P(r_1 = n | r_2 = y)$.

library(tidyverse)

# Probabilities
# Initial probabilities (e.g. 80% say yes initially)
# (Tests should have the same outcome if p0$y = p0$n = .5?)
p0 = list(y = .8, n = .2)

# Changes (eg. y2n is the probabilty of changing from yes to no)
p = list(
  y2y = .9, y2n = .1,
  n2y = .4, n2n = .6
)
prob_table = matrix(t(p), nrow = 2, 
                    dimnames = list(To = c('Y', 'N'), 
                                    From = c('Y', 'N')))
prob_table

##    From
## To  Y   N  
##   Y 0.9 0.4
##   N 0.1 0.6


N = 100
counts = N * c(
  p0$y * p$y2y, p0$y * p$y2n,
  p0$n * p$n2y, p0$n * p$n2n
) 
obs_table = matrix(t(counts), nrow = 2, 
                   dimnames = list(t2=c('Yes', 'No'),
                                   t1=c('Yes', 'No')))
obs_table

##      t1
## t2    Yes No
##   Yes  72  8
##   No    8 12


# Logistic regression tests
# if P(t2 = y | t1 = n) = P(t2 = n | t1 = y)
# In this example...
#   P(t2 = y | t1 = n) = .1
#   P(t2 = n | t1 = y) = .4

data = epitools::expand.table(obs_table) %>%
  mutate(has_changed = ifelse(t1 == t2, 1, 0))
m = glm(has_changed ~ t1, data = data, family = binomial)
summary(m)

##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   2.1972     0.3726   5.896 3.72e-09 ***
## t1No         -1.7918     0.5892  -3.041  0.00236 ** 


# McNemar tests
# if P(t2 = y & t1 = n) = P(t2 = n & t1 = y)
# In this example...
#  P(t2 = y & t1 = n) = 8/100
#  P(t2 = n & t1 = y) = 8/100

mcnemar.test(obs_table)
## 
##  McNemar's Chi-squared test
## 
## data:  obs_table
## McNemar's chi-squared = 0, df = 1, p-value = 1

Original response

In the model sex_changed ~ 1, the null hypothesis is that the log-odds of sex changing are equal to zero, corresponding to a probability of 0.5. This isn't what you want.

You want the model sex_changed ~ 1 + pre_intervention. A significant effect of pre_intervention here would mean that the probability of change for cases that start out as male is signigantly different from the probability for those that start out as female.

Note, however, that logistic regression has a perfect separation issue when all of the outcomes are the same in one of the groups, so will not give sensible results for your example data. I hope you have more than 10 cases!


data = data.frame(pre = rep(c('M', 'F'), each = 5))
data$change = c(0,0,0,0,1,0,0,0,0,0)

m = glm(change ~ 1 + pre, data = data, family = binomial)
summary(m)
Call:
glm(formula = change ~ 1 + pre, family = binomial, data = data)

Deviance Residuals: 
     Min        1Q    Median        3Q       Max  
-1.35373  -0.00008  -0.00008   0.75806   1.01077  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)
(Intercept)   -19.57    4809.34  -0.004    0.997
preM           19.97    4809.34   0.004    0.997

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 12.2173  on 9  degrees of freedom
Residual deviance:  6.7301  on 8  degrees of freedom
AIC: 10.73

Number of Fisher Scoring iterations: 18
```
$\endgroup$
0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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