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(pixel~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.