I have a dataset with 0/1 response data for a set of test questions. There are approx 7.5k items that have been randomly allocated to participants in an online quiz. Just over 10,000 participants (set to increase to around 25,000) have responded 50 questions each. The majority of cells in the data set are therefore blank.
How do you run a IRT model in this environment? Please help.
A very simple model for this would be a fixed or random effects logistic regression with random effects for both participant and question on the logit of a correct response. The fixed or random effects for participants and questions represent how skilled participants are and how easy questions are (you can make this more complex, if e.g. questions are in different topics or some other structure exists that you want to take into account, or if you want to somehow capture that certain people are better at certain types of questions).
The key bit is to fit this type of model in such a way that you get estimates of the effects for each question (and participant, but that's not important for what you want to do). E.g. a simple fixed effects model would be something like glm(correct ~ participant + question, family=binomial(link='logit')) after making sure that participant and question are factor variables (i.e. not being treated as continuous numbers). That may be a model you can fit despite the sparse nature of the data, but you might also run into trouble (e.g. if for one question everyone got it right or wrong, or if one person got everything right or wrong, which with so many questions and people seems almost certain - plus there could be more complicated examples of complete separation).
An approach to get around issues of too sparse data would be to fit a random effects model in such a way that you get estimates of the random effects. Doing so is easy by default when you fit this model in a Bayesian way (a fixed effects model with the Firth penalty is another idea). For example, in R one could use the brms package e.g. like this: brmfit1 = brm(formula = correct ~ (1|participant+question), family=bernoulli()) and then you can do (using also the tidyverse and tidybayes packages) brmfit1 %>% tidybayes::spread_draws(Intercept, r_question[participant,RE_is_on]) %>% mutate(easiness_score = Intercept + r_question).
