Dear R users, I'm struggling with applying a user-defined link function in lmer. For analyzing data of a 2AFC psychophysical experiment, I would like to model my binary data with a logistic function with a lower limit at 0.5 instead of 0. In a previous question this has been described as a halflogit function. To do so I wrote my own link function and would like to submit it to lmer, however this results in error messages. Below I've listed some of my attempts to solve this problem based on previous postings on this list. I've also added the corresponding problems I encountered with each of these attempts.
1) I've used the link function as specified at https://stat.ethz.ch/pipermail/r-help/2007-August/138436.html # define link function halflogit=function(){ half.logit=function(mu) qlogis(2*mu-1) half.logit.inv=function(eta) .5*plogis(eta)+.5 half.logit.deriv=function(eta) .5*(exp(eta/2)+exp(-eta/2))^-2 half.logit.inv.indicator=function(eta) TRUE half.logit.indicator=function(mu) mu>.5 & mu<1 link <- "half.logit" structure(list(linkfun = half.logit, linkinv = half.logit.inv, mu.eta = half.logit.deriv, validmu = half.logit.indicator ,valideta = half.logit.inv.indicator, name = link), class = "link-glm") } # submit halflogit function to lmer rcNrPSrcNxPSxP<-lmer(Correct~-1+cNoise*f_Pos+f_Shape_con+f_Shape_con:f_Pos+(-1+cNoise*f_Pos|f_Subject), data=dens_data, family=binomial(halflogit())) #resulting error message Error in if (!(validmu(mu) && valideta(eta))) stop("cannot find valid starting values: please specify some", : missing value where TRUE/FALSE needed In addition: Warning message: In qlogis(p, location, scale, lower.tail, log.p) : NaNs produced 2) Also making the variable global seems not to help (see also http://permalink.gmane.org/gmane.comp.lang.r.lme4.devel/4501) 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") } lmer(Correct~-1+cNoise*f_Pos+f_Shape_con+f_Shape_con:f_Pos+(-1+cNoise*f_Pos|f_Subject), data=dens_data, family=binomial(mafc.logit(2))) Error in famType(glmFit$family) : unknown link: 'mafc.logit(2)' 3) The link function does work with glm so the problem is in lmer, I guess, not in the newly defined link function. Unfortunately I would like to model random effect, so I prefer to use the lmer function. 4) I've also tried with other versions of the package: lme4a and lme4.0, but lme4a seems not to work under the current R version, and lme4.0 gives the same errors as lme4. 5) I've tried to transform the data myself with my own function and fit a linear model on the transformed data. However transforming results in Inf, -Inf and NaN for proportion of 1, 0.5 and <0.5 respectively and a linear model cannot handle Inf, -Inf and NaN. In addition, I'm not sure if this is an appropriate way to do my analyses because the variance on the data would be estimated as an additional parameter in a linear model while in logistic regression the variance is prop(1-prop) and is not estimated as an additional parameter. 6) Last I tried using glmmPQL but this does not reach convergence and I'm not sure if it's suitable for my analyses. I would appreciate any help on applying a halflogit function to lmer. Kind regards, Kathleen --------------------------------------------------------------------------- Kathleen Vancleef PhD student Laboratory of Experimental Psychology University of Leuven Tiensestraat 102 box 3711 B-3000 Leuven Belgium Tel: +32 (0)16/32.62.83 Email: kathleen.vancl...@ppw.kuleuven.be http://www.gestaltrevision.be ______________________________________________ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
