Analysis for ordinal categorical outcome

I am not sure how to evaluate my data correctly, since the outcome (dependent variable) is on an ordinal categorical scale and there are two independent variables.

To simplify my experiment:

Each subject is assigned to one of three different Times (morning, noon or afternoon) and receives one Treatment (A,B or C). Afterwards a test is performed and the outcome is an ordinal categorical Score (1, 2 or 3).

The data can be summarized in three tables (one for each Time):

Time = Morning

Score
Treatment  "1"     "2"     "3"
A        19      10      20
B        6       7       5
C        17      16      11

Time = Noon

Score
Treatment  "1"     "2"     "3"
A        23      14      12
B        15      33      15
C        14      15      19

Time = Afternoon

Score
Treatment  "1"     "2"     "3"
A        11      15      12
B        17      18      15
C        9       15      19


So the research question is: which treatment dependent on the time results in the best scores?

If the outcome would be some continuous normal distributed data, I would apply pair-wise ANOVAs and post hoc t-tests.

However in the given situation I have no experience and googling did not help me, that is why I am asking here for help.

Here it is stated that for tables larger than 2x2, in the case that the chi-squared test is significant, post hoc chi-squared tests of the format 2x2 can be applied.

My question is the following: Do I summarize the data in one big table, something like:

Morning/A
Morning/B
Morning/C
Noon/A
Noon/B
Noon/C
...


applying the chi-squared test and then following by pair-wise tests (like Morning/A vs Morning/B)? Or which statistical approach is the appropriate with post hoc tests to answer the research question?

P.S. The analysis is done with R, so if anyone has a tutorial or a link with an example that would be really helpful.

Edit 1)

I decide to edit my question and append my current approach, as it is still a part of the current question, if I should open a new question please let me know.

Taking the answer of @(Frank Harrell) into account I tried using "Ordinal Logistic Regression" (OLR).

My hypothesis is that:

Treatment C results in higher scores for Noon and Afternoon.

However reading the course Material from Prof. Harrell, I have the "idea" that the hypothesis should be:

Treatment interacts with Time.

Question E1) Am I correct about the adjustment of the hypothesis for (OLR)?

So this is what I tried(I should say that my data has two more score values than in my example above, it was designed exemplary, in my case there are the following scores 0,1,2,3 and 4 ):

polr(formula = Outcome ~ Treatment + Time , data = Data, Hess = T)


getting the following results

                          Value       Std. Error   t value    p value
TreatmentB               -0.26373775  1.0059962 -0.26216576   0.793
TreatmentC               -0.62342836  0.8070730 -0.77245595   0.440
TimeNoon                  0.08236084  0.8257454  0.09974121   0.921
TimeAfterNoon             0.74163895  0.9282499  0.79896475   0.424
TreatmentB:TimeNoon      -0.53242839  1.1327578 -0.47002844   0.638
TreatmentC:TimeNoon       0.19885488  0.9406062  0.21141141   0.833
TreatmentB:TimeAfterNoon -0.56413612  1.2321034 -0.45786428   0.647
TreatmentC:TimeAfterNoon -1.32190703  1.0447585 -1.26527520   0.206
0|1                      -3.96290826  0.8499135 -4.66271949   0.000
1|2                      -2.60173952  0.7747232 -3.35828277   0.001
2|3                      -0.56198118  0.7460030 -0.75332292   0.451
3|4                       2.46049479  0.7878248  3.12314968   0.002


Questions E2)

As the interaction are all not significant, $$H_0$$ can not be rejected? So the conclusion here is: The outcome does not depend on the interaction between Treatmentand Time?

Let us assume the interactions would have been significant (TreatmentB:TimeNoon and TreatmentC:TimeNoon) and $$H_0$$ is rejected:

My conclusion would be that the best scores are reach for Treatment B and C during Noon, is this assumption correct? In the second step I would like to make a statement which of the two treatments is the better one. Is it enough to check whether on of both is significant or which approach is the correct one here? Do I need further post hoc test to make a statement like this?

For ordinary response variables Y, use a method that takes advantage of the ordinality, to increase power by reducing the number of parameters in the model. Ordinary $$\chi^2$$ tests for contingency (frequency) tables involve a lot of parameters because they consider Y as polytomous (multinomial; categorical).
To handle the repeated-within-subject nature of the data, you can incorporate random subject effects. So combining a proportional odds ordinal logistic model with random effects is advised, e.g., R package ordinal. Even better (fewer approximations) would be a Bayesian mixed effects ordinal model.
• Thank you Prof. Harrell for your fast response. It seems I haven't made it clear, the study was not designed in a repeated-within-subject nature, each participant was only tested for one Time with one Treatment. Nevertheless, thanks for guiding me into the direction of ordinal and Bayesian models. Do you have any links where I can see how to apply this approaches for hypotheses-testing, especially in the case with two IVs? Commented May 28, 2019 at 13:39
• It's looking reasonable. One suggestion: statisticians are trying to get away from significance testing with thresholds for 'significance'. So I would fit the pre-specified model and make contrasts to get comparisons of interest, taking into account interacting factors. Even if you don't have within-subject correlations, the course notes case study on longitutinal modeling will show how time-related contrasts are made. You could use the R rms package lrm or orm functions for fitting and its contrast function (full internal name contrast.rms) for getting any contrasts(odds ratios) Commented Jun 2, 2019 at 11:32