initialize expression in 'quasi' (PR#8486)

My apologies, I thought it was going as a text attachement.  Looks like
Microsoft second guessing me again.

The problem is the 'mustart' value (if you are going to base it on the
variance structure and not the link as well) should not allow values to
touch the point 
where the variance becomes zero.  The fudge that is in there now obeys
this
'rule' for all but the mu(1-mu).  It is a bit too strong for the
'constant'
case, but I guess we are stuck with that now for back compatibility
reasons.

Here is my suggestion:

##############(cut
here)##################################################

quasi <- function (link = "identity", variance = "constant")  {
    linktemp <- substitute(link)
    if (is.expression(linktemp) || is.call(linktemp)) 
        linktemp <- link
    else if (!is.character(linktemp)) 
        linktemp <- deparse(linktemp)
    if (is.character(linktemp)) 
        stats <- make.link(linktemp)
    else stats <- linktemp
    variancetemp <- substitute(variance)
    if (!is.character(variancetemp)) {
        variancetemp <- deparse(variancetemp)
        if (linktemp == "variance") 
            variancetemp <- eval(variance)
    }
    initialize <- NULL
#######change
    switch(variancetemp,
    constant = {
        variance <- function(mu) rep.int(1, length(mu))
        dev.resids <- function(y, mu, wt) wt * ((y - mu)^2)
        validmu <- function(mu) TRUE
    },
    "mu(1-mu)" = {
        variance <- function(mu) mu * (1 - mu)
        validmu <- function(mu) all(mu > 0) && all(mu < 1)
        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))))
        initialize <- expression({
#######change
          n <- rep.int(1, nobs)
#######change
          mustart <- pmax(0.001, pmin(0.999, y))	#######change
        })
#######change
    },
    mu = {
        variance <- function(mu) mu
        validmu <- function(mu) all(mu > 0)
        dev.resids <- function(y, mu, wt) 2 * wt * (y * log(ifelse(y == 
            0, 1, y/mu)) - (y - mu))
    },
    "mu^2" = {
        variance <- function(mu) mu^2
        validmu <- function(mu) all(mu > 0)
        dev.resids <- function(y, mu, wt) pmax(-2 * wt * (log(ifelse(y
== 
            0, 1, y)/mu) - (y - mu)/mu), 0)
    },
    "mu^3" = {
        variance <- function(mu) mu^3
        validmu <- function(mu) all(mu > 0)
        dev.resids <- function(y, mu, wt) wt * ((y - mu)^2)/(y * 
            mu^2)
    },
    stop(gettextf("'variance' \"%s\" is invalid: possible values are
\"mu(1-mu)\", \"mu\", \"mu^2\", \"mu^3\" and \"constant\"",
        variancetemp), domain = NA))
    if(is.null(initialize))
#######change
    initialize <- expression({
        n <- rep.int(1, nobs)
        mustart <- y + 0.1 * (y == 0)
    })
    aic <- function(y, n, mu, wt, dev) NA
    structure(list(family = "quasi", link = linktemp, linkfun =
stats$linkfun, 
        linkinv = stats$linkinv, variance = variance, dev.resids =
dev.resids, 
        aic = aic, mu.eta = stats$mu.eta, initialize = initialize, 
        validmu = validmu, valideta = stats$valideta, varfun =
variancetemp), 
        class = "family")
}

###############(cut here)##############################

Bill Venables, 
CMIS, CSIRO Laboratories, 
PO Box 120, Cleveland, Qld. 4163 
AUSTRALIA 
Office Phone (email preferred): +61 7 3826 7251 
Fax (if absolutely necessary):    +61 7 3826 7304 
Mobile (rarely used):                +61 4 1963 4642 
Home Phone:                          +61 7 3286 7700 
mailto:Bill.Venables at csiro.au 
http://www.cmis.csiro.au/bill.venables/ 

-----Original Message-----
From: Prof Brian Ripley [mailto:ripley at stats.ox.ac.uk] 
Sent: Monday, 16 January 2006 6:28 PM
To: Venables, Bill (CMIS, Cleveland)
Cc: r-devel at stat.math.ethz.ch
Subject: Re: [Rd] initialize expression in 'quasi' (PR#8486)

Bill,

As you see, R-bugs does not work with encoded (binary) attachments. 
Could you please send this again to R-devel with (if possible) text 
attachments?

The source code says

# 0.1 fudge here matches poisson: S has 1/6.
     initialize <- expression({ n <- rep.int(1, nobs); mustart <- y +
0.1 * (y == 0)})

although I believe S-PLUS has diverged here.

Brian