Friday, 15 July 2011

r - How do I use a custom link function in glm? -



r - How do I use a custom link function in glm? -

i don't want utilize standard log link in glm poisson regression, since have zeros. consider next code:

foo = 0:10 bar = 2 * foo glm(bar ~ foo, family = poisson(link = "identity"))

i error:

error: no valid set of coefficients has been found: please supply starting values

i'm not means. "identity" link function think (i.e. doesn't transform info @ all)? error mean , how can resolve it?

you can reply if start somewhere other default (0,0) starting point. start parameter vector containing intercept , slope of response, on scale of link function. problem r reporting typically calculated (negative) log-likelihood becomes infinite starting values. can check yourself: -sum(dpois(bar,0+0*foo,log=true)) inf (because setting poisson 0 mean, non-zero response).

however, isn't finish explanation, because starting points (0,2) starting negative-log-likelihood finite (-sum(dpois(bar,0+2*foo,log=true)) 20), same error occurs -- 1 have dig in deeper see what's matter, can imagine instance poisson mean of 0 not allowed @ in code. log-likelihood of poisson (a constant plus) x*log(lambda)-lambda: though works out ok if lambda , x both zero, that's not obvious in math. in particular, if @ poisson()$validmu, function glm uses found whether set of calculated means poisson ok, you'll see definition function (mu) { all(mu > 0) }. (it possible modify allow 0 values mu, plenty problem you'd need reason -- tried it, , there's problem because variances calculated equal zero. in short, easier through custom maximum likelihood estimator (e.g. bbmle::mle2()) hack glm ...)

however, starting point there no 0 estimates of poisson mean works out fine, although there plenty of warnings:

glm(bar ~ foo, family = poisson(link = "identity"), start=c(1,0))

however: want point out you're misunderstanding purpose of link function. it's fine have zeros in response variable of poisson regression, standard log link. glm model poisson regression y ~ poisson(exp(a+b*x)), not log(y) = + b*x. latter bad if y=0, former ok. glm(bar ~ foo, family = poisson()) works fine.

in general, non-canonical link functions bit of pain: they're need (although you've said i'm not convinced true in case), tend fussier , harder fit canonical links.

one final note: refer want "non-canonical" or "non-standard" link; custom link function, me, 1 wasn't provided family() command in r, had write link function (e.g. see http://rpubs.com/bbolker/4082 )

r glm

No comments:

Post a Comment