erreur lors de l'obtention des prédictions d'un objet lme

8

J'essaie d'obtenir des prédictions pour les observations d'un objet lme. C'est censé être assez simple. Pourtant, comme je reçois différents types d'erreurs pour différents essais, il me semble que je manque quelque chose. Mon modèle est le suivant:

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))

Il fonctionne bien, s'adapte bien aux données et les résultats ont du sens. Maintenant, pour obtenir des prédictions, j'ai essayé ce qui suit:

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

Si j'exclus, disons, la France, et que j'inclus uniquement les pays dans lesquels cluster = "OneTransition, Middle Income", je reçois une erreur différente

# 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

Dans cet exemple, le problème est dû à cluster = "One Transition, Middle Income" tout le temps.

Je ne comprends pas pourquoi c'est un problème. Si je veux que Predict () fonctionne, je dois inclure toutes les variables du modèle, non? Évidemment, les données d'entrée dans l'appel du modèle n'incluront pas de facteur réglé aux mêmes valeurs pour tous les cas. Pourtant, si je veux obtenir des prédictions uniquement pour un sous-ensemble de données ou pour de nouvelles observations, je ne serai intéressé que dans les cas où un facteur est toujours réglé pour être le même. Est-ce que ça fait du sens? Comment puis-je obtenir des prédictions dans ce cas?

Antonio Pedro Ramos
la source
Je soupçonne que c'est parce que lorsque vous voyez deux facteurs, l'un le sous-ensemble de l'autre, R voit deux facteurs non liés. Juste une idée: essayez de démarrer une nouvelle session R, de taper options(stringsAsFactors = FALSE), puis d'exécuter votre code. Cela empêcherait votre original test.predd'avoir ses propres facteurs.
Matt Parker
Merci Matt, mais ça ne marche toujours pas. Je suis en fait un casse-tête. Ce doit être une sorte d'erreur.
Antonio Pedro Ramos
Juste un travail autour, par exemple un facteur de 3 niveaux A, B, C, vous pouvez faire une prédiction pour le niveau A avec des données de 100 A, 1 B et 1 C.
Verbal

Réponses:

8

Merci d'avoir fourni les données pour que je puisse effectuer des diagnostics. En fait, c'est un bug épique de predict.lme. Vous factorsavez plus de niveaux dans vos données initiales (par exemple, vous avez plus de 4 pays) que dans vos nouvelles données. Une ligne de code provoque spécifiquement l'élimination des niveaux inutilisés, vous vous retrouvez avec des matrices de différentes dimensions, d'oùnon-conformable arguments

J'ai supprimé cette ligne et mis le code ici .

Dans R, vous pouvez faire

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

Cela enregistre une nouvelle fonction predict.lmequi sera invoquée à la place de celle du package nlmeet vous pouvez exécuter votre code. Au moins, cela a fonctionné pour moi.

Avertissement: Le code publié et la méthode ne sont ni un remplacement ni une véritable correction de bogue du package. La fonction corrigée n'a pas été testée au-delà de sa capacité à exécuter le bit de code de l'OP.

gui11aume
la source
En fait, c'est le cas. J'ai country.xdans les deux. Le jury est toujours dehors.
Antonio Pedro Ramos
Ouais ... c'est vrai. Désolé pour ça. Il semble que certains de vos types de données ne soient pas les mêmes dans votre entrée initiale et dans vos nouvelles données. Ce sera très difficile de se passer des données. Si elle n'est pas confidentielle, et pas trop grande, pourriez-vous sauvegarder la session R et la mettre quelque part (ou me l'envoyer par mail)?
gui11aume
Merci beaucoup. Avez-vous un e-mail pour que je vous envoie un exemple de code et de données?
Antonio Pedro Ramos
Juste une question rapide: cette version fonctionne bien pour levels=1mais pas pour level=0?
Antonio Pedro Ramos