Analiza wariancji z powtórzeniami w R: kiedy użyć a kiedy wybrać model mieszany

Cel analizy danych z powtórzeniami: czego naprawdę szuka analityk

Osoba, która sięga po analizę wariancji z powtórzeniami w R, zwykle ma jedno z dwóch pytań: czy średnia zmienia się w czasie lub pomiędzy warunkami w tej samej grupie jednostek, oraz czy kształt tej zmiany różni się między grupami. Różnica względem klasycznej ANOVA z niezależnymi próbami polega na tym, że te same osoby, próbki lub działki doświadczalne pojawiają się w wielu wierszach zbioru danych i nie można ich traktować jako niezależnych obserwacji.

Kluczowy dylemat brzmi: czy wystarczy klasyczna ANOVA z powtórzeniami, czy należy od razu przejść do modelu mieszanego (np. w pakiecie lme4)? Odpowiedź zależy od struktury danych, liczby poziomów czasu/warunku, obecności braków danych i tego, jak poważnie naruszane są założenia sferyczności oraz równych odstępów czasowych.

Schematy danych z powtórzeniami: co faktycznie analizujemy

Jednostka losowa kontra jednostka obserwacji

W powtórzonych pomiarach często myli się jednostkę losowania z jednostką obserwacji. Jednostką losową jest typowo uczestnik badania, działka w eksperymencie polowym albo próbka laboratoryjna. Jednostką obserwacji jest natomiast pojedynczy pomiar, np. wynik w czasie T1, T2, T3. W klasycznej ANOVA z powtórzeniami każdy uczestnik pojawia się w wielu obserwacjach, ale źródłem losowości jest wciąż ta sama osoba.

Ważny wniosek: pomiary w obrębie jednego uczestnika nie są niezależne. Analiza, która zignoruje ten fakt i potraktuje wszystkie obserwacje jako niezależne w klasycznej jednoczynnikowej ANOVA, niemal zawsze przeszacuje istotność efektów (zawyża moc, zaniża błędy standardowe).

Typowe układy: czas, warunki, krzyżowanie czynników

Najczęstsze schematy danych z powtórzeniami w praktyce to:

• Pomiar w czasie: ciśnienie krwi mierzone u tych samych pacjentów przed terapią, po 4 tygodniach i po 8 tygodniach.
• Kilka warunków w tym samym dniu: czas reakcji na bodźce wizualne vs słuchowe u tych samych uczestników w jednej sesji.
• Połączenie czasu i warunku: dwa typy interwencji (A, B), pomiar w czterech punktach czasowych, ci sami uczestnicy przypisani do jednej interwencji, ale w ramach niej badani wielokrotnie.
• Struktury krzyżowe: każdy uczestnik testuje kilka bodźców, a każdy bodziec jest oceniany przez wielu uczestników – wtedy często pojawiają się dwa efekty losowe (uczestnik, bodziec).

Każdy z tych układów można formalnie opisać zarówno przy użyciu ANOVA z powtórzeniami, jak i modeli mieszanych. Różnica polega na tym, jak elastycznie można modelować wariancję pomiędzy jednostkami oraz kowariancję pomiędzy pomiarami w czasie.

Identyfikacja jednostki losowej: uczestnik, próbka, działka

Efekt losowy nie musi oznaczać wyłącznie osoby. W eksperymentach rolniczych jednostką losową może być działka, w badaniach laboratoryjnych – próbka biologiczna mierzona w kilku warunkach. Kluczowe jest pytanie: co zostało wylosowane lub przydzielone do warunków jako całość. Ta jednostka powinna mieć swój identyfikator (np. id) i być później używana jako efekt losowy w modelu mieszanym lub jako czynnik w komponencie Error() w ANOVA z powtórzeniami.

Brak świadomej identyfikacji jednostki losowej jest typową pułapką – prowadzi do błędów typu pseudoreplikacja: traktowania wielu pomiarów tej samej działki jako niezależnych replikacji. Modele mieszane zmuszają do jawnego zdefiniowania tej jednostki, ANOVA z powtórzeniami – jeśli w ogóle jest poprawnie zapisana – również, ale składnia aov() sprzyja pominięciom.

Praktyczny przykład: reakcje na bodźce w kilku sesjach

Wyobraźmy sobie badanie, w którym 30 uczestników wykonuje zadanie reakcji na bodźce w trzech sesjach (np. tydzień po tygodniu). W każdej sesji mierzony jest średni czas reakcji na bodźce wizualne i słuchowe. Dane mogą mieć strukturę:

• id – identyfikator uczestnika,
• sesja – 1, 2, 3,
• modalność – „wizualna”, „słuchowa”,
• rt – średni czas reakcji.

Pytanie badawcze: czy czas reakcji zmienia się w sesjach i czy zmiana ta jest inna dla bodźców wizualnych niż słuchowych? To klasyczny przypadek powtórzonych pomiarów – każdy uczestnik dostarcza wiele obserwacji, a różnice między uczestnikami (szybsi vs wolniejsi) powinny być potraktowane jako efekt losowy.

Klasyczna ANOVA z powtórzeniami – założenia, o których rzadko się mówi

Sferyczność, homogeniczność i normalność – praktyczne znaczenie

ANOVA z powtórzeniami opiera się na kilku założeniach, które w praktyce są często jedynie wspomniane, ale rzadko sprawdzane dogłębnie.

Sferyczność to założenie, że różnice między kolejnymi poziomami czynnika wewnątrzosobniczego mają tę samą wariancję. Mówiąc intuicyjnie: kowariancja pomiędzy pomiarami w czasie jest taka sama dla każdej pary punktów czasowych – niezależnie od odległości między nimi. W efekcie macierz kowariancji przyjmuje bardzo specyficzną, „kulistą” strukturę.

Homogeniczność wariancji odnosi się do równości wariancji pomiędzy grupami (np. między różnymi warunkami międzyosobniczymi). Przy powtórzeniach dochodzi do tego założenie o równych wariancjach reszt w obrębie każdej kombinacji czynników.

Normalność reszt oznacza, że rozkład błędów w modelu jest w przybliżeniu normalny. Przy wystarczająco dużej liczbie jednostek losowych centralne twierdzenie graniczne zwykle łagodzi to założenie, ale przy małych próbach odstępstwa mogą być istotne.

Wspólnym problemem jest to, że klasyczna ANOVA z powtórzeniami zakłada dość sztywną strukturę kowariancji, której często nie da się obronić przy danych z realnych badań (zwłaszcza, gdy punkty czasowe są nierówno rozłożone lub efekt czasu jest wyraźnie nieliniowy).

Pełne dane: brak braków, stały harmonogram

Standardowa repeated measures ANOVA w R (np. z aov() czy ezANOVA()) zakłada, że:

• każda jednostka losowa (uczestnik) ma kompletny zestaw pomiarów – brak braków danych,
• harmonogram pomiarów jest taki sam dla wszystkich – te same punkty czasu, te same warunki,
• próba jest w dużej mierze zbalansowana, bez drastycznie różnej liczby jednostek w grupach.

Gdy pojawiają się braki danych (uczestnik nie przyszedł na sesję, pomiar nie udał się technicznie), prosta ANOVA z powtórzeniami zaczyna się sypać: wiele implementacji domyślnie wyklucza całego uczestnika, jeśli brakuje choć jednego pomiaru. Z perspektywy efektywnego wykorzystania danych jest to poważna strata, zwłaszcza przy małej próbie.

Niezależność między jednostkami i „upakowana” zależność wewnątrz osoby

Klasyczna ANOVA z powtórzeniami zakłada, że obserwacje są niezależne między jednostkami (jedna osoba nie wpływa na drugą), a zależność wewnątrz jednostki jest w pełni „opakowana” w czynniku id (osoba) oraz strukturze Error(). Jednak ta zależność jest modelowana w bardzo uproszczony sposób – brak elastycznego modelowania tego, jak zmienia się wariancja i kowariancja w czasie.

To ograniczenie staje się szczególnie dotkliwe, gdy:

• odstępy czasowe między pomiarami są nierówne,
• efekt czasu jest nieliniowy (np. szybki spadek, później plateau),
• zmiany w czasie są bardzo zróżnicowane między osobami (różne „trajektorie”).

W takich scenariuszach klasyczna ANOVA z powtórzeniami zaczyna mocno rozmijać się z rzeczywistą strukturą danych. To typowy moment, w którym model mieszany przestaje być opcją „premium”, a staje się po prostu poprawnym narzędziem.

Kiedy ANOVA z powtórzeniami jest nieadekwatna

Sytuacje, gdy ANOVA z powtórzeniami zazwyczaj nie jest właściwym wyborem:

• silnie nierówne odstępy czasu (np. 1 dzień, 7 dni, 30 dni),
• dużo brakujących danych, przy czym nie są one całkowicie losowe (MNAR/MCAR wątpliwe),
• duża liczba poziomów czasu (np. 10–20 pomiarów na osobę),
• złożone układy czynników – kilka czynników wewnątrz i kilka międzyosobniczych,
• nierówne liczby obserwacji na jednostkę (np. różna liczba prób w sesji dla uczestników).

W takich przypadkach upieranie się przy ANOVA z powtórzeniami zwykle oznacza redukcję danych, stosowanie ciężkich korekt i łatanie problemów, które modele mieszane rozwiązują wprost, przez bardziej realistyczne modelowanie losowych efektów i kowariancji.

Dane szerokie vs długie i przygotowanie do analizy w R

Format wide vs long dla powtórzonych pomiarów

W wielu badaniach dane zbierane są w formacie szerokim (wide): każda kolumna to osobny pomiar w czasie lub warunek, a każdy wiersz to jedna osoba. Przykładowy układ:

• id,
• rt_T1, rt_T2, rt_T3,
• grupa.

Do większości analiz w R (szczególnie przy modelach mieszanych) znacznie wygodniejszy jest format długi (long): każdy wiersz to jeden pomiar, a informacje o czasie i warunku przechowywane są w osobnych kolumnach.

Transformacja za pomocą tidyr::pivot_longer

Przykładowa transformacja z formatu wide na long dla danych z trzema punktami czasu:

library(tidyr)
library(dplyr)

dane_long <- dane_wide |> 
  pivot_longer(
    cols = starts_with("rt_T"),
    names_to = "czas",
    values_to = "rt"
  ) |> 
  mutate(
    czas = factor(czas, levels = c("rt_T1", "rt_T2", "rt_T3"),
                  labels = c("T1", "T2", "T3"))
  )

Tak przygotowany zbiór jest od razu gotowy do użycia zarówno w aov() z komponentem Error(), jak i w lmer() z pakietu lme4. Warto na tym etapie świadomie ustalić:

• czy czas ma być traktowany jako czynnik (factor) czy zmienna ciągła,
• czy inne zmienne (np. grupa) mają być zapisane jako factor czy numeric.

Sprawdzanie poprawności struktury: liczba wierszy, poziomy ID

Po transformacji do formatu long dobrze jest sprawdzić kilka prostych rzeczy, zanim rozpocznie się anova z powtórzeniami w R lub modelowanie mieszane:

• liczbę unikalnych id:

length(unique(dane_long$id))

• liczbę wierszy – czy zgadza się z oczekiwaniem (np. liczba_id × liczba_pomiarów):

nrow(dane_long)

• rozłożenie liczby pomiarów na osobę:

dane_long |> 
  count(id) |> 
  arrange(n)

Jeśli liczba pomiarów mocno różni się między osobami, klasyczna ANOVA z powtórzeniami zaczyna być problematyczna. Dla modeli mieszanych brak balansu jest znacznie mniej kłopotliwy – choć nadal wymaga ostrożności w interpretacji.

Czas jako factor czy jako liczba

W anova z powtórzeniami czas zwykle występuje jako czynnik (np. czas = factor(T1, T2, T3)), co pozwala testować ogólny efekt czasu bez przyjmowania konkretnego kształtu trendu. W modelach mieszanych można podejść na dwa sposoby:

• traktować czas jako factor – bardziej „ANOVA-podobne” podejście,
• traktować czas jako zmienną liczbową – wygodne dla prostych trendów liniowych, z możliwością rozszerzenia o wielomiany lub splajny.

ANOVA z powtórzeniami w R: praktyczne wzorce formuł

Jednoczynnikowa ANOVA z powtórzeniami: wzorzec z Error(id/czas)

Najprostsza sytuacja: jeden czynnik wewnątrzosobniczy (np. czas), bez dodatkowych grup międzyosobniczych. Typowy kod w R:

model_aov <- aov(
  rt ~ czas + Error(id/czas),
  data = dane_long
)

summary(model_aov)

Struktura Error(id/czas) oznacza, że dla każdego uczestnika modelowany jest osobny poziom bazowy (id) oraz zagnieżdżony w nim czynnik czas. To przekłada się na klasyczną strukturę ANOVA z powtórzeniami, w której:

• id traktowane jest jako efekt losowy,
• czas jest efektem wewnątrzosobniczym testowanym z właściwym błędem.

Jeśli wszystkie założenia (sferyczność, kompletny schemat, brak braków) są rozsądnie spełnione, taki model jest szybki, zrozumiały i zaskakująco trudny do zastąpienia czymś prostszym.

Dodanie czynnika międzyosobniczego: interakcja grupa × czas

Częsty układ to jeden czynnik międzyosobniczy (np. grupa: kontrolna vs eksperymentalna) i jeden czynnik wewnątrzosobniczy (czas). Formuła ANOVA z powtórzeniami w R:

model_aov2 <- aov(
  rt ~ grupa * czas + Error(id/czas),
  data = dane_long
)

summary(model_aov2)

Tutaj grupa jest testowana na poziomie międzyosobniczym, a czas i interakcja grupa:czas – na poziomie wewnątrzosobniczym, z odpowiednimi składnikami błędu tworzonymi przez Error(). Ten zapis jest poprawny, o ile każda osoba należy tylko do jednej grupy, a grupy są względnie zbalansowane.

Popularna rada, którą często się powtarza, to „zawsze użyj aov() z Error() do powtórzonych pomiarów”. Problemy zaczynają się, gdy do gry wchodzą braki danych, nierówna liczba sesji albo bardziej złożona struktura czasu – wtedy taka formuła formalnie zadziała, ale wyniki mogą być mocno mylące.

Rozszerzony schemat: kilka czynników wewnątrz i między osobami

W badaniach poznawczych czy medycznych nie kończy się zwykle na jednym czynniku. Przykładowy układ:

• czas – czynnik wewnątrzosobniczy,
• warunek – drugi czynnik wewnątrzosobniczy (np. typ bodźca),
• grupa – czynnik międzyosobniczy.

Możliwy zapis w ANOVA z powtórzeniami:

model_aov3 <- aov(
  rt ~ grupa * czas * warunek + Error(id/(czas * warunek)),
  data = dane_long
)

summary(model_aov3)

Error(id/(czas * warunek)) definiuje strukturę zagnieżdżenia obu czynników wewnątrzosobniczych w osobie. To już punkt, w którym wiele osób gubi się w interpretacji poszczególnych tabel w summary(). Jeśli do tego dochodzą braki danych lub niestandardowy harmonogram, znacznie mniej ryzykowną drogą staje się model mieszany, mimo że formuła wydaje się na pierwszy rzut oka bardziej skomplikowana.

Przejście do modeli mieszanych: te same dane, inna filozofia

Podstawowy model mieszany: losowy punkt przecięcia dla osoby

Najprostszy odpowiednik jednoczynnikowej ANOVA z powtórzeniami w modelu mieszanym (przy czasie jako czynniku) wygląda tak:

library(lme4)

model_lmm <- lmer(
  rt ~ czas + (1 | id),
  data = dane_long
)

summary(model_lmm)

Część (1 | id) oznacza losowy intercept dla każdej osoby – różne poziomy bazowe reakcji. To odpowiednik traktowania id jako efektu losowego w ANOVA, ale już bez wymogu sferyczności i kompletności danych. Jeśli brakuje pojedynczych pomiarów, model nadal korzysta z pozostałych, zamiast wyrzucać całą osobę.

Efekt losowy nachylenia: różne trajektorie w czasie

Gdy interesuje nie tylko różny poziom bazowy, ale też różny trend w czasie pomiędzy osobami, struktura losowa powinna się rozszerzyć:

model_lmm_rs <- lmer(
  rt ~ czas + (czas | id),
  data = dane_long
)

Zapis (czas | id) oznacza, że dla każdej osoby estymowany jest:

• losowy punkt przecięcia (intercept),
• losowe nachylenie dla czasu (slope),
• oraz kowariancja między nimi.

To podejście jest w praktyce przeciwieństwem sztywnego założenia sferyczności: pozwala, by trajektorie osób różniły się kształtem i tempem, zamiast wtłaczać wszystkich w tę samą strukturę kowariancji. Popularna rada „dodaj losowe nachylenia zawsze, gdy masz czas” jest sensowna, dopóki model da się oszacować stabilnie; przy bardzo małej liczbie osób lub niewielkiej liczbie pomiarów bywa to nierealne.

Czas jako liczba: prosty trend w modelu mieszanym

Przy wielu punktach czasu przełączanie się z factor na reprezentację liczbową często daje czytelniejszy model i bardziej informatywne parametry:

dane_long <- dane_long |> 
  mutate(
    czas_num = as.numeric(czas)  # np. 1, 2, 3, ...
  )

model_lmm_lin <- lmer(
  rt ~ czas_num + (czas_num | id),
  data = dane_long
)

Nachylenie przy czas_num to wtedy średnia zmiana reakcji przy przesunięciu o jednostkę czasu. Przy równych odstępach taka parametryzacja jest przejrzysta. Gdy odstępy są nierówne (np. dni = 0, 3, 10, 30), sensowniejsze jest wpisanie czas_num jako rzeczywistego czasu (np. liczba dni od pierwszego pomiaru), zamiast używania sztucznych indeksów 1, 2, 3, 4.

Łączenie czynników wewnątrz i między osobami w lmer()

Odpowiednik schematu z grupa × czas w modelu mieszanym:

model_lmm2 <- lmer(
  rt ~ grupa * czas + (czas | id),
  data = dane_long
)

Tutaj:

• grupa – efekt międzyosobniczy,
• czas – efekt wewnątrzosobniczy,
• grupa:czas – różnice w zmianie w czasie między grupami,
• (czas | id) – indywidualne trajektorie, na które nakładają się efekty grupowe.

To jeden z momentów, w którym klasyczna ANOVA z powtórzeniami traci przewagę: gdy pojawiają się braki danych, nierówne odstępy i zróżnicowane trajektorie, taki model mieszany lepiej oddaje strukturę, przy mniejszej liczbie „awaryjnych” korekt.

Jak wybierać strukturę losową: minimalizm kontra „maksymalny” model

Koncepcja „maksymalnego” modelu losowego i jej ograniczenia

Popularna w literaturze psycholingwistycznej rada brzmi: „używaj maksymalnego możliwego modelu losowego, który jest zgodny z planem badania”. W praktyce przekłada się to na dodanie losowych nachyleń dla wszystkich istotnych czynników wewnątrzosobniczych:

model_max <- lmer(
  rt ~ grupa * czas * warunek +
    (czas * warunek | id),
  data = dane_long
)

Taka struktura bywa bardzo atrakcyjna teoretycznie, ale ma dwie praktyczne pułapki:

• konwergencja modelu – przy małej liczbie osób lub krótkich seriach pomiarów algorytm estymacji może się nie dogadać,
• przeparametryzowanie – macierz kowariancji efektów losowych jest trudna do oszacowania, co kończy się osobliwymi rozwiązaniami (wariancje ~ 0, kowariancje graniczne).

Jeśli model sygnalizuje problemy z konwergencją, a wyniki dla losowych efektów wyglądają nienaturalnie, dokładanie kolejnych nachyleń tylko po to, by „było maksymalnie”, mija się z celem.

Praktyczny kompromis: od prostego do bardziej złożonego

Rozsądniejsza strategia w wielu projektach to budowanie struktury losowej stopniowo. Przykładowy schemat:

1. Zacząć od modelu z samym losowym interceptem: (1 | id).
2. Dodać najważniejszy czynnik wewnątrzosobniczy jako losowe nachylenie: (czas | id).
3. Jeśli są kolejne czynniki w obrębie osoby (np. warunek), rozważyć (czas + warunek | id), ale tylko przy wystarczającej liczbie obserwacji.

Na każdym etapie warto spojrzeć nie tylko na p-wartości, ale też na oszacowane wariancje efektów losowych. Jeżeli wariancja losowego nachylenia jest bardzo bliska zeru, model sygnalizuje, że ten poziom złożoności nie wnosi wiele ponad prostszy wariant.

Gdy danych jest mało: kiedy minimalny model jest uczciwszy

W małych próbach (np. kilkanaście osób, 3–4 punkty w czasie) forsowanie skomplikowanej struktury losowej bywa intelektualnie kuszące, ale statystycznie wątpliwe. Minimalny, stabilny model z (1 | id) i prostą strukturą stałych efektów może być wtedy bardziej uczciwy niż wymuszony „maksymalny” model, który:

• nie konwerguje,
• ma osobliwe rozwiązania,
• trudno go zinterpretować i zreplikować.

Zamiast ślepo trzymać się zasady „maksymalnej” struktury losowej, lepiej dopasować poziom złożoności do liczby jednostek losowych, liczby powtórzeń i jakości danych.

Modele mieszane a brakujące dane i nierówne odstępy czasu

Mechanizm braków: kiedy lmer() faktycznie „ratuje” dane

Jedną z głównych zalet modeli mieszanych jest łatwiejsze obchodzenie się z brakującymi danymi na poziomie obserwacji. Przykład z praktyki: badanie z czterema sesjami, część uczestników wypadła po trzeciej. W ANOVA z powtórzeniami zwykle trzeba usunąć tych uczestników całkowicie lub przeprowadzać imputacje. W modelu mieszanym:

• osoba z 3/4 sesji nadal wnosi informację o poziomie i trendzie,
• przy założeniu braków MAR (brak zależy od obserwowanych, ale nie nieobserwowanych zmiennych) estymatory pozostają nieobciążone.

To nie jest magiczne „uzupełnianie” braków, tylko wykorzystanie pełnej informacji, którą da się oszacować na podstawie struktury modelu.

Nierówne odstępy czasu: przewaga modelowania czasu liczbowego

Przy schematach typu: pomiar 1 dzień po interwencji, potem 7, potem 30 dni, traktowanie czasu jako czystego factor w ANOVA ignoruje informację o rzeczywistych odstępach. Model mieszany, w którym czas zapisany jest np. jako liczba dni od pierwszego pomiaru, pozwala:

• szacować tempo zmiany per jednostkę czasu,
• dopuszczać różne trajektorie między osobami (losowe nachylenia),
• korzystać z nierównomiernej siatki pomiarów.

dane_long <- dane_long |> 
  mutate(
    dni = as.numeric(dni_od_startu)  # np. 1, 7, 30
  )

model_lmm_time <- lmer(
  rt ~ dni + (dni | id),
  data = dane_long
)

Klasyczna ANOVA z powtórzeniami nie ma tu dobrego wyjścia: wymagałaby albo zignorowania informacji o rzeczywistych odstępach, albo sztucznego przerabiania osi czasu (co i tak nie rozwiązuje problemu struktury kowariancji).

Nieregularna liczba pomiarów na osobę

W eksperymentach behawioralnych część uczestników kończy wszystkie próby, inni rezygnują w połowie, a jeszcze inni mają przerwane sesje. W ANOVA z powtórzeniami taka nieregularność jest kłopotliwa, bo:

• wielu uczestników trzeba wykluczyć,
• schemat analizy przestaje być „czysty” i porównywalny między grupami.

Model mieszany z (1 | id) lub (czas | id) radzi sobie z tym znacznie lepiej, o ile:

• mechanizm braków nie jest skrajnie zależny od nieobserwowanych wyników (np. osoby z bardzo złymi wynikami systematycznie odpadają),
• liczba dostępnych pomiarów na osobę jest wystarczająca do oszacowania trajektorii.

Porównywanie ANOVA z powtórzeniami i modeli mieszanych na tych samych danych

Ten sam schemat, dwa narzędzia: jak wygląda praktyczne porównanie

Przy prostym schemacie: kilka pomiarów w czasie, dwa poziomy czynnika międzyosobniczego i żadnych braków, klasyczna ANOVA i model mieszany dadzą zwykle podobne wnioski co do efektów głównych i interakcji. Różnice ujawniają się, kiedy:

• naruszone są założenia ANOVA (np. sferyczność),
• wprowadzamy kowariaty,
• chcemy modelować różnice indywidualne w trajektoriach.

Prosty test porównawczy można wykonać na tych samych danych, licząc ANOVA z powtórzeniami (np. przez aov_car() z pakietu afex) i model mieszany (lmer()), a następnie sprawdzając zgodność kierunku i rzędu wielkości efektów.

library(afex)
library(lme4)

# ANOVA z powtórzeniami
aov_res <- aov_car(
  rt ~ grupa * czas + Error(id/czas),
  data = dane_long,
  factorize = FALSE
)

# Model mieszany
lmm_res <- lmer(
  rt ~ grupa * czas + (czas | id),
  data = dane_long
)

Jeśli efekty są stabilne, główne wnioski (np. „grupy różnią się trendem w czasie”) powinny się zgadzać. Gdy tak nie jest, pierwszym podejrzanym jest zwykle struktura kowariancji i korekty w ANOVA, a nie sam fakt użycia lmer().

Kiedy ANOVA i lmer() dają różne odpowiedzi

Typowa sytuacja rozjazdu pojawia się, gdy ANOVA używa korekty Greenhouse–Geissera (silne naruszenie sferyczności), a model mieszany ma bogatą strukturę losową (np. losowe nachylenia). Objawy:

• p-wartości dla efektu czas i grupa:czas są w ANOVA znacznie „łagodniejsze” lub odwrotnie – dużo ostrzejsze,
• średnie z ANOVA sugerują „gładką” trajektorię, podczas gdy predykcje z lmer() pokazują dużą heterogeniczność indywidualną.

Jeżeli różnice biorą się głównie z korekt stopni swobody w ANOVA, a model mieszany wydaje się dobrze dopasowany (brak ostrzeżeń, sensowne wariancje losowe), wnioskowanie z modelu mieszanego bywa bardziej wiarygodne, bo nie opiera się na uproszczeniu struktury kowariancji.

Redukcja modelu mieszanego do „odpowiednika ANOVA”

Często pada rada: „użyj lmer() zamiast ANOVA, ale ustaw bardzo prostą strukturę losową, żeby wyniki były bardziej zbliżone do ANOVA”. Ten zabieg ma sens tylko częściowo. Przykład „minimalnego” modelu odpowiadającego ANOVA z jednym czynnikiem wewnątrz:

model_min <- lmer(
  rt ~ grupa * czas + (1 | id),
  data = dane_long
)

Takie podejście działa jako przybliżenie, ale:

• ignoruje różnice indywidualne w trendzie (brak (czas | id)),
• wymusza jednorodną strukturę kowariancji wewnątrz osoby (różne od ANOVA z sferycznością),
• może zaniżać wariancję błędu dla efektu czas, jeśli rzeczywiste trajektorie mocno się różnią.

Jeżeli plan badania i liczba obserwacji pozwalają, bardziej uczciwe jest dołożenie losowego nachylenia dla czasu i zaakceptowanie faktu, że wyniki nie będą idealnym „lustrem” klasycznej ANOVA.

Jak testować efekty w modelach mieszanych: F, t, LRT i inne wynalazki

Popularne rekomendacje i gdzie się rozjeżdżają

Wokół testowania w modelach mieszanych narosło sporo „prostych” rad: używać LRT (likelihood ratio test), patrzeć na t > |2|, stosować korekty Kenwarda–Rogera’a albo Satterthwaite’a. Każda z tych metod bywa przydatna, ale żadna nie jest uniwersalnym złotym standardem.

LRT: testowanie przez porównanie z modelem zredukowanym

Porównywanie modeli przez LRT ma sens, jeśli efekty testowane wchodzą hierarchicznie w strukturę modelu. Przykład: test interakcji grupa:czas:

model_full <- lmer(rt ~ grupa * czas + (czas | id), data = dane_long)
model_red  <- lmer(rt ~ grupa + czas + (czas | id), data = dane_long)

anova(model_red, model_full)

Ta strategia jest spójna przy dużych próbach, ale przy małej liczbie jednostek losowych p-wartości z LRT są często zbyt „agresywne” (zaniżone), bo rozkład graniczny χ² jest tylko przybliżeniem. W eksperymentach z kilkunastoma osobami takie testy łatwo przeceniają istotność interakcji.

Korekty Satterthwaite’a i Kenwarda–Rogera’a

Podejście implementowane m.in. w pakiecie lmerTest (Satterthwaite) czy pbkrtest (Kenward–Roger) polega na przybliżaniu rozkładu statystyk t i F z „efektywnymi” stopniami swobody. Schemat użycia:

library(lmerTest)

model_lmm <- lmer(
  rt ~ grupa * czas + (czas | id),
  data = dane_long
)

summary(model_lmm)  # zawiera p-wartości z przybliżeniem Satterthwaite’a

Zalety:

• większa zgodność z ANOVA w prostych schematach,
• testy oparte na standardowych rozkładach (t, F),
• praktycznie użyteczne przy średnich licznościach osób.

Słaby punkt: przy bardzo skomplikowanej strukturze losowej i niewielu jednostkach losowych oszacowanie efektywnych stopni swobody staje się niepewne, a wyniki potrafią być wrażliwe na drobne zmiany modelu.

Reguła „|t| > 2”: kiedy ma sens, a kiedy jest myląca

Często słychać poradę: „przy dużej liczbie stopni swobody test t można przybliżyć N(0,1), więc |t| > 2 ~ p < .05”. To przybliżenie bywa sensowne przy naprawdę dużej liczbie obserwacji i jednostek losowych, ale w typowym badaniu psychologicznym (np. 30–60 osób) jest to raczej szybka heurystyka niż solidne wnioskowanie.

Zdarza się też odwrotna skrajność: całkowite odrzucenie tej reguły. Rozsądne użycie: traktować |t| > 2 jako wskaźnik efektów wartych dalszego badania (szczególnie przy eksploracji), a nie ostateczne kryterium „istotności”.

Dobór struktury stałych efektów: kiedy komplikować, a kiedy upraszczać

Interakcje wyższego rzędu: czy zawsze „trzeba” je dodać

Kanon mówi: „jeśli masz trzy czynniki, zawsze uwzględnij pełną interakcję trójczynnikową, bo inaczej model jest błędny”. To podejście ma sens przy dużych próbach i precyzyjnej hipotezie o takiej interakcji. W wielu projektach eksperymentalnych jest jednak odwrotnie: interakcja trójczynnikowa jest mgliście zdefiniowana, a próba za mała, by ją wiarygodnie oszacować.

Model mieszany z pełną interakcją:

model_full <- lmer(
  rt ~ grupa * czas * warunek + (czas * warunek | id),
  data = dane_long
)

może szybko stać się nieidentyfikowalny, a nawet jeśli się zbiega, parametry interakcji wyższego rzędu mają ogromne błędy standardowe. W takiej sytuacji sensowniejsze jest sformułowanie konkretnych kontrastów (np. różnica między grupami w zmianie między t1 i t2 w określonym warunku), zamiast mechanicznego testowania wszystkich interakcji.

Kontrasty a „gołe” poziomy czynników

Zamiast interpretować główne efekty w obecności interakcji na poziomie kodowania dummy, wygodniej bywa użyć kontrastów odzwierciedlających konkretne pytania. Przykład: interesuje nas zmiana od t1 do t3 w każdej grupie osobno.

library(emmeans)

emm <- emmeans(model_lmm2, ~ grupa | czas)

contrast(emm, "revpairwise")  # porównania między grupami w każdym czasie

# lub konkretne kontrasty czasu wewnątrz grupy
emm_time <- emmeans(model_lmm2, ~ czas | grupa)
contrast(emm_time, "consec")  # kolejne różnice w czasie w każdej grupie

Podobną logikę da się zastosować bez ANOVA: nie trzeba budować coraz bardziej skomplikowanych stałych efektów, jeśli sensowne pytania można ująć w postaci kontrastów na już dopasowanym modelu.

Kowariaty: kiedy dodanie zmiennej dodatkowej ma sens

Częsta rada: „dodaj IQ, wiek lub zmienną wyjściową jako kowariatę, żeby wyczyścić wariancję”. W modelach mieszanych, szczególnie z czasem jako czynnikiem, sprawa jest subtelniejsza. Przykładowo:

• dodanie wyniku wyjściowego jako kowariaty przy jednoczesnym modelowaniu interceptu może prowadzić do kolinearności,
• kowariata może mieć różny efekt w czasie (np. wiek wpływa mocniej na późniejsze pomiary), co sugeruje interakcję kowariata × czas.

model_cov <- lmer(
  rt ~ grupa * czas + wiek + (czas | id),
  data = dane_long
)

Jeśli wiek jest istotnie skorelowany z czasem lub warunkiem (np. starsi częściej trafiają do jednej grupy), prosty dodatek jako kowariaty nie usuwa problemu planu quasi-eksperymentalnego. W takich przypadkach lepiej od początku traktować model bardziej opisowo niż przyczynowo, a efekty „oczyszczone” interpretować ostrożnie.

Implementacja w R: od prostych modeli ANOVA do elastycznych modeli mieszanych

afex jako pomost między światem ANOVA i lmer()

Pakiet afex bywa praktycznym kompromisem: pozwala zadeklarować projekt w stylu ANOVA, a w tle dopasować model mieszany i wyciągać estymaty metodami lme4. Przykład:

library(afex)

aov_mixed_res <- aov_car(
  rt ~ grupa * czas + Error(id/czas),
  data = dane_long,
  factorize = TRUE,
  include_aov = TRUE  # umożliwia porównanie z klasycznym aov
)

# ten sam projekt jako model mieszany z afex
mixed_res <- mixed(
  rt ~ grupa * czas,
  data = dane_long,
  id = "id",
  within = "czas",
  between = "grupa",
  method = "KR"  # Kenward-Roger
)

Takie podejście jest wygodne przy przechodzeniu z myślenia „ANOVA-owego” do mieszanego: zachowujemy znany zapis efektów stałych, a jednocześnie otrzymujemy poprawniejszą obsługę kowariancji i braków danych.

lme4 dla osób, które chcą pełnej kontroli

Gdy znajomość modeli mieszanych rośnie, naturalne jest przejście na bezpośrednie użycie lmer(). Pozwala to:

• definiować dowolne struktury losowe (w tym zagnieżdżenia, skrzyżowania),
• modelować nieliniowe funkcje czasu (np. polinom, splajny),
• łączyć kilka poziomów zagnieżdżenia (np. pomiary w zadaniach w osobach).

# zagnieżdżenie prób w zadaniu i zadania w osobie
model_nested <- lmer(
  rt ~ grupa * czas +
    (czas | id) +
    (1 | id:zadanie),
  data = dane_long
)

Struktura (1 | id:zadanie) oddaje dodatkową wariancję między zadaniami w obrębie osoby, którą ANOVA z powtórzeniami musiałaby traktować jako część błędu resztowego.

nlme i jawne struktury korelacji

Tam, gdzie szczególnie zależy na modelowaniu konkretnej struktury korelacji (np. autoregresja AR(1) w długich szeregach czasowych), pakiet nlme oferuje jawne zadawanie korelacji między obserwacjami:

library(nlme)

model_ar1 <- lme(
  rt ~ dni + grupa,
  random = ~ 1 | id,
  correlation = corAR1(form = ~ dni | id),
  data = dane_long
)

To podejście jest bliższe „modelom powtarzanych pomiarów” znanym z pakietów typu SPSS czy SAS. Gdy główną troską jest specyficzna forma korelacji, a nie bogata struktura losowych nachyleń, taki model bywa bardziej adekwatny niż klasyczny lmer() z domyślną, nieustrukturyzowaną macierzą kowariancji efektów losowych.

Typowe pułapki interpretacyjne przy przejściu z ANOVA na modele mieszane

Interpretacja efektów głównych w obecności losowych nachyleń

W ANOVA wiele osób jest przyzwyczajonych, że „główny efekt czasu” to różnica średnich między poziomami, uśredniona po grupach. W modelu mieszanym z (czas | id) ten sam współczynnik opisuje średni trend populacyjny, ale „na tle” losowych trajektorii jednostek.