error al obtener predicciones de un objeto lme

8

Estoy tratando de obtener predicciones para observaciones de un objeto lme. Se supone que esto es bastante sencillo. Sin embargo, dado que recibo diferentes tipos de errores para diferentes ensayos, me parece que me falta algo. Mi modelo es el siguiente:

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

Funciona bien, se ajusta bien a los datos y los resultados tienen sentido. Ahora para obtener predicciones he intentado lo siguiente:

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 excluyo, digamos, Francia, y solo incluyo países en los que cluster = "OneTransition, Middle Income" obtengo un error diferente

# 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

En este ejemplo, el problema se debe a cluster = "One Transition, Middle Income" todo el tiempo.

No entiendo por qué esto es un problema. Si quiero que predic () funcione, tengo que incluir todas las variables del modelo, ¿verdad? Obviamente, los datos de entrada en la llamada del modelo no incluirán el factor establecido en los mismos valores para todos los casos. Sin embargo, si quiero obtener predicciones solo para un subconjunto de datos, o para nuevas observaciones, podría interesarme solo en los casos en que algún factor siempre sea el mismo. ¿Tiene sentido? ¿Cómo puedo obtener predicciones en ese caso?

Antonio Pedro Ramos
fuente
Sospecho que esto se debe a que cuando ve dos factores, uno el subconjunto del otro, R ve dos factores no relacionados. Solo una corazonada: intente iniciar una nueva sesión de R, escribir options(stringsAsFactors = FALSE)y luego ejecutar su código. Eso evitaría que su original test.predtenga sus propios factores.
Matt Parker,
Gracias Matt, pero aún no funciona. En realidad soy un rompecabezas. Debe ser algún tipo de error.
Antonio Pedro Ramos
Solo una solución, por ejemplo, un factor de 3 niveles A, B, C puede hacer una predicción para el nivel A con datos de 100 A, 1 B y 1 C.
Verbal

Respuestas:

8

Gracias por proporcionar los datos para poder realizar algunos diagnósticos. En realidad, este es un error épico de predict.lme. Su factorstiene más niveles en su inicial de datos (por ejemplo, si tiene más de 4 países) que en sus nuevos datos. Una línea de código específicamente hace que los niveles no utilizados se descarten para que termines con matrices de diferentes dimensiones, de donde elnon-conformable arguments

Quité esa línea y puse el código aquí .

En R puedes hacer

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

Esto registra una nueva función predict.lmeque se invocará en lugar de la del paquete nlmey puede ejecutar su código. Al menos funcionó para mí.

Advertencia: El código publicado y el método no son un reemplazo ni una corrección real de errores del paquete. La función parcheada no se ha probado más allá de su capacidad para ejecutar el bit de código del OP.

gui11aume
fuente
En realidad lo hace. Tengo country.xen ambos. El jurado aún está deliberando.
Antonio Pedro Ramos
Sí ... eso es correcto. Lo siento por eso. Parece que algunos de sus tipos de datos no son los mismos en su entrada inicial y sus nuevos datos. Esto será muy difícil de hacer sin los datos. Si no es confidencial y no es demasiado grande, ¿podría guardar la sesión R y ponerla en algún lugar (o enviármela por correo)?
gui11aume
Muchas gracias. ¿Tiene un correo electrónico para que le envíe un código de muestra y datos?
Antonio Pedro Ramos
Solo una pregunta rápida: ¿esta versión funciona bien levels=1pero no para level=0?
Antonio Pedro Ramos