-----BEGIN PGP SIGNED MESSAGE-----
Hash: SHA1
On 15-03-23 12:55 PM, Thierry Onkelinx wrote:
> Dear Ben,
>
> Last week I was struggling with incorporating lme4 into a package.
> I traced the problem and made a reproducible example (
> https://github.com/ThierryO/testlme4). It looks very simular to
> the problem you describe.
>
> The 'tests' directory contains the reproducible examples. confint()
> of a model as returned by a function fails. It even fails when I
> try to calculate the confint() inside the same function as the
> glmer() call (see the fit_model_ci function).
>
> Best regards,
>
> Thierry
Ugh. I can get this to work if I also try searching up the call
stack, as follows (within update.merMod). This feels like "code
smell" to me though -- i.e., if I have to hack this hard I must be
doing something wrong/misunderstanding how the problem *should* be done.
if (evaluate) {
ff <- environment(formula(object))
pf <- parent.frame() ## save parent frame in case we need it
sf <- sys.frames()[[1]]
tryCatch(eval(call, env=ff),
error=function(e) {
tryCatch(eval(call, env=sf),
error=function(e) {
eval(call, pf)
})
})
} else call
Here is an adapted even-more-minimal version of your code, which
seems to work with the version of update.merMod I just pushed to
github, but fails for glm():
## https://github.com/ThierryO/testlme4/blob/master/R/fit_model_ci.R
fit_model_ci <- function(formula, dataset, mfun=glmer){
model <- mfun(
formula = formula,
data = dataset,
family = "poisson"
)
ci <- confint(model)
return(list(model = model, confint = ci))
}
library("lme4")
set.seed(101)
dd <- data.frame(f=factor(rep(1:10,each=100)),
y=rpois(2,1000))
fit_model_ci(y~(1|f),dataset=dd)
fit_model_ci(y~(1|f),dataset=dd,mfun=glm)
>
>
> ir. Thierry Onkelinx Instituut voor natuur- en bosonderzoek /
> Research Institute for Nature and Forest team Biometrie &
> Kwaliteitszorg / team Biometrics & Quality Assurance Kliniekstraat
> 25 1070 Anderlecht Belgium
>
> To call in the statistician after the experiment is done may be no
> more than asking him to perform a post-mortem examination: he may
> be able to say what the experiment died of. ~ Sir Ronald Aylmer
> Fisher The plural of anecdote is not data. ~ Roger Brinner The
> combination of some data and an aching desire for an answer does
> not ensure that a reasonable answer can be extracted from a given
> body of data. ~ John Tukey
>
> 2015-03-22 17:45 GMT+01:00 Ben Bolker <bbol...@gmail.com>:
>
> WARNING: this is long. Sorry I couldn't find a way to compress
> it.
>
> Is there a reasonable way to design an update method so that it's
> robust to a variety of reasonable use cases of generating calls or
> data inside or outside a function? Is it even possible? Should I
> just tell users "don't do that"?
>
> * `update.default()` uses `eval(call, parent.frame())`; this fails
> when the call depends on objects that were defined in a different
> environment (e.g., when the data are generated and the model
> initially fitted within a function scope)
>
> * an alternative is to store the original environment in which the
> fitting is done in the environment of the formula and use
> `eval(call, env=environment(formula(object)))`; this fails if the
> user tries to update the model originally fitted outside a
> function with data modified within a function ...
>
> * I think I've got a hack that works below, which first tries in
> the environment of the formula and falls back to the parent frame
> if that fails, but I wonder if I'm missing something much simpler
> ..
>
> Thoughts? My understanding of environments and frames is still,
> after all these years, not what it should be ...
>
> I've thought of some other workarounds, none entirely
> satisfactory:
>
> * force evaluation of all elements in the original call * printing
> components of the call can get ugly (can save the original call
> before evaluating) * large objects in the call get duplicated *
> don't use `eval(call)` for updates; instead try to store
> everything internally * this works OK but has the same drawback of
> potentially storing large extra copies * we could try to use the
> model frame (which is stored already), but there are issues with
> this (the basis of a whole separate rant) because the model frame
> stores something in between predictor variables and input
> variables. For example
>
> d <- data.frame(y=1:10,x=runif(10))
> names(model.frame(lm(y~log(x),data=d))) ## "y" "log(x)"
>
> So if we wanted to do something like update to "y ~ sqrt(x)", it
> wouldn't work ...
>
> ================== update.envformula <- function(object,...) {
> extras <- match.call(expand.dots = FALSE)$... call <-
> getCall(object) for (i in names(extras)) { existing <-
> !is.na(match(names(extras), names(call))) for (a in
> names(extras)[existing]) call[[a]] <- extras[[a]] if
> (any(!existing)) { call <- c(as.list(call), extras[!existing]) call
> <- as.call(call) } } eval(call, env=environment(formula(object)))
> ## enclos=parent.frame() doesn't help }
>
> update.both <- function(object,...) { extras <-
> match.call(expand.dots = FALSE)$... call <- getCall(object) for (i
> in names(extras)) { existing <- !is.na(match(names(extras),
> names(call))) for (a in names(extras)[existing]) call[[a]] <-
> extras[[a]] if (any(!existing)) { call <- c(as.list(call),
> extras[!existing]) call <- as.call(call) } } pf <- parent.frame()
> ## save parent frame in case we need it tryCatch(eval(call,
> env=environment(formula(object))), error=function(e) { eval(call,
> pf) }) }
>
> ### TEST CASES
>
> set.seed(101) d <- data.frame(x=1:10,y=rnorm(10)) m1 <-
> lm(y~x,data=d)
>
> ##' define data within function, return fitted model f1 <-
> function() { d2 <- d lm(y~x,data=d2) return(lm(y~x,data=d2)) } ##'
> define (and modify) data within function, try to update ##' model
> fitted elsewhere f2 <- function(m) { d2 <- d; d2[1] <- d2[1]+0 ##
> force copy update.default(m,data=d2) } ##' define (and modify) data
> within function, try to update ##' model fitted elsewhere (use
> envformula) f3 <- function(m) { d2 <- d; d2[1] <- d2[1]+0 ## force
> copy update.envformula(m,data=d2) }
>
> ##' hack: first try the formula, then the parent frame ##' if that
> doesn't work for any reason f4 <- function(m) { d2 <- d; d2[1] <-
> d2[1]+0 ## force copy update.both(m,data=d2) }
>
> ## Case 1: fit within function m2 <- f1() try(update.default(m2))
> ## default: object 'd2' not found m3A <- update.envformula(m2) ##
> envformula: works m3B <- update.both(m2) ## works
>
> ## Case 2: update within function m4A <- f2(m1) ## default: works
> try(f3(m1)) ## envformula: object 'd2' not found m4B <- f4(m1)
> ## works
>
>>
>> ______________________________________________
>> R-devel@r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-devel
>>
>
-----BEGIN PGP SIGNATURE-----
Version: GnuPG v1.4.11 (GNU/Linux)
iQEcBAEBAgAGBQJVEfl/AAoJEOCV5YRblxUHWdgH/AqLAhDqKV8aRg6jnX9rO96D
nwzqv0ClMIxVr2dzD4eSQTL2caWZnXVkws+lg9N7bc4BaWplcYxLNRBw5M8zHOPJ
E7JlhG3EecvmeAEt9OY0/q6I0D6vdoEjcH7wzzuyLLIqllu9OskxURi/azMs0XRo
tiN+oG5aOKsMYsEGjtiWySRDzhJh2TM40A1HHjAViqpxZcqilAZ6RiNEFe1t1JY0
IvDI8yesSuHnKtgAiqk9ivGw4BCCGoBSIHB3GrJIi11j06iYKw0ugVHIlKYO8cqf
AYTvEX2sSxsJgKWYTiG/1dr/kiFTntTDji03zRLVUdPKIZATJMczv+KB+0bpoVY=
=Z34K
-----END PGP SIGNATURE-----
______________________________________________
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
