Pokračujeme v našem malém úvodním seriálu o kognitivní psychologii ilustrovaném na příkladu klasického Posnerova výzkumu skryté zrakové pozornosti. Již jsme se dozvěděli něco o psychologickém konceptu pozornosti, jak se dá takový kognitivistický konstrukt objektivně zkoumat, jak se přistupuje k designu konkrétního experimentu a sběru dat, a v tomto příspěvku si získaná data základně zanalyzujeme.

Možná že někomu bude tento příspěvek připadat suchý, ale takto skutečně vypadá podstatná část vědecké práce (a není nic krásnějšího, než si moci hrát s výsledky své práce – a nalézt tam vážně to, co jsme nalézt chtěli!:-)). Snad to bude pro někoho prakticky užitečné; kéž bych tohle vše například znal v době, kdy jsem dělal na své diplomové práci! Kroky zpracování dat budu vysvětlovat slovně i ilustrovat pomocí jazyka R, což je jeden z nejpoužívanějších statistických nástrojů vůbec, a navíc, je zadarmo (doporučuji k práci grafické prostředí R Studio). Pokud chcete pracovat s daty, naučit se R se rozhodně vyplatí. Excel ani SPSS vám takovou jednoduchost a sílu nenabídnou.

Začínáme od surových dat

Formát, v jakém se nám data z experimentu uložila, závisí na použitém experimentálním programu. Většinou ale mají formu dlouhé tabulky (ve formátu txt, csv, xls..), kde každý řádek je jeden trial a jednotlivé sloupce obsahují nějaké informace o tomto trialu. Např. jeho číslo, datum a čas, číslo probanda, reakční čas, správnost odpovědi, všechny experimentální faktory/podmínky (zda zrovna šlo o validní trial, jaké bylo SOA..) a případně množství dalších informací.

Dva sloupce obsahují závislé proměnné (DV, dependent variables): RT (reakční čas, reaction time) a ACC (správnost odpovědi, accuracy). To jsou hodnoty (measures), které nás zajímají; hlavní výsledky experimentu. Další dva sloupce obsahují nezávislé proměnné (IV, independent variables): Validity (zda nápověď ukazovala směrem k cílovému stimulu) a SOA (časová prodleva mezi nápovědí a cílovým stimulem). To jsou faktory, které jsme v experimentu manipulovali. Zajímá nás, jak jsou DV ovlivněny jednotlivými IV. Jak to zjistíme? Především, nejde jenom teď rozdělit trialy do skupin dle IV a spočítat třeba průměrný reakční čas v každé z nich a ty pak porovnat. Je potřeba tohle rozdělení dle IV a spočítání průměrů (či mediánů apod.) nejprve udělat pro každého participanta. Tyto individuální průměry pak užijeme ke statistickému porovnání, jak si ukážeme dále.

Nejprve se můžeme podívat na DV správnost odpovědi. Pro každého probanda, každou IV (Validity a SOA) si spočítáme počet jedniček (správná odpověď) a nul (nesprávná odpověď) a z těch pak spočítáme procento úspěšnosti. Anebo můžeme využít toho, že máme tyto hodnoty šikovně zaznamenané právě jako nuly a jedničky a spočítáme jednoduše průměrnou hodnotu ACC (pro každou kombinaci Subject x Validity x SOA) a rovnou získáme procento úspěšnosti, např. 0.87 pro 87%! V eRku to spočítáme velmi snadno:

AccuracyRates <- aggregate(ACC ~ Subject + Validity + SOA, data = TrialData, FUN = mean)

Vymyšlený výsledek v proměnné AccuracyRates:

Data trochu vyčistíme..

Závislou proměnnou Accuracy se ale nebudu nyní více zabývat a přeskočíme k počítání obvykle zajímavějších reakčních časů. Zajímají nás pravděpodobně reakční časy jen v trialech, kde byla podána správná odpověď, takže vyloučíme řádky, kde jsou nuly v ACC.

TrialDataCorrect <- subset(TrialData, ACC!=0)

Někdy můžeme vyloučit i některé jiné trialy, pokud si pro to stanovíme kritéria před provedením analýzy (zde závisí na poctivosti výzkumníka, aby si pak nevymýšlel kritéria post hoc, aby mu to hezky vyšlo, a neskončil pak v některém z kruhů vědeckého pekla). Např. někdy se vylučují reakční časy, které se příliš odlišují od “běžných” reakčních časů – jsou nereálně rychlé (třeba 120 ms, proband asi začal mačkat tlačítko ještě před zobrazením cíle a zrovna se náhodou strefil) či podezřele pomalé (např. 1800 ms, proband asi zaspal) – a které nazýváme outliers. Není na škodu podívat se na histogram či boxplot reakčních časů, v těchto grafech lze zjevné outliery spatřit. Jejich vyloučení je vhodné dělat zvlášť pro každou kombinací Subjekt x jednotlivé IV (protože pochopitelně každý člověk má sám o sobě jinak rychlé reakce a v každé experimentální podmínce je ten samý člověk také bude mít odlišné). Kritériem vyloučení může být např. vzdálenost více než 3 směrodatných odchylek od průměru. Metod je vícero a někdy se i zastává názor, že je lepší se tomuto odstraňování vyhnout, protože to zvyšuje výzkumníkovo stupně volnosti v provádění analýzy (a použít např. medián místo průměru). Většinou se ale nějaká metoda používá. Data, která budu dále analyzovat, skutečně obsahovala výrazné outliery (reakční čas i několik sekund), níže je pročistím jen na základě hrubých limitů (po prohlédnutí histogramů, viz dále).

 

Je možné vyloučit nějen trialy, ale rovnou vyloučit některé participanty. Např. Pokud si a priori (!) stanovíme kritérium, že pokud má člověk moc velkou chybovost (např. více než 30 %), nebo nějak poznáme, že úlohu nepochopil, data od něj asi nebudou moc kvalitní a bude lepší obětovat všechna data od tohoto účastníka.

Spočítáme hlavní výsledky…

Následně můžeme spočítat průměrné reakční časy ve všech kombinacích IV pro každého probanda.

RTData.Long <- aggregate(RT ~ Subject + Validity + SOA, data = TrialDataCorrect, FUN = mean)

Vysledek bude vypadat podobně jako výše u AccuracyRates. Tomuto získanému formátu dat se říká long format. Ve sloupcích máme jednotlivé IV a v samostatném sloupci DV (reakční čas), na řádcích pak všechny existující kombinace IV (a subjektů). Dlouhý formát je vhodný k provádění mnoha druhů analýz a obecně k další práci s daty.

…na základě reálných dat

Obrázek s takovými daty neukážu, protože od teď končíme s vymyšlenými daty a přecházíme k reálným experimentálním datům nasbíraným na mnoha světových univerzitách! K nasbírání dat byla možná použita tato online implementace Posnerovy úlohy (která má IMHO značné mouchy, např. absence inter-trial intervalu, proto si raději zkuste moji verzi;-)), možná ale nějaká úplně jiná verze na lokálních počítačích těch univerzit. Nasbíraná data lze stáhnout zde. Pokud si chcete zkusit analýzu na svém počítači, vyberte zde formát Comma Separated, jako experiment vyberte Covert Attention, pomocí CTRL+klikání si vyberte datasety (až dvacet najednou) a zvolte Download Data. Já jsem stáhl data od 389 participantů z mnoha experimentů z různých pracovišť, což je na běžné poměry experimentátorské praxe opravdu hodně dat. Náhled náhodně zvolených dat z této databáze:

Tomuto formátu, kdy každý řádek je jeden subjekt/jedno měření/jeden případ (case) a jednotlivé sloupce zahycují veškerá data spojená s tímto případem se říká wide format, široký formát, a také se hodí k některým účelům. Dobrá, data máme stažená v několika souborech s příponou .csv ve společné složce, jak je v eRku zpracujeme?

## Balicky s uzitecnymi funkcemi
library(plyr)     #Funkce na transformace dat
library(reshape2) #Dalsi funkce na transformace dat, ktere by asi zvladl i plyr
library(ggplot2)  #Tvorba peknych grafu
library(ez)       #Jednoducha statistika

## Nacist soubory s daty (jsou v "sirokem formatu", kazdy radek je jeden participant a vsechna jeho mereni)
datapath <- "d:/Dokumenty/R/Posner/"             #Kde jsou soubory s daty?
Files <- list.files(datapath, pattern = ".csv")  #Ziskame seznam souboru v te slozce s priponou csv
Data.Wide <- data.frame()                        #Vytvorime prazdny datovy ramec
for (i in 1:length(Files)){                      #Cyklime skrze seznam souboru s daty
    name <- paste(datapath, Files[i], sep="")    #Ziskame celou cestu k souboru
    x <- read.table(name, header = TRUE, sep = ",")  #Nacteme cely soubor jako tabulku (datovy ramec) do pomocne promenne
    Data.Wide <- rbind(Data.Wide, x)             #Pridame data do hlavniho datoveho ramce
}

## Nektere sloupce nepotrebujeme
Data.Wide$ClassID <- NULL
Data.Wide$Gender <- NULL
Data.Wide$Age <- NULL
Data.Wide$DateTaken <- NULL

## Data jsou bohuzel v sirokem formatu, pro analyzy je potreba prevest data na "dlouhy format"
## Nasleduje trocha magie, ktere zatim netreba rozumet
Data.Long <- melt(Data.Wide, id.vars="UserID", value.name="RT")
Data.Long$Validity <- gsub("[0-9]+", "\\1",Data.Long$variable)
Data.Long$SOA <- gsub("[^0-9]+", "\\1",Data.Long$variable)

Data.Long$RT <- Data.Long$RT * 1000  #Prevedeme reakcni casy ze sekund na milisekundy, je to hezci:-)
Data.Long$variable <- NULL

# Vycistime zjevne outliers
Outliers <- unique(subset(Data.Long, RT < 100 | RT > 1500)$UserID) # Kteri participanti meli RT mimo limity?
Data.Long <- subset(Data.Long, !(UserID %in% Outliers))            # Zbavime se vsech dat od nich

Po tomhle všem máme data v dlouhém formátu: průměrný reakční čas pro všechny participanty a kombinace IV, jak jsme to měli výše ukázáno na vymyšlených datech. Reálná, co jsem právě přetransformoval, vypadají takto:

Výsledky si zobrazíme v grafu

Jeden z nejužitečnějších pohledů na data je histogram – jak časté je zastoupení různě rychlých reakčních časů (napříč participanty)?

ggplot(Data.Long, aes(x=RT, fill = Validity)) +         # Mapovani dat
  geom_histogram(colour = NA, size=1, binwidth = 50) +  # Vykreslime histogram
  facet_wrap(Validity ~ SOA, ncol=5) +                  # pro kazdou kombinaci Validity a SOA
  theme_bw()                                            # Graficky styl

Špičatý tvar ukazuje na tzv. normální rozdělení, možná blíž k log-normálnímu, kvůli dlouhému “ocasu”. Už zde nejsou outlieři, ti byli vyloučeni dle výše uvedeného kritéria. Zajímá nás, jak hodně se od sebe reakční časy liší napříč faktory SOA a Validitou a to z těchto grafů není snadno patrné. Proto si data připravíme ve více zagregované formě (už nepotřebujeme data pro jednotlivé subjekty – můžeme je zprůměrnit dohromady) a vypočteme si několik dalších užitečných sloupců:

## Data si pripravime do formatu k zobrazeni v grafu, nyni jiz zagregujeme i subjekty
Data.Graph <- ddply(Data.Long, ~ SOA + Validity, summarize,  #Pro vsechny kombinace SOA a Validity..
                    Mean = mean(RT),                         #Prumerny reakcni cas
                    SD = sd(RT),                             #Smerodatna odchylka
                    N = length(RT),                          #Pocet zaznamu (=pocet participantu)
                    SEM = SD / sqrt(N),                      #Stredni chyba prumeru
                    CI95 = SEM * abs(qt((1-0.95)/2, N-1)),   #95% interval spolehlivosti (není upraven pro within-subject design)
                    LowerCI = Mean - CI95,                   #jeho spodni hranice
                    UpperCI = Mean + CI95)                   #a horni hranice

Data upravená pro graf vypadají takto:

Takto vykreslíme graf:

dodge <- position_dodge(width=0.2)  # Aby se nam cary moc neprekryvaly, tak je trochu posuneme
ggplot(data=Data.Graph, aes(x=SOA, y=Mean, group=Validity, colour=Validity, shape=Validity)) +
  geom_line(position=dodge) +       # Vykreslime cary
  geom_point(position=dodge) +      # Vykreslime body
  geom_errorbar(aes(ymax=LowerCI, ymin=UpperCI), position=dodge, width=0.25) +  # Vykreslime chybove usecky
  theme_bw()

A takto vypadá výsledek:

Interpretace grafu reakčních časů

Z grafu vidíme několik věcí. Validita nápovědi má zjevně vliv na průměrný reakční čas: Validní nápověď zrychluje reakce oproti neutrální nápovědi v průměru (přes všechny SOA) o nějakých 40 ms, což je v takovýchto úlohách relativně hodně. Největší (50 ms) je tento efekt při 100 ms. Invalidní nápověď má zkraje negativní dopad okolo 25 ms, ale po krátkém čase je zmenšený. Vidíme také, že se reakční čas snižuje nezávisle na validitě minimálně do SOA 400 ms. Ideální příležitost poslechnout Matúšův komentář ukázat si, jak si (podezření na) vztah mezi RT a SOA modelovat jednoduchou lineární regresí (viz níže).

Někoho by možná napadlo statisticky testovat rozdíly mezi kombinacemi faktorů “od oka”, pomocí překryvu chybových úseček. Z těch však nejde přímo vyčíst, zda jsou rozdíly mezi průměry statisticky významné. (Což je stále hlavní kritérium v (psychologické) vědě, a ač se proti tomu zvedají silné kritické hlasy, v tomto příspěvku je budu ignorovat:-) ). Toho bychom docílili úpravou intervalů spolehlivosti pro within-subject design (všechny faktory byly manipulovány “uvnitř subjektů” = každý člověk byl vystaven všem jejich kombinacím, takže je možné eliminovat variabilitu “uvnitř subjektů”), upravením magickým faktorem (2^0.5)/2 a případnou korekcí pro mnohočetná porovnání. Bez oné korekce pro multiple comparisons graf vypadá takto:

Pro srovnání, takto vypadala data získaná Posnerem (Posner et al. 1978):

Trocha statistiky

Protože graf je jen hrubý nástroj pro vyhodnocení výsledků, spočítáme si základní inferenční statistiku (zobecnění ze získaného vzorku dat). To by bylo, aby nám na takovém množství lidí všechno nevyšlo statisticky významně! Varování: styl “nasbíráme hodně dat a pak na to pustíme všechny možné statistické metody, ono něco vyleze” není zcela korektní přístup:-)

Následující tedy berte jen jako ukázku jednoduché a možná lehce explorativní statistiky. I když je pravda, že jsme do toho šli s nějakou tou konkrétní hypotézou: Např. očekáváme, že validní nápovědi ovlivňují efektivitu zpracování informací z “napovězené” (cued) lokace.

Lineární regrese

Začneme nejjednodušším přístupem: lineární regresí. Nejprve se tedy podíváme na naše z grafu plynoucí podezření, že RT lineárně klesá s rostoucí SOA:

LinModel <- lm(RT ~ SOA, data=Data.Long) # Linearni model RT = a + b * SOA
summary(LinModel) # Vypiseme ziskane koeficienty modelu

Coefficients:
             Estimate Std. Error t value Pr(>|t|)
(Intercept) 535.71438    5.74202  93.297  < 2e-16 ***  # koeficient a
SOA          -0.04116    0.01393  -2.956  0.00313 **   # koeficient b

# Model je: RT = 535.71438 - 0.04116 * SOA

Regresní model nám říká, že na každých dalších 100 ms SOA poklesne RT v průměru o 4.1 milisekundy, z počáteční hodnoty 535.7 ms. Můžeme model rozšířit a přidat kategorickou proměnou Validity:

LinModel <- lm(RT ~ SOA + Validity, data=Data.Long)
summary(LinModel)

Coefficients:
              Estimate Std. Error t value Pr(>|t|)
(Intercept) 547.237818   4.173787 131.113  < 2e-16 ***  # koeficient a
SOA          -0.028025   0.007267  -3.857 0.000116 ***  # koeficient b
ValidityN   -12.965352   5.040351  -2.572 0.010127 *    # tolik odecteme, pokud misto nevalidni napovedi mame neutralni
ValidityV   -54.215794   5.029908 -10.779  < 2e-16 ***  # tolik odecteme, pokud misto nevalidni napovedi mame validni

# V modelu byla zvolena jako zakladni uroven faktoru Validita In (nevalidni napoved) a pokud chceme znat hodnotu RT
# pro jine typy napovedi, mame na to nalezene dalsi dva koeficienty

Analýza rozptylu (ANOVA)

Dostáváme se k oné kontroverzní, avšak masově používané praktice statistického testování vůči nulové hypotéze (NHST). Provedeme tzv. analýzu rozptylu. Zjednodušeně, chceme zkusit vyvrátit nulovou hypotézu(y), že se průměrné reakční časy neliší mezi jednotlivými úrovni faktorů SOA a Validita.

## Test hlavnich efektu a interakce faktoru (Analyza rozptylu). Funkce je z balicku "ez".
ezANOVA(
  data = Data.Long,           # Zdroj dat
  dv = RT,                    # Zavisla promenna
  wid = UserID,               # Identifikace participanta
  within = .(Validity, SOA),  # Nezavisle promenne, within-subject
  type = 3)                   # Metoda sums of squares jako treba v SPSS

        Effect DFn  DFd          F            p p<.05          ges
      Validity   2  732 89.1502675 2.235731e-35     * 0.0236833456  # Hlavni efekt Validity
           SOA   4 1464  7.9209830 2.576118e-06     * 0.0034150247  # Hlavni efekt SOA
  Validity:SOA   8 2928  0.7249572 6.696396e-01       0.0006146399  # Interakce mezi SOA a Validitou

Faktory Validity i SOA mají aditivní efekty (main effect) na hodnotu reakčního času (p-hodnoty hluboce pod arbitrární magickou hranicí 0.05), což nám už ukázala lineární regrese. Není mezi nimi však významná interakce. Na oba faktory se teď můžeme podívat ve větším detailu: Zamítli jsme nulovou hypotézu o absenci vlivu validity nápovědi. Kde konkrétně se skrývají nalezené efekty? Poskytuje validní nápověď benefit oproti neutrální nápovědi, a přináší nevalidní nápověď dodatečné “náklady” oproti neutrální?

Follow-up t-testy (mnohočetné porovnání)

## Prumerne RT pro jednotlive urovne faktoru Validity
aggregate(data=Data.Long, RT ~ Validity, FUN=mean)
  Validity       RT
        In 533.6928
         N 520.0644
         V 479.1248

## Mnohocetne porovnani (s korekci) vyznamnosti rozdilu mezi temito prumery pomoci t-testu
pairwise.t.test(Data.Long$RT, Data.Long$Validity, paired=TRUE, p.adj="holm")
  In     N
N 0.0024 -
V <2e-16 <2e-16

Všechna porovnání jsou významně odlišná (p<0.05 = efekty bychom nejspíše našli i u dalších členů zkoumané populace). Nevalidní napověď tedy má negativní vliv oproti žádné nápovědi a validní má zcela zjevně pozitivní vliv na rychlost reakce na stimuly. Nyní SOA:

aggregate(data=Data.Long, RT ~ SOA, FUN=mean)
  SOA       RT
    0 526.0204
  100 514.6172
  200 508.5478
  400 503.4818
  800 502.1360

pairwise.t.test(Data.Long$RT, Data.Long$SOA, paired=TRUE, p.adj="holm")
    0       100    200    400
100 0.1077  -      -      -
200 0.0019  0.6949 -      -
400 0.0001  0.1220 0.6949 -
800 1.9e-05 0.0721 0.6949 0.7810

SOA 0 ms se odlišuje od SOA 200, 400 a 800 ms. Vypadá to tedy, že 200 ms mezi nápovědí a zobrazením cíle je dostatečný čas pro zpracování této informace, což vcelku odpovídá i Posnerovo výsledkům. Absence interakce oproti Posnerovým datům je ale trochu zvláštní: pro nulovou SOA bychom skutečně čekali, že se efekt validity nebude projevovat. V našem případě to může být dáno špatným designem onoho online experimentu (pokud byl skutečně použit k získání těchto dat). I proto si pamatujte:

Nejspíš toho už bylo až až, takže případnou další interpretaci výsledků a třeba i nějaké real-world důsledky a aplikace tohoto výzkumu necháme na jindy..

Reference:

Posner, M. I., Nissen, M. J., and Ogden, W. C. (1978). “Attended and unattended processing modes: the role of set for spatial location,” in Modes of Perceiving and Processing Information, eds H. L. Pick and N. J. Saltzman (Hillsdale, NJ: Lawrence Erlbaum Associates).