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

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.

• See the proportional odds model implementation of the rank difference test in hbiostat.org/bbr/nonpar Jan 28 at 14:37

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