Wiele porównań w modelu efektów mieszanych

31

Próbuję analizować niektóre dane przy użyciu modelu efektu mieszanego. Zebrane przeze mnie dane przedstawiają masę niektórych młodych zwierząt o różnym genotypie w czasie.

Korzystam z zaproponowanego tutaj podejścia: https://gribblelab.wordpress.com/2009/03/09/repeated-measures-anova-using-r/

W szczególności używam rozwiązania nr 2

Więc mam coś takiego

require(nlme)
model <- lme(weight ~ time * Genotype, random = ~1|Animal/time, 
         data=weights)    
av <- anova(model)

Teraz chciałbym mieć kilka porównań. Za pomocą multcompmogę zrobić:

require(multcomp)
comp.geno <- glht(model, linfct=mcp(Genotype="Tukey"))
print(summary(comp.geno))

I oczywiście mogłem zrobić to samo z czasem.

Mam dwa pytania:

  1. Jak mogę mcpzobaczyć interakcję między czasem a genotypem?
  2. Po uruchomieniu glhtpojawia się następujące ostrzeżenie:

    covariate interactions found -- default contrast might be inappropriate

    Co to znaczy? Czy mogę to bezpiecznie zignorować? Lub co powinienem zrobić, aby tego uniknąć?

EDYCJA: Znalazłem ten plik PDF, który mówi:

Ponieważ w tym przypadku niemożliwe jest automatyczne określenie parametrów będących przedmiotem zainteresowania, mcp () w multcomp domyślnie generuje porównania tylko dla głównych efektów, ignorując zmienne towarzyszące i interakcje . Od wersji 1.1-2 można określić średnią dla warunków interakcji i zmiennych towarzyszących, używając odpowiednio argumentów: interakcji_średnia = PRAWDA i zmienna_obiektywna = PRAWDA, natomiast wersje starsze niż 1.0-0 są automatycznie uśredniane względem warunków interakcji. Sugerujemy jednak użytkownikom, aby zapisali ręcznie zestaw kontrastów, których chcą.Należy to zrobić, gdy istnieją wątpliwości co do tego, co mierzą domyślne kontrasty, co zwykle dzieje się w modelach z warunkami interakcji wyższego rzędu. Odnosimy się do Hsu (1996), rozdział ~ 7 i Searle (1971), rozdział ~ 7.3, w celu dalszych dyskusji i przykładów na ten temat.

Nie mam dostępu do tych książek, ale może ktoś tu ma?

Nico
źródło
Spójrz na InvivoStat ( invivostat.co.uk ), powinien wykonać typ analizy, której szukasz

Odpowiedzi:

29

Jeśli timei Genotypesą predyktorami jakościowymi, i wydają się zainteresowane porównywaniem wszystkich par czas / genotyp, możesz po prostu utworzyć jedną zmienną interakcji i użyć na niej kontrastów Tukeya:

weights$TimeGeno <- interaction(weigths$Time, weights$Geno)
model <- lme(weight ~ TimeGeno, random = ~1|Animal/time, data=weights) 
comp.timegeno <- glht(model, linfct=mcp(TimeGeno="Tukey")) 

Jeśli interesują Cię inne kontrasty, możesz skorzystać z faktu, że linfctargument może przyjąć macierz współczynników dla kontrastów - w ten sposób możesz ustawić dokładnie takie porównania, jakie chcesz.

EDYTOWAĆ

W uwagach pojawia się obawa, że ​​model wyposażony w TimeGenopredyktor różni się od modelu oryginalnego wyposażonego w Time * Genotypepredyktor. Tak nie jest , modele są równoważne. Jedyną różnicą jest parametryzacja stałych efektów, która jest skonfigurowana w celu ułatwienia korzystania z glhtfunkcji.

Użyłem jednego z wbudowanych zestawów danych (ma dietę zamiast genotypu), aby wykazać, że oba podejścia mają takie samo prawdopodobieństwo, przewidywane wartości itp .:

> # extract a subset of a built-in dataset for the example
> data(BodyWeight)
> ex <- as.data.frame(subset(BodyWeight, Time %in% c(1, 22, 44)))
> ex$Time <- factor(ex$Time)
> 
> #create interaction variable
> ex$TimeDiet <- interaction(ex$Time, ex$Diet)
    > 
    > model1 <- lme(weight ~ Time * Diet, random = ~1|Rat/Time,  data=ex)    
    > model2 <- lme(weight ~ TimeDiet, random = ~1|Rat/Time, data=ex)    
    > 
    > # the degrees of freedom, AIC, BIC, log-likelihood are all the same 
    > anova(model1, model2)
           Model df      AIC      BIC    logLik
    model1     1 12 367.4266 387.3893 -171.7133
    model2     2 12 367.4266 387.3893 -171.7133
    Warning message:
    In anova.lme(model1, model2) :
      fitted objects with different fixed effects. REML comparisons are not meaningful.
    > 
    > # the second model collapses the main and interaction effects of the first model
    > anova(model1)
                numDF denDF   F-value p-value
    (Intercept)     1    26 1719.5059  <.0001
    Time            2    26   28.9986  <.0001
    Diet            2    13   85.3659  <.0001
    Time:Diet       4    26    1.7610  0.1671
    > anova(model2)
                numDF denDF   F-value p-value
    (Intercept)     1    24 1719.5059  <.0001
    TimeDiet        8    24   29.4716  <.0001
    > 
    > # they give the same predicted values
    > newdata <- expand.grid(Time=levels(ex$Time), Diet=levels(ex$Diet))
    > newdata$TimeDiet <- interaction(newdata$Time, newdata$Diet)
> newdata$pred1 <- predict(model1, newdata=newdata, level=0)
    > newdata$pred2 <- predict(model2, newdata=newdata, level=0)
> newdata
  Time Diet TimeDiet   pred1   pred2
1    1    1      1.1 250.625 250.625
2   22    1     22.1 261.875 261.875
3   44    1     44.1 267.250 267.250
4    1    2      1.2 453.750 453.750
5   22    2     22.2 475.000 475.000
6   44    2     44.2 488.750 488.750
7    1    3      1.3 508.750 508.750
8   22    3     22.3 518.250 518.250
9   44    3     44.3 530.000 530.000

Jedyną różnicą jest to, że hipotezy są łatwe do przetestowania. Na przykład w pierwszym modelu łatwo jest sprawdzić, czy dwa predyktory oddziałują na siebie, w drugim modelu nie ma na to wyraźnego testu. Z drugiej strony, łączny efekt dwóch predyktorów jest łatwy do przetestowania w drugim modelu, ale nie w pierwszym. Pozostałe hipotezy są testowalne, ich konfiguracja jest po prostu więcej pracy.

Aniko
źródło
3
glhtużywa stopni swobody podanych w modelu Lme. Nie jestem pewien, czy te stopnie swobody są odpowiednie ...?
Stéphane Laurent,
2
Jestem również ciekawy, jak najlepiej to zrobić. Podejście to daje jednak efekty z innego modelu - takiego, który zasadniczo testuje tylko interakcję. Drugi model w ogóle nie obejmuje dwóch potencjalnych głównych efektów. Nie wydaje się to odpowiednią metodą sprawdzania efektów w pierwszym modelu.
Marcus Morrisey
@Aniko, myślałem o połączeniu 2 zmiennych kategorialnych w jedną, tak jak właśnie to zrobiłeś, ale zawahałem się, ponieważ tylko jedna ze zmiennych jest w obrębie podmiotu, a druga jest pomiędzy. Czy możesz potwierdzić, że to nie pasuje? Zauważyłem, że w przykładzie, Animal/timektóry trzymasz, który teraz nie jest jednym z czynników. Czy to naprawdę understand?
toto_tico
@toto_tico, edytowałem odpowiedź, aby pokazać, że drugi model jest równoważny pierwszemu.
Aniko
3
@toto_tico, dałem ci powtarzalny przykład. Dlaczego nie spróbujesz all.equal(resid(model1), resid(model2))przekonać się, że są takie same, zanim zgadniesz inaczej? Jedynie parametryzacja efektów stałych jest inna. TimeDietnie jest czystym pojęciem interakcji i nie jest równoważne Time:Diet, ale raczej Time + Diet + Time:Diet.
Aniko