Re: [R] New Family object for GLM models...

From: Thomas Lumley <tlumley_at_u.washington.edu>
Date: Wed 15 Jun 2005 - 00:41:20 EST

On Tue, 14 Jun 2005, Abatih E. wrote:
> Dear R-Users, I wish to create a new family object based on the Binomial
> family. The only difference will be with the link function. Thus instead
> if using the 'logit(u)' link function, i plan to use '-log(i-u)'. So
> far, i have tried to write the function following that of the Binomial
> and Negative Binomial families. The major problem i have here is with
> the definition of the initial values. The definitions i have attempted
> so far don't yield convergence when implemented in a model.

Are you sure the problem is with initial values? Your link function is capable of producing mu<0 from large enough negative eta. If the MLE is on the boundary of the parameter space it will not solve the likelihood equations and so won't be found by iterative reweighted least squares. It might be found by step-halving in glm, but that isn't guaranteed.

If the mle is in the interior of the parameter space then in my experience the starting values aren't terribly important (though my experience is with the log-binomial rather than -log(1-mu) link).

         -thomas

> I still
> can't figure out how the initial values are defined. I don't know if it
> is based on some property of the fitted model or the domain and /or
> range of the link function. I wish someone could help me provide a
> solution to this problem. I have appended the function i wrote.
>
>
> Add.haz<-function () {
>
>
>
> env <- new.env(parent = .GlobalEnv)
>
> assign(".nziek", nziek, envir = env)
>
> famname<-"Addhaza"
>
> link="addlink"
>
> linkfun<-function(mu) -log(1-mu)
>
> linkinv<-function(eta) 1-exp(-eta)
>
> variance<-function(mu) mu*(1-mu)
>
> validmu<-function(mu) all(mu > 0) && all(mu < 1)
>
> mu.eta<-function(mu) 1/(1-mu)
>
> dev.resids <- function(y, mu, wt) 2 * wt * (y * log(ifelse(y ==
>
> 0, 1, y/mu)) + (1 - y) * log(ifelse(y == 1, 1, (1 - y)/(1 -
>
> mu))))
>
> aic <- function(y, n, mu, wt, dev) {
>
> m <- if (any(n > 1))
>
> n
>
> else wt
>
> -2 * sum(ifelse(m > 0, (wt/m), 0) * dbinom(round(m *
>
> y), round(m), mu, log = TRUE))
>
> }
>
>
>
> initialize<- expression({
>
> n<-rep(1,nobs)
>
> if (any(y < 0 | y > 1)) stop("y values must be 0 <= y <= 1")
>
> mustart <- (weights * y + 0.5)/(weights + 1)
>
> m <- weights * y
>
> if (any(abs(m - round(m)) > 0.001)) warning("non-integer #successes in a binomial glm!")
>
>
>
>
>
> })
>
>
>
> environment(variance) <- environment(validmu) <- environment(dev.resids) <- environment(aic) <- env
>
>
>
> structure(list(family = famname, link = link, linkfun = linkfun,
>
> linkinv = linkinv, variance = variance, dev.resids = dev.resids,
>
> aic = aic, mu.eta =mu.eta, initialize = initialize,
>
> validmu = validmu), class = "family")
> }
>
> Thank you for your kind attention.
> Emmanuel
>
>
> EMMANUEL NJI ABATIH
> Rues des deux Eglises 140
> 1210 St Josse-Ten-Node, Bruxelles, BELGIUM
> MOBLIE: 0032486958988(anytime)
> Fix: 0032 2642 5038(8am to 6pm)
> http://www.iph.fgov.be/epidemio/epien
>
> ---------------------------------
>
>
> [[alternative HTML version deleted]]
>
> ______________________________________________
> R-help@stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
>

Thomas Lumley			Assoc. Professor, Biostatistics
tlumley@u.washington.edu	University of Washington, Seattle

______________________________________________
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html Received on Wed Jun 15 00:54:47 2005

This archive was generated by hypermail 2.1.8 : Fri 03 Mar 2006 - 03:32:36 EST