Como acelerar o cálculo dos efeitos fixos em um GLMM?

Estou fazendo um estudo de simulação que requer estimativas de bootstrapping obtidas de um modelo misto linear generalizado (na verdade, o produto de duas estimativas para efeitos fixos, um de um GLMM e outro de um LMM). Para fazer o estudo, seria necessário cerca de 1000 simulações com 1000 ou 1500 replicações de autoinicialização por vez. Isso leva uma quantidade significativa de tempo no meu computador (muitos dias).

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

Para ser mais específico, tenho sujeitos que são medidos repetidamente de três maneiras, dando origem às variáveis X, M e Y, onde X e M são contínuos e Y é binário. Temos duas equações de regressão que Y é a variável contínua latente subjacente para e os erros não são iid. A estatística que queremos inicializar é . Assim, cada replicação de autoinicialização requer o ajuste de um LMM e um GLMM. Meu código R é (usando lme4)

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$

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

Percebo que obtenho a mesma estimativa para se apenas o ajustar como modelo linear, economizando tempo, mas o mesmo truque não funciona para . $\alpha_1$ $\beta_2$

Eu só preciso comprar um computador mais rápido? :)

r mixed-model BR
fonte

@BR qual é o gargalo aqui? Basicamente, o que está demorando Rprof.

suncoolsu

Uma maneira é simplesmente ignorar a parte "mista" do GLMM. Apenas caber uma GLM comum, as estimativas não vai mudar muito, mas seus erros padrão provavelmente irá

probabilityislogic

@probabilityislogic. Além da sua observação, também acho que se a resposta diferiria muito depende do tamanho do grupo e do comportamento individual do grupo. Como Gelman e Hill dizem: os resultados do modelo de efeitos mistos estariam entre pool e sem pool. (. Obv isto é para modelos hierárquicos bayesianos, mas modelos mistos são uma maneira clássica de fazer o mesmo.)

suncoolsu

@probabilityislogic: Isso funciona para o LMM, mas parece falhar no GLMM (o que significa que eu executei modelos com e sem o M extra nos mesmos dados e acabei com resultados significativamente diferentes). A menos, é claro, que haja um erro na implementação do glmer.

@ suncoolsu: o gargalo é a estimativa do GLMM, que pode demorar alguns segundos (especialmente com vários efeitos aleatórios). Mas faça isso 1000 * 1000 vezes, e são 280 horas de computação. A instalação de um GLM leva cerca de 1/100 do tempo.

Respostas:

Deve ajudar a especificar valores iniciais, embora seja difícil saber quanto. Ao fazer a simulação e a inicialização, você deve conhecer os valores "verdadeiros" ou as estimativas sem inicialização ou ambas. Tente usar aqueles na start =opção de glmer.

Você também pode considerar se as tolerâncias para declarar convergência são mais rigorosas do que você precisa. Não estou claro como alterá-los da lme4documentação.

uma parada
fonte

Duas outras possibilidades também devem ser consideradas antes de comprar um novo computador.

Computação paralela - a inicialização é fácil de executar em paralelo. Se o seu computador for razoavelmente novo, você provavelmente possui quatro núcleos. Dê uma olhada na biblioteca multicore em R.
A computação em nuvem também é uma possibilidade e razoavelmente barata. Tenho colegas que usaram a nuvem amazon para executar scripts R. Eles descobriram que era bastante rentável.

csgillespie
fonte

Obrigado pela resposta. De alguma forma, eu negligenciei o fato de ter dois núcleos (meu computador não é muito novo). Eu deveria ter olhado multicore há muito tempo.

Poderia ser um computador mais rápido. Mas aqui está um truque que pode 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$

a = \dots

$a=\dots$

b = 0

$b=0$

d o s = 1, \dots, S

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

b_{e s t} = l m e r (Y_{s} \dots)

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

b = b + \frac{1}{s} (b_{e s t} - b)

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

e n d

$end$

r e t u r n (a * b)

$return(a*b)$

Deixe-me saber se preciso explicar algo um pouco mais claro

probabilityislogic
fonte

Obrigado pela resposta, vou demorar um pouco para digeri-la (e eu já tenho planos para minha noite de sábado). É diferente o suficiente para que não esteja claro se ela dá a mesma resposta que a abordagem GLMM, mas preciso pensar mais sobre isso.