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

[Thomas Lumley]

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 at 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 at u.washington.edu	University of Washington, Seattle