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