Open MichaelChirico opened 4 years ago
MichaelChirico
commented
4 years ago Bumping this. Here is a slightly more severe case where neither coef() nor summary() work ...
library(nlme) data("cbpp",package="lme4") fm6 <- nlme::lmList(incidence∼ period | herd, data=cbpp) try(coef(fm6)) ## coef does not work here try(summary(fm6))
sapply(fm6,"[[","coefficients")
shows that the coefficients are NULL for one case ...
MichaelChirico
commented
4 years ago Bumping this. Here is a slightly more severe case where neither coef() nor summary() work ...
library(nlme) data("cbpp",package="lme4") fm6 <- nlme::lmList(incidence∼ period | herd, data=cbpp) try(coef(fm6)) ## coef does not work here try(summary(fm6))
sapply(fm6,"[[","coefficients")
shows that the coefficients are NULL for one case ...
github-actions[bot]
commented
4 years ago fm6 <- nlme::lmList(incidence∼ period | herd, data=cbpp)
Produces a warning for me in R 4.0.2.
I get
Warning message:
1 error caught in contrasts<-(*tmp*, value = contr.funs[1 + isOF[nn]]): contrasts can be applied only to factors with 2 or more levels
That says I think the documentation is confusing. This is what it says about the returned object:
Value
a list of lm objects with as many components as the number of groups defined by the grouping factor. Generic functions such as coef, fixed.effects, lme, pairs, plot, predict, random.effects, summary, and update have methods that can be applied to an lmList object.
and summary(fm6[[1]])
Returns without a problem (as do the other elements on the list. I also noticed that there is no summary.lmList() .
So I think either we need to make the documentation clearly say that summary() must be applied to each element in the list, or else create a summary.lmList function that applies summary.lm to all of the elements.
github-actions[bot]
commented
4 years ago Well, now that I read the StackOverflow link more carefully I see that there does seem to be a summary that works for a list without models that are not estimable.
github-actions[bot]
commented
4 years ago lmList() runs a separate model for each combination and puts them into a list of results. Because of this, when it encounters the NA it creates a model and results that simply do not have the variable with the NA. Using the original example ...
In terms of the summary function the "template" is (Intercept) day I(day^2) (Intercept) 0.266243029 -0.082171573 0.0047497082 day -0.082171573 0.043476465 -0.0029275587 I(day^2) 0.004749708 -0.002927559 0.0002126677
But lst[[10] (or lst[["9]])
9 : (Intercept) day (Intercept) 2.500 -0.3750 day -0.375 0.0625
Which does not match the template and hence gets subscript out of bounds.
"Assigned" usually means that a patch would be welcome so I will make some suggestions about what that might include.
It seems to me that there are a few options here.
Update the documentation to say that summary won't work in the presence of missing data, but that seems not serious in that the model itself runs and is fine.
Make a better error message, but again that seems not useful since the model is there it is that summary is not handling it well.
Change the summary function to cope with the situation where the dimensions of the list element are not the same as that for the template. I believe this could entail adding the missing column and filling it with NAs, as one solution.
nlme::summary.lmList can fail when some of the groups are too small for the model being fitted and hence give different-size output vcov matrices.
Over on StackOverflow someone is trying to get this code to work:
library("nlme") fm2Pixel.lis<-lmList(pixe∼day+I(day^2)|Dog, Pixel) summary(fm2Pixel.lis)
(results: Error in
[<-(*tmp*, use, use, ii, value = lst[[ii]]) : subscript out of bounds)http://stackoverflow.com/questions/32680306/plotting-individual-confidence-intervals-for-the-coefficients-in-the-lmlist-fit
The problem occurs because one Dog (#9) in the data set only has data from two days, so the quadratic term is automatically dropped from the corresponding lm() fit. (See with(Pixel,table(Dog,day)) )
nlme::summary.lmList seems to have some reasonably fancy code that tries to address this, but I haven't been able to get to the bottom of the problem yet.
METADATA