Co to jest struktura R struktura G w glmm?

16

MCMCglmmOstatnio korzystam z pakietu. Jestem zdezorientowany tym, co w dokumentacji nazywane jest strukturą R i strukturą G. Wydaje się, że odnoszą się one do efektów losowych - w szczególności określają parametry wcześniejszego rozkładu na nich, ale dyskusja w dokumentacji wydaje się zakładać, że czytelnik wie, jakie są te warunki. Na przykład:

opcjonalna lista wcześniejszych specyfikacji zawierających 3 możliwe elementy: R (struktura R) G (struktura G) i B (efekty stałe) ............ Priory dla struktur wariancji (R i G ) to listy z oczekiwanymi (ko) wariancjami (V) i parametrem stopnia przekonania (nu) dla odwrotnego Wishart

... zabrane stąd .

EDYCJA: Pamiętaj, że napisałem resztę pytania po komentarzach Stephane'a.

Ktoś może rzucić światło na R Struktura i G-konstrukcji są, w ramach prostego modelu składniki wariancji w którym czynnikiem liniowych

β0+e0ij+u0j
o e0ijN(0,σ0e2) i u0jN(0,σ0u2)

Zrobiłem następujący przykład z niektórymi danymi, które pochodzą 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 ***

Na podstawie komentarzy Stephane'a myślę, że struktura G jest dla . Ale komentarze mówią również, że struktura R jest dla σ 2 0 e, ale wydaje się, że nie pojawia się to na wyjściu.σ0u2σ0e2lme4

Zauważ, że wyniki z lme4/glmer()są zgodne z oboma przykładami z MCMC MCMCglmm.

Czy więc struktura R dla i dlaczego nie pojawia się na wyjściu ?σ0e2lme4/glmer()

Joe King
źródło
1
Z terminologią SAS (ale jest to prawdopodobnie bardziej powszechna terminologia), macierz G jest macierzą wariancji efektów losowych, a macierz R jest macierzą wariancji „warunków błędów” (w twoim przypadku być może jest to szacunkowa wartość resztkowa wariancja ?)σ0e2
Stéphane Laurent
@ StéphaneLaurent dziękuję. Zastanawiałem się, czy można to oszacować ale kiedy po raz pierwszy dowiedziałem się o uogólnionym modelu liniowym, pamiętam, że σ 2 0 e nie jest oszacowane - obliczane jest tylko „odchylenie” (jak w przypadku ). Może czegoś mi brakuje?σ0e2σ0e2lme4
Joe King
1
może sens resztkowej wariancji nie jest jasny, gdy rodzina dystrybucyjna nie jest gaussowska
Stéphane Laurent
1
@ Stéphane Laurent Tak! Proszę zobaczyć mój komentarz do odpowiedzi Michaela przed chwilą - aby uzyskać wynik binarny, należy go naprawić (jak w moich modelach w moim OP)
Joe King
1
Gdy masz model ME / wielopoziomowy, istnieje kilka wariantów. Wyobraź sobie najprostszy przypadek: . Występuje wariancja w punktach przechwytywania bi i w błędzie ε i . G jest często stosowane do macierzy var-covar efektów losowych (w tym przypadku skalara, σ 2 b ) i R iYi=β0+β1X+bi+εibiεiGσb2Ri jest do macierzy var-covar resztkowych wariancji εipo uwzględnieniu stałych i losowych efektów tego klastra. Zwykle jest pomyślany jako macierz diagonalna . Ponadto oba dystanse są uważane za wielowymiarowe normalne w / średnia = 0. σ2
gung - Przywróć Monikę

Odpowiedzi:

8

Wolałbym zamieścić moje komentarze poniżej jako komentarz, ale to nie wystarczy. To są raczej pytania niż odpowiedź (po prostu na @gung, nie czuję się wystarczająco silny na ten temat).

Mam wrażenie, że MCMCglmm nie implementuje „prawdziwego” Bayesowskiego glmm. Prawdziwy model bayesowski opisano w rozdziale 2 tego artykułu . Podobnie jak w modelu częstym, jeden ma a parametr rozproszenia ϕ 1 jest wymagany dodatkowo oprócz stałych parametrów βg(E(yu))=Xβ+Zuϕ1β i wariancji „G” efekt losowy .u

Ale zgodnie z tą winietą MCMCglmm model zaimplementowany w MCMCglmm jest określony przez i nie obejmuje parametru dyspersji ϕ 1 . Nie jest podobny do klasycznego modelu częstokroć.g(E(yu,e))=Xβ+Zu+eϕ1

Dlatego nie zdziwiłbym się, że nie ma analogii z brokatem.σe

Przepraszam za te szorstkie komentarze, po prostu rzuciłem okiem na to.

Stéphane Laurent
źródło
σeglmerMCMCglmmMCMCglmmϕ1
Przepraszam, moje słowa nie były całkowicie odpowiednie. MCMCglmm jest naprawdę Bayesowski, ale nie do końca implementuje klasyczny glmm (tak myślę). Ponadto należy zdawać sobie sprawę z tego, że trudno jest ustawić priorytety, które prowadzą do wnioskowania na temat składników wariancji zbliżonych do wnioskowania częstych.
Stéphane Laurent,
Dzięki jeszcze raz. Podczas moich badań odkryłem, że mogę użyć domyślnego rozkładu odwrotnego życzenia dla składników wariancji MCMCglmmprzy użyciu różnych parametrów, a 95% wiarygodne przedziały zawsze zawierają wartość wariancji dla oszacowania efektów losowych przez, glmerwięc czułem, że było to uzasadnione , ale jak mam zinterpretować ten przypadek, co może nie być typowe, w wyniku czego MCMCglmminterwały nie są bardzo wrażliwe na wybór wcześniejszego? Może powinienem zadać nowe pytanie na ten temat?
Joe King,
σe=0σe0
σeσeσuglmer
11

Rσe21 (również prawdziwe dla modelu probit). W przypadku danych zliczania (np. Z rozkładem Poissona) nie naprawia się tego, a to zasadniczo automatycznie szacuje parametr nadmiernej dyspersji.

Gsol.

Ostatnia uwaga, ponieważ wariancja rezydualna nie jest ustalona na zero, oszacowania nie będą pasować do tych z glmer. Musisz je przeskalować. Oto mały przykład (nieużywanie efektów losowych, ale uogólnia). Zwróć uwagę, w jaki sposób wariancja struktury R jest ustalona na 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)

Oto stała przeskalowania dla rodziny dwumianowej:

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

Teraz podziel przez to rozwiązanie i uzyskaj tylne tryby

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

Co powinno być dość zbliżone do tego, co otrzymujemy glm

summary(glm(vs ~mpg, data = mtcars, family = binomial))
Jozuego
źródło
Czy wiesz, jak określić heteroskedastyczność na poziomie pierwszym w MCMCglmm? Czy to jest struktura R? Jaka jest zatem składnia?
Maxim.K.
@Joshua, czy możesz wyjaśnić „stałą przeskalowania dla rodziny dwumianowej”? PS: W przypadku nasion 123otrzymuję (z korektą) m2wartości -8.164i 0.421; oraz z glmwartości -8.833i 0.430.
Qaswed
Stałą przeskalowania można znaleźć w Diggle i in. glin. ( amazon.de/Analysis-Longitudinal-Oxford-Statistic-Science/dp/… ) - zgodnie z cran.r-project.org/web/packages/MCMCglmm/vignettes/… eq. 2.14 na stronie 47.
Qaswed