Cormac,
thanks for the detailed information on this. Three comments:
(a) The "horrible failure" that encomptest() produces in your second issue
is for me:
R> encomptest(models12[[1]],models1[[1]]) ##horrible failure!
Error in encomptest(models12[[1]], models1[[1]]) :
models were not all fitted to the same size of dataset
which seems to be a reasonable (and correct) error message.
(b) In the first issue: Both lmList() and encomptest() are intended to
make life with linear models somewhat easier or more convenient. But
apparently they don't cooperate very well because both assume some
standard setting. I don't see an easy way how I could change encomptest()
to cooperate generally with lmList(). And my impression is that the same
is true for lmList. But I don't know that implementation well enough to
say whether your fix will cover all cases.
(c) With such specific requests about CRAN packages, it is typically
better to contact the package maintainers directly or at least cc them in
the mail to the list (as the posting guide suggests).
Best,
Z
On Tue, 19 Jun 2012, Cormac Long wrote:
Hello R-Help,
-----------------------------------------------------------------------------------------------------------------------------------------
Issues (there are 2):
1) Possible bug when using lmtest::encomptest() with a linear model
created using nlme::lmList()
2) Possible modification to lmtest::encomptest() to fix confusing fail
when models provided are, in fact, nested.
I have had a look around the web and in R-devel, but not been able to find
any references to these
issues.
-----------------------------------------------------------------------------------------------------------------------------------------
-----------------------------------------------------------------------------------------------------------------------------------------
Issue 1:
-----------------------------------------------------------------------------------------------------------------------------------------
-----------------------------------------------------------------------------------------------------------------------------------------
Description:
----------------
encomptest fails with error
Error in is.data.frame(data) : object 'dat' not found
when applied to a pair of linear models generated via the nlme::lmList()
function
Reproducibility:
---------------------
#######################################################################
##Demonstration of encomptest breaking when models are created via the
##nlme::lmList() function.
##Required libraries:
library(nlme)
library(lmtest)
library(plyr)
##Create a noisy slope (spoofing real-world data):
fakeData<-c(rnorm(10))+seq(0,5,length=10)
##Create a data frame of data to test. The columns
##model1 and model2 contain my basis functions
modelFrame<-as.data.frame(
list(
srcData=fakeData,
model1=seq(0,1,length=10),
model2=c(rep(0,5),rep(1,5))
)
)
##Create all data to be fitted
allData<-ldply(1:2,function(x){modelFrame$modelFactor<-x;return(modelFrame)})
##Apply models by factor:
models1<-lmList(srcData ~ model1 | modelFactor,data=allData)
models2<-lmList(srcData ~ model2 | modelFactor,data=allData)
##Attempt to apply encomptest - it fails!
encomptest(models1[[1]],models2[[1]])
##Now, create two models by hand
modelFrame1<-modelFrame
modelFrame2<-modelFrame
models1<-list(lm(srcData ~ model1 ,data=allData))
models2<-list(lm(srcData ~ model2 ,data=allData))
##Attempt to apply encomptest - it works!
encomptest(models1[[1]],models2[[1]])
##End demonstration:
#######################################################################
Traceback and investigation:
---------------------------------------
This error is raised at the line containing the code
fm <- update(fm1, fm)
in the implementation of lmtest::encomptest(). This appears to be due to
the implementation
of nlme::lmList.formula(). The offending line is the evaluation
val <- lapply(split(data, groups), function(dat, form, na.action)
{lm(formula = form, data = dat, na.action = na.action)}, form = object,
na.action = na.action)
of nlme::lmList.formula(). The resulting model has call
lm(formula = form, data = dat, na.action = na.action)
This call is parsed in stats::update.default and the resulting call is
evaluated, looking for the
symbol 'dat'. However, the only place that the data frame 'data' is
referred to as 'dat' that I can
see is is in the inline function
function(dat, form, na.action) {lm(formula = form, data = dat,
na.action = na.action)}
as previously observed in nlme::lmList.formula()
Suggested fix:
-------------------
Replace
function(dat, form, na.action) {lm(formula = form, data = dat,
na.action = na.action)}
with
function(data, form, na.action) {lm(formula = form, data = data,
na.action = na.action)}
-----------------------------------------------------------------------------------------------------------------------------------------
-----------------------------------------------------------------------------------------------------------------------------------------
Issue 2:
-----------------------------------------------------------------------------------------------------------------------------------------
-----------------------------------------------------------------------------------------------------------------------------------------
Description:
----------------
encomptest fails with error
Error in parse(text = x) : <text>:2:0: unexpected end of input
1: . ~ . +
^
when applied to a pair of nested linear models
Reproducibility:
----------------------
#######################################################################
##Demonstration of encomptest breaking when models are, in fact, nested
##Required libraries:
library(lmtest)
##Create models:
models12<-lmList(srcData ~ model1 + model2 | modelFactor,data=allData)
##Function to test if models are nested
isNested<-function(model1, model2){
##Test if model1 is nested in model2:
m1labs<-attr(terms(model1), "term.labels")
m2labs<-attr(terms(model2), "term.labels")
rv<-m1labs[!(m1labs %in% m2labs)]
return(ifelse(length(rv),FALSE,TRUE))
}
##Show models are nested:
isNested(models1[[1]],models12[[1]]) ##See, models are nested!
##Force encomptest to break:
encomptest(models12[[1]],models1[[1]]) ##horrible failure!
##End demonstration:
#######################################################################
Traceback and investigation:
---------------------------------------
This happens because model1 is contained in model12. This error is raised
at the
line containing the code
fm <- as.formula(paste(". ~ . +", fm))
in the encompassing test implementation in lmtest.
Suggested fix:
-------------------
This could be avoided if additional error checking is added into
encomptest(): check if
provided models are nested and called a halt via stop() if they are.
-----------------------------------------------------------------------------------------------------------------------------------------
-----------------------------------------------------------------------------------------------------------------------------------------
Many thanks (in advance) for your help!
Best wishes,
Dr. Cormac Long.
[[alternative HTML version deleted]]
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
