Wiele porównań w modelu efektów mieszanych

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?

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.