conjugateprior
conjugateprior

Reputation: 398

R model fitting routine fails to run in a function

I'm using Firth and Turner's BradleyTerry2 package for paired comparisons but have run into a mysterious problem using the main fitting function BTm. Here is a minimal data setup from their own example:

data(citations, package = "BradleyTerry2")
citations.sf <- countsToBinomial(citations)
names(citations.sf)[1:2] <- c("journal1", "journal2")

So at the console the following works:

citeModel <- BTm(cbind(win1, win2), journal1, journal2, data = citations.sf)

but the following does not work

f1 <- function(x){ BTm(cbind(win1, win2), journal1, journal2, data=x) }
f1(citations.sf)

while this (statistically nonsensical but) structurally similar linear model example does work, just as I would expect:

f2 <- function(x){ lm(log(win1/win2) ~ journal1, data=x) }
f2(citations.sf)

The error with f1 is "Error in eval(substitute(expr), data, enclos = parent.frame()): invalid 'envir' argument". But this is not telling me anything I can understand.

Thoughts?

Upvotes: 1

Views: 998

Answers (1)

Heather Turner
Heather Turner

Reputation: 3314

Thanks for linking to this post on the existing bug report [#2294] objects not found when BTm not called in global environment.

When setting up the data internally, BTm looks for objects in the environment of the formula and then in the global environment. If a formula is not specified to BTm, then the formula is defined internally and, as currently implemented, inherits the environment constructed when BTm is called, rather than the environment from which BTm was called.

Until I get round to fixing this, there is a simple work-around - always specify a formula when using BTm inside a function, so that the environment of the formula is that created when your function is called. E.g.

> f1 <- function(x){ BTm(cbind(win1, win2), journal1, journal2, ~.., data=x) }
> f1(citations.sf)
Bradley Terry model fit by glm.fit 

Call:  BTm(outcome = cbind(win1, win2), player1 = journal1, player2 = journal2, 
    formula = ~.., data = x)

Coefficients:
..Comm Statist          ..JASA        ..JRSS-B  
       -2.9491         -0.4796          0.2690  

Degrees of Freedom: 6 Total (i.e. Null);  3 Residual
Null Deviance:  1925 
Residual Deviance: 4.293    AIC: 46.39 

Upvotes: 2

Related Questions