[Rd] initialize expression in 'quasi' (PR#8486)
Prof Brian Ripley
ripley at stats.ox.ac.uk
Mon Jan 16 09:28:19 CET 2006
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_Australia=
> .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"
>
> InF1YXNpIiA8LQ0KZnVuY3Rpb24gKGxpbmsgPSAiaWRlbnRpdHkiLCB2YXJpYW5jZSA9ICJjb25z
> dGFudCIpIA0Kew0KICAgIGxpbmt0ZW1wIDwtIHN1YnN0aXR1dGUobGluaykNCiAgICBpZiAoaXMu
> ZXhwcmVzc2lvbihsaW5rdGVtcCkgfHwgaXMuY2FsbChsaW5rdGVtcCkpIA0KICAgICAgICBsaW5r
> dGVtcCA8LSBsaW5rDQogICAgZWxzZSBpZiAoIWlzLmNoYXJhY3RlcihsaW5rdGVtcCkpIA0KICAg
> ICAgICBsaW5rdGVtcCA8LSBkZXBhcnNlKGxpbmt0ZW1wKQ0KICAgIGlmIChpcy5jaGFyYWN0ZXIo
> bGlua3RlbXApKSANCiAgICAgICAgc3RhdHMgPC0gbWFrZS5saW5rKGxpbmt0ZW1wKQ0KICAgIGVs
> c2Ugc3RhdHMgPC0gbGlua3RlbXANCiAgICB2YXJpYW5jZXRlbXAgPC0gc3Vic3RpdHV0ZSh2YXJp
> YW5jZSkNCiAgICBpZiAoIWlzLmNoYXJhY3Rlcih2YXJpYW5jZXRlbXApKSB7DQogICAgICAgIHZh
> cmlhbmNldGVtcCA8LSBkZXBhcnNlKHZhcmlhbmNldGVtcCkNCiAgICAgICAgaWYgKGxpbmt0ZW1w
> ID09ICJ2YXJpYW5jZSIpIA0KICAgICAgICAgICAgdmFyaWFuY2V0ZW1wIDwtIGV2YWwodmFyaWFu
> Y2UpDQogICAgfQ0KICAgIGluaXRpYWxpemUgPC0gTlVMTCAgICAgICAgICAgICAgICAgICAgICAg
> ICAgICAgICAgICAgICAgICAgICAgICAgIyMjIyBjaGFuZ2UNCiAgICBzd2l0Y2godmFyaWFuY2V0
> ZW1wLCBjb25zdGFudCA9IHsNCiAgICAgICAgdmFyaWFuY2UgPC0gZnVuY3Rpb24obXUpIHJlcC5p
> bnQoMSwgbGVuZ3RoKG11KSkNCiAgICAgICAgZGV2LnJlc2lkcyA8LSBmdW5jdGlvbih5LCBtdSwg
> d3QpIHd0ICogKCh5IC0gbXUpXjIpDQogICAgICAgIHZhbGlkbXUgPC0gZnVuY3Rpb24obXUpIFRS
> VUUNCiAgICB9LCAibXUoMS1tdSkiID0gew0KICAgICAgICB2YXJpYW5jZSA8LSBmdW5jdGlvbiht
> dSkgbXUgKiAoMSAtIG11KQ0KICAgICAgICB2YWxpZG11IDwtIGZ1bmN0aW9uKG11KSBhbGwobXUg
> PiAwKSAmJiBhbGwobXUgPCAxKQ0KICAgICAgICBkZXYucmVzaWRzIDwtIGZ1bmN0aW9uKHksIG11
> LCB3dCkgMiAqIHd0ICogKHkgKiBsb2coaWZlbHNlKHkgPT0gDQogICAgICAgICAgICAwLCAxLCB5
> L211KSkgKyAoMSAtIHkpICogbG9nKGlmZWxzZSh5ID09IDEsIDEsICgxIC0gDQogICAgICAgICAg
> ICB5KS8oMSAtIG11KSkpKQ0KICAgICAgICBpbml0aWFsaXplIDwtIGV4cHJlc3Npb24oeyAgICAg
> ICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIyMjIyBjaGFuZ2UNCiAgICAgICAgICBuIDwt
> IHJlcC5pbnQoMSwgbm9icykgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICMj
> IyMgY2hhbmdlDQogICAgICAgICAgbXVzdGFydCA8LSBwbWF4KDAuMDAxLCBwbWluKDAuOTk5LCB5
> KSkgICAgICAgICAgICAgICAgICAgICAjIyMjIGNoYW5nZQ0KICAgICAgICB9KSAgICAgICAgICAg
> ICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIyMjIyBjaGFu
> Z2UNCiAgICB9LCBtdSA9IHsNCiAgICAgICAgdmFyaWFuY2UgPC0gZnVuY3Rpb24obXUpIG11DQog
> ICAgICAgIHZhbGlkbXUgPC0gZnVuY3Rpb24obXUpIGFsbChtdSA+IDApDQogICAgICAgIGRldi5y
> ZXNpZHMgPC0gZnVuY3Rpb24oeSwgbXUsIHd0KSAyICogd3QgKiAoeSAqIGxvZyhpZmVsc2UoeSA9
> PSANCiAgICAgICAgICAgIDAsIDEsIHkvbXUpKSAtICh5IC0gbXUpKQ0KICAgIH0sICJtdV4yIiA9
> IHsNCiAgICAgICAgdmFyaWFuY2UgPC0gZnVuY3Rpb24obXUpIG11XjINCiAgICAgICAgdmFsaWRt
> dSA8LSBmdW5jdGlvbihtdSkgYWxsKG11ID4gMCkNCiAgICAgICAgZGV2LnJlc2lkcyA8LSBmdW5j
> dGlvbih5LCBtdSwgd3QpIHBtYXgoLTIgKiB3dCAqIChsb2coaWZlbHNlKHkgPT0gDQogICAgICAg
> ICAgICAwLCAxLCB5KS9tdSkgLSAoeSAtIG11KS9tdSksIDApDQogICAgfSwgIm11XjMiID0gew0K
> ICAgICAgICB2YXJpYW5jZSA8LSBmdW5jdGlvbihtdSkgbXVeMw0KICAgICAgICB2YWxpZG11IDwt
> IGZ1bmN0aW9uKG11KSBhbGwobXUgPiAwKQ0KICAgICAgICBkZXYucmVzaWRzIDwtIGZ1bmN0aW9u
> KHksIG11LCB3dCkgd3QgKiAoKHkgLSBtdSleMikvKHkgKiANCiAgICAgICAgICAgIG11XjIpDQog
> ICAgfSwgc3RvcChnZXR0ZXh0ZigiJ3ZhcmlhbmNlJyBcIiVzXCIgaXMgaW52YWxpZDogcG9zc2li
> bGUgdmFsdWVzIGFyZSBcIm11KDEtbXUpXCIsIFwibXVcIiwgXCJtdV4yXCIsIFwibXVeM1wiIGFu
> ZCBcImNvbnN0YW50XCIiLCANCiAgICAgICAgdmFyaWFuY2V0ZW1wKSwgZG9tYWluID0gTkEpKQ0K
> ICAgIGlmKGlzLm51bGwoaW5pdGlhbGl6ZSkpICAgICAgICAgICAgICAgICAgICAgICAgICAgICAg
> ICAgICAgICAgICAgIyMjIyBjaGFuZ2UNCiAgICBpbml0aWFsaXplIDwtIGV4cHJlc3Npb24oew0K
> ICAgICAgICBuIDwtIHJlcC5pbnQoMSwgbm9icykNCiAgICAgICAgbXVzdGFydCA8LSB5ICsgMC4x
> ICogKHkgPT0gMCkNCiAgICB9KQ0KICAgIGFpYyA8LSBmdW5jdGlvbih5LCBuLCBtdSwgd3QsIGRl
> dikgTkENCiAgICBzdHJ1Y3R1cmUobGlzdChmYW1pbHkgPSAicXVhc2kiLCBsaW5rID0gbGlua3Rl
> bXAsIGxpbmtmdW4gPSBzdGF0cyRsaW5rZnVuLCANCiAgICAgICAgbGlua2ludiA9IHN0YXRzJGxp
> bmtpbnYsIHZhcmlhbmNlID0gdmFyaWFuY2UsIGRldi5yZXNpZHMgPSBkZXYucmVzaWRzLCANCiAg
> ICAgICAgYWljID0gYWljLCBtdS5ldGEgPSBzdGF0cyRtdS5ldGEsIGluaXRpYWxpemUgPSBpbml0
> aWFsaXplLCANCiAgICAgICAgdmFsaWRtdSA9IHZhbGlkbXUsIHZhbGlkZXRhID0gc3RhdHMkdmFs
> aWRldGEsIHZhcmZ1biA9IHZhcmlhbmNldGVtcCksIA0KICAgICAgICBjbGFzcyA9ICJmYW1pbHki
> 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"
>
> c2V0LnNlZWQoNjY2KQ0KZGF0IDwtIGRhdGEuZnJhbWUoeCA9IHJlcCgoLTEwKToxMCwgZWFjaCA9
> IDUpLCB3ID0gcmVwKDE6NSwgMjEpKQ0KZGF0IDwtIHRyYW5zZm9ybShkYXQsIHkgPSByYmlub20o
> eCwgc2l6ZSA9IHcsIHByb2IgPSBwY2F1Y2h5KDEgKyAyKngpKSkNCg0KbW9kRml0IDwtIGdsbSh5
> L3cgfiB4LCBxdWFzaShsaW5rID0gY2F1Y2hpdCwgdmFyaWFuY2UgPSAibXUoMS1tdSkiKSwNCiAg
> ICAgICAgICAgICAgZGF0LCB3ZWlnaHRzID0gdywgdHJhY2UgPSBUKQ0KDQptb2RGaXQgPC0gZ2xt
> KHkvdyB+IHgsIHF1YXNpKGxpbmsgPSBjYXVjaGl0LCB2YXJpYW5jZSA9ICJtdSgxLW11KSIpLA0K
> ICAgICAgICAgICAgICBkYXQsIHdlaWdodHMgPSB3LCB0cmFjZSA9IFQsIG11c3RhcnQgPSBwbWF4
> KDAuMDAxLCBwbWluKDAuOTk5LCB5L3cpKSkNCiMgcXVhc2kgPC0gcXVhc2kubmV3DQoNCg0KbW9k
> Rml0IDwtIGdsbSh5L3cgfiB4LCBxdWFzaShsaW5rID0gY2F1Y2hpdCwgdmFyaWFuY2UgPSAibXUo
> 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