[Rd] initialize expression in 'quasi' (PR#8486)
Bill.Venables@csiro.au
Bill.Venables at csiro.au
Mon Jan 16 23:55:08 CET 2006
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
On Sun, 15 Jan 2006 Bill.Venables at csiro.au wrote:
> This is a multi-part message in MIME format.
>
> ------_=_NextPart_001_01C61962.212F35AB
> Content-Type: text/plain;
> charset="us-ascii"
> Content-Transfer-Encoding: quoted-printable
>
> This is not so much a bug as an infelicity in the code that can easily
> be fixed.
>
> The initialize expression in the quasi family function is, (uniformly
> for all links and all variance functions):
>
>
> initialize <- expression({
> n <- rep.int(1, nobs)
> mustart <- y + 0.1 * (y =3D=3D 0)
> })
>
> This is inappropriate (and often fails) for variance function
> "mu(1-mu)". Here is a short demo to show it:
> #################################################
>
> set.seed(666)
> dat <- data.frame(x =3D rep((-10):10, each =3D 5), w =3D rep(1:5, 21))
> dat <- transform(dat, y =3D rbinom(x, size =3D w, prob =3D pcauchy(1 +
=
> 2*x)))
>
> modFit <- glm(y/w ~ x, quasi(link =3D cauchit, variance =3D
"mu(1-mu)"),
> dat, weights =3D w, trace =3D T)
>
> Deviance =3D 309.2785 Iterations - 1=20
> Deviance =3D 3257.501 Iterations - 2=20
> Deviance =3D 1043.455 Iterations - 3=20
> ..
> Deviance =3D 1733.824 Iterations - 24=20
> Deviance =3D 1665.487 Iterations - 25=20
> Warning message:
> algorithm did not converge in: glm.fit(x =3D X, ...
> #################################################
>
> A comprehensive fix would involve tying the initialize expression to
> both the link and the variance function, but that would involve
changing
> make.link(), which is probably not a good idea for other reasons
(though
> ultimately that might be the way to go). There are at least three
> possible work-arounds:
>
> 1) use quasibinomial for this kind of model (but that's not possible
> here and an unnecessary complication if you are transferring S-PLUS
> code, which is how I found the problem) but quasibinomial could be
> extended to take this link as well, of course.
>
> 2) warn people in the help information that with quasi() they should
> always give the algorithm a bit more help and supply an appropriate
> mustart. This works fine (even if the coefficient estimates are a
bit
> ropey!):
>
> #################################################
> modFit <- glm(y/w ~ x, quasi(link =3D cauchit, variance =3D
"mu(1-mu)"),
> dat, weights =3D w, trace =3D T, mustart =3D pmax(0.001, pmin(0.999,
=
> y/w)))
>
> Deviance =3D 218.9552 Iterations - 1=20
> Deviance =3D 123.2773 Iterations - 2=20
> Deviance =3D 86.13804 Iterations - 3=20
> Deviance =3D 74.23746 Iterations - 4=20
> Deviance =3D 72.03787 Iterations - 5=20
> Deviance =3D 71.89566 Iterations - 6=20
> Deviance =3D 71.89395 Iterations - 7=20
> Deviance =3D 71.89395 Iterations - 8=20
> Deviance =3D 71.89395 Iterations - 9=20
>> =20
> #################################################
>
> 3) change quasi() slightly to cover this case, at least, in a better
> way.
>
> I have included a minimally modified version of quasi that I think
> achieves this and the demo shown above.=20
>
> With the changed version the performance, in this case, is identical
to
> what you get above when mustart is supplied. The changes cannot
affect
> performance with any other variance functions and with this variance
> function should only make things better, but it just _might_ make
things
> work worse in extreme and unusual cases. I have not found one,
though.
>
> Bill Venables.=20
>
>
>
> --please do not edit the information below--
>
> Version:
> platform =3D i386-pc-mingw32
> arch =3D i386
> os =3D mingw32
> system =3D i386, mingw32
> status =3D=20
> major =3D 2
> minor =3D 2.1
> year =3D 2005
> month =3D 12
> day =3D 20
> svn rev =3D 36812
> language =3D R
>
> Windows XP Professional (build 2600) Service Pack 2.0
>
> Locale:
>
LC_COLLATE=3DEnglish_Australia.1252;LC_CTYPE=3DEnglish_Australia.1252;LC
_=
> MON
>
ETARY=3DEnglish_Australia.1252;LC_NUMERIC=3DC;LC_TIME=3DEnglish_Australi
a=
> .1252
>
> Search Path:
> .GlobalEnv, package:mgcv, package:AusNew, package:MASS, package:RODBC,
> package:lattice, package:splines, package:methods, package:stats,
> package:graphics, package:grDevices, package:datasets,
package:RBigData,
> package:mvbutils, mvb.session.info, package:tools, package:utils,
> package:RBigData, package:RUtilities, package:RBigLibrary,
> package:g.data, Autoloads, package:base
>
> Bill Venables,=20
> CMIS, CSIRO Laboratories,=20
> PO Box 120, Cleveland, Qld. 4163=20
> AUSTRALIA=20
> Office Phone (email preferred): +61 7 3826 7251=20
> Fax (if absolutely necessary): +61 7 3826 7304=20
> Mobile (rarely used): +61 4 1963 4642=20
> Home Phone: +61 7 3286 7700=20
> mailto:Bill.Venables at csiro.au=20
> http://www.cmis.csiro.au/bill.venables/=20
>
>
>
> ------_=_NextPart_001_01C61962.212F35AB
> Content-Type: application/octet-stream;
> name="_quasi.R"
> Content-Transfer-Encoding: base64
> Content-Description: _quasi.R
> Content-Disposition: attachment;
> filename="_quasi.R"
>
>
InF1YXNpIiA8LQ0KZnVuY3Rpb24gKGxpbmsgPSAiaWRlbnRpdHkiLCB2YXJpYW5jZSA9ICJj
b25z
>
dGFudCIpIA0Kew0KICAgIGxpbmt0ZW1wIDwtIHN1YnN0aXR1dGUobGluaykNCiAgICBpZiAo
aXMu
>
ZXhwcmVzc2lvbihsaW5rdGVtcCkgfHwgaXMuY2FsbChsaW5rdGVtcCkpIA0KICAgICAgICBs
aW5r
>
dGVtcCA8LSBsaW5rDQogICAgZWxzZSBpZiAoIWlzLmNoYXJhY3RlcihsaW5rdGVtcCkpIA0K
ICAg
>
ICAgICBsaW5rdGVtcCA8LSBkZXBhcnNlKGxpbmt0ZW1wKQ0KICAgIGlmIChpcy5jaGFyYWN0
ZXIo
>
bGlua3RlbXApKSANCiAgICAgICAgc3RhdHMgPC0gbWFrZS5saW5rKGxpbmt0ZW1wKQ0KICAg
IGVs
>
c2Ugc3RhdHMgPC0gbGlua3RlbXANCiAgICB2YXJpYW5jZXRlbXAgPC0gc3Vic3RpdHV0ZSh2
YXJp
>
YW5jZSkNCiAgICBpZiAoIWlzLmNoYXJhY3Rlcih2YXJpYW5jZXRlbXApKSB7DQogICAgICAg
IHZh
>
cmlhbmNldGVtcCA8LSBkZXBhcnNlKHZhcmlhbmNldGVtcCkNCiAgICAgICAgaWYgKGxpbmt0
ZW1w
>
ID09ICJ2YXJpYW5jZSIpIA0KICAgICAgICAgICAgdmFyaWFuY2V0ZW1wIDwtIGV2YWwodmFy
aWFu
>
Y2UpDQogICAgfQ0KICAgIGluaXRpYWxpemUgPC0gTlVMTCAgICAgICAgICAgICAgICAgICAg
ICAg
>
ICAgICAgICAgICAgICAgICAgICAgICAgIyMjIyBjaGFuZ2UNCiAgICBzd2l0Y2godmFyaWFu
Y2V0
>
ZW1wLCBjb25zdGFudCA9IHsNCiAgICAgICAgdmFyaWFuY2UgPC0gZnVuY3Rpb24obXUpIHJl
cC5p
>
bnQoMSwgbGVuZ3RoKG11KSkNCiAgICAgICAgZGV2LnJlc2lkcyA8LSBmdW5jdGlvbih5LCBt
dSwg
>
d3QpIHd0ICogKCh5IC0gbXUpXjIpDQogICAgICAgIHZhbGlkbXUgPC0gZnVuY3Rpb24obXUp
IFRS
>
VUUNCiAgICB9LCAibXUoMS1tdSkiID0gew0KICAgICAgICB2YXJpYW5jZSA8LSBmdW5jdGlv
biht
>
dSkgbXUgKiAoMSAtIG11KQ0KICAgICAgICB2YWxpZG11IDwtIGZ1bmN0aW9uKG11KSBhbGwo
bXUg
>
PiAwKSAmJiBhbGwobXUgPCAxKQ0KICAgICAgICBkZXYucmVzaWRzIDwtIGZ1bmN0aW9uKHks
IG11
>
LCB3dCkgMiAqIHd0ICogKHkgKiBsb2coaWZlbHNlKHkgPT0gDQogICAgICAgICAgICAwLCAx
LCB5
>
L211KSkgKyAoMSAtIHkpICogbG9nKGlmZWxzZSh5ID09IDEsIDEsICgxIC0gDQogICAgICAg
ICAg
>
ICB5KS8oMSAtIG11KSkpKQ0KICAgICAgICBpbml0aWFsaXplIDwtIGV4cHJlc3Npb24oeyAg
ICAg
>
ICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIyMjIyBjaGFuZ2UNCiAgICAgICAgICBu
IDwt
>
IHJlcC5pbnQoMSwgbm9icykgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAg
ICMj
>
IyMgY2hhbmdlDQogICAgICAgICAgbXVzdGFydCA8LSBwbWF4KDAuMDAxLCBwbWluKDAuOTk5
LCB5
>
KSkgICAgICAgICAgICAgICAgICAgICAjIyMjIGNoYW5nZQ0KICAgICAgICB9KSAgICAgICAg
ICAg
>
ICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIyMjIyBj
aGFu
>
Z2UNCiAgICB9LCBtdSA9IHsNCiAgICAgICAgdmFyaWFuY2UgPC0gZnVuY3Rpb24obXUpIG11
DQog
>
ICAgICAgIHZhbGlkbXUgPC0gZnVuY3Rpb24obXUpIGFsbChtdSA+IDApDQogICAgICAgIGRl
di5y
>
ZXNpZHMgPC0gZnVuY3Rpb24oeSwgbXUsIHd0KSAyICogd3QgKiAoeSAqIGxvZyhpZmVsc2Uo
eSA9
>
PSANCiAgICAgICAgICAgIDAsIDEsIHkvbXUpKSAtICh5IC0gbXUpKQ0KICAgIH0sICJtdV4y
IiA9
>
IHsNCiAgICAgICAgdmFyaWFuY2UgPC0gZnVuY3Rpb24obXUpIG11XjINCiAgICAgICAgdmFs
aWRt
>
dSA8LSBmdW5jdGlvbihtdSkgYWxsKG11ID4gMCkNCiAgICAgICAgZGV2LnJlc2lkcyA8LSBm
dW5j
>
dGlvbih5LCBtdSwgd3QpIHBtYXgoLTIgKiB3dCAqIChsb2coaWZlbHNlKHkgPT0gDQogICAg
ICAg
>
ICAgICAwLCAxLCB5KS9tdSkgLSAoeSAtIG11KS9tdSksIDApDQogICAgfSwgIm11XjMiID0g
ew0K
>
ICAgICAgICB2YXJpYW5jZSA8LSBmdW5jdGlvbihtdSkgbXVeMw0KICAgICAgICB2YWxpZG11
IDwt
>
IGZ1bmN0aW9uKG11KSBhbGwobXUgPiAwKQ0KICAgICAgICBkZXYucmVzaWRzIDwtIGZ1bmN0
aW9u
>
KHksIG11LCB3dCkgd3QgKiAoKHkgLSBtdSleMikvKHkgKiANCiAgICAgICAgICAgIG11XjIp
DQog
>
ICAgfSwgc3RvcChnZXR0ZXh0ZigiJ3ZhcmlhbmNlJyBcIiVzXCIgaXMgaW52YWxpZDogcG9z
c2li
>
bGUgdmFsdWVzIGFyZSBcIm11KDEtbXUpXCIsIFwibXVcIiwgXCJtdV4yXCIsIFwibXVeM1wi
IGFu
>
ZCBcImNvbnN0YW50XCIiLCANCiAgICAgICAgdmFyaWFuY2V0ZW1wKSwgZG9tYWluID0gTkEp
KQ0K
>
ICAgIGlmKGlzLm51bGwoaW5pdGlhbGl6ZSkpICAgICAgICAgICAgICAgICAgICAgICAgICAg
ICAg
>
ICAgICAgICAgICAgIyMjIyBjaGFuZ2UNCiAgICBpbml0aWFsaXplIDwtIGV4cHJlc3Npb24o
ew0K
>
ICAgICAgICBuIDwtIHJlcC5pbnQoMSwgbm9icykNCiAgICAgICAgbXVzdGFydCA8LSB5ICsg
MC4x
>
ICogKHkgPT0gMCkNCiAgICB9KQ0KICAgIGFpYyA8LSBmdW5jdGlvbih5LCBuLCBtdSwgd3Qs
IGRl
>
dikgTkENCiAgICBzdHJ1Y3R1cmUobGlzdChmYW1pbHkgPSAicXVhc2kiLCBsaW5rID0gbGlu
a3Rl
>
bXAsIGxpbmtmdW4gPSBzdGF0cyRsaW5rZnVuLCANCiAgICAgICAgbGlua2ludiA9IHN0YXRz
JGxp
>
bmtpbnYsIHZhcmlhbmNlID0gdmFyaWFuY2UsIGRldi5yZXNpZHMgPSBkZXYucmVzaWRzLCAN
CiAg
>
ICAgICAgYWljID0gYWljLCBtdS5ldGEgPSBzdGF0cyRtdS5ldGEsIGluaXRpYWxpemUgPSBp
bml0
>
aWFsaXplLCANCiAgICAgICAgdmFsaWRtdSA9IHZhbGlkbXUsIHZhbGlkZXRhID0gc3RhdHMk
dmFs
>
aWRldGEsIHZhcmZ1biA9IHZhcmlhbmNldGVtcCksIA0KICAgICAgICBjbGFzcyA9ICJmYW1p
bHki
> KQ0KfQ0K
>
> ------_=_NextPart_001_01C61962.212F35AB
> Content-Type: application/octet-stream;
> name="demo_quasi.R"
> Content-Transfer-Encoding: base64
> Content-Description: demo_quasi.R
> Content-Disposition: attachment;
> filename="demo_quasi.R"
>
>
c2V0LnNlZWQoNjY2KQ0KZGF0IDwtIGRhdGEuZnJhbWUoeCA9IHJlcCgoLTEwKToxMCwgZWFj
aCA9
>
IDUpLCB3ID0gcmVwKDE6NSwgMjEpKQ0KZGF0IDwtIHRyYW5zZm9ybShkYXQsIHkgPSByYmlu
b20o
>
eCwgc2l6ZSA9IHcsIHByb2IgPSBwY2F1Y2h5KDEgKyAyKngpKSkNCg0KbW9kRml0IDwtIGds
bSh5
>
L3cgfiB4LCBxdWFzaShsaW5rID0gY2F1Y2hpdCwgdmFyaWFuY2UgPSAibXUoMS1tdSkiKSwN
CiAg
>
ICAgICAgICAgICAgZGF0LCB3ZWlnaHRzID0gdywgdHJhY2UgPSBUKQ0KDQptb2RGaXQgPC0g
Z2xt
>
KHkvdyB+IHgsIHF1YXNpKGxpbmsgPSBjYXVjaGl0LCB2YXJpYW5jZSA9ICJtdSgxLW11KSIp
LA0K
>
ICAgICAgICAgICAgICBkYXQsIHdlaWdodHMgPSB3LCB0cmFjZSA9IFQsIG11c3RhcnQgPSBw
bWF4
>
KDAuMDAxLCBwbWluKDAuOTk5LCB5L3cpKSkNCiMgcXVhc2kgPC0gcXVhc2kubmV3DQoNCg0K
bW9k
>
Rml0IDwtIGdsbSh5L3cgfiB4LCBxdWFzaShsaW5rID0gY2F1Y2hpdCwgdmFyaWFuY2UgPSAi
bXUo
>
MS1tdSkiKSwNCiAgICAgICAgICAgICAgZGF0LCB3ZWlnaHRzID0gdywgdHJhY2UgPSBUKQ0K
>
> ------_=_NextPart_001_01C61962.212F35AB--
>
> ______________________________________________
> R-devel at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-devel
>
>
--
Brian D. Ripley, ripley at stats.ox.ac.uk
Professor of Applied Statistics, http://www.stats.ox.ac.uk/~ripley/
University of Oxford, Tel: +44 1865 272861 (self)
1 South Parks Road, +44 1865 272866 (PA)
Oxford OX1 3TG, UK Fax: +44 1865 272595
More information about the R-devel
mailing list