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?
fuente
options(stringsAsFactors = FALSE)
y luego ejecutar su código. Eso evitaría que su originaltest.pred
tenga sus propios factores.Respuestas:
Gracias por proporcionar los datos para poder realizar algunos diagnósticos. En realidad, este es un error épico de
predict.lme
. Sufactors
tiene 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
Esto registra una nueva función
predict.lme
que se invocará en lugar de la del paquetenlme
y 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.
fuente
country.x
en ambos. El jurado aún está deliberando.levels=1
pero no paralevel=0
?