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

Ken Knoblauch ken.knoblauch at inserm.fr
Fri Oct 1 16:45:20 CEST 2010

```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)
>
> 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")
mu <- pmax(mu, 1/.m + .Machine\$double.eps)
qlogis((.m * mu - 1)/(.m - 1))
}
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 = "")
}

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

```