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 at 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.
>