¿Cómo puedo acelerar el cálculo de los efectos fijos en un GLMM?

Estoy haciendo un estudio de simulación que requiere estimaciones de arranque obtenidas de un modelo mixto lineal generalizado (en realidad, el producto de dos estimaciones para efectos fijos, una de un GLMM y otra de un LMM). Hacer bien el estudio requeriría alrededor de 1000 simulaciones con 1000 o 1500 repeticiones de arranque cada vez. Esto lleva una cantidad significativa de tiempo en mi computadora (muchos días).

How can I speed up the computation of these fixed effects?

Para ser más específico, tengo sujetos que se miden repetidamente de tres maneras, dando lugar a las variables X, M e Y, donde X e M son continuos e Y es binario. Tenemos dos ecuaciones de regresión donde Y es la variable continua latente subyacente para y los errores no son iid. La estadística que queremos arrancar es

M = α_{0} + α_{1} X + ϵ_{1}

$M=\alpha_0+\alpha_1X+\epsilon_1$

Y^{*} = β_{0} + β_{1} X + β_{2} M + ϵ_{2}

$Y^*=\beta_0+\beta_1X+\beta_2M+\epsilon_2$

^{*}

$^*$

Y

$Y$

α_{1} β_{2}

$\alpha_1\beta_2$ . Por lo tanto, cada replicación de arranque requiere ajustar un LMM y un GLMM. Mi código R es (usando lme4)

    stat=function(dat){
        a=fixef(lmer(M~X+(X|person),data=dat))["X"]
        b=fixef(glmer(Y~X+M+(X+M|person),data=dat,family="binomial"))["M"]
        return(a*b)
    }

Me doy cuenta de que obtengo la misma estimación para si solo lo ajusto como un modelo lineal, por lo que ahorra algo de tiempo, pero el mismo truco no funciona para . $\alpha_1$ $\beta_2$

¿Solo necesito comprar una computadora más rápida? :)

r mixed-model BR
fuente

@BR ¿cuál es el cuello de la botella aquí? Básicamente lo que está tomando tiempo Rprof.

suncoolsu

Una forma es ignorar la parte "mixta" del GLMM. Sólo una ajustarse GLM ordinaria, las estimaciones no cambiarán mucho, pero sus errores estándar probable que

probabilityislogic

@probabilityislogic. Además de su comentario, también creo que si la respuesta sería muy diferente depende del tamaño del grupo y del comportamiento individual en el grupo. Como dicen Gelman y Hill: los resultados del modelo de efectos mixtos estarían entre agrupar y no agrupar. (Obv. Esto es para los modelos jerárquicos bayesianos, pero los modelos mixtos son una forma clásica de hacer lo mismo.)

suncoolsu

@probabilityislogic: Eso funciona para LMM, pero parece fallar para GLMM (lo que significa que ejecuté modelos con y sin la M adicional en los mismos datos y terminé con resultados significativamente diferentes). A menos, por supuesto, que haya un error en la implementación de glmer.

@suncoolsu: el cuello de botella es la estimación del GLMM, que puede tardar un par de segundos (especialmente con varios efectos aleatorios). Pero haz eso 1000 * 1000 veces, y eso es 280 horas de cálculo. Ajustar un GLM toma aproximadamente 1/100 del tiempo.

Respuestas:

Debería ayudar a especificar valores iniciales, aunque es difícil saber cuánto. Mientras realiza simulación y bootstrapping, debe conocer los valores "verdaderos" o las estimaciones sin bootstrap o ambas. Intente usar aquellos en la start =opción de glmer.

También podría considerar investigar si las tolerancias para declarar la convergencia son más estrictas de lo que necesita. Sin lme4embargo, no tengo claro cómo modificarlos de la documentación.

una parada
fuente

Considere también otras dos posibilidades, antes de comprar una computadora nueva.

Computación en paralelo: el arranque es fácil de ejecutar en paralelo. Si su computadora es razonablemente nueva, probablemente tenga cuatro núcleos. Echa un vistazo a la biblioteca multinúcleo en R.
La computación en la nube también es una posibilidad y es razonablemente barata. Tengo colegas que han usado la nube de Amazon para ejecutar scripts R. Descubrieron que era bastante rentable.

csgillespie
fuente

Gracias por la respuesta. De alguna manera, he pasado por alto el hecho de que tengo dos núcleos (mi computadora no es muy nueva). Debería haber mirado multinúcleo hace mucho tiempo.

Posiblemente podría ser una computadora más rápida. Pero aquí hay un truco que puede funcionar.

$Y^*$ $Y$ $Y^*$

$g(.)$ $Y=1$ $Y^*$ $g(z)=log \Big(\frac{z}{1-z}\Big)$

$Y\rightarrow Y\sim Bernoulli(p)$ $p\sim Beta(Y_{obs}+\frac{1}{2},1-Y_{obs}+\frac{1}{2})$ $Y's$ $1$ $n_i$ $Y_i$

$p_{sim}$ $Y_{sim}=g(p_{sim})$ $Y_{sim}$ $N\times S$ $N$ $S$

$gmler()$ $lmer()$ $Y$ $lmer()$ $b$

una = ...

$a=\dots$

si = 0 0

$b=0$

re o s = 1, ..., S

$do \ s=1,\dots,S$

{si}_{mi s t} = l metro mi r (Y_{s} ...)

$b_{est}=lmer(Y_s\dots)$

si = si + \frac{1}{s} ({si}_{mi s t} - si)

$b=b+\frac{1}{s}(b_{est}-b)$

mi norte re

$end$

r mi t tu r norte (una * si)

$return(a*b)$

Avíseme si necesito explicar algo un poco más claro

probabilidadislogica
fuente

Gracias por la respuesta, me llevará un poco digerirlo (y ya tengo planes para mi sábado por la noche). Es lo suficientemente diferente como para que no me quede claro si da la misma respuesta que el enfoque GLMM, pero necesito pensarlo más.