initialize expression in 'quasi' (PR#8486)

classic Classic list List threaded Threaded
3 messages Options
Reply | Threaded
Open this post in threaded view
|

initialize expression in 'quasi' (PR#8486)

Bill.Venables
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
Reply | Threaded
Open this post in threaded view
|

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

Brian Ripley
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 [hidden email] 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:[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
>
>

--
Brian D. Ripley,                  [hidden email]
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

______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
Reply | Threaded
Open this post in threaded view
|

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

Bill.Venables
In reply to this post by Bill.Venables
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:[hidden email]
http://www.cmis.csiro.au/bill.venables/ 



-----Original Message-----
From: Prof Brian Ripley [mailto:[hidden email]]
Sent: Monday, 16 January 2006 6:28 PM
To: Venables, Bill (CMIS, Cleveland)
Cc: [hidden email]
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 [hidden email] 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:[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"
>
>
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--
>
> ______________________________________________
> [hidden email] mailing list
> https://stat.ethz.ch/mailman/listinfo/r-devel
>
>

--
Brian D. Ripley,                  [hidden email]
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

______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel