2
$\begingroup$

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?

Improve this question
$\endgroup$
6
  • 1
    $\begingroup$ 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. $\endgroup$
    – Momo
    Commented Sep 4, 2014 at 16:30
  • 2
    $\begingroup$ 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. $\endgroup$
    – Jake Westfall
    Commented Sep 4, 2014 at 16:56
  • $\begingroup$ @Momo, I guess you're right. Sorry for not mentioning. I thought of GRM with constrained slope parameter. $\endgroup$
    – KH Kim
    Commented Sep 5, 2014 at 1:14
  • $\begingroup$ @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 $\endgroup$
    – KH Kim
    Commented Sep 5, 2014 at 1:15
  • $\begingroup$ @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...??? $\endgroup$
    – KH Kim
    Commented Sep 5, 2014 at 1:23

1 Answer 1

Reset to default
1
$\begingroup$

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

Improve this answer
$\endgroup$

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.