Relationship between statistical models: mcnemar test vs logistic regression

Question

I always find useful and reliable answers in this forum and I hope this pattern will not change this time.

I was spending some time to check whether an intervention program had any effect on an outcome. Unfortunately, all examples I got from stats book did not reply to this question, once the majority of books create "simple" examples very distant from real data or empirical questions.

Let's imagine that the goal of qn intervention is to increase an "Assertive" Response and decrease all other responses.

To make it clear, let's say 60 participants were submitted to "psychotherapy" vs "Psychodynamics" and we have a "control" group. The participant's attitude was evaluated before and after the treatment using an ordinal category (assertive, passive and neutral).

The question is very simple: did the intervention work?

I know the Mcnemar test is a right option for this kind of question, however, Mcnemar does not include covariates (such as sex). In addition, Mcnemar does not "adjust" for different baselines. Some participants started with a "passive" attitude, whereas others started with an "assertive" attitude.

I have a "feeling" that an ordered logistic regression (considering random effects) can deal with this question, but I'm struggling to figure out how to model that using R.

I really appreciate any help and you can reproduce this code in R to clarify my question:

The closest we have in this forum is here.

Please inform if this question is duplicated or very hard to understand.\ Thanks

library(tidyverse)
set.seed(123)
ds <- data.frame(ID=seq.int(nrow(ds)),
                 resp_t2=rep(c("agressive","assertive","neutral"),times=c(5,40,15)),
                 resp_t1=rep(c("agressive","assertive","neutral")),
                 sex=c("male","female"),
                 group=c("Psychotherapy","Psychodynamics","Control"))

mcnemar.test(ds$resp_t2, ds$resp_t1)

Answer 1

Because your dependent variable is ordinal in nature, ordinal regression is probably your best bet. Luckily the ordinal package in R allows for fixed effects and mixed effects models.

If I understand correctly, you want a data frame with variables Individual, Time, Sex, Group, Response. Your model might look something like this: Response ~ Time + Group + Time:Group + (1|Individual).

I think the effect of Time:Group will tell you what you want to know, though you may want to look at the effect of Time|Group. For this, you might look at the emmeans package.

Edit:

Below is some code for ordinal regression in R. I'm not saying that is is necessarily to the correct model for your situation. Also, please be sure to read the documentation for the ordinal package, as there are certain assumptions that must be met that aren't addressed here.

I changed the data.

Install packages.

if(!require(ordinal)){install.packages("ordinal")}
if(!require(multcompView)){install.packages("multcompView")}
if(!require(emmeans)){install.packages("emmeans")}

Create data.

set.seed(122)

ds <- data.frame(Individual = rep(1:60, 2),
                 Time = rep(1:2, each=60),
                 Response=c(sample(c("agressive","assertive","neutral"),60, replace=TRUE),
                            sample(c("agressive","assertive","neutral",
                                     "neutral","neutral","neutral"),
                                     20, replace=TRUE),
                            sample(c("agressive","assertive","assertive",
                                     "assertive","assertive","neutral"),
                                     20, replace=TRUE),
                            sample(c("agressive","agressive","agressive",
                                     "agressive","assertive","neutral"),
                                     20, replace=TRUE)),
                 Sex=c("male","female"),
                 Group=rep(c("Psychodynamics","Control","Psychotherapyl"), each=20))

ds$Individual = factor(ds$Individual)

ds$Time = factor(ds$Time)

ds$Response = factor(ds$Response, ordered=TRUE,
                 levels=c("neutral","assertive","agressive"))

Ordinal regression and anova-like table.

library(ordinal)

model = clmm(Response ~ Time + Group + Time:Group + (1|Individual), data=ds)

library(emmeans)

joint_tests(model)

   ### model term df1 df2 F.ratio p.value
   ### Time         1 Inf   0.441  0.5065
   ### Group        2 Inf   2.871  0.0567
   ### Time:Group   2 Inf   6.817  0.0011

Post-hoc comparisons. In this case for Time within each Group.

marginal = emmeans(model, ~ Time|Group)

pairs(marginal)

   ### Group = Control:
   ###  contrast   estimate        SE  df z.ratio p.value
   ###  1 - 2     0.1040226 0.5990195 Inf   0.174  0.8621
   ### 
   ### Group = Psychodynamics:
   ###  contrast   estimate        SE  df z.ratio p.value
   ###  1 - 2     2.1633167 0.6972041 Inf   3.103  0.0019
   ### 
   ### Group = Psychotherapyl:
   ###  contrast   estimate        SE  df z.ratio p.value
   ###  1 - 2    -1.5462085 0.6548526 Inf  -2.361  0.0182

Instead, you could look at the EM means as a compact letter display for the interaction. This can be useful for the purposes of plotting the results, though the emmeans themselves may not be easy for your audience to interpret.

marginal = emmeans(model, ~ Time:Group)

CLD(marginal, Letters=letters)

   ###  Time Group                emmean        SE  df  asymp.LCL  asymp.UCL .group
   ###   2    Psychodynamics -1.788468998 0.5554225 Inf -2.8770770 -0.6998610  a    
   ###   1    Psychotherapyl -0.279180862 0.4459297 Inf -1.1531871  0.5948253  ab   
   ###   2    Control        -0.105695451 0.4132884 Inf -0.9157258  0.7043349  ab   
   ###   1    Control        -0.001672842 0.4710520 Inf -0.9249178  0.9215721  ab   
   ###   1    Psychodynamics  0.374847718 0.4375040 Inf -0.4826443  1.2323397   b   
   ###   2    Psychotherapyl  1.267027620 0.4995453 Inf  0.2879369  2.2461183   b   
   ###  
   ###  Confidence level used: 0.95 
   ###  P value adjustment: tukey method for comparing a family of 6 estimates 
   ###  significance level used: alpha = 0.05 

As an alternative, you could look at the simple means, with the ordinal categories coded as 1, 2, 3.

if(!require(FSA)){install.packages("FSA")}

library(FSA)

Summarize(as.numeric(Response) ~ Time + Group, data=ds)

###   Time          Group  n mean        sd min Q1 median   Q3 max
### 1    1        Control 20 2.00 0.9176629   1  1      2 3.00   3
### 2    2        Control 20 1.95 0.6048053   1  2      2 2.00   3
### 3    1 Psychodynamics 20 2.15 0.7451598   1  2      2 3.00   3
### 4    2 Psychodynamics 20 1.40 0.6805570   1  1      1 2.00   3
### 5    1 Psychotherapyl 20 1.90 0.7880689   1  1      2 2.25   3
### 6    2 Psychotherapyl 20 2.45 0.7591547   1  2      3 3.00   3