O que são estruturas G de estrutura R em um glmm?

Eu tenho usado o MCMCglmmpacote recentemente. Estou confuso com o que é referido na documentação como estrutura R e estrutura G. Eles parecem estar relacionados aos efeitos aleatórios - em particular especificando os parâmetros para a distribuição anterior sobre eles, mas a discussão na documentação parece assumir que o leitor sabe quais são esses termos. Por exemplo:

lista opcional de especificações anteriores com 3 elementos possíveis: R (estrutura R) G (estrutura G) e B (efeitos fixos) ............ Os anteriores para as estruturas de variação (R e G ) são listas com as (co) variações esperadas (V) e o parâmetro do grau de crença (nu) para o inverso-Wishart

... tirado daqui .

EDIT: Observe que reescrevi o restante da pergunta após os comentários de Stephane.

Alguém pode esclarecer o que são a estrutura R e a estrutura G, no contexto de um modelo de componentes de variação simples em que o preditor linear é com e

β_{0 0} + e_{0 0 Eu j} + {você}_{0 0 j}

$\beta_0 + e_{0ij} + u_{0j}$

e_{0 i j} \sim N (0, σ_{0 e}^{2})

$e_{0ij} \sim N(0,\sigma_{0e}^2)$

u_{0 j} \sim N (0, σ_{0 u}^{2})

$u_{0j} \sim N(0,\sigma_{0u}^2)$

Fiz o exemplo a seguir com alguns dados que acompanham MCMCglmm

> require(MCMCglmm)
> require(lme4)
> data(PlodiaRB)
> prior1 = list(R = list(V = 1, fix=1), G = list(G1 = list(V = 1, nu = 0.002)))
> m1 <- MCMCglmm(Pupated ~1, random = ~FSfamily, family = "categorical", 
+ data = PlodiaRB, prior = prior1, verbose = FALSE)
> summary(m1)


 G-structure:  ~FSfamily

         post.mean l-95% CI u-95% CI eff.samp
FSfamily    0.8529   0.2951    1.455      160

 R-structure:  ~units

      post.mean l-95% CI u-95% CI eff.samp
units         1        1        1        0

 Location effects: Pupated ~ 1 

            post.mean l-95% CI u-95% CI eff.samp  pMCMC    
(Intercept)   -1.1630  -1.4558  -0.8119    463.1 <0.001 ***
---

> prior2 = list(R = list(V = 1, nu = 0), G = list(G1 = list(V = 1, nu = 0.002)))
> m2 <- MCMCglmm(Pupated ~1, random = ~FSfamily, family = "categorical", 
+ data = PlodiaRB, prior = prior2, verbose = FALSE)
> summary(m2)


 G-structure:  ~FSfamily

         post.mean l-95% CI u-95% CI eff.samp
FSfamily    0.8325   0.3101    1.438    79.25

 R-structure:  ~units

      post.mean l-95% CI u-95% CI eff.samp
units    0.7212  0.04808    2.427    3.125

 Location effects: Pupated ~ 1 

            post.mean l-95% CI u-95% CI eff.samp  pMCMC    
(Intercept)   -1.1042  -1.5191  -0.7078    20.99 <0.001 ***
---

> m2 <- glmer(Pupated ~ 1+ (1|FSfamily), family="binomial",data=PlodiaRB)
> summary(m2)
Generalized linear mixed model fit by the Laplace approximation 
Formula: Pupated ~ 1 + (1 | FSfamily) 
   Data: PlodiaRB 
  AIC  BIC logLik deviance
 1020 1029   -508     1016
Random effects:
 Groups   Name        Variance Std.Dev.
 FSfamily (Intercept) 0.56023  0.74849 
Number of obs: 874, groups: FSfamily, 49

Fixed effects:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  -0.9861     0.1344  -7.336  2.2e-13 ***

Portanto, com base nos comentários de Stephane, acho que a estrutura G é para . Mas os comentários também dizem que a estrutura R é para mas isso não parece aparecer na saída. $\sigma_{0u}^2$ $\sigma_{0e}^2$ lme4

Observe que os resultados de lme4/glmer()são consistentes com os dois exemplos do MCMC MCMCglmm.

Então, é a estrutura R para e por que isso não aparece na saída ? $\sigma_{0e}^2$ lme4/glmer()

r bayesian mixed-model lme4-nlme Joe King
fonte

Com a terminologia SAS (mas é possivelmente uma terminologia mais comum), a matriz G é a matriz de variação dos efeitos aleatórios e a matriz R é a matriz de variação dos "termos de erros" (no seu caso, talvez seja o resíduo estimado variance ?)

σ_{0 e}^{2}

$\sigma_{0e}^2$

Stéphane Laurent

@ StéphaneLaurent obrigado. Eu me perguntei se poderia ser estimado mas quando eu aprendi sobre o modelo linear generalizado, lembro que não é estimado - apenas "desvio" é calculado (como em ). Talvez esteja faltando alguma coisa?

σ_{0 e}^{2}

$\sigma_{0e}^2$

σ_{0 e}^{2}

$\sigma_{0e}^2$ lme4

31812 Joe King

talvez o sentido da variação residual não está claro quando a família distribuição não é o Gaussian

Stéphane Laurent

@ Stéphane Laurent Sim! Por favor, veja o meu comentário para a resposta de Michael um minuto atrás - para o resultado binário, que deve ser fixado (como em meus modelos no meu OP)

Joe King

Quando você possui um modelo ME / Multinível, existem várias variações. Imagine o caso mais simples: . Há variação nas interceptações e no termo de erro . é freqüentemente usado para a matriz var-covar dos efeitos aleatórios (neste caso, um escalar, ) & é para a matriz var-covar das variâncias residuais após contabilizar a aleatória fixa e desse cluster efeitos Geralmente é concebido como uma matriz diagonal de 's. Além disso, ambos os discos são considerados normais multivariados com média = 0.

Y_{i} = β_{0} + β_{1} X + b_{i} + ε_{i}

$Y_i=\beta_0+\beta_1X+b_i+\varepsilon_i$

b_{i}

$b_i$

ε_{i}

$\varepsilon_i$

G

$G$

σ_{b}^{2}

$\sigma^2_b$

R_{i}

$R_i$

ε_{i}

$\varepsilon_i$

σ^{2}

$\sigma^2$

- Reinstate Monica

Respostas:

Eu preferiria postar meus comentários abaixo como um comentário, mas isso não seria suficiente. Essas são perguntas, e não uma resposta (semelhante a @gung, não me sinto forte o suficiente sobre o assunto).

Estou com a impressão de que o MCMCglmm não implementa um glmm bayesiano "verdadeiro". O verdadeiro modelo bayesiano é descrito na seção 2 deste artigo . Da mesma forma que o modelo freqüentista, um tem e há um requisito prévio no parâmetro de dispersão , além dos parâmetros fixos e da variação "G" do efeito aleatório . $g(E(y \mid u)) = X\beta + Zu$ $\phi_1$ $\beta$ $u$

Porém, de acordo com essa vinheta do MCMCglmm , o modelo implementado no MCMCglmm é dado por , e não envolve o parâmetro de dispersão . Não é semelhante ao modelo freqüentista clássico. $g(E(y \mid u,e)) = X\beta + Zu + e$ $\phi_1$

Portanto, não ficaria surpreso que não exista um análogo de com o glmer. $\sigma_e$

Por favor, peça desculpas por esses comentários grosseiros, apenas dei uma olhada rápida sobre isso.

Stéphane Laurent
fonte

Obrigado. Esse tópico deveria ser difícil, porque eu estou achando bastante difícil? Acho que estou satisfeito com o significado da estrutura R e G agora. Eu ainda estou confuso sobre a falta de com e estou muito curioso sobre o seu comentário de que não é verdadeiramente Bayesian. Sinceramente, não posso dizer que entendi todo o artigo ao qual você vinculou e também estou lutando com partes da vinheta, mas apenas pela perspectiva do meu exemplo, acredito que o parâmetro de dispersão deve ser constante (porque o exemplo é binomial). O que estou perdendo ?

σ_{e}

$\sigma_e$ glmerMCMCglmmMCMCglmm

ϕ_{1}

$\phi_1$

Joe King

Desculpe, minhas palavras não foram totalmente apropriadas. O MCMCglmm é verdadeiramente bayesiano, mas não implementa exatamente o glmm clássico (eu acho). Além disso, você deve estar ciente de que é difícil definir antecedentes que produzam uma inferência nos componentes de variação próximos à inferência freqüentista.

Stéphane Laurent

Obrigado novamente. Nos meus estudos, descobri que posso usar a distribuição inversa-wishart padrão para componentes de variação MCMCglmmusando vários parâmetros, e os intervalos de 95% credíveis sempre contêm o valor de variação para os efeitos aleatórios estimados, glmerentão achei que isso era razoável. , mas como devo interpretar esse caso, que pode não ser típico, onde o resultado é que os MCMCglmmintervalos não são muito sensíveis à escolha do anterior? Talvez eu deva fazer uma nova pergunta sobre isso?

31812 Joe King

Talvez você tenha um grande tamanho de amostra? Em relação à sua pergunta inicial, tenho a impressão de que, pelo menos para o caso binomial, o modelo glmer é equivalente ao modelo MCMCglmm com . O que acontece se você definir um prior em altamente concentrado em ?

σ_{e} = 0

$\sigma_e=0$

σ_{e}

$\sigma_e$

0

$0$

Stéphane Laurent

Sim, tenho um tamanho de amostra bastante grande: 50.000 observações em 225 clusters (meus próprios dados, não o exemplo da minha pergunta). Quando defino um prior muito concentrado próximo de zero em , definindo V = 0,01 e nu = 100, obtenho 0,25 (IC: 0,16, 0,29) para e 0,53 (0,38, 0,73) para . Quando defino um anterior menos informativo, com V = 10 e nu = 0,01, obtenho 0,18 (0,12, 0,23) e 0,49 (0,34, 0,63), respectivamente. Isso se compara a 0,51 de . Eu até tentei um plano inadequado anterior, que deu 0,10 (0,08, 0,13) e 0,47 (0,25, 0,68).

σ_{e}

$\sigma_e$

σ_{e}

$\sigma_e$

σ_{u}

$\sigma_u$ glmer

31712 Joe King

Estou atrasado para o jogo, mas algumas notas. A estrutura é a estrutura residual. No seu caso, a "estrutura" possui apenas um único elemento (mas isso não precisa ser o caso). Para a variável de resposta gaussiana, a variação residual, é normalmente estimada. Para resultados binários, é mantido constante. Devido à forma como o MCMCglmm está configurado, você não pode corrigi-lo em zero, mas é relativamente padrão corrigi-lo em (também verdadeiro para um modelo probit). Para dados de contagem (por exemplo, com uma distribuição de poisson), você não o corrige e isso estima automaticamente um parâmetro de superdispersão essencialmente. $\mathbf{R}$ $\sigma^{2}_{e}$ $1$

A estrutura é a estrutura de efeitos aleatórios. Novamente no seu caso, apenas uma interceptação aleatória, mas se você tivesse vários efeitos aleatórios, eles formariam uma matriz de variância-covariância, . $\mathbf{G}$ $\mathbf{G}$

Uma observação final, como a variação residual não é fixada em zero, as estimativas não serão iguais às de glmer. Você precisa revalorizá-los. Aqui está um pequeno exemplo (sem usar efeitos aleatórios, mas generaliza). Observe como a variação da estrutura R é fixada em 1.

# example showing how close the match is to ML without separation
m2 <- MCMCglmm(vs ~ mpg, data = mtcars, family = "categorical",
  prior = list(
    B = list(mu = c(0, 0), V = diag(2) * 1e10),
    R = list(V = 1, fix = 1)),
  nitt = 1e6, thin = 500, burnin = 10000)
summary(m2)

Aqui está a constante de redimensionamento da família binomial:

k <- ((16*sqrt(3))/(15*pi))^2

Agora divida a solução por ela e obtenha os modos posteriores

posterior.mode(m2$Sol/(sqrt(1 + k)))

O que deve ser bastante próximo do que obtemos glm

summary(glm(vs ~mpg, data = mtcars, family = binomial))

Joshua
fonte

você saberia como especificar a heterocedasticidade no nível um no MCMCglmm? Essa é a estrutura R? Qual é a sintaxe então?

precisa saber é o seguinte

@ Josué, você pode explicar a "constante de redimensionamento para a família binomial"? PS: Para sementes 123, recebo (com a correção) dos m2valores -8.164e 0.421; e a partir glmdos valores -8.833e 0.430.

precisa saber é o seguinte

A constante de redimensionamento pode ser encontrada em Diggle et. al. ( amazon.de/Analysis-Longitudinal-Oxford-Statistical-Science/dp/… ) - de acordo com cran.r-project.org/web/packages/MCMCglmm/vignettes/… eq. 2,14 na página 47.

Qaswed 11/09/16