Wie man Varianzkomponenten mit lmer für Modelle mit zufälligen Effekten schätzt und mit lme-Ergebnissen vergleicht

14

Ich habe ein Experiment durchgeführt, bei dem ich verschiedene Familien aus zwei verschiedenen Bevölkerungsgruppen großgezogen habe. Jede Familie erhielt eine von zwei Behandlungen. Nach dem Experiment habe ich mehrere Merkmale an jedem Individuum gemessen. Um die Wirkung einer Behandlung oder einer Quelle sowie deren Wechselwirkung zu testen, verwendete ich ein lineares Mischeffektmodell mit der Familie als Zufallsfaktor, d. H

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

Soweit so gut. Jetzt muss ich die relativen Varianzkomponenten berechnen, dh den Prozentsatz der Variation, der durch die Behandlung oder die Quelle sowie die Wechselwirkung erklärt wird.

Ohne einen zufälligen Effekt könnte ich leicht die Quadratsummen (SS) verwenden, um die Varianz zu berechnen, die durch jeden Faktor erklärt wird. Da es für ein gemischtes Modell (mit ML-Schätzung) keine SS gibt, dachte ich, ich könnte Treatment und Source auch als Zufallseffekte verwenden, um die Varianz zu schätzen, d. H

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

In einigen Fällen konvergiert lme jedoch nicht. Daher habe ich lmer aus dem Paket lme4 verwendet:

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

Wo ich die Varianzen aus dem Modell mit der Zusammenfassungsfunktion extrahiere:

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

Ich erhalte die gleichen Werte wie bei der VarCorr-Funktion. Ich benutze dann diese Werte, um den tatsächlichen Prozentsatz der Variation zu berechnen, wobei die Summe als Gesamtvariation genommen wird.

Ich habe Probleme mit der Interpretation der Ergebnisse aus dem ersten Modell (mit Behandlung und Quelle als festen Effekten) und dem Zufallsmodell zur Schätzung der Varianzkomponenten (mit Behandlung und Quelle als zufälligen Effekten). Ich finde in den meisten Fällen, dass der Prozentsatz der durch jeden Faktor erklärten Varianz nicht der Signifikanz des festen Effekts entspricht.

Zum Beispiel für das Merkmal HD deutet der erste Abschnitt auf eine Tendenz zur Wechselwirkung sowie auf eine Bedeutung für die Behandlung hin. Ich stelle fest, dass die Behandlung nach einem Rückwärtsverfahren eine beinahe signifikante Tendenz aufweist. Bei der Schätzung der Varianzkomponenten stelle ich jedoch fest, dass die Quelle die höchste Varianz aufweist und bis zu 26,7% der Gesamtvarianz ausmacht.

Die lme:

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

Und der 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

Daher ist meine Frage, ist es richtig, was ich tue? Oder sollte ich eine andere Methode verwenden, um den Betrag der Varianz zu schätzen, der durch jeden Faktor erklärt wird (dh Behandlung, Quelle und deren Wechselwirkung). Wären die Effektgrößen beispielsweise ein geeigneterer Weg?

KL_STKBK
quelle
Der Behandlungsfaktor hat 40x so viele Freiheitsgrade wie der Quellfaktor (Pseudoreplikation?). Dies verringert zweifellos den P-Wert der Behandlung.

Antworten:

1

Eine gebräuchliche Methode, um den relativen Beitrag jedes Faktors zu einem Modell zu bestimmen, besteht darin, den Faktor zu entfernen und die relative Wahrscheinlichkeit mit einem Chi-Quadrat-Test zu vergleichen:

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

Da die Art und Weise, wie Wahrscheinlichkeiten zwischen Funktionen berechnet werden, leicht unterschiedlich sein kann, vergleiche ich sie normalerweise nur zwischen derselben Methode.

Bill Denney
quelle
1
sollte es nicht sein 1-pchisq(logLik(model1) - logLik(model2), 1)?
user81411