Latent Variable constraints and Graded response item effects

[Mark White]

Hello

I'm having some trouble identifying whether it is possible to do a few things with the mirt package.  I have a data set where raters scored teachers on a 1-7 scale.  One thing that I would like to do for each model is adjust for a rater effect on each item (starting with mean shift of thresholds).  With a binary outcome, I appear to be able to do this through providing an equation that interacts RaterID with items in the fixed parameter.  Is there a way to do this for polytomous items?  If I treat each rater as a separate item, is there a way to constrain the differences across parameters to be a fixed value?  My second question has to do with restraints on latent variables.  Is it possible to constrain the means or variances of latent variables to be equal to each other with their value estimated?

Thank you,

Mark

[Phil Chalmers]

Hi Mark,

On Fri, Sep 25, 2015 at 2:08 PM, Mark White <mark.c....@gmail.com> wrote:

Hello

I'm having some trouble identifying whether it is possible to do a few things with the mirt package.  I have a data set where raters scored teachers on a 1-7 scale.  One thing that I would like to do for each model is adjust for a rater effect on each item (starting with mean shift of thresholds).  With a binary outcome, I appear to be able to do this through providing an equation that interacts RaterID with items in the fixed parameter.  Is there a way to do this for polytomous items?  If I treat each rater as a separate item, is there a way to constrain the differences across parameters to be a fixed value?  

I'm having trouble conceptualizing this, so I can't recommend anything at the moment. Could you create a dummy dataset/source code to work with so I can see what you are trying to setup?

 

My second question has to do with restraints on latent variables.  Is it possible to constrain the means or variances of latent variables to be equal to each other with their value estimated?

Yes, but I haven't included syntax support for this yet. You can however do it by explicitly creating lists of constraints and passing them to the constrain argument. It requires that the parameters numbers are known, which you can obtain through the data.frame of starting values by passing pars = 'values'. There's some examples of this spattered throughout the package, especially in ?mirt and ?multipleGroup. HTH.

Phil

 

Thank you,

Mark

--
You received this message because you are subscribed to the Google Groups "mirt-package" group.
To unsubscribe from this group and stop receiving emails from it, send an email to mirt-package...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[Mark White]

Thank you for your prompt reply.  I see now the examples around constraining latent variances. 

I'm not the greatest at creating random data, but below should do the job using graded response model.  It does require arm package for invlogit function.  The most relevant part to my question is the r matrix which includes a rater "severity" effect for each of 5 raters that varies over items.  This leads to a different intercept for each unique rater by item combination.   The conceptual idea being that some raters are more severe than others and this severity is a function of the specific item.  

Mark

**************

require(arm)

set.seed(1)

N=500 # number cases; N/R unique persons (fully crossed person by rater)

R=5 # number raters

I=10 # number items

head(dat<-data.frame(id=rep(1:(N/R),each=R),rater=rep(1:R,N/R),theta=rep(rnorm(N/R),each = R))) # each person scored by all raters

a <- runif(I,.5,.8)

d1<- runif(I,-2,-1.5)

d2<- d1+runif(I,1,1.4)

d3<- d2+runif(I,1,1.4)

# Generate Random rater effect for each rater (unique for each item)

r <-matrix(rnorm(R*I,0,0.1),R,I)

for (i in 1:I) {

  # Generate cumulative probabilities of graded response model (note intercept due to rater and item)

  dat[,paste0("CumProb",i,".1")] <- with(dat,invlogit(a[i]*theta-d1[i]-r[rater,i]))

  dat[,paste0("CumProb",i,".2")] <- with(dat,invlogit(a[i]*theta-d2[i]-r[rater,i]))

  dat[,paste0("CumProb",i,".3")] <- with(dat,invlogit(a[i]*theta-d3[i]-r[rater,i]))

  # Generate Break Points for in each category by adding width to previous value

  dat[,paste0("Prob",i,".0")] <- 1             -dat$CumProb1.1

  dat[,paste0("Prob",i,".1")] <- dat[,paste0("Prob",i,".0")]+(dat$CumProb1.1-dat$CumProb1.2)

  dat[,paste0("Prob",i,".2")] <- dat[,paste0("Prob",i,".1")]+(dat$CumProb1.2-dat$CumProb1.3)

  # Random variable to select category

  dat$R01 <- runif(N,0,1)

  # Create final Y value

  dat[,paste0("Y",i)] <- ifelse(dat$R01<=dat[,paste0("Prob",i,".0")],0,

                         ifelse(dat$R01<=dat[,paste0("Prob",i,".1")],1,

                         ifelse(dat$R01<=dat[,paste0("Prob",i,".2")],2,3)))

                                                     

}

dat<-dat[,c("id","rater","theta","Y1","Y2","Y3","Y4","Y5","Y6","Y7","Y8","Y9","Y10")]

head(dat)

summary(sapply(dat,factor))

[Phil Chalmers]

Perfect, I see what you are saying now. Unfortunately there doesn't seem to be an easy way to do this in the package, mostly because mixedmirt() doesn't support rating scale items for polytomous data (yet). It's on the todo list, but it's implementation is rather cumbersome. 

You could approach the problem slightly differently though through the mirt() function, and it would just require that the data be reorganized into the wide form (where each subject is on only one row). Continuing with your example code, here are three models which increase in the number of constraints imposed on the model (the last one seems to be what you are looking for):

library(reshape2)

mlt <- melt(dat, id.vars = c('id', 'rater', 'theta'))

wide <- dcast(mlt, id + theta ~ rater + variable)

items <- wide[,-c(1:2)]

apply(items, 2, function(x) length(unique(x)))

# extremely loose rating scale model for each item

mod <- mirt(items, 1, 'grsm') 

coef(mod, simplify=TRUE)

# more constrained, accounting for specific rating effects

model <- 'F=1-50

          FIXED = (1-10, c)

          START = (1-10, c, 0.0)

          CONSTRAIN = (11-20, c), (21-30, c), (31-40, c), (41-50, c)'

mod2 <- mirt(items, model, 'grsm')

coef(mod2, simplify=TRUE)

# even more constrained (equal slopes)

# cat(rowSums(expand.grid(c(1,11,21,31,41), 0:9)), sep=',')

model2 <- "F=1-50

           FIXED = (1-10, c)

           START = (1-10, c, 0.0)

           CONSTRAIN = (11-20, c), (21-30, c), (31-40, c), (41-50, c),

             (1,11,21,31,41,a1),(2,12,22,32,42,a1),(3,13,23,33,43,a1),(4,14,24,34,44,a1),

             (5,15,25,35,45,a1),(6,16,26,36,46,a1),(7,17,27,37,47,a1),(8,18,28,38,48,a1),

             (9,19,29,39,49,a1),(10,20,30,40,50,a1)"

                    

mod3 <- mirt(items, model2, 'grsm')

coef(mod3, simplify=TRUE)

anova(mod, mod2)

anova(mod2, mod3)

anova(mod, mod3)

Cheers.

Phil

[Mark White]

This is a useful way to look at increasing levels of restrictions, thank you.  The grsm model though is a bit more restrictive than I was hoping to use.  Namely the d1-d3 parameters are constant across all the items.  Is there a way to remove this restriction across items, but retain it across raters?  I see the grsm.block argument, which, based on my reading of the documentation, creates blocks of variables that are restricted to have the same d1-d3 values, but it does not seem to be functional. 

mod <- mirt(items, 1, 'grsm') 

mod4 <- mirt(items, 1, 'grsm',grsm.block=rep(1:10,5)) # set all Y1 variables in block 1, Y2 in block 2....

anova(mod,mod4)

# anova shows the models are equivalent

[Phil Chalmers]

Yes this should work, but unfortunately you've stumbled upon a bug that was recently patched. If you install the dev version on Github this will function correctly.

Phil

[Mark White]

Thanks so much!

[KK Sasa]

Hi Phil,

I found this post is useful to model the rater effect. However, it seems, taking a look at manual, only grsm has the single 'c' parameter, not in gpcm/nominal model. I prefer PCM because it belongs to Rasch model family. The marginal scores of the data matrix would be sufficient statistiscs for person, item, and rater. Would it be possible to do the same analysis on PCM? Many thanks!

[Phil Chalmers]

On Fri, Dec 18, 2015 at 4:44 AM, KK Sasa <genw...@gmail.com> wrote:

Hi Phil,

I found this post is useful to model the rater effect. However, it seems, taking a look at manual, only grsm has the single 'c' parameter, not in gpcm/nominal model. I prefer gpcm model because it belongs to Rasch model family.

Yes, and this is the code which is commented out on the dev version. Pull requests are welcome as I'm busy with many other things at this time of year.

 

The marginal scores of the data matrix would be sufficient statistiscs for person, item, and rater. Would it be possible to do the same analysis on gpcm model? Many thanks!

Yes of course, the reasoning is exactly the same. Though the 'sufficient statistic' estimation is something that only relates to CML, which I don't care much for anyway and which isn't supported in mirt (likely never will be). MML and CML are essentially equivalent in all respects (after all, they are just different ways to achieve ML estimation), and I think the literature has largely blown the whole sufficient statistics aspect out of proportion. There really is no benefit of CML over MML, and in fact MML has some nicer properties than CML in larger tests due to the divide-and-conquer nature of the EM.

Phil

[KK Sasa]

Hi Phil,

Thank you for your prompt reply.

I tried a grsm with nominal model (2-dimensional model). 1-25 items follow grsm and 21-25 items follow nominal model. I tried the following syntax:

load("test")
cmod <- mirt(s$data, model=s$cmodel, pars=s$sv, itemtype=s$itemtype,SE=TRUE,SE.type='complete')

Error: Constraint applied to fixed parameter(s) 85 91 97 103 109  but should only be applied to
                 estimated parameters. Please fix!

Those parameters seem being "d0" of nominal item which shall be fixed to 0... The errors looks odd.


Phil Chalmers於 2015年12月18日星期五 UTC+8下午11時34分46秒寫道：

[Phil Chalmers]

Looks like this was a bug, should be patched now (https://github.com/philchalmers/mirt/commit/a70d2b2e19e02e141bca55ded9f80762f7ed46e2). This model is somewhat pointless though, given that all the items have dichotomous responses and therefore nominal/rating-scale models are not at all needed....I'm not sure if that was deliberate or not, but just stating it in case anyone else were to stumble upon this post. Cheers.

Phil

[KK Sasa]

Thank you, Phil.

Actually, I found this bug as fitting my original model. For testing, I reduced it to a simpler model and also found the bug. Following is my original model. I constrain the 'a1' and 'a2' of nominal model to be summed to 1 in each nominal item (item 21 to 25). And, item 1 to 10 for rater1, item11 to 20 for rater2. So the threshold 'd1' are constrained to be equal (e.g., item1 and item11, item2 and item12, item3 and item 13 ...). But, I got error again:

~~~
require("Rsolnp")
eqfun <- function(para, optim_args) {   
    a1 = para[31] + para[32]
    a2 = para[34] + para[35]
    a3 = para[37] + para[38]
    a4 = para[40] + para[41]
    a5 = para[43] + para[44]
    d1 = para[1] - para[11]
    d2 = para[2] - para[13]; d3 = para[3] - para[15]
    d4 = para[4] - para[17]; d5 = para[5] - para[19]
    d6 = para[6] - para[21]; d7 = para[7] - para[23]
    d8 = para[8] - para[25]; d9 = para[9] - para[27]
    d10 = para[10] - para[29];
    return(c(a1,a2,a3,a4,a5,d1,d2,d3,d4,d5,d6,d7,d8,d9,d10))
}

aa = s$sv[s$sv$est, ]
solnp_args <- list(eqfun=eqfun, eqB=c(rep(1,5),rep(0,10)),LB= aa$lbound[1:(dim(aa)[1]-2)], control=list(rho=1))

cmod <- mirt(s$data, model=s$cmodel, pars=s$sv, itemtype=s$itemtype,SE=TRUE,SE.type='complete',optimizer = 'solnp', solnp_args=solnp_args)
~~~
Error: Error :
solnp-->error: LB length not equal to parameter length



Phil Chalmers於 2015年12月22日星期二 UTC+8上午11時00分46秒寫道：

[Phil Chalmers]

I see. This is a different kind of error that's slightly harder to solve. Currently the parameter vector is actually smaller than you think because linear constraints are being tracked internally for the grsm types (hence the error). So, considerably fewer unique parameters are being tracked, all redundant intercepts in this case.

The easiest solution would be to remove all constraints when these special optimizers are used, however that means that users would have to supply the rating scale model constraints as well, which tend to be cumbersome and are largely automatic (and could cause optimizers like this to choke due to the higher dimensional search....though don't quote me on that). 

Perhaps I'll add a technical argument to disable all internal constraints so that it's more clear what is going on. There's already and option like this for mixedmirt(), so might as well make it available elsewhere. Cheers.

Phil

[KK Sasa]

Thank you for pointing out of this, Phil. The mixedmirt is more flexible, but seems having an another issue of not allowing solnp :)

cmod <- mixedmirt(s$data, model=s$cmodel, pars=s$sv, itemtype=s$itemtype,SE=TRUE,SE.type='complete',internal_constraints=FALSE,optimizer = 'solnp', solnp_args=solnp_args)

Error: solnp only supported for optimization with EM algorithm

as well as

Error: alabama only supported for optimization with EM algorithm

Phil Chalmers於 2015年12月22日星期二 UTC+8下午12時24分03秒寫道：

[Phil Chalmers]

It's now possible to disable the internal constraints for rating scale models by passing technical = list(internal_constraints = FALSE). Cheers.

Phil