Solved – How to fit Graded Response Model with lme4::glmer

glmmitem-response-theoryr

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).

=== ADDED

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.)

for 1-5 category(ordinal) answer.

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?

Best Answer

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
head(ltm::Science)

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 
quadrature approximation with 31 quadrature points 

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

 link  threshold nobs logLik   AIC     niter   max.grad
 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 
quadrature approximation with 31 quadrature points 

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

Related Solutions

Longitudinal Item Response – How to Apply Longitudinal Item Response Theory Models in R

As a precursor, the IRT approach to this problem is very demanding computationally due to the higher dimensionality. It may be worthwhile to look into structural equation modeling (SEM) alternatives using the WLSMV estimator for ordinal data since I imagine less issues will exist. Plus, including external covariates is much easier within that framework. Both approaches I describe here are also possible in SEM.

There are two ways that I know of which you can estimate unidimensional longitudinal IRT models that are not Rasch in nature. The the first approach requires a unique latent factor for each time block and a specific residual variation term for each item. A different approach, similar to what one would find in the SEM literature, is via a latent growth curve model whereby only a fixed number of factors are estimated (three if the relationship over time is believed to be linear). Fixed loadings are used in this approach, so computationally it may be much more stable due to the reduced number of estimated parameters, so I would tend to prefer the growth curve model for both the smaller dimensionality and fewer estimated parameters.

The idea for both approaches is to set up latent time factors indicating how person level $\theta$ values change over each test administration, and constrain the influence of their loadings across time as well so that their hyper parameters can be estimated (i.e., the latent mean and covariances). Item constraints must also be imposed across time to remain invariable so that the person differences are only captured in the hyper parameters. Since this approach can require a huge number of integration dimensions, so you'll need to use something like the dimensional reduction algorithm which is available in mirt under the bfactor() function.

Instead of going through a worked example here, which would take a lot of code, I'll simply point to a worked versions of these analyses. A word of warning though, these are very computationally demanding and may take more than an hour to converge on your computer since you have 4 dimensions of integration in the first case and 3 dimensions in the second. Or, if you don't have much RAM you could experience issues when increase the number of quadpts.

Data simulation script: https://github.com/philchalmers/mirt/blob/gh-pages/data-scripts/Longitudinal-IRT.R

Analysis output: http://philchalmers.github.io/mirt/html/Longitudinal-IRT.html

In the first example, if you save the factor scores by using fscores() you'll obtain estimates for each time point regarding how individual $\theta$ values are changing. In the second example, using the linear growth curve approach, the first column of the factor scores will represent the initial $\theta$ estimates while the second column will indicate the slope/change occurring on average over time. In the example, I set up a constant mean change of .5, so the values in fscores() should all be around 0.5 for each individual. Both analyses give roughly the same conclusions but are somewhat different approaches to the problem. However, if you are familiar with longitudinal models in SEM then these should be fairly natural to interpret.

Solved – Implementing nominal model from Polytomous Item Response Theory

The $a$ and $b$ parameters likely don't function the way you would think. Specifically, the $a$s are not really discrimination parameters in the sense of the usual IRT dichotomous/graded response models, but instead represent the parameter ordering. Under Bock's (1972) original formulation, both sets of parameters were constrained to sum to 1 for proper identification and interpretation, so some of the $a$s will necessarily have to be negative as well as the $b$s. Lower values of $a$ represent lower categories (where the rank of the values indicate the empirical ordering), while lower values of $b$ represent the relative 'easiness' of the category to select.

What you've presented here is two response category that have exactly the same ordering (neither one is theoretically higher than the other) since both the $a$ and $b$ are the same. Try changing the value of the parameters to see what happens, but try to do so under the $\sum_k a_k = 1$ and $\sum_k b_k = 1$ constraints (not required for plots of course, but good to think about to see how real data would behave).

For a quick interactive graphical representation already available, you could always check out the shiny interface that ships with the mirt package. Simply put the following in R:

library(mirt)
itemplot(shiny=TRUE)

Then select the nominal itemtype from the drop down menu, and click the checkbox to use the traditional IRT metric. Play around after that with the parameter values sliders to see what happens to the output plots. Hope that helps.

References

Bock, R. D. Estimating item parameters and latent ability when the responses are scored in two or more nominal categories Psychometrika, 1972, 37, 29-51.

Best Answer

Related Solutions

Longitudinal Item Response – How to Apply Longitudinal Item Response Theory Models in R

Solved – Implementing nominal model from Polytomous Item Response Theory

Related Question