Dear Phil,
I have a question concerning the simulation of data according the graded response model and the estimation of the M2 statistic for testing the model fit.
I performed 10 000 repetitions (looping with for) of the following calculation:
#-------------------------------------------------------------------------------------------------------------------------------------------------------
# 1. Simulate data according to the graded response model
library(mirt)
ol <- matrix(rlnorm(15, .2, .3))
ci <- matrix(rnorm(15*3, 0, 2), 15)
ci <- t(apply(ci, 1, sort, decreasing=TRUE))
poly.grm <- simdata(ol, ci, 200, itemtype = 'graded')
# 2. dichotomization of the data [1 & 2 = 0] & [3 & 4 = 1]
dich.grm <- ifelse(poly.grm < 3,0,1)
# 3. estimation of the graded response model for the polytomous and dichotomized data
est.grm.poly <- mirt(poly.grm, model=1, itemtype="graded")
est.grm.dich <- mirt(dich.grm, model=1, itemtype="graded")
# 4. calculation of the M2 statistic for testing the model fit
# p-values which are smaller than alpha are coded as 1, else 0
alpha <- 0.05
m.grm.poly <- M2(est.grm.poly)
p.grm.poly <- m.grm.poly[,3]
alpha.grm.poly <- ifelse(p.grm.poly < alpha,1,0)
m.grm.dich <- M2(est.grm.dich)
p.grm.dich <- m.grm.dich[,3]
alpha.grm.dich <- ifelse(p.grm.dich < alpha,1,0)
#-------------------------------------------------------------------------------------------------------------------------------------------------------
The results were surprisingly. The M2 statistic (for testing the model fit) had significant p-values in 2129 from 10000 (0.2129) cases when calculating the polytomous model,
but only 570 from 10000 (0.057) significant p-values when calculating the dichotomous model.
How can that be? Have I made a mistake in the calculation?
Thank you very much for your Help.
Best regards Isabell
--
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.
Many thanks for the prompt reply.
You are right, I suppressed the warnings in my simulation. So I didn't get the 'EM cycles terminated after 500 iterations' -warning (thanks for the hint).
Would it be better to simulate the data as follows:
ol <- matrix(rlnorm(15, .2, .3))
diffs <- t(apply(matrix(runif(40*4, .3, 1), 40), 1, cumsum))
diffs <- -(diffs - rowMeans(diffs))
d <- diffs + rnorm(40)
poly.grm <- simdata(ol, d, 400, itemtype = 'graded')
and perhaps also to increase the sampling size to 400 or maybe even 800?
One final question i would like to ask is: Are there any other functions for testing the modelfit of the Graded Response Model?
Best regards I.B
PS Your packages is great.
Many thanks for the prompt reply.
You are right, I suppressed the warnings in my simulation. So I didn't get the 'EM cycles terminated after 500 iterations' -warning (thanks for the hint).
Would it be better to simulate the data as follows:
ol <- matrix(rlnorm(15, .2, .3))
diffs <- t(apply(matrix(runif(40*4, .3, 1), 40), 1, cumsum))
diffs <- -(diffs - rowMeans(diffs))
d <- diffs + rnorm(40)
poly.grm <- simdata(ol, d, 400, itemtype = 'graded')
and perhaps also to increase the sampling size to 400 or maybe even 800?
One final question i would like to ask is: Are there any other functions for testing the modelfit of the Graded Response Model?
Best regards I.B
PS Your packages is great.
On Sat, Mar 21, 2015 at 1:09 PM, ibaldauf <isabell...@gmail.com> wrote:Many thanks for the prompt reply.
You are right, I suppressed the warnings in my simulation. So I didn't get the 'EM cycles terminated after 500 iterations' -warning (thanks for the hint).
Would it be better to simulate the data as follows:
ol <- matrix(rlnorm(15, .2, .3))
diffs <- t(apply(matrix(runif(40*4, .3, 1), 40), 1, cumsum))
diffs <- -(diffs - rowMeans(diffs))
d <- diffs + rnorm(40)
poly.grm <- simdata(ol, d, 400, itemtype = 'graded')
and perhaps also to increase the sampling size to 400 or maybe even 800?
Yes the sample size will also help. It's somewhat of an unwritten rule that you should have at least 5 * #parameters when fitting your model, but of course it depends on a lot of other factors in model to make that definitive. Your last attempt was a little on the lower side though, with 15 * 5 = 75 parameters, with out N = 200.One final question i would like to ask is: Are there any other functions for testing the modelfit of the Graded Response Model?
There also is one based on collapsing factor scores to sum scores, which you can get with fscores(mod, method = 'EAPsum'). They are just the usual X2 and G2 statistics, which can be calculated from the observed and expected returned scores (I don't return them in the function, but probably should for simulation reasons). Cheers.