[R-sig-ME] Logistic modelling with guessing parameter

[Ben Bolker]

-----BEGIN PGP SIGNED MESSAGE-----
Hash: SHA1

  PS I had a closer look and most of your individual groups seem to
show complete separation ... this is generally going to be a problem ...

making audio_name a *random* effect seems to work OK.  You could also
look into doing some shrinkage/penalization via blmer ...


library(plyr)
sumvals <- ddply(data,c("audio_name","accuracy"),
      summarise,
      responses=mean(responses))

library(ggplot2); theme_set(theme_bw())
ggplot(data,aes(accuracy,responses,colour=audio_name))+
    stat_sum()+geom_smooth(method="glm",
                           family = binomial(mafc.logit(2)),
                           se=FALSE)+
        geom_point(data=sumvals,alpha=0.25)+
        geom_line(data=sumvals,alpha=0.25)
m1 <- lme4::lmList(response~accuracy|audio_name,
                   family = binomial(mafc.logit(2)),
                   data=data)

debug(lme4::lmList)
ff2 <- function(dat, formula=responses~accuracy) {
    data <- as.data.frame(dat)
    glm(formula=formula, family=binomial(mafc.logit(2)), data)
}
ss <- split(data,data$audio_name)
## 4: OK
sapply(ss,function(d) {
           !any(coef(ff2(d))>1e6)
       })
ff2(ss[[6]])
mod1 <- glmer(responses ~ accuracy + (1|p_ID),
              family = binomial(mafc.logit(2)), data=data)
mod2 <- update(mod1,. ~. + audio_name)
mod3 <- update(mod1,. ~. + (1|audio_name))
summary(mod1)



On 15-05-07 01:13 PM, Peter Harrison wrote:
> Hello everyone,
> 
> I'm currently modelling response data for a musical error detection
> experiment. Participants' responses can either be correct (1) or
> incorrect (0); the difficulty of the question is predicted by a
> number of factors, including accuracy (accurate performance ->
> errors are difficult to detect) and musical piece ("audio_name").
> 
> The model I'm fitting is similar to a 3-PL IRT model, with a
> constrained guessing parameter of 0.5 (i.e. the most difficult
> questions should be answered with a 50% success rate). However,
> instead of the normal conception of item difficulty, probability of
> correct response is given by a linear combination of
> difficulty-related predictors (such as accuracy and audio_name).
> Therefore:
> 
> probability of correct response = 0.5 + 0.5 / (1 +
> exp(predictors))
> 
> where predictors = e.g. a*accuracy + b*audio_name + c (but with
> audio_name dummy coded, with different b for each level of
> audio_name).
> 
> I've been using a script kindly posted by Ken Knoblauch back in
> 2010
> (https://stat.ethz.ch/pipermail/r-sig-mixed-models/2010q4/004531.html)
> to achieve this using lme4, simply changing a few terms to keep it
> up-to-date with current package versions (see bottom of message for
> the script). It seems to work fine with "accuracy" as a fixed
> effect and "p_ID" as a random effect. However, I want to add
> "audio_name" as a fixed effect, i.e. responses ~ accuracy +
> audio_name + (1|p_ID). But when I do this, I get the following
> error message: (maxstephalfit) PIRLS step-halvings failed to reduce
> deviance in pwrssUpdate.
> 
> I'm not sure what this error message means, but have been assuming
> it is talking about a convergence problem. I've tried various
> things to fix it, including
> control=glmerControl(optimizer="bobyqa"),
> control=glmerControl(optimizer="Nelder_Mead"), and nAgQ=3, but none
> of these have worked.
> 
> Would anyone be able to help me with this error message, and what
> to do about it?
> 
> Thanks so much! Peter
> 
> 
> Key to data file:
> 
> accuracy = accuracy of the music clip (higher accuracy -> question
> is more difficult) = continuous variable audio_name = musical piece
> played = categorical variable (27 levels, coded with text strings) 
> p_ID = participant ID = categorical variable (229 levels, coded
> 1:229) behind_beat = whether clip was behind the beat or not =
> dichotomous (1 = TRUE, 0 = FALSE)
> 
> Script:
> 
> library(lme4) data <-
> read.csv("https://dl.dropboxusercontent.com/s/szq2e2sxwfsuxo6/data.csv",
> header = TRUE)
> 
> mafc.logit <- function (.m = 2) { .m <- as.integer(.m) if (.m < 2) 
> stop(".m must be an integer > 1") linkfun <- function(mu) { mu <-
> pmax(mu, 1/.m + .Machine$double.eps) qlogis((.m * mu - 1)/(.m -
> 1)) } linkinv <- function(eta) { 1/.m + (.m - 1)/.m *
> binomial()$linkinv(eta) } mu.eta <- function(eta) ((.m - 1)/.m) *
> binomial()$mu.eta(eta) valideta <- function(eta) TRUE link <-
> paste("mafc.logit(", .m, ")", sep = "") structure(list(linkfun =
> linkfun, linkinv = linkinv, mu.eta = mu.eta, valideta = valideta,
> name = link), class = "link-glm") }
> 
> mod1 <- glmer(responses ~ accuracy + (1|p_ID), family =
> binomial(mafc.logit(2)), data=data) summary(mod1)
> 
> [[alternative HTML version deleted]]
> 
> _______________________________________________ 
> R-sig-mixed-models at r-project.org mailing list 
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
> 

-----BEGIN PGP SIGNATURE-----
Version: GnuPG v1.4.11 (GNU/Linux)

iQEcBAEBAgAGBQJVS78SAAoJEOCV5YRblxUHt2QIAK5g3BvoQxKO3RLO9wwOFiU8
eIz77IpkDDMt+uec6dLp30I1/HlwVUuW2zlpK3mkEy4kAWF7YU8m4BtztvvLrDPX
mHfSH6oEgCqQHRbQwvH1FyqKt2Th7JW3pdfAsDfNKa0ctIwnK7qNmJJSKEG/jWtm
Z2ERNit8U2kHY0eEPhBwjPMBjxRqsR2gBncckqq0KntWw8l8QqBFJM3+qhPzLPB0
ftkA2tJfDkIZ0dqPTHqdeGyzW1IK84Pr8HObDiRBUMOgmS0UNhBXQWYutjCz6pwE
zuFzQoXvtvVBVVVPTgO8lOacX4IEgydYFncKh3L4wM1AdmJx73Wg8wi4eppXtuI=
=vRQM
-----END PGP SIGNATURE-----