Jak oszacować komponenty wariancji za pomocą lmera dla modeli z efektami losowymi i porównać je z wynikami lme

Przeprowadziłem eksperyment, w którym wychowałem różne rodziny pochodzące z dwóch różnych populacji źródłowych. Każdej rodzinie przydzielono jeden z dwóch zabiegów. Po eksperymencie zmierzyłem kilka cech u każdej osoby. Aby przetestować wpływ leczenia lub źródła, a także ich interakcji, zastosowałem liniowy model efektu mieszanego z rodziną jako czynnikiem losowym, tj.

lme(fixed=Trait~Treatment*Source,random=~1|Family,method="ML")

jak dotąd tak dobrze, teraz muszę obliczyć względne wariancje składowe, tj. procent zmienności, który jest wyjaśniony przez leczenie lub źródło, a także interakcję.

Bez efektu losowego mógłbym łatwo użyć sum kwadratów (SS) do obliczenia wariancji wyjaśnionej przez każdy czynnik. Ale dla modelu mieszanego (z oszacowaniem ML) nie ma SS, dlatego pomyślałem, że mógłbym użyć Leczenia i Źródła jako efektów losowych również do oszacowania wariancji, tj.

lme(fixed=Trait~1,random=~(Treatment*Source)|Family, method="REML")

Jednak w niektórych przypadkach lme nie jest zbieżne, dlatego użyłem lmera z pakietu lme4:

lmer(Trait~1+(Treatment*Source|Family),data=DATA)

Gdzie wyodrębniam wariancje z modelu za pomocą funkcji podsumowania:

model<-lmer(Trait~1+(Treatment*Source|Family),data=regrexpdat)
results<-VarCorr(model)
variances<-results[,3]

Otrzymuję te same wartości, co w przypadku funkcji VarCorr. Następnie używam tych wartości do obliczenia rzeczywistego procentu wariancji, przyjmując sumę jako całkowitą zmienność.

Walczę z interpretacją wyników z początkowego modelu Lme (z traktowaniem i źródłem jako ustalonymi efektami) oraz z modelem losowym do oszacowania składników wariancji (z traktowaniem i źródłem jako efektem losowym). Uważam w większości przypadków, że procent wariancji wyjaśniony przez każdy czynnik nie odpowiada znaczeniu ustalonego efektu.

Na przykład dla cechy HD, Początkowa wartość sugeruje tendencję do interakcji, a także znaczenie dla leczenia. Stosując procedurę wsteczną, stwierdzam, że leczenie ma tendencję zbliżoną do znaczącej. Jednak oceniając komponenty wariancji, stwierdzam, że Źródło ma największą wariancję, stanowiącą do 26,7% całkowitej wariancji.

Theme:

anova(lme(fixed=HD~as.factor(Treatment)*as.factor(Source),random=~1|as.factor(Family),method="ML",data=test),type="m")
                                      numDF denDF  F-value p-value
(Intercept)                                1   426 0.044523  0.8330
as.factor(Treatment)                       1   426 5.935189  0.0153
as.factor(Source)                          1    11 0.042662  0.8401
as.factor(Treatment):as.factor(Source)     1   426 3.754112  0.0533

A lmer:

summary(lmer(HD~1+(as.factor(Treatment)*as.factor(Source)|Family),data=regrexpdat))
Linear mixed model fit by REML 
Formula: HD ~ 1 + (as.factor(Treatment) * as.factor(Source) | Family) 
   Data: regrexpdat 
    AIC    BIC logLik deviance REMLdev
 -103.5 -54.43  63.75   -132.5  -127.5
Random effects:
 Groups   Name                                      Variance  Std.Dev. Corr                 
 Family   (Intercept)                               0.0113276 0.106431                      
          as.factor(Treatment)                      0.0063710 0.079819  0.405               
          as.factor(Source)                         0.0235294 0.153393 -0.134 -0.157        
          as.factor(Treatment)L:as.factor(Source)   0.0076353 0.087380 -0.578 -0.589 -0.585 
 Residual                                           0.0394610 0.198648                      
Number of obs: 441, groups: Family, 13

Fixed effects:
            Estimate Std. Error t value
(Intercept) -0.02740    0.03237  -0.846

Stąd moje pytanie, czy to, co robię, jest prawidłowe? A może powinienem użyć innego sposobu oszacowania wariancji wyjaśnionej przez każdy czynnik (tj. Leczenie, Źródło i ich interakcja). Na przykład, czy rozmiary efektów byłyby bardziej odpowiednie?

Jednym z powszechnych sposobów określania względnego udziału każdego czynnika w modelu jest usunięcie współczynnika i porównanie względnego prawdopodobieństwa z czymś w rodzaju testu chi-kwadrat:

pchisq(logLik(model1) - logLik(model2), 1)

Ponieważ sposób obliczania prawdopodobieństw między funkcjami może być nieco inny, zazwyczaj porównuję je tylko dla tej samej metody.