Multiple nonlinear constraints

[KK Sasa]

Dear Phil,

1.
I run an example in mirt manual but got error:

dat <- expand.table(LSAT6)
model <- mirt.model('
Theta = 1-5
MEAN = Theta
COV = Theta*Theta')
sv <- mirt(dat, model, itemtype = 'Rasch', pars = 'values')
eqfun <- function(p, optim_args) {
    sum(p[1:5])
     }
solnp_args <- list(eqfun=eqfun, eqB=0)
mod <- mirt(dat, model, itemtype = 'Rasch', optimizer = 'solnp', solnp_args=solnp_args)

> Error: expecting a single value

2. And, I still cannot understand why the "p" in eqfun would automatically be equal to the item difficulty? Didn't see the relation.

3. At the same time, what if I want to add more constraints, such as the the "prod(slope)=1", could you give me an example to do this?

Many thanks.

[Phil Chalmers]

Hi KK,

On Mon, Apr 6, 2015 at 5:55 AM, KK Sasa <genw...@gmail.com> wrote:

Dear Phil,

1.
I run an example in mirt manual but got error:

dat <- expand.table(LSAT6)
model <- mirt.model('
Theta = 1-5
MEAN = Theta
COV = Theta*Theta')
sv <- mirt(dat, model, itemtype = 'Rasch', pars = 'values')
eqfun <- function(p, optim_args) {
    sum(p[1:5])
     }
solnp_args <- list(eqfun=eqfun, eqB=0)
mod <- mirt(dat, model, itemtype = 'Rasch', optimizer = 'solnp', solnp_args=solnp_args)

> Error: expecting a single value

The code runs fine to me (e.g., http://philchalmers.github.io/mirt/html/mirt.html). What is your sessionInfo()? Perhaps some packages are out of date.

 

2. And, I still cannot understand why the "p" in eqfun would automatically be equal to the item difficulty? Didn't see the relation.

Pass pars = 'values' to obtain the starting values, and index the rows where est = TRUE. The values in this reduced data.frame forms the values that make up the vector that is p, so the first 5 values in this case is equal to the item difficulty. If you used the 2PL model instead, it would be every other free parameter, because slopes would be estimated too.

 

3. At the same time, what if I want to add more constraints, such as the the "prod(slope)=1", could you give me an example to do this?

```

dat <- expand.table(LSAT6)

model <- mirt.model('

Theta = 1-5

COV = Theta*Theta')

sv <- mirt(dat, model, itemtype = '2PL', pars = 'values')

sv[sv$est,]

eqfun <- function(p, optim_args) {

    prod(p[seq(1,10,2)])

}

solnp_args <- list(eqfun=eqfun, eqB=1)

mod <- mirt(dat, model, itemtype = '2PL', optimizer = 'solnp', solnp_args=solnp_args)

coef(mod, simplify=TRUE)

prod(coef(mod, simplify=TRUE)$items[,1])

```

Phil

 

Many thanks.

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

[KK Sasa]

My infor is as follows.

> sessionInfo()
R version 3.0.3 (2014-03-06)
Platform: i386-w64-mingw32/i386 (32-bit)

locale:
[1] LC_COLLATE=English_United States.1252
[2] LC_CTYPE=English_United States.1252  
[3] LC_MONETARY=English_United States.1252
[4] LC_NUMERIC=C                         
[5] LC_TIME=English_United States.1252   

attached base packages:
[1] stats4    stats     graphics  grDevices utils     datasets  methods 
[8] base    

other attached packages:
[1] mirt_1.9         lattice_0.20-27  rpart_4.1-5    

loaded via a namespace (and not attached):
 [1] codetools_0.2-8       foreach_1.4.2         GPArotation_2014.11-1
 [4] grid_3.0.3            iterators_1.0.7       pkgXMLBuilder_1.0   
 [7] Rcpp_0.11.4           revoIpe_1.0           sfsmisc_1.0-27      
[10] tools_3.0.3           XML_3.98-1.1        
>

Phil Chalmers於 2015年4月6日星期一 UTC+8下午7時13分10秒寫道：

[Phil Chalmers]

Have you installed the Rsolnp package? mirt doesn't require installing it as a dependency (it's more a plug-in). http://cran.r-project.org/web/packages/Rsolnp/index.html

Phil 

[KK Sasa]

Ohoh, it's my mistake missing it. It works now. Thank you, Phil. Now I understand the idea behind the mirt scene.

Phil Chalmers於 2015年4月6日星期一 UTC+8下午7時32分34秒寫道：

[KK Sasa]

Hi Phil,

Here comes another problem as follows.

Phil Chalmers於 2015年4月6日星期一 UTC+8下午7時13分10秒寫道：

3. At the same time, what if I want to add more constraints, such as the the "prod(slope)=1", could you give me an example to do this?

```

dat <- expand.table(LSAT6)

model <- mirt.model('

Theta = 1-5

COV = Theta*Theta')

sv <- mirt(dat, model, itemtype = '2PL', pars = 'values')

sv[sv$est,]

eqfun <- function(p, optim_args) {

    prod(p[seq(1,10,2)])

}

solnp_args <- list(eqfun=eqfun, eqB=1)

mod <- mirt(dat, model, itemtype = '2PL', optimizer = 'solnp', solnp_args=solnp_args)

coef(mod, simplify=TRUE)

prod(coef(mod, simplify=TRUE)$items[,1])

I try whether the variance could be constrained, but it seems failed. I modifed

eqfun <- function(p, optim_args) {
    z1 = prod(p[seq(1,10,2)])
    z2 = p[11]
    return(c(z1,z2))
}
solnp_args <- list(eqfun=eqfun, eqB=c(1,1))

Error in Mstep(pars = pars, est = est, longpars = longpars, ngroups = ngroups,  :
  Error in svd(x) : infinite or missing values in 'x'

Seems it doesn't work for GroupPars parameters?

 

[Phil Chalmers]

That's correct, GroupPars is not included in the constraint optimization (it's pretty rare to constrain these anyway...and what you are doing here is pretty unconventional. Just set the group pars to a constant in the syntax rather than this function). 

It may be helpful for you to put a browser() in your eqfun() to get an idea of what's going on:

eqfun <- function(p, optim_args) {

    browser()

    prod(p[seq(1,10,2)])

}

Passing that, you'll go into R's debugging mode, where you can simply type p in the console to see what you are working with in the optimizer. Use Q to exit the browser when you are done inspecting.

Phil

[KK Sasa]

Hi Phil,

Actually, I will compare a IRT model with its another equivalent model. For the equivalent model, the GroupPars are required to be constrained nonlinear with the slope parameters. I would like to test if they could get similar estimates numerically. For this, might I know would you plan to include GroupPars to optimization? Just a suggestion. Thank you.

Phil Chalmers於 2015年4月6日星期一 UTC+8下午9時43分32秒寫道：

[Phil Chalmers]

On Mon, Apr 6, 2015 at 10:06 AM, KK Sasa <genw...@gmail.com> wrote:

Hi Phil,

Actually, I will compare a IRT model with its another equivalent model. For the equivalent model, the GroupPars are required to be constrained nonlinear with the slope parameters. I would like to test if they could get similar estimates numerically. For this, might I know would you plan to include GroupPars to optimization? Just a suggestion. Thank you.

That would require changing quite a bit around internally, so it's not likely. What kind of parameterization are you trying to achieve exactly? 

[KK Sasa]

A two dimensional model with latent traits theta and gamma. Constrain their correlation and var(gamma) as

var(gamma) = 1 - cor*cor,

and

theta's slope = cor * gamma's slope (on gamma's scale)

[Phil Chalmers]

I see, yes that definitely isn't possible right now. I'll consider adding a more flexible objective function to be optimized sometime to make this possible. Is there any urgency for this or are you just playing around at the moment?

Phil

[KK Sasa]

Not in urgency. I don't like to be pushy. Please take your time. Many thanks for creating so flexible package, Phil.

Phil Chalmers於 2015年4月7日星期二 UTC+8下午10時12分01秒寫道：

[Phil Chalmers]

Just as an update on this thread: it's now possible to include the hyper-parameters in the objective function directly. So these non-linear constraints should be possible. Cheers.

Phil