[R-sig-ME] Logistic and nonlinear mixed models: Accounting for guessing probability

[Ken Knoblauch]

Robert Miller <miller at ...> writes:

> 
> Hello everyone,
> 
> Recently i tried to predict the discrimination probability of a chemosignal
> by its concentration and an experimental manipulation factor (term:
> concentration*x + test*b + concentration*test*c + d) with nested factor
> "manipulation" within "participants". For statistical analysis i needed to
> incorporate a fixed guessing probability into my model (similiar to a 3-PL
> IRT model) resulting in the following equation:
> 
> P(correct) = 0.33 + 0.67*(exp(term)/(1 + exp(term)))
> 
> As i found no way to do so via the glmer()-function of the lme4-package, i
> tried to use nlmer() but unfortunately even the simplest analysis with just
> the concentration factor and intercept resulted in cryptic error messages.
> 
> Syntax:
> library(lme4)
> rawdata <- read.csv2("http://dl.dropbox.com/u/7147679/AND_data.csv")
> 
> mod1 <- glmer(Correct ~ log(Concentrat) * Test + (Test|Code), family =
> binomial, data=rawdata) #works fine but is inappropriate

Here is a temporary solution that Doug Bates encouraged me to share
with you to get you going, and until we come up with a more elegant/safer
method.  It requires the lme4a package as suggested and a modified link
from the psyphy package that I'll include below.  In essence, we solve
the problem of the visibility of m by making it global, which is not the
most elegant or safest of solutions, but it seems to work for now.

library(lme4a)

mafc.logit <- function (.m = 2) 
{
    .m <<- as.integer(.m)  # .m goes global
    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 * .Call("logit_linkinv", eta, PACKAGE = "stats")
    }
    mu.eta <- function(eta) ((.m - 1)/.m) * .Call("logit_mu_eta", 
        eta, PACKAGE = "stats")
    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")
}

rawdata <- read.csv2("http://dl.dropbox.com/u/7147679/AND_data.csv")
mod1 <- glmer(Correct ~ log(Concentrat) * Test + (Test|Code), family =
     binomial(mafc.logit(3)), data=rawdata)
summary(mod1)
Generalized linear mixed model fit by maximum likelihood ['summary.mer']
 Family: binomial 
Formula: Correct ~ log(Concentrat) * Test + (Test | Code) 
   Data: rawdata 
      AIC       BIC    logLik  deviance 
 277.0939  301.4584 -131.5469  263.0939 

Random effects:
 Groups Name        Variance Std.Dev. Corr 
 Code   (Intercept) 112.6    10.61         
        Test        103.0    10.15    0.994
Number of obs: 240, groups: Code, 10

Fixed effects:
                     Estimate Std. Error z value
(Intercept)           -26.942     10.715  -2.514
log(Concentrat)         5.041      1.876   2.687
Test                   -1.918      6.931  -0.277
log(Concentrat):Test    2.136      1.304   1.638

Correlation of Fixed Effects:
            (Intr) lg(Cn) Test  
lg(Cncntrt) -0.934              
Test        -0.136  0.269       
lg(Cncnt):T -0.087  0.097 -0.739

range(fitted(mod1))
[1] 0.3333333 1.0000000

HTH

Ken

> Right now I'm quite desperate and would appreciate any help.
> Thank you
> Robert Miller
> 
> 

-- 
Ken Knoblauch
Inserm U846
Stem-cell and Brain Research Institute
Department of Integrative Neurosciences
18 avenue du Doyen Lépine
69500 Bron
France
tel: +33 (0)4 72 91 34 77
fax: +33 (0)4 72 91 34 61
portable: +33 (0)6 84 10 64 10
http://www.sbri.fr/members/kenneth-knoblauch.html