Calculate (quasi) AIC for mixed-effect baseline-category (multinomial logit) model in R

Question

I am doing a discrete choice experiment where respondents are presented with different patient profiles, and for each profile, respondents need to choose one (out of three) treatment options. An example of the dataset could be downloaded at https://github.com/trangtph/DCE_test/blob/main/dce_toy.csv.

To investigate the effect of patient characteristics on treatment choice, I am running a mixed-effect baseline-category model with respondent ID as random effect. I use function mblogit in package mclogit.

The full model is:

mtnom_full <- mblogit(choice ~ response*resist_prob + ambu_regimen*resist_prob + 
                      expo_his*resist_prob, random = ~ 1|id, data=dce_toy, 
                      method ="PQL", estimator = "ML")

Model output:

I need to perform model selection to choose what interaction effects stay in the model (the main effects are always kept in the model). I would like to use AIC. However in the manual of mclogit it says that the random-effect model is estimated using penalized quasi-likelihood techniques based on Laplace approximation. I read that for quasi-likelihood model the AIC is not applicable, and Quasi-AIC should be used.

I have tried to calculate QAIC for the mblogit model using following code but it returns NA:

MuMIn::QAIC(mtnom_full, chat= deviance(mtnom_full)/df.residual(mtnom_full))

Warning: 'chat' given is < 1, increased to 1 [1] NA

Yet using the conventional AIC still gives an output

stats::AIC(mtnom_full)

[1] 531.3431

My questions:

1. Is the AIC values for that model + nested models reliable to use for model selection?
2. If the QAIC should be used instead, how to calculate QAIC for mblogit model in R?
3. Is there any other model selection approaches for this issue? Thanks a lot!