# How to fit Graded Response Model with lme4::glmer

Thanks to Rijmen et al.(2003), we can fit GRM to the data with lme4::glmer.

I think Rasch model is straightforward, with data.frame with columns like this

 response  person  item
0         1      1
0         1      2
1         1      3
...
1         2      1
0         2      2


we can fit Rasch model like this

 glmer(response ~ -1 + item + (1|person), data=   , family="binomial")


But how about GRM? The data would be like this

 response  person  item
2         1      1
4         1      2
3         1      3
...
1         2      1
4         2      2
...


For a Likert scale (1 to 5). I thought converting the data like this

 response  person item  category
1          1    1       2
0          1    1       3
1          1    2       4
0          1    2       5


Because for person1, item1, the response is 2, which means that for response 2, it's yes and for response 3, it's no.

The model would be

 response ~ item:category + (1|person)


But I am not quite sure this is the right way to do...

Note: person, item, category variables are all factors

According to De Boeck et al. (2011), GRM cannot be fitted with lmer which is rather in contrast to Rijmen et al(2003).

Now I think I am pretty sure it will work, at least for GRM with no slope parameter.

Data should be coded like this.

response  person item  category
0          1    1       1
1          1    1       2
1          1    1       3
1          1    1       4
1          1    1       5  (which is always true so should be omitted.)


Main benefit of using GLMM for IRT model is you can put other covariates (person, item, person-item) into the model.

And for GRM, you can set the difference between the ordinal response is the same, which can't be handled by ordinary GRM function, for example, ltm::grm. (Oh, I see ordinal::clmm can handle this, but I doubt it can be useful for a model like this)

  response ~ item + (1 + category|person)


or this

  response ~ item + (-1 + category|item) + (1|person)


in this case, category is integer and would be better if coded as -2, -1, 0, 1, 2.

References

Rijmen, F., Tuerlinckx, F., De Boeck, P., & Kuppens, P. (2003). A nonlinear mixed model framework for item response theory. Psychological methods, 8(2), 185.

De Boeck, P., Bakker, M., Zwitser, R., Nivard, M., Hofman, A., Tuerlinckx, F., & Partchev, I. (2011). The estimation of item response models with the lmer function from the lme4 package in R. Journal of Statistical Software, 39(12), 1-28.

====

Here's my source.

library(ltm)
#Science[c(1,3,4,7)]
Sci.df <- Science[c(1,3,4,7)] # Comfort, Work, Future, Benefit
Sci.df$id = 1:nrow(Sci.df) Sci.long <- reshape(Sci.df, varying=colnames(Sci.df[-5]), v.names="Response", timevar="item", idvar=c("id"), direction="long") Sci.long$id <- as.factor(Sci.long$id) Sci.long$item <- as.factor(Sci.long$item) library(ordinal) Sci.long.clmm <- clmm(Response ~ (1|id)+item, data=Sci.long, threshold="flexible", nAGQ=-21) summary(Sci.long.clmm) Positive1=as.integer(Sci.long$Response)<=1
Positive2=as.integer(Sci.long$Response)<=2 Positive3=as.integer(Sci.long$Response)<=3

Sci.long.sep1=Sci.long
Sci.long.sep1$Response=1; Sci.long.sep1$Positive=Positive1

Sci.long.sep2=Sci.long
Sci.long.sep2$Response=2; Sci.long.sep2$Positive=Positive2

Sci.long.sep3=Sci.long
Sci.long.sep3$Response=3; Sci.long.sep3$Positive=Positive3

Sci.long.sep = rbind(Sci.long.sep1, Sci.long.sep2, Sci.long.sep3)

Sci.long.sep$Response=as.factor(Sci.long.sep$Response)

Sci.long.sep.glmm <- glmer(Positive ~ -1 + Response + item + (1|id), data=Sci.long.sep, family=binomial,
nAGQ=21, control=glmerControl(optimizer="optimx",
optCtrl=list(method="nlminb"), check.conv.grad= .makeCC("warning", tol = 1e-4, relTol = NULL) ))
summary(Sci.long.sep.glmm)


I tried my best to make it same for clmm and glmer... but the log likelihood is different.

logLik = -1730.6 for glmer logLik = -1633.5 for clmm

and the parameters r not the same but similar.

Does anyone know why the log likehoods are different?

• I'm pretty sure you can't fit GRM in glmer, as the gl stands for generalized linear and the GRM is a generalized nonlinear model. This means the structural part of the model is not a linear predictor but is of the form $a*(b+c)$ with $a,b,c$ unknown. Linear would be something like $b+c$ (e.g., a Rasch model). I therefore agree with de Boeck. – Momo Sep 4 '14 at 16:30
• 1st comment: As far as I know, Rijmen et al. do not say that one can fit a GRM using lme4::glmer(), and indeed I don't see how one could--please provide a page reference in Rijmen et al. if you disagree. 2nd comment: I think you can fit GRM in ordinal::clmm(), although I haven't personally looked into it closely. 3rd comment: I don't know what your actual question is. Please edit your question to very clearly state what it is you want to know or are confused about. – Jake Westfall Sep 4 '14 at 16:56
• @Momo, I guess you're right. Sorry for not mentioning. I thought of GRM with constrained slope parameter. – KH Kim Sep 5 '14 at 1:14
• @Jake Westfall, As far as I figure, he is alluding that Polytomous Data can be modeled by GLMM(Generalized Linear Mixed Model, I am not referring to any particular package here). See page 191 of Rijmen et al.(2003) $L_{nij}=x_{nij} beta + z'_{nij} theta _n$ So everything's fine here. Link function and linear predictors – KH Kim Sep 5 '14 at 1:15
• @Jake Westfall, I guess I am wrong about the way of implementing my data. For each category, the code should be 0 and 1 and the only category the subject selected should be 1. and the logit function should be different, that is cumulative logit function. So I guess I have to redefine link function...??? – KH Kim Sep 5 '14 at 1:23

For anyone who might be curious about the same thing, clmm can fit GRM with slope 1... but I don't think clmm can allow different slopes for different items.

Here's the code.

library(ltm) # For data Science
library(ordinal)

# Science to Long format

Sci <- ltm::Science;
Sci$id <- 1:nrow(Sci) Sci.long <- reshape(Sci, varying = list(1:7), v.names = "resp", idvar="id", direction="long", timevar="item") Sci.long$item <- factor(Sci.long$item, labels=colnames(ltm::Science)) head(Sci.long) Sci.clmm <- clmm(resp ~ item + (1|id), data = Sci.long, link = "logit", threshold = "flexible", nAGQ=-31, control = list(method = "nlminb", trace = 1, # trace : 1=outer op, -1=inner op. maxIter = 1000, gradTol = 1e-5, maxLineIter = 1000, useMatrix = FALSE, innerCtrl = "giveError")) Sci.clmm2 <- clmm2(location = resp ~ item , data = Sci.long, random = factor(id), link = "logistic", threshold = "flexible", nAGQ=-31, control = clmm2.control(method = "nlminb", trace = 1, # trace : 1=outer op, -1=inner op. maxIter = 1000, gradTol = 1e-5, maxLineIter = 1000)) library(mirt) model = "F1 = 1 - 7 COV = F1 * F1 START = (1-7, a1, 1) FIXED = (1-7, a1) CONSTRAIN = (1-7, d1), (1-7, d2), (1-7, d3) " #as.numeric(as.character(ltm::Science)) Science2 <- sapply(ltm::Science, as.numeric) Sci.mirt <- mirt(Science2, mirt.model(model), itemtype="grsm", TOL=1e-5, quadpts = 31) # COMPARE The Result from mirt and clmm/clmm2 coef(Sci.mirt) Sci.clmm Sci.clmm2  So the results are identical, one thing to give heed to is that mixed model allows free random effect variance, but IRT model assumes variance of theta = 1. > coef(Sci.mirt)$Comfort
a1    d1    d2     d3 c
par  1 3.642 1.642 -0.888 0

$Environment a1 d1 d2 d3 c par 1 3.642 1.642 -0.888 -0.292 ...$GroupPars
MEAN_1 COV_11
par      0  0.679

> Sci.clmm
Cumulative Link Mixed Model fitted with the Gauss-Hermite

formula: resp ~ item + (1 | id)
data:    Sci.long

logit flexible  2744 -3072.83 6165.66 1119(3) 1.04e-03

Random effects:
Groups Name        Variance Std.Dev.
id     (Intercept) 0.6788   0.8239
Number of groups:  id 392

Coefficients:
itemEnvironment        itemWork      itemFuture  itemTechnology            itemIndustry     itemBenefit
-0.2922         -0.9476         -0.2901         -0.2036          0.4632         -0.6788

Thresholds:
strongly disagree|disagree             disagree|agree       agree|strongly     agree
-3.6419                    -1.6420                         0.8886
> Sci.clmm2
Cumulative Link Mixed Model fitted with the Gauss-Hermite

Call:
clmm2(location = resp ~ item, random = factor(id), data = Sci.long,
link = "logistic", control = clmm2.control(method = "nlminb",
trace = 1, maxIter = 1000, gradTol = 1e-05, maxLineIter = 1000),
nAGQ = -31, threshold = "flexible")

Random effects:
Var   Std.Dev
factor(id) 0.6788056 0.8238966

Location coefficients:
itemEnvironment        itemWork      itemFuture  itemTechnology        itemIndustry     itemBenefit
-0.2922101      -0.9476163      -0.2901400      -0.2035440           0.4631815      -0.6788248

No Scale coefficients

Threshold coefficients:
strongly disagree|disagree             disagree|agree       agree|strongly agree
-3.6419330                 -1.6420389                      0.8885574

log-likelihood: -3072.828
AIC: 6165.656


And one way to use glmer to fit GRM is mention as "An approximate ML estimation method for CLMs" in https://cran.r-project.org/web/packages/ordinal/vignettes/clm_intro.pdf