error in getting predictions from a lme object

Question

I am trying to get predictions for observations from a lme object. This is supposed to be quite straightforward. Yet, since I am get different types of errors for different trials, it seems to me I am missing something. My model is the following:

model <- lme(log(child_mortality) ~ as.factor(cluster)*time +
         my.new.time.one.transition.low.and.middle + ttd +
         maternal_educ+ log(IHME_id_gdppc) + hiv_prev-1,
         merged0,na.action=na.omit,method="ML",weights=varPower(form=~time),
         random= ~ time| country.x,
         correlation=corAR1(form = ~ time),
         control=lmeControl(msMaxIter = 200, msVerbose = TRUE))

It runs fine, fits the data well and the results make sense. Now to get predictions I have tried the following:

test.pred <- data.frame(time=c(10,10,10,10),country.x=c("Poland","Brazil",
            "Argentina","France"),    
             my.new.time.one.transition.low.and.middle=c(1,1,1,0),
             ttd=c(0,0,0,0),maternal_educ=c(10,10,10,10),
             IHME_id_gdppc=c(log(5000),log(8000),log(8000),log(15000)),   
             hiv_prev=c(.005,.005,.005,.005), 
             cluster=c("One Transition, Middle Income","One Transition,   
             Middle Income","One Transition, Middle Income","Democracy, 
             High Income"))
>
> predict(model,test.pred,level=0)


Error in X %*% fixef(object) : non-conformable arguments

If I exclude, say, France, and only include countries in which cluster="OneTransition, Middle Income" then I get a different error

# create a toy data set
test.pred0 <-
    expand.grid(time=20:29,country.x=c("Poland","Brazil","Argentina"))
z0 <-as.data.frame(cbind(my.new.time.one.transition.low.and.middle = 
                         c(0,0,0,0,0,0,1,2,3,4), ttd=c(0,0,0,0,0,0,1,0,0,0),
                         maternal_educ=seq(from=10.0, to=12.0, length.out=10),
                         IHME_id_gdppc=log(seq(from=5000, to=8000, length.out=10)),
                         hiv_prev=rep(.005,10),
                         cluster=rep("One Transition, Middle Income",10)))

z <- rbind(z0,z0,z0)
test.pred <- cbind(test.pred0,z)
# check
head(test.pred)
>  time country.x my.new.time.one.transition.low.and.middle ttd
> maternal_educ    IHME_id_gdppc hiv_prev
> 1   20    Poland                                         0   0
>   10 8.51719319141624    0.005
> 2   21    Poland                                         0   0
> 10.2222222222222 8.58173171255381    0.005
> 3   22    Poland                                         0   0
> 10.4444444444444 8.64235633437024    0.005
> 4   23    Poland                                         0   0
> 10.6666666666667 8.69951474821019    0.005
> 5   24    Poland                                         0   0
> 10.8888888888889 8.75358196948047    0.005
> 6   25    Poland                                         0   0
> 11.1111111111111 8.80487526386802    0.005
>                         cluster
> 1 One Transition, Middle Income
> 2 One Transition, Middle Income
> 3 One Transition, Middle Income
> 4 One Transition, Middle Income
> 5 One Transition, Middle Income
> 6 One Transition, Middle Income

# run the predictions
predict(model,test.pred,level=0)
> Error in `contrasts<-`(`*tmp*`, value = contr.funs[1 + isOF[nn]]) :
>   contrasts can be applied only to factors with 2 or more levels

In this example, the problem is due to cluster="One Transition, Middle Income" all the time.

I don't understand why this is a problem. If I want to get predict() to work I have to include all variables from the model, right? Obviously, the input data in the model's call will not include factor set to the same values for all cases. Yet, if I want to get predictions just for subset of the data, or for new observations, I may be interested in only in cases where some factor is always set to be the same. Does it make sense? How can I get predictions in that case?

Answer 1

Thanks for providing the data so that I could perform some diagnostics. Actually, this is an epic bug of predict.lme. Your factors have more levels in your initial data (for example you have more than 4 countries) than in your new data. A line of code specifically causes the unused levels to be discarded so you end up with matrices of different dimensions, whence the non-conformable arguments

I removed that line and put the code here.

In R you can do

library(nlme)
source("http://lab.thegrandlocus.com/static/code/predict.lme_patched.txt")

This registers a new function predict.lme that will be invoked instead of the one from the package nlme and you can run your code. At least it worked for me.

Warning: The posted code and the method are neither a replacement nor a real bug fix of the package. The patched function has not been tested beyond its ability to run the bit of code of the OP.