8
$\begingroup$

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?

$\endgroup$
3
  • $\begingroup$ I suspect this is because where you see two factors, one the subset of the other, R sees two unrelated factors. Just a hunch: try starting a fresh R session, typing options(stringsAsFactors = FALSE), and then running your code. That would prevent your original test.pred from having its own factors. $\endgroup$ May 30 '12 at 22:28
  • $\begingroup$ Thank you Matt, but is still not working. I am actually puzzle. It must be some kind of mistake. $\endgroup$ Jun 5 '12 at 2:41
  • $\begingroup$ Just a work around, e.g. a factor of 3 levels A, B, C you can make prediction for level A with a data of 100 A, 1 B and 1 C. $\endgroup$
    – Verbal
    Oct 15 '15 at 23:40
8
$\begingroup$

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.

$\endgroup$
4
  • $\begingroup$ Actually, it does. I have country.x in both. The jury is still out. $\endgroup$ Jun 1 '12 at 15:50
  • $\begingroup$ Yep... that's right. Sorry about that. It seems that some of your data types are not the same in your initial input and your new data. This will be very hard to do without the data. If it is not confidential, and not too big, could you save the R session and put it somewhere (or send it to me by mail)? $\endgroup$
    – gui11aume
    Jun 1 '12 at 18:03
  • $\begingroup$ Thank you very much. Do you have an email for me to send you sample code and data? $\endgroup$ Jun 5 '12 at 2:35
  • $\begingroup$ Just a quick question: this version works fine for levels=1 but not for level=0? $\endgroup$ Jun 11 '12 at 23:08

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.