Which test should be used for an ordinal dependent variable with pretest-posttest and two factors?

Question

Consider a pretest-posttest design, where the dependent variable (y) is measured by 3 ordinal values 1,2,3. There are also two factors: One factor (say A) consists of a control and an experimental group. The other factor (say B) consists of two age groups. The total sample size is 60, where each possible combination of the factors A and B groups has 15 subjects (The total possible combinations is 4, so 4*15=60). Now we want to test the following cases:

1. The effect of A on y,
2. The effect of A and B on y.

If y was continuous following a normal distribution, we could use analysis of covariance or ANCOVA (one-way ANCOVA for the 1st case and two-way ANCOVA for the 2nd case). But here y is not normal, nor can be transformed to a normal variable, since it has only 3 ordinal values. So what statistical test do you suggest for these 2 cases?

It would be also appreciated if you could answer this in case of y being a dichotomous variable.

Answer 1

In addition to @Frank Harrell's suggestion, you can also use proportional odds logistic regression.

set.seed(18478893)

X <- data.frame(A = rep(c("Control", "Test"), each = 30),
                B = rep(c("Age1", "Age2"), each = 15, times = 2),
                p1 = rep(c(.5, .4, .3, .2), each = 15),
                p2 = rep(c(.1, .05, .2, .1), each = 15))
X$p3 <- 1 - X$p1 - X$p2

X$y <- apply(X, 1, function(x) {
  sample(1:3, size = 1, replace = FALSE, prob = c(x["p1"], x["p2"], x["p3"]))
})

X$y <- ordered(X$y)
str(X)
#> 'data.frame':    60 obs. of  6 variables:
#>  $ A : chr  "Control" "Control" "Control" "Control" ...
#>  $ B : chr  "Age1" "Age1" "Age1" "Age1" ...
#>  $ p1: num  0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 ...
#>  $ p2: num  0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 ...
#>  $ p3: num  0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 ...
#>  $ y : Ord.factor w/ 3 levels "1"<"2"<"3": 1 1 3 3 1 1 1 3 2 1 ...

mod1 <- MASS::polr(y ~ A*B, data = X, method = "logistic")
summary(mod1)
#> 
#> Re-fitting to get Hessian
#> Call:
#> MASS::polr(formula = y ~ A * B, data = X, method = "logistic")
#> 
#> Coefficients:
#>               Value Std. Error t value
#> ATest        0.6729     0.6847  0.9829
#> BAge2        2.5095     0.9282  2.7035
#> ATest:BAge2 -1.4428     1.1764 -1.2265
#> 
#> Intercepts:
#>     Value   Std. Error t value
#> 1|2 -0.0930  0.5242    -0.1775
#> 2|3  0.7207  0.5356     1.3458
#> 
#> Residual Deviance: 101.3168 
#> AIC: 111.3168
confint(mod1, parm = "ATest:BAge2")
#> Waiting for profiling to be done...
#> 
#> Re-fitting to get Hessian
#>      2.5 %     97.5 % 
#> -3.9071558  0.8110326

# you can take the model fitting process from here

# if y is binary

X$y2 <- sapply(X$p1, function(x) rbinom(1, 1, x))

mod2 <- glm(y2 ~ A*B, data = X, family = binomial("logit"))
summary(mod2)
#> 
#> Call:
#> glm(formula = y2 ~ A * B, family = binomial("logit"), data = X)
#> 
#> Deviance Residuals: 
#>     Min       1Q   Median       3Q      Max  
#> -0.9005  -0.9005  -0.7876   1.4823   1.6259  
#> 
#> Coefficients:
#>               Estimate Std. Error z value Pr(>|z|)
#> (Intercept) -6.931e-01  5.477e-01  -1.266    0.206
#> ATest       -6.053e-17  7.746e-01   0.000    1.000
#> BAge2       -3.185e-01  8.006e-01  -0.398    0.691
#> ATest:BAge2  0.000e+00  1.132e+00   0.000    1.000
#> 
#> (Dispersion parameter for binomial family taken to be 1)
#> 
#>     Null deviance: 73.304  on 59  degrees of freedom
#> Residual deviance: 72.986  on 56  degrees of freedom
#> AIC: 80.986
#> 
#> Number of Fisher Scoring iterations: 4

^{Created on 2024-01-28 with reprex v2.0.2}