# initialize expression in 'quasi' (PR#8486)

3 messages
Open this post in threaded view
|

## initialize expression in 'quasi' (PR#8486)

 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:[hidden email]=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-- ______________________________________________ [hidden email] mailing list https://stat.ethz.ch/mailman/listinfo/r-devel