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.
