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

Question

Estoy haciendo un estudio de simulación que requiere de arranque de estimaciones obtenidas a partir de un lineales generalizados mixtos modelo (en realidad, el producto de dos estimaciones por efectos fijos, uno de un GLMM y uno de un LMM). Para hacer el estudio también requeriría alrededor de 1000 simulaciones con 1000 o 1500 bootstrap repeticiones cada vez. Esto toma una cantidad significativa de tiempo en mi equipo (muchos días).

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

Para ser más específicos, tengo los sujetos que se mide repetidamente en tres maneras, dando lugar a las variables X, M, e y, donde X y M son continuos y y es binario. Tenemos dos ecuaciones de regresión $$M=\alpha_0+\alpha_1X+\epsilon_1$$ $$Y^*=\beta_0+\beta_1X+\beta_2M+\epsilon_2$$ donde Y$^*$ es el trasfondo latente variable continua para $Y$ y los errores no son iid.
La estadística queremos bootstrap es $\alpha_1\beta_2$. Por lo tanto, cada bootstrap de replicación requiere el ajuste de 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 tengo la misma estimación para $\alpha_1$ si se me acaba de encajar como un modelo lineal, por lo que ahorra algo de tiempo, pero el mismo truco no funciona para $\beta_2$.

¿Necesito comprar un equipo más rápido? :)

Answer 1

Debe ayudar a especificar los valores partidos, aunque es difícil saber cuánto. Como haces la simulación Bootstrap, usted debe saber los valores 'true', las estimaciones no arranque o ambos. Trate de usarlos de la start = opción de glmer .

Usted también podría considerar buscar en si las tolerancias para declarar convergencia son más estrictas de lo necesario. No tengo claro cómo alterar desde el lme4 documentación aunque.

Answer 2

Dos otras posibilidades también considerar, antes de comprar un equipo nuevo.

1. La computación paralela - bootstrapping es fácil de ejecutar en paralelo. Si su ordenador está razonablemente nuevo, usted probablemente tiene cuatro núcleos. Echa un vistazo a la de los núcleos de la biblioteca en R.
2. El Cloud computing es también una posibilidad y razonablemente barato. Tengo colegas que han utilizado la nube de amazon para la ejecución de secuencias de comandos de R. Encontraron que era muy rentable.

Answer 3

Que podría ser un equipo más rápido. Pero aquí hay un truco que puede funcionar.

Generar un simultation de $Y^*$, pero sólo condicional en $Y$, a continuación, sólo hacer OLS o LMM en la simulación de la $Y^*$ valores.

Suponiendo que su función de enlace es $g(.)$. esto explica la forma en que se obtiene de la probabilidad de $Y=1$ $Y^*$ valor, y es más probable que la función logística $g(z)=log \Big(\frac{z}{1-z}\Big)$.

Así que si usted asume una bernouli distribución de muestreo para $Y\rightarrow Y\sim Bernoulli(p)$, y, a continuación, utilizar la jeffreys previo para la probabilidad de obtener una beta posterior para $p\sim Beta(Y_{obs}+\frac{1}{2},1-Y_{obs}+\frac{1}{2})$. La simulación de este debe ser igual que la iluminación, y si no, entonces usted necesita un equipo más rápido. Además, las muestras son independientes, por lo que no tienen necesidad de controlar la "convergencia" de los diagnósticos como en MCMC, y probablemente no necesita tantas muestras - 100 puede funcionar bien para tu caso. Si usted tiene binomio $Y's$, a continuación, sólo tiene que sustituir el $1$ arriba en la parte posterior con $n_i$, el número de ensayos de la binomial para cada una de las $Y_i$.

Así que usted tiene un conjunto de valores simulados $p_{sim}$. A continuación, aplicar la función de enlace para cada uno de estos valores, para obtener el $Y_{sim}=g(p_{sim})$. Ajuste un LMM a $Y_{sim}$, que es probablemente más rápido que el GLMM programa. Básicamente, puede ignorar el original de los valores binarios (pero no eliminar por ellos!!!), y sólo trabajar con la simulación de "matrix" ($N\times S$donde $N$ es el tamaño de la muestra, y $S$ es el número de simulaciones).

Así, en su programa, yo sustituiría a la $gmler()$ función con la $lmer()$ función, y $Y$ con un solo simultation, entonces Usted podría crear una especie de bucle que se aplica el $lmer()$ función para cada simulación, y luego toma el promedio como la estimación de $b$. Algo como $$a=\dots$$ $$b=0$$ $$do \ s=1,\dots,S$$ $$b_{est}=lmer(Y_s\dots)$$ $$b=b+\frac{1}{s}(b_{est}-b)$$ $$end$$ $$return(a*b)$$

Déjame saber si tengo que explicar las cosas un poco más claras