[R] half-logit and glm (again)
Richard D. Morey
moreyr at missouri.edu
Fri Aug 10 20:37:25 CEST 2007
I know this has been dealt with before on this list, but the previous
messages lacked detail, and I haven't figured it out yet.
The model is:
\x_{ij} = \mu + \alpha_i + \beta_j
\alpha is a random effect (subjects), and \beta is a fixed effect
(condition).
I have a link function:
p_{ij} = .5 + .5( 1 / (1 + exp{ -x_{ij} } ) )
Which is simply a logistic transformed to be between .5 and 1.
The data y_{ij} ~ Binomial( p_{ij}, N_{ij} )
I've generated data using this model, and I'd like to fit it. My data is
a data frame with 3 columns, "response" (0/1), "subject" (a factor), and
"condition" (another factor).
Here is my link definition:
#############################
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")
}
binomial(halflogit())
Family: binomial
Link function: half.logit
#############################
I based this off the help for the family() function.
So I try to call glmmPQL (based on other R-help posts, this is the
easiest to use?)
#################
glmmPQL(response ~ condition, random = ~ 1|subject, family =
binomial(halflogit()), data = dat)
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:
NaNs produced in: qlogis(p, location, scale, lower.tail, log.p)
#################
It looks like I've misdefined something and it is going outside the
specified domains for the functions. I can't find any reference to
starting starting values in help for glmmPQL() or lme().
If anyone has any working code where they've done a user defined link
function, it would be greatly appreciated.
Thanks,
Richard
More information about the R-help
mailing list