[Rd] initialize expression in 'quasi' (PR#8486)

Sun Jan 15 00:28:50 CET 2006

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--