I am using glmer from lm4 in R to run an explanatory IRT analysis. The aim is to model person parameters with items nested in persons and to control for age (the test scale has 10 items; all persons respond to every item). I want to run a cross validation and I have a training dataset (700 persons distributed over an age range of 7 years, 10 items). The test dataset (cross validation) comprises 84000 persons. I would like to set up the model based on the training data, retrieve the person parameters (while controlling for age) and the to apply it to the test data.

The data is in long format with an ID specifying the person, an age variable and the response in each item. Since I use simulated data, I know the latent ability of each person and the original item parameter. The code for the modelling the training data:

          # explanatory IRT modelling
          glmer.model <- glmer(response ~ -1 + age + item + (1 | ID), family=binomial("logit"), data=long.train, control=control)
          glmer.item.params.train <- fixef(glmer.model)
          glmer.person.params.train <- coef(glmer.model)$ID[, 1]

Both, the item parameters and the person parameters of the simulated data are well replicated (though not perfectly, given the limited reliability of the simulated scale). The original item parameters and those modelled with glmer correlate with r = .997. Great!

My problem lies in applying the glmer model to the test data. The test data has the same structure. I am trying to predict the responses for each item and aggregate over the persons to obtain the person parameters with the following code, but it completely fails to produce plausible results:

          long.cross$predicted.prob <- predict(glmer.model, newdata=long.cross, re.form=NA, type="response")
          person.abil <- aggregate(predicted.prob ~ ID, data=long.cross, FUN=mean)

How do I get the fixed effects for the predicted data (-> predicted person parameters)?

I tried to vary the predict call and how the data is aggregated over the single responses.

  • $\begingroup$ are you looking for the ranef function? Did you have a look at De Boeck et al. 2011. The Estimation of Item Response Models with the lmer Function from the lme4 Package in R? $\endgroup$
    – Tom
    Commented Aug 27, 2023 at 20:17
  • $\begingroup$ Many thanks for the suggestion. Yes, I read the de Broek et al. paper. The problem is to apply a model to new data and to extract the fixed effects. ranef extracts the random effects. $\endgroup$ Commented Aug 28, 2023 at 6:19
  • $\begingroup$ can you add example data (maybe use some of existing packages) $\endgroup$
    – Tom
    Commented Aug 28, 2023 at 7:10
  • $\begingroup$ Many thanks! I have prepared an OSF project and uploaded simulated training and cross validation data. It is available via osf.io/ahu6w I think, my last reply is flawed. I am not looking at the fixed effect, but at the intercepts - at least I think so - of the person parameters via "coef". The goal is to get the person parameter, with an controlled age effect. $\endgroup$ Commented Aug 29, 2023 at 8:10

1 Answer 1


Trying to answer your question, I learnt several things; so you might excuse a lengthy answer. Also, there might be some shortcut I'm not aware of. However, at the moment I don't know of any standard method that does what you're looking for. There might be in the future. I did come up with a prediction that works in your case, though.

Note that in your case the ranef (ranef(glmer.model)$ID[, 1]) and the coef function (coef(glmer.model)$ID[, 1]) produce the same result, which is the coditional mode of the random effects (see JSS paper on Fitting Linear Mixed-Effects Models Using lme4, section 5.1).

> all.equal(
+   ranef(glmer.model)$ID[, 1], 
+     coef(glmer.model)$ID[, 1])
[1] TRUE

In specialized IRT software you'll most often find a similar method that extracts person parameters (e.g., ltm::factor.scores, similar in mirt and TAM). Most often this method is able to predict person parameters given a response model and new responses. However in lme4 that does not seem to be the case.

Lessons learnt: the different predict options

Given the de Boeck Paper or his handout or this stats.stackexchange discussion we find that under the present model, given some predictors $X$,

$$\textrm{logit}\left(P\left(Y_{pi}=1\right)\right) \sim \sum_{k=0}^{K}b_{k}X_{ik}+\theta_{p}+\varepsilon_{pi}.$$

From the help file for the predict method for objects of class mermod we finde several options

  • re.form: (formula, NULL, or NA) specify which random effects to condition on when predicting. If NULL, include all random effects; if NA or ~0, include no random effects.

Here, we can only choose between all (re.form = NULL) or no random effect (re.form = NA), as we only have one random effect, i.e. whether or not we include the prediction of $\theta_{p}$ in the total prediction.

  • random.only: (logical) ignore fixed effects, making predictions only using random effects?

I.e. whether or not we include the prediction of $\sum_{k=0}^{K}b_{k}X_{ik}$ in the total prediction.

  • type: character string - either "link", the default, or "response" indicating the type of prediction object returned.

I.e., whether we predict $\textrm{logit}\left(P\left(Y_{pi}=1\right)\right)$ (with type="link"; i.e. the z-scores from your data frame) or $P\left(Y_{pi}=1\right)$ (with type="response").

  • allow.new.levels: logical if new levels (or NA values) in newdata are allowed. If FALSE (default), such new values in newdata will trigger an error; if TRUE, then the prediction will use the unconditional (population-level) values for data with previously unobserved levels (or NAs).

I.e., in case of (re.form = NULL) and new observations. However, since the population model for a single random effect (with mean 0) is the fixed part of the model, there is no difference for the test data in whether you use re.form = NULL (and allow.new.levels = TRUE) or re.form = NA. E.g. you would find that

> #
> summary(predict(object = glmer.model, long.cross, random.only = TRUE,
+                 allow.new.levels = TRUE, type = "link"))
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
      0       0       0       0       0       0 

You can check that behavior with predict(object = glmer.model, long.train,...) and also find that

> all.equal(aggregate(
+   predict(object = glmer.model, long.train, random.only = TRUE, type = "link"), 
+   by = long.train["ID"], FUN = mean)$x, 
+      coef(glmer.model)$ID[, 1])
[1] TRUE


> all.equal(
+   predict(object = glmer.model, long.train, random.only = TRUE, type = "link"), 
+     (predict(object = glmer.model, long.train, re.form = NULL, type = "link") - 
+        predict(object = glmer.model, long.train, re.form = NA, type = "link"))
+ )
[1] TRUE

and that, predicting manually,

> long.train$predicted.prob <- predict(object = glmer.model, long.train, 
+                                      re.form = NA, type = "response")
> long.train$itemdiff <- fixef(glmer.model)[paste0("item", long.train$item)]
> long.train$agediff <- long.train$age * fixef(glmer.model)["age"]
> long.train$predicted.prob2 <- long.train$agediff + long.train$itemdiff
> all.equal(
+   car::logit(long.train$predicted.prob), 
+     long.train$predicted.prob2
+ )
[1] TRUE

Using a trick to get the required predictions

If, in typical IRT software, there is no method to extract person parameters given a response model and new responses, there might be some tricks to get your data in the correct shape in order to force the method to (only) estimate the conditional mode of the random effects (or any other person parameter). One would be to estimate a new model with constrained item parameters; e.g. like this. Another would be to use imputation techniques like the mice.impute.ml.lmer function. Unfortunately I didn't get either to work.

However, further analysing the predicted z-scores from the training data, we find a perfect relationship between age and the predicted z-scores given a specific raw score (which is no surprise, since in Rasch-type IRT models the raw score is a sufficient statistic for the person parameter).

> long.train$predicted.z <- round(predict(object = glmer.model, long.train, 
+                                         random.only = TRUE, type = "link"), 8)
> trafo <- aggregate(long.train["response"], 
+                    long.train[c("ID", "age", "predicted.z")], sum)
> plot(trafo[, c("age", "predicted.z")], pch = 19, col = factor(trafo$response))

enter image description here

That is, we can use this relationship between age together with raw score and the predicted z-scores to generate a transformation table

> cross.data$data$predicted.z <- NA
> for(rr in sort(unique(trafo$response))){
+   md.lm <- lm(predicted.z~age, trafo[trafo$response == rr, ])
+   ind <- cross.data$data$raw == rr
+   cross.data$data$predicted.z[ind] <- 
+     predict(md.lm, cross.data$data[ind, "age", drop = FALSE])
+ }
Warning message:
In predict.lm(md.lm, cross.data$data[ind, "age", drop = FALSE]) :
  prediction from rank-deficient fit; attr(*, "non-estim") has doubtful cases
> cor(cross.data$data$z, cross.data$data$predicted.z)
[1] 0.745833

Or, similar

> md.lm3 <- lm(predicted.z~age*as.factor(response), trafo)
> cross.data$data$predicted.z3 <- 
+   predict(md.lm3, setNames(cross.data$data[, c("raw", "age")],
+                            c("response", "age")))
Warning message:
In predict.lm(md.lm3, setNames(cross.data$data[, c("raw", "age")],  :
  prediction from rank-deficient fit; attr(*, "non-estim") has doubtful cases
> cor(cross.data$data$z, cross.data$data$predicted.z3)
[1] 0.7460221

The warning message means that this approach works best if there are at least two cases of different age within each raw score. If we trust that the relationship is the same (that the lines are parallel), we can use

> md.lm2 <- lm(predicted.z~age+as.factor(response), trafo)
> cross.data$data$predicted.z2 <- 
+ predict(md.lm2, setNames(cross.data$data[, c("raw", "age")],
+                          c("response", "age")))
> cor(cross.data$data$z, cross.data$data$predicted.z2)
[1] 0.7427774
  • 1
    $\begingroup$ Tom, you are a genius. Many thanks for the huge amount of time and expertise, you spent on this question. I understand your solution and I think, it pretty much touches the role, IRT models play in psychometrics - not only from a data analysis perspective, but as well in individual diagnostics after the test has been constructed. If this is ok for you, I would like to contact you via mail. There could be potential for further discussion and maybe participation in a methodology paper. Our focus is to unite explanatory IRT and regression-based norming, BTW. $\endgroup$ Commented Aug 31, 2023 at 19:42
  • $\begingroup$ If you like my response, kindly mark it as an answer. I'll send you an e-mail regarding the other matter. $\endgroup$
    – Tom
    Commented Sep 4, 2023 at 7:19

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.