Logistic and nonlinear mixed models: Accounting for guessing probability

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