Error in glm(..., family=quasi(..., variance=list(...)))

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

Error in glm(..., family=quasi(..., variance=list(...)))

Wollschlaeger, Daniel
In a glm() call using a quasi() family, one may define a custom variance function in the form of a "list containing components varfun, validmu, dev.resids, initialize and name" (quoting the help page for family). In trying to do so, I run into the following issue that I have not seen discussed previously:

x  <- runif(1000, min=0, max=1)
y  <- x + rnorm(1000, mean=0, sd=1)*x^(3/4)
vf <- function(mu) { abs(mu)^(3/4) }
vm <- function(mu) { rep(TRUE, length(mu)) }
dr <- function(y, mu, wt) { (y-mu)^2 }
it <- expression({ n <- rep.int(1, nobs); mustart <- y })
glm(y ~ x, family=quasi(link="identity", variance=list(varfun=vf, validmu=vm, dev.resids=dr, initialize=it, name="custom")))

This gives "Error in switch(vtemp, constant = { : EXPR must be a length 1 vector" from line 576 in file family.R (https://github.com/wch/r-source/blob/5a156a0865362bb8381dcd69ac335f5174a4f60c/src/library/stats/R/family.R#L576). I believe this is due to line 573 "vtemp <- substitute(variance)" and 574 "if (!is.character(vtemp)) vtemp <- deparse(vtemp)" where vtemp becomes a length 2 character vector because, by default, deparse() breaks lines at width.cutoff = 60L characters. In stepping through quasi() during debug, setting vtemp <- (vtemp, collapse=" ") on line 576 avoids the error.

A workaround from https://tolstoy.newcastle.edu.au/R/help/05/06/6795.html appears to be to define one's own complete quasi2() function with the desired variance function pre-stored.

Is this known/expected, or should I file a bug?

Many thanks and best regards

Daniel

> sessionInfo()
R version 3.6.0 (2019-04-26)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 10 x64 (build 14393)

Matrix products: default

locale:
[1] LC_COLLATE=German_Germany.1252  LC_CTYPE=German_Germany.1252  
[3] LC_MONETARY=German_Germany.1252 LC_NUMERIC=C                  
[5] LC_TIME=German_Germany.1252    

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base    

loaded via a namespace (and not attached):
[1] compiler_3.6.0

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

Re: Error in glm(..., family=quasi(..., variance=list(...)))

Martin Maechler
>>>>> Wollschlaeger, Daniel
>>>>>     on Fri, 26 Apr 2019 15:13:36 +0000 writes:

    > In a glm() call using a quasi() family, one may define a custom variance function in the form of a "list containing components varfun, validmu, dev.resids, initialize and name" (quoting the help page for family). In trying to do so, I run into the following issue that I have not seen discussed previously:

    > x  <- runif(1000, min=0, max=1)
    > y  <- x + rnorm(1000, mean=0, sd=1)*x^(3/4)
    > vf <- function(mu) { abs(mu)^(3/4) }
    > vm <- function(mu) { rep(TRUE, length(mu)) }
    > dr <- function(y, mu, wt) { (y-mu)^2 }
    > it <- expression({ n <- rep.int(1, nobs); mustart <- y })
    > glm(y ~ x, family=quasi(link="identity", variance=list(varfun=vf, validmu=vm, dev.resids=dr, initialize=it, name="custom")))

    > This gives "Error in switch(vtemp, constant = { : EXPR must be a length 1 vector"
 > from line 576 in file family.R
 > (https://github.com/wch/r-source/blob/5a156a0865362bb8381dcd69ac335f5174a4f60c/src/library/stats/R/family.R#L576).
 > I believe this is due to line 573 "vtemp <- substitute(variance)" and 574 "if (!is.character(vtemp)) vtemp <- deparse(vtemp)" where vtemp becomes a length 2 character vector because, by default, deparse() breaks lines at width.cutoff = 60L characters. In stepping through quasi() during debug, setting vtemp <- (vtemp, collapse=" ") on line 576 avoids the error.

but really in this case, neither substitute() nor deparse()
should be called !

    > A workaround from https://tolstoy.newcastle.edu.au/R/help/05/06/6795.html appears to be to define one's own complete quasi2() function with the desired variance function pre-stored.

    > Is this known/expected, or should I file a bug?

and you *have* filed a bug report.  Thank you!
--> https://bugs.r-project.org/bugzilla/show_bug.cgi?id=17560

Note that the above R-help "workaround" goes back to 2005, and I
find it amazing this has never been taken up.

I've committed fix to this bug avoiding the misleading
substitute() and deparse() altogether in this case.

The plan is to port the bug fix also to "R 3.6.0 patched" so the
problem would be solved in the next released version of R.


    > Many thanks and best regards
    > Daniel

    >> sessionInfo()
    > R version 3.6.0 (2019-04-26)
    ..........

Thank you again!
Martin

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