Escrevendo a equação matemática para um modelo de efeitos mistos multinível

A questão CV

Estou tentando fornecer (a) representações matemáticas detalhadas e concisas de um modelo de efeitos mistos. Estou usando o lme4pacote em R. Qual é a representação matemática correta para o meu modelo?

---

Dados, Questão Científica e Código R

Meu conjunto de dados consiste em espécies em diferentes regiões. Estou testando se a prevalência de uma espécie muda no tempo que leva à extinção (extinções não são necessariamente permanentes; podem recolonizar), ou após uma colonização.

lmer(prevalence ~ time + time:type + (1 + time + type:time | reg) + (1 + time + type:time | reg:spp))

• Prevalência é a proporção de estratos ocupados por uma espécie em uma região / ano
• O tempo é uma variável contínua que indica o tempo para extinção ou colonização; é sempre positivo
• Tipo é uma variável categórica com dois níveis. Esses dois níveis são "-" e "+". Quando o tipo é -, é uma colonização (nível padrão). Quando o tipo é +, é uma extinção.
• Reg é uma variável categórica com nove níveis, indicando a região
• Spp é uma variável categórica; o número de níveis varia entre as regiões e varia entre 48 e 144.

Em palavras: variável de resposta é prevalência (proporção de estratos ocupados). Os efeitos fixos incluíram 1) e interceptar, 2) tempo do evento e 3) a interação entre o tempo do evento e o tipo de evento (colonização ou extinção). Cada um desses três efeitos fixos variava aleatoriamente entre as regiões. Dentro de uma região, cada um dos efeitos variou aleatoriamente entre as espécies.

Estou tentando descobrir como escrever a equação matemática para o modelo. Eu acho que entendo o que está acontecendo no código R (embora eu tenha certeza que tenho algumas lacunas de conhecimento, e espero que escrever a expressão matemática formal melhore minha compreensão).

Eu pesquisei bastante na web e nesses fóruns. Encontrei toneladas de informações úteis, com certeza (e talvez eu vincule algumas delas em uma edição desta pergunta). No entanto, eu não consegui encontrar a "Rosetta Stone" do código R traduzida para matemática (estou mais à vontade com o código) que realmente me ajudaria a confirmar que tenho essas equações corretas. Na verdade, eu sei que já existem algumas lacunas, mas vamos chegar a isso.

---

Minha tentativa

A forma básica de um modelo de efeitos mistos, em notação matricial, é (pelo meu entendimento):

$$Y = X \beta + Z \gamma + \epsilon$$

$$X = \begin{bmatrix} 1 & \Delta t & \Delta t_{+} \\ \vdots & \vdots & \vdots \\ 1 & \Delta t_n & \Delta t_{+,n} \end{bmatrix}$$ Z= [ 1 I ( r 1 ) Δ t I ( r 1 ) Δ t + I ( r 1 ) … 1 I ( r 9 ) Δ t I ( r 9 ) Δ t + I ( r 9 ) ⋮ ⋮ ⋮ ⋱ $$\beta^{'} = \begin{bmatrix} \beta_0 & \beta_1 & \beta_2 \end{bmatrix}$$ $$Z = \begin{bmatrix} 1 I(r_1) & \Delta t I(r_1) & \Delta t_{+} I(r_1) & \dots & 1 I(r_9) & \Delta t I(r_9) & \Delta t_{+} I(r_9) \\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots & \vdots \\ 1 I(r_{1,n}) & \Delta t_n I(r_{1,n}) & \Delta t_{+,n} I(r_{1,n}) & \dots & 1 I(r_{9,n}) & \Delta t I(r_{9,n}) & \Delta t_{+,n} I(r_{9,n}) \\ \end{bmatrix}$$ ε~N(0,Σ $$\gamma^{'} = \begin{bmatrix} \gamma_{0,1} & \gamma_{1,1} &\gamma_{2,1} & \dots & \gamma_{0,9} & \gamma_{1,9} &\gamma_{2,9} \end{bmatrix}$$ $$\epsilon \sim \mathcal{N}(0,\Sigma)$$

• Δ t Δ t $X$ é a matriz de design para os efeitos fixos, é o tempo após a colonização ( ) e é o tempo após a extinção ( )$\Delta t$time$\Delta t_{+}$time:type
• $Z$ é a matriz de design dos efeitos aleatórios (nível 1?), I () é a função indicadora que fornece 1 se a amostra pertence à região designada e 0 caso contrário, r é indexado para indicar uma das nove regiões.
• $\beta$ e contêm parâmetros$\gamma$
• $\epsilon$ são erros; Não tenho muita certeza de como explicar , embora perceba que uma dessas matrizes de variância / covariância expressará covariâncias entre inclinações e interceptações, por exemplo,$\Sigma$

Supondo que as coisas até agora estejam corretas, isso significa que eu sou bom no nível superior. No entanto, explicar a variação específica da espécie nos parâmetros, que está aninhada em cada região, me surpreendeu ainda mais.

Mas fiz uma rachadura em algo que talvez faça sentido ...

Cada um dos parâmetros em é derivado de uma combinação linear de preditores e parâmetros específicos da espécie em uma região. Para cada região, existem 3 linhas de, correspondentes às 3 variáveis ​​preditoras. Cada pode ser expresso individualmente como$\gamma$$\gamma$

• $\gamma_{p,r} = U_{p,r} b_{p,r} + \eta_{p,r}$
  • onde é uma matriz de desenho específico para a região e preditor , é um 1 por S matriz de parâmetros para a região (riqueza na região = , por exemplo, 48 ou 144), e é uma matriz de termos de erro$U_{p,r}$$r$$p$$b_{p,r}$$S$$\eta_{p,r}$

Especificamente, para uma determinada região, cada um dos seria:$\gamma_{p,r}$

$$\gamma_{0,r} = U_{0,r} b_{0,r} + \eta_{0,r}$$ $$\gamma_{0,r} = \begin{bmatrix} 1 I(s_1) \dots 1 I(s_S) \end{bmatrix} + \begin{bmatrix} b_{0,1}\\ \vdots \\ b_{0,S} \end{bmatrix} + \eta_{0,r}$$ $$\gamma_{1,r} = U_{1,r} b_{1,r} + \eta_{1,r}$$ $$\gamma_{1,r} = \begin{bmatrix} \Delta t I(s_1) \dots \Delta t I(s_S) \end{bmatrix} + \begin{bmatrix} b_{1,1}\\ \vdots \\ b_{1,S} \end{bmatrix} + \eta_{1,r}$$ $$\gamma_{2,r} = U_{2,r} b_{2,r} + \eta_{2,r}$$ $$\gamma_{2,r} = \begin{bmatrix} \Delta t_+ I(s_1) \dots \Delta t_+ I(s_S) \end{bmatrix} + \begin{bmatrix} b_{2,1}\\ \vdots \\ b_{2,S} \end{bmatrix} + \eta_{2,r}$$

Isso seria repetido para cada região. Em seguida, , como . Embora, talvez em vez de , exista outra letra, como , que é comumente usada.$\eta \sim \mathcal{N}(0,\Sigma_{\eta})$$\epsilon$$\Sigma$$G$

---

Edit: outras perguntas / respostas que foram úteis

• Este Q / A foi bom, mas não escreveu as coisas na forma de matriz completa

Se eu entendi o código corretamente, por que não simplesmente escrever algo como

$$y_{i} = \Big(\alpha + \nu_{j[i]}^{(\alpha)} + \eta_{k[i]}^{(\alpha)}\Big) + \Big(\beta + \nu_{j[i]}^{(\beta)} + \eta_{k[i]}^{(\beta)}\Big)T_{i} + \Big(\delta + \nu_{j[i]}^{(\delta)} + \eta_{k[i]}^{(\delta)}\Big)(T_{i} * Z_{i}) + \epsilon_i$$ $$\begin{aligned} \Big[\nu_{j}^{(\alpha)}, \nu_j^{(\beta)}, \nu_j^{(\delta)}\Big] &\sim \text{Multi-Normal}(\mathbf 0, \boldsymbol \Sigma_\nu) \\ \Big[\eta_{j}^{(\alpha)}, \eta_j^{(\beta)}, \eta_j^{(\delta)}\Big] &\sim \text{Multi-Normal}(\mathbf 0, \boldsymbol \Sigma_\eta)\\ \epsilon_i & \sim \text{Normal}(0, \sigma_\epsilon) \end{aligned}$$ $$y_{i} = \alpha_{j[i],k[i]} + \beta_{j[i],k[i]}T_{i} + \delta_{j[i],k[i]}(T_i * Z_i) + \epsilon_i$$ $$\begin{aligned} \alpha_{j[i],k[i]} &= \alpha + \nu_{j}^{(\alpha)} + \eta_{k}^{(\alpha)} \\ \beta_{j[i],k[i]}&=\beta + \nu_{j}^{(\beta)} + \eta_{k}^{(\beta)}\\ \delta_{j[i],k[i]}&=\delta + \nu_{j}^{(\delta)} + \eta_{k}^{(\delta)}\\ \end{aligned}$$