How to model and plot a distribution of choice frequency (categorical variable) split by 2 levels of within-subject

Question

I measured the difference in people's habits between this time last year and now. Specifically, I want to explore the difference in group size when people go to camping: do people tend to go with fewer people this Fall (covid era) compared to Fall of last year?

I go ahead and survey 100 people at random. Unfortunately, my sample is not representative in terms of gender and age. I got 20% females and 80% males, which then means that it doesn't reflect the typical 50-50 split in the population. It's also not representative in terms of age, as the age in my sample is 40, but in the population it's about 50. Therefore, I would like that my analysis will account for correction of the gender and age variables as well.

Data

library(tidyverse)
library(magrittr)
library(truncnorm)


set.seed(13)
df <- tibble(id = rep(seq(1, 50), each = 2),
       is_male = rep(sample(c(0, 1), prob = c(0.2, 0.8), size = 50, replace = T), each = 2),
       age = round(rtruncnorm(n = 100, a=20, b = 80, mean= 25, sd=25.09)),
       tp = rep(c('this_year', "last_year"), time = 50),
       n_of_ppl = sample(c("1_5", "6_10", "11_15", "16_20", "20_plus"), size = 100, replace = TRUE))

> df
## # A tibble: 100 x 5
##       id is_male   age tp        n_of_ppl
##    <int>   <dbl> <dbl> <chr>     <chr>   
##  1     1       1    29 this_year 6_10    
##  2     1       1    25 last_year 11_15   
##  3     2       1    46 this_year 20_plus 
##  4     2       1    28 last_year 11_15   
##  5     3       1    65 this_year 1_5     
##  6     3       1    39 last_year 20_plus 
##  7     4       1    66 this_year 1_5     
##  8     4       1    32 last_year 11_15   
##  9     5       0    60 this_year 6_10    
## 10     5       0    32 last_year 11_15   
## # ... with 90 more rows

If I had just wanted to see the difference in distributions between this year and last year, regardless of within-subjectness or gender or age, I could have done the following

df %<>% mutate(across(n_of_ppl, fct_relevel, "1_5", "6_10", "11_15", "16_20", "20_plus")) 
  

help <- function(count, group) {
  count / tapply(count, group, sum)[group]
}


df %>%
  ggplot(data = ., aes(x = n_of_ppl, y = ..prop.., group = tp)) + 
  geom_histogram(aes(y = help(..count.., ..group..)), stat = "count" , fill = "darkblue") +
  geom_text(aes(label = scales::percent(help(..count.., ..group..), accuracy = 1),
                y = help(..count.., ..group..) ), stat= "count", vjust = -.5, color = "darkblue") +
  scale_y_continuous(labels = scales::percent) +
  facet_grid(~ tp)

But I believe that my data should be modeled more accurately by taking into account the within-subject nature of the data. At the same time, I want to correct for gender and age representativeness, probably using a prediction function (predict?), where is_male = 0.5 and age = 50. How can I get a side-by-side bar plot like the one above, only that it would take into account both within-subjectness and age & gender correction?