Motywacja

Niewątpliwie rozwój nauk eksperymentlanych idzie w parze z rozwojem metod ilościowych oraz wizualizacji danych. W naukach takich jak ekonomia, psychologia czy socjologia, coraz częściej stosuje się tego rodzaju metody przy prowadzeniu badań. Celem tego projektu jest stworzenie odpowiedniej pomocy do tworzenia raportów przy badaniach psychologicznych. Wykorzystywana baza pytań pochodzi z badania eksperymentalnego. Projekt ten pozwoli każdemu, kto prowadzi i analizuje badanie eksperymentalne do napisania końcowego raportu. Odpowiednie części będą odpowiedzialne za: wizualizację danych, testowanie hipotez statystycznych oraz metodę bootstrap. Projekt ten umożliwi znaczne szybsze przygotowanie i analizę potrzebnych danych w projektach psychologicznych.

Wczytanie pakietów

Zacznijmy od wczytania odpowiednich pakietów:

# W tym miejscu zadeklarowane są wszystkie potrzebne biblioteki.
biblioteki = c("readr", "tidyverse", "mosaic", "knitr", "ggthemes", "hrbrthemes", "gdata", "viridis", "tm","tidytext", "broom", "stringr", "sjPlot", "sjmisc" , "wesanderson", "hunspell")  

# W tym miejscu wykonywana jest pętla, która sprawdza czy dany pakiet jest zainstalowany oraz wczytany. 
for(i in biblioteki){if(!require(i,character.only = TRUE)) install.packages(i)}
for(i in biblioteki){if(!require(i,character.only = TRUE)) library(i,character.only = TRUE) }

Analiza sentymentu – czy dobrze sformuowaliśmy pytania?

Zanim przejdziemy do analizy danych, powinniśmy sprawdzić czy nasze pytania z ankiety różnią się. Analizowane pytania pochodzą z eksperymentu psychologicnzego, wykonanego podczas zajęć z psychologii eksperymentalnej. Sam eksperyment polegał na sprawdzeniu, tego jak ludzie podejmują decyzję. W tym celu, wraz z zespołem badawczym, sprawdzaliśmy jak ludzie oceniają decyzję i ich konsekwencje. Manipulacja eksperymentalna polegała na przedstawieniu jednej grupie konsekwencji pozytywnych, a drugiej grupie konsekwencji negatywnej. Musimy zatem sprawdzić, czy odpowiednie pytania różniły się wyłącznie wartościami możliwych strat i zysków, a ich pytania nie były nacechowane pozytywnie/negatywnie. Wykorzystamy do tego text mining i analizę sentymentu. Sprawdzimy na koniec, czy pytania różnią się od siebie pod względem emocjonalnym. Jeśli nie, to uznamy, że są one dobrze sformułowane.

Dane do pobrania są dostępne pod linkiem: https://github.com/wojszymcz/danePsychologia

# Wczytujemy pliki tekstowe, w których znajdują się treści pytań przygotowanych do umieszczenia w ankiecie 
pozytywne <- readLines("Pytania_Pozytywne.txt", encoding = "UTF-8") 
negatywne <- readLines("Pytania_negatywne.txt", encoding = "UTF-8")

# Usuwamy puste linie z tekstu
pozytywne <- pozytywne[pozytywne != '']
negatywne <- negatywne[negatywne != '']

# Spójrzmy przykładowo, jak wygląda jeden z tych tekstów
head(pozytywne)
## [1] "Gracz gieldowy zaobserwowal tendencje zmian cen akcji spólki A. Posiada on obecnie pewna liczbe akcji tej spólki. Wedlug jego analizy, jesli sprzeda je w tym momencie to uda mu sie uzyskac stope zwrotu na poziomie 20%. Istnieje jednak ryzyko, oszacowane prawdopodobienstwem na 25%, ze sprzedaz akcji spowoduje strate na poziomie 30%. Gracz decyduje sie sprzedac akcje. Udaje mu sie osiagnac oczekiwany zysk."
## [2] "Gracz gieldowy zaobserwowal tendencje spadkowe spólki B. Posiada on obecnie pewna liczbe akcji tej spólki. Wedlug jego analizy, jesli sprzeda je w tym momencie to uda mu sie uniknac straty. Istnieje jednak prawdopodobienstwo oszacowane na 35%, ze akcje niedlugo z powrotem zaczna rosnac. Gracz decyduje sie zatrzymac akcje. Akcje rosna na wartosci i gracz zyskuje 15% wartosci swojego udzialu."              
## [3] "Gracz gieldowy postanowil czesc swoich akcji zamienic na akcje spólki C, która miala niedlugo wprowadzic na rynek innowacyjny system operacyjny. Ryzyko starty w zwiazku z nowymi doniesieniami o wadliwosci systemu okreslil na 80% i postanowil sprzedac akcje tracac 5% ich wartosci. Plotki okazaly sie jednak klamstwem i spólka odniosla sukces na rynku."
# Konwersja na data.frame
pozytwyne_df <- data_frame(line = 1:3, tekst = pozytywne)
negatywne_df <- data_frame(line = 1:3, tekst = negatywne)

# Na początku musimy zaimportować listę polskich stopwords, czyli słów, które nie wnoszą nic nowego do tekstu, są to np. spójniki:
pl_stop_words <- readLines("polish_stopwords.txt", encoding = "UTF-8") 
pl_stop_words <- str_split(pl_stop_words, ', ', simplify = T)
pl_stop_words <- drop(pl_stop_words)

# Data frame ze stop wordsami
stopwords_pl <- data_frame(word = pl_stop_words)

# Spójrzmy na przykłady 
head(stopwords_pl)
## # A tibble: 6 x 1
##   word      
##   <chr>     
## 1 a         
## 2 aby       
## 3 ach       
## 4 acz       
## 5 aczkolwiek
## 6 aj
# Pozbywamy się stop_wordsów z obu tekstów
pozytywne_df3 <- pozytwyne_df %>% 
  unnest_tokens(word, tekst) %>% 
  anti_join(stopwords_pl)

negatywne_df3 <- negatywne_df %>% 
  unnest_tokens(word, tekst) %>% 
  anti_join(stopwords_pl)

Sys.setenv(DICPATH = './Polski')

# Policzenie lematyzacji, czyli zamiany słów użytych w tekście na ich podstawowe formy np. jeżdzący -> jeździć, uległa -> uległy. Korzystamy z biblioteki hunspell oraz wykorzystujemy słownik języka polskiegi "pl_PL"
pozytywne_df3 <- pozytywne_df3 %>%
  mutate(word2 = hunspell_stem(word, dict = dictionary('pl_PL')))
negatywne_df3 <- negatywne_df3 %>%
  mutate(word2 = hunspell_stem(word, dict = dictionary('pl_PL')))

# Uproszczenie wyniku
pozytywne_df3$word3 <- unlist(lapply(pozytywne_df3$word2, function(x) x[1]))
negatywne_df3$word3 <- unlist(lapply(negatywne_df3$word2, function(x) x[1]))


# Wczytanie słownika sentymentu, korzystamy z pliku opisującego sentyment
pl_words_sentiment <- read.csv("pl_sent.csv", fileEncoding = 'UTF-8',
                               stringsAsFactors = F)

# Połączenie słów
slowa_sentiment_pozytywne <- inner_join(pozytywne_df3,
                              pl_words_sentiment,
                              by = c("word3" = "word"))

slowa_sentiment_negatywne <- inner_join(pozytywne_df3,
                              pl_words_sentiment,
                              by = c("word3" = "word"))

# Spójrzmy na przykłady 
head(slowa_sentiment_pozytywne)
## # A tibble: 6 x 16
##    line word     word2     word3 category mean.Happiness mean.Anger mean.Sadness
##   <int> <chr>    <list>    <chr> <chr>             <dbl>      <dbl>        <dbl>
## 1     1 zmian    <chr [1]> zmia~ F                  3.19       2.08         2   
## 2     1 cen      <chr [1]> cena  U                  1.81       3.38         2.88
## 3     1 akcji    <chr [1]> akcja U                  3.85       1.69         1.5 
## 4     1 liczbe   <chr [1]> licz~ U                  2.93       1.78         1.67
## 5     1 akcji    <chr [1]> akcja U                  3.85       1.69         1.5 
## 6     1 poziomie <chr [1]> pozi~ U                  2.74       1.52         1.52
## # ... with 8 more variables: mean.Fear <dbl>, mean.Disgust <dbl>,
## #   distance.to.H <dbl>, distance.to.A <dbl>, distance.to.S <dbl>,
## #   distance.to.F <dbl>, distance.to.D <dbl>, distance.to.N <dbl>
head(slowa_sentiment_negatywne)
## # A tibble: 6 x 16
##    line word     word2     word3 category mean.Happiness mean.Anger mean.Sadness
##   <int> <chr>    <list>    <chr> <chr>             <dbl>      <dbl>        <dbl>
## 1     1 zmian    <chr [1]> zmia~ F                  3.19       2.08         2   
## 2     1 cen      <chr [1]> cena  U                  1.81       3.38         2.88
## 3     1 akcji    <chr [1]> akcja U                  3.85       1.69         1.5 
## 4     1 liczbe   <chr [1]> licz~ U                  2.93       1.78         1.67
## 5     1 akcji    <chr [1]> akcja U                  3.85       1.69         1.5 
## 6     1 poziomie <chr [1]> pozi~ U                  2.74       1.52         1.52
## # ... with 8 more variables: mean.Fear <dbl>, mean.Disgust <dbl>,
## #   distance.to.H <dbl>, distance.to.A <dbl>, distance.to.S <dbl>,
## #   distance.to.F <dbl>, distance.to.D <dbl>, distance.to.N <dbl>
#Sprawdzamy warunek, czy pytania różnią się średniami opisującymi emocje
mean(slowa_sentiment_negatywne$mean.Happiness)==mean(slowa_sentiment_pozytywne$mean.Happiness)
## [1] TRUE
mean(slowa_sentiment_negatywne$mean.Anger)==mean(slowa_sentiment_pozytywne$mean.Anger)
## [1] TRUE
mean(slowa_sentiment_negatywne$mean.Disgust)==mean(slowa_sentiment_pozytywne$mean.Disgust)
## [1] TRUE
mean(slowa_sentiment_negatywne$mean.Fear)==mean(slowa_sentiment_pozytywne$mean.Fear)
## [1] TRUE
mean(slowa_sentiment_negatywne$mean.Sadness)==mean(slowa_sentiment_pozytywne$mean.Sadness)
## [1] TRUE

Analiza danych

Standardem studenckich badań (ale nie tylko!) jest korzystanie z platform takich jak google forms. Przy ich pomocy możemy uzyskać zbiór danych takich jak tabela poniżej. Specjalnie nic nie zostało zmienione, tak aby krok po kroku przejść przez analizę danych na podstawię bazy ściągniętej z google forms. Baza została przedstawiona poniżej DO POBRANIA

PSYCH_DANE <- read_csv("PSYCH_DANE.csv")

#Opis bazy 
#P1_E-P3_E -> Pytania z grupy eksperymentalnej
#WIEK -> Grupa wiekowa badanych
#WIEDZA -> Wiedza finansowa, podzielona na kategorie: brak, hobbystyczna, z wykszałcenia
#GRUPA -> Wybrana grupa, eksperymentalna lub kontrolna
#P1_K-P3_K -> Pytania z grupy badawczej
#PLEC -> Płeć Badanych 

PSYCH_DANE
## # A tibble: 47 x 10
##     P1_E  P2_E  P3_E PLEC      WIEK  WIEDZA GRUPA  P1_K  P2_K  P3_K
##    <dbl> <dbl> <dbl> <chr>     <chr>  <dbl> <chr> <dbl> <dbl> <dbl>
##  1     3     6     6 Mezczyzna 18-25      2 E        NA    NA    NA
##  2     6     3     3 Mezczyzna 18-25      3 E        NA    NA    NA
##  3    NA    NA    NA Kobieta   18-25      2 K         5     5     7
##  4     6     2     2 Kobieta   18-25      3 E        NA    NA    NA
##  5     7     4     4 Mezczyzna 18-25      1 E        NA    NA    NA
##  6    NA    NA    NA Mezczyzna 18-25      3 K         5     5     2
##  7     5     6     6 Kobieta   18-25      3 E        NA    NA    NA
##  8    NA    NA    NA Kobieta   18-25      3 K         6     5     6
##  9     6     5     5 Kobieta   18-25      2 E        NA    NA    NA
## 10    NA    NA    NA Mezczyzna 18-25      3 K         5     3     6
## # ... with 37 more rows

Interesujące nas dane, czyli odpowiedzi na pytania możemy teraz zwizualizować, aby uzyskać pierwsze informacje na temat naszej grupy badanej.

# Tworzymy nową bazę z samymi pytaniami, aby móc zwizualizować je porównawczo
dane <- data.frame(P1_E=PSYCH_DANE$P1_E,P2_E=PSYCH_DANE$P2_E,P3_E=PSYCH_DANE$P3_E,P1_K=PSYCH_DANE$P1_K,P2_K=PSYCH_DANE$P2_K,P3_K=PSYCH_DANE$P3_K)

#Zmieniamy bazę tak, aby można było ją zwizualizowac
dane %>% 
  gather(key = items, value = odpowiedz, na.rm=TRUE) %>% # Baza będzie pogrupowana według pytania i odpowiedzi na nie
  mutate(odpowiedz = factor(odpowiedz),
       items = factor(items), na.rm=TRUE) -> data2 # Zapisujemy w bazie data2

# Tworzymy wykres słupkowy - rozkład odpowiedzi w każdej z prób
ggplot(data2, aes(x = items)) +
  geom_bar(aes(fill = odpowiedz), position = "fill") +
  coord_flip()+
  theme_ipsum()+
  labs(y="udział badanych w próbie", x="pytania")+
  scale_x_discrete(limits = rev(levels(data2$items))) -> p2
p2 + scale_fill_manual(values = wes_palette("BottleRocket1")) #Dodajemy estetykę z palety z filmów Wesa Andersona

# Możemy również wykorzytsać pakiet do wizualizacji skali Likerta
mydf <- find_var(dane, pattern = "P", out = "df") # Parametr p odpowiada za parametr poszukiwany w bazie -  my poszukujemy pytań. 
plot_likert(
  mydf, # Baza danych
  grid.range = c(1.2, 1.4), # Podziałka siatki na wykresie 
  expand.grid = FALSE,
  values = "sum.outside", # Wyświetlenie etykiet danych poza wykresem
  show.prc.sign = TRUE, # Czy pokazywać znak procentowy
  geom.colors = wes_palette("IsleofDogs1"),# Kolor z palety wes_palette
  title = "Rozkład odpowiedzi w grupach" # Tytuł 
  
)

# Spróbujemy zwizualizować jeszcze informacje o wieku 

 ggplot(PSYCH_DANE)+ 
     geom_bar(aes(x=WIEK, fill= WIEK))+
     scale_fill_manual(values = wes_palette("Rushmore1"))+
     theme_ipsum()+
     ylab("liczba wystąpień")+
     ggtitle("Rozkład wieku w grupie badanej")

#spróbujmy opisać jeszcze płeć osób badanych

  ggplot(subset(PSYCH_DANE, !is.na(PLEC)))+ 
     geom_bar(aes(x=PLEC, fill= PLEC), na.rm = TRUE)+
     scale_fill_manual(values = wes_palette("GrandBudapest2"))+
     theme_ipsum()+
     ylab("liczba wystąpień")+
     ggtitle("Rozkład płci w grupie badanej")

Pobrane dane często pobrane są w niewygodnej do wizualizacji formie. Dobrym przykładem, jest używana w psychologii skala Likerta. Odpowiedzi od “Zdecydowanie się zgadzam” do “Zdecydowanie się nie zgadzam” zapisywane są w formacie 1-5, lub 1-7, co jest bardzo niewygodne do wizualizacji. W naszej bazie ten problem wystąpił również w przypadku wykształcenia.

# Jak wyglądają dane bez zmiany 

ggplot(PSYCH_DANE)+
  geom_bar(aes(x=WIEDZA, fill=WIEDZA))+
  theme_ipsum()+
     ylab("liczba wystąpień")+
     ggtitle("Rozkład wiedzy finansowej w grupie badanej")

# Nic nie mówią nam informacje na temat 1,2,3. 
# zmienimy zatem strukturę danych. 
# na potrzeby zadania stworzymy nową bazę

wyksztalcenie <- c("Brak", "Hobbystyczne", "Z wykształcenia")
dane_wyksztalcenie <- PSYCH_DANE

for(i in 1:length(wyksztalcenie)){
dane_wyksztalcenie$WIEDZA[dane_wyksztalcenie$WIEDZA==i]<- wyksztalcenie[i]
}

#Zwizualizujemy ponownie dane 

ggplot(dane_wyksztalcenie)+
  geom_bar(aes(x=WIEDZA, fill=WIEDZA))+
  theme_ipsum()+
  scale_fill_manual(values=wes_palette("Darjeeling2"))+
  ylab("liczba wystąpień")+
  ggtitle("Rozkład wiedzy finansowej w grupie badanej")

Taki wykres jest o wiele bardziej informatywny i potencjalnie przydatny przy pisaniu raportu

Możemy też napotkać sytuację, kiedy niekoniecznie chcemy zmienić wartości w naszej bazie. W przypadku, gdy będziemy chcieli liczyć statystyki opisowe, a jednocześnie przedstawić na wykresie czytelne dane nie chcielibyśmy zmieniać samych danych. Aby to zrobić możemy edytować samą legendę.

# Posłużymy się rozkładem odpowiedzi 
# Wczytujemy wektor z treścią odpowiedzi w kolejności wartości 1-7
mylevels<-c("niepoprawna i niewybaczalna", "niepoprawna, biorąc wszystko pod uwagę", "niepoprawna, ale nie jest nierozsądna", "równie dobra jak decyzja przeciwna",  "poprawna, ale przeciwna decyzja też byłaby rozsądna","poprawna, biorąc wszystko pod uwagę","w oczywisty sposób poprawna, inna decyzja byłaby niedopuszczalna") 

ggplot(data2, aes(x = items)) + # Wczytujemy dane
  geom_bar(aes(fill = odpowiedz), position = "fill") +
  coord_flip()+ # Odwracamy
  theme_ipsum()+ # Wczytujemy motyw
  labs(y="udział badanych w próbie", x="pytania")+ # Oznaczenie osi
  scale_x_discrete(limits = rev(levels(data2$items))) -> p2
p2 + scale_fill_manual(values = wes_palette( "Rushmore1",7, type = "continuous"), labels=mylevels) # Do labels, oznaczeń legendy, wczytujemy wektor odpowiedzi.

Czasem poza grupą badaną interesuje nas też charakterystyka grup spełniających wybrane kryteria, na przykład chcemy zobaczyć jak wygląda rozkład wieku, ale tylko wśród kobiet.

#zobaczymy jak wygłądał rozkład wieku wśród kobiet w naszej grupie badanej

PSYCH_DANE %>% #wybieramy pozdzbiór 
  filter(PLEC=="Kobieta") %>% #filtrujemy
  ggplot()+ #resztę wykonujemy analogicznie do przypadków wcześniejszych 
  geom_bar(aes(x=WIEK, fill= WIEK))+
  scale_fill_manual(values = wes_palette("Chevalier1"))+
  theme_ipsum()+
  ylab("liczba wystąpień")+
  ggtitle("Rozkład wieku wśród kobiet z grupy badanej")

Skoro już zrozumieliśmy jak wygląda nasza grupa badawcza, to możemy teraz sprawdzić odpowiednie statystyki opisowe i zwizualizować rozkład odpowiedzi.

# Wybieramy  odpowiednie pytania
pytania_pozytywne <- c(PSYCH_DANE$P1_E, PSYCH_DANE$P2_E, PSYCH_DANE$P3_E) 
pytania_pozytywne <- na.omit(pytania_pozytywne) # Pomijamy puste wartości
pytania_pozytywneDF <- data.frame(pytania_pozytywne) # Konwertujemy na data frame

# Analogicznie dla drugiej grupy
pytania_negatywne <- c(PSYCH_DANE$P1_K, PSYCH_DANE$P2_K, PSYCH_DANE$P3_K)
pytania_negatywne <- na.omit(pytania_negatywne)
pytania_negatywneDF <- data.frame(pytania_negatywne)

# Tworzymy statystyczne podsumowanie naszej próby
pytania_pozytywne_pod <- pytania_pozytywneDF %>%
  summarize(sample_size = n(), # Liczba wystąpień
    mean = mean(pytania_pozytywne), # Średnie wynik - kontrowersyjne (skala Likerta)
    sd = sd(pytania_pozytywne), # Odchylenie standardowe
    minimum = min(pytania_pozytywne), # Wartość najmniejsza
    lower_quartile = quantile(pytania_pozytywne, 0.25), # Pierwszy kwantyl
    median = median(pytania_pozytywne), # Mediana
    upper_quartile = quantile(pytania_pozytywne, 0.75), # Trzeci kwantyl
    max = max(pytania_pozytywne)) # Wartość największa
kable(pytania_pozytywne_pod) # Podsumowanie
sample_size mean sd minimum lower_quartile median upper_quartile max
75 4.653333 1.278654 2 4 5 6 7 
# Analogiczne podsumowanie dla drugiej grupy
pytania_negatywne_pod <- pytania_negatywneDF %>%
  summarize(sample_size = n(),
    mean = mean(pytania_negatywne),
    sd = sd(pytania_negatywne),
    minimum = min(pytania_negatywne),
    lower_quartile = quantile(pytania_negatywne, 0.25),
    median = median(pytania_negatywne),
    upper_quartile = quantile(pytania_negatywne, 0.75),
    max = max(pytania_negatywne))
kable(pytania_negatywne_pod)
sample_size mean sd minimum lower_quartile median upper_quartile max
66 4.681818 1.618728 1 3 5 6 7 
# Tworzymy histogram z rozkładem odpowiedzi
pytania_pozytywneDF %>% 
  ggplot(aes(x = pytania_pozytywne)) +# Wczytujemy baze
  geom_histogram(binwidth = 1, fill = "#1DACE8", color="#1C366B")+ #wybieramy szerokość kolumn w histogramie, u nas jeden - skala dyskretna z zakresu 1-7, tym razem wczytamy kolor oraz obramowanie, żeby pokazać jak to zrobić
  geom_vline(aes(xintercept=mean(pytania_pozytywne)),color="black",linetype="dashed")+ # Dodajemy wartość średnią w formie linii przerywanej
  labs(title="Rozkład odpowiedzi w pytaniach pozytwynych",x="Ocena decyzji", y = "Liczba wystąpień")+ # Oznaczenia osi 
  theme_ipsum()+ # Dodajemy motyw
  xlim(1,7) # Dodajemy ograniczeni osii OX

# Analogicznie w przypadku drugiej grupy

pytania_negatywneDF %>% ggplot(aes(x = pytania_negatywne)) +
  geom_histogram(binwidth = 1, fill = "#1DACE8", color="#1C366B")+geom_vline(aes(xintercept=mean(pytania_negatywne)), color="black",linetype="dashed")+
  labs(title="Rozkład odpowiedzi w pytaniach negatywnych",x="Ocena decyzji", y = "Liczba wystąpień")+ 
  theme_ipsum()+
  xlim(1,7)

# Tworzymy wspólną bazę do wizualizacji rozkładów

pytania <- data.frame(
  type = c( rep("pozytywne", length(pytania_pozytywne)), rep("negatywne", length(pytania_negatywne)) ), # Powtarzamy w polu "type" słowo pozytywne tyle razy ile jest danych w wektorze, a następnie analogicznie dla negatywnych słow. W ten sposób otrzymujemy data frame z wynikami dla obu grup
  value = c(pytania_pozytywne, pytania_negatywne )
)

# Wizualizujemy analogicznie jak wyżej
ggplot(pytania, aes(x = value, fill = type)) +                    
  geom_histogram(position = "identity", alpha = 0.7, bins = 7)+
  scale_fill_manual(values=c("#F8DF4F", "#541F12")) +
  theme_ipsum() +
  labs(fill="") +
  xlab("Ocena decyzji")+
  ylab("Liczba wystąpień")

Weryfikacja hipotez statystycznych

Spróbujemy teraz zweryfikować hipotezę na poziomie a=0.1. H0: średnie w obu grupach są takie same H1: średnia w grupie z odpowiedziami pozytywnymi jest wyższa niż w grupie z odpowiedziami negatywnymi. 

# Wybieramy test t-studenta, i dwa zbiory danych. Zaznaczamy alternatywę, że pierwszy zbiór danych jest większy od drugiego
t<-t.test(pytania_pozytywne, pytania_negatywne, alternative ="g" )
#Patrzymy na dane
t 
## 
##  Welch Two Sample t-test
## 
## data:  pytania_pozytywne and pytania_negatywne
## t = -0.11486, df = 123.32, p-value = 0.5456
## alternative hypothesis: true difference in means is greater than 0
## 95 percent confidence interval:
##  -0.4394848        Inf
## sample estimates:
## mean of x mean of y 
##  4.653333  4.681818
#patrzymy na wartość p-value
t$p.value 
## [1] 0.5456291
#widzimy, że wielkość próby jest niewielka, co może zaważać na jakości estymacji
length(pytania_negatywne)
## [1] 66
length(pytania_pozytywne)
## [1] 75

Jak widzimy, nie mamy podstaw do odrzucenia H0, ponieważ otrzymane p-value jest zdecydowanie większe niż przyjęty poziom istotności. Widzimy jednak, że w każdej z grup otrzymaliśmy bardzo mało odpowiedzi na pytania. Spróbujemy zatem zweryfikować hipotezę za pomocą metody bootstrap.

Bootstrap

Meteoda bootstrap polega na wielokrotnym losowaniu ze zwracaniem z próby. Polega ona na stworzeniu nowego rozkładu na podstawie otrzymanych zmiennych. Metoda ta, może być bardzo pomocna, kiedy mamy do czynienia ze skalą likerta, czyli z ograniczonym zbiorem danych do 7 możliwych wyników. Jeśli nasza próba byłaby reprezentatywna, to ostatecznie poprzez replikację wyników, moglibyśmy otrzymać bardziej wiarygodną estymację.

Dokonamy teraz bootstrapu, oraz ponowanie zwizualizujemy dane na temat rokzładów.

# Dokonujemy bootsrapów
pytania_pozytywne_boot<- resample(pytania_pozytywne,replace = TRUE,  10000) 
# Losujemy ze zwraceniem 10,000 razy pytania pozytywne
pytania_negatywne_boot <- resample(pytania_negatywne, replace=TRUE, 10000)

# Losujemy ze zwraceniem 10,000 razy pytania negatywne
pytania_boot <- data.frame(pytania_pozytywne_boot, pytania_negatywne_boot)

# Podsumowujemy statystyki dla "bootowanych" danych, analogicznie jak we wcześniejszych przypadkach
pytania_pozytywne_boot_pod <- pytania_boot %>%
  summarize(sample_size = n(),
    mean = mean(pytania_pozytywne_boot),
    sd = sd(pytania_pozytywne),
    minimum = min(pytania_pozytywne_boot),
    lower_quartile = quantile(pytania_pozytywne_boot, 0.25),
    median = median(pytania_pozytywne_boot),
    upper_quartile = quantile(pytania_pozytywne_boot, 0.75),
    max = max(pytania_pozytywne_boot))
kable(pytania_pozytywne_boot_pod)
sample_size mean sd minimum lower_quartile median upper_quartile max
10000 4.6542 1.278654 2 4 5 6 7 
# Podsumowujemy statystyki dla "bootowanych" danych, analogicznie jak we wcześniejszych przypadkach
pytania_negatywne_boot_pod <- pytania_boot %>%
  summarize(sample_size = n(),
    mean = mean(pytania_negatywne_boot),
    sd = sd(pytania_negatywne_boot),
    minimum = min(pytania_negatywne_boot),
    lower_quartile = quantile(pytania_negatywne_boot, 0.25),
    median = median(pytania_negatywne_boot),
    upper_quartile = quantile(pytania_negatywne_boot, 0.75),
    max = max(pytania_negatywne_boot))
kable(pytania_negatywne_boot_pod)
sample_size mean sd minimum lower_quartile median upper_quartile max
10000 4.6806 1.60276 1 3 5 6 7

Widzimy, że zbiory danych nie różnią się znacznie od wcześniejszych zbiorów. Zobaczmy jak wyglądają na jednym wykresie.

# Analogiczna wizualizacja jak w przypadku nie-bootwanych danych 

pytania_boot %>% ggplot(aes(x = pytania_pozytywne_boot)) +
  geom_histogram(binwidth = 1, fill = "#A35E60", color="#541F12")+geom_vline(aes(xintercept=mean(pytania_pozytywne_boot)), color="black",linetype="dashed")+
  labs(title="Rozkład odpowiedzi w pytaniach pozytwynych",subtitle="Metoda bootstrap",x="Ocena decyzji", y = "Liczba wystąpień")+ 
  theme_ipsum()+
  xlim(1,7)

pytania_boot %>% ggplot(aes(x = pytania_negatywne_boot)) +
  geom_histogram(binwidth = 1, fill = "#A35E60", color="#541F12")+geom_vline(aes(xintercept=mean(pytania_negatywne_boot)), color="black",linetype="dashed")+
  labs(title="Rozkład odpowiedzi w pytaniach negatywnych",subtitle="Metoda bootstrap",x="Ocena decyzji", y = "Liczba wystąpień")+ 
  theme_ipsum()+
  xlim(1,7)

Jak widzimy, rozkłady te wyglądają bardzo podobnie do tych przed dokonaniem resamplingu. Trzeba pamiętać, że metoda bootstrap nie dostarcza nam żadnych nowych informacji na temat próby, a pozwala jedynie otrzymać lepsze wartości estymacji. Sprawdzimy w takim razie na podstawie testu ANOVA, czy rozkłady przed resamplingiem, a po resamplingu można uznać za takie same. Porównamy najpierw rozkłady przy użyciu wykresu pudełkowego.

# Tworzymy data frame, gdzie będą zarówno pytania zbootowane, jak i oryginalne.
pytania_porownanie <- data.frame(type = c( rep("pozytywne", length(pytania_pozytywne)), rep("negatywne", length(pytania_negatywne)), rep("pozytywne_boot", length(pytania_pozytywne_boot)), rep("negatywne_boot", length(pytania_negatywne_boot))),
  value = c(pytania_pozytywne, pytania_negatywne, pytania_pozytywne_boot, pytania_negatywne_boot )
) 

#Wizualizujemy efekt bootsrapu
pytania_porownanie %>% 
  ggplot( aes(x=type, y=value, fill=type)) + # Dane
  geom_boxplot() + # Typ wykresu, tym razem wykres pudełkowy - dobry do porównania rozkładów
  scale_fill_stata() + # Wybieramy kolory ze staty
  geom_jitter(color="black", size=0.4, alpha=0.3) + # Wybieramy kolory i przezroczystość punktów 
  theme_ipsum() + # Wczytujemy motyw
  theme(
    legend.position="none",
    plot.title = element_text(size=11)
  ) +
  ggtitle("Wykres pudełkowy - rokzłady odpowiedzi")+ # Tytuł wykresu 
  xlab("") + # Opis osi OX
  ylab("ocena decyzji") # Opis osi OY

Annova_results <- aov(value ~ type, data= pytania_porownanie) # Test ANOVA, sprawdzamy zależność wyników od typu zbioru 
summary(Annova_results) # Patrzymy na wyniki testu ANOVA
##                Df Sum Sq Mean Sq F value Pr(>F)
## type            3      4   1.171   0.558  0.643
## Residuals   20137  42281   2.100

Widzimy, że poziom p-value jest wysoki. Oznacza to, że nie mamy podstaw do odrzucenia hipotezy zerowej o równości średniej. Widzimy, że bootsrap dał nam lepsze oszacowanie w tym teście. Sprawdzimy jeszcze jak wygląda to w teście T-studenta.

t1 <- t.test(pytania_pozytywne_boot, pytania_negatywne_boot, alternative = "greater", var.equal = FALSE)
t1$p.value
## [1] 0.9011672

Widzimy, że p-value jest wysokie, więc brak podstaw do odrzucenia hipotezy zerowej o równości średnich. Co więcej otrzymane p-value jest znacznie większe, niż to otrzymane przed bootstrapem. Widzimy, że metoda ta, znacznie zwiększyła jakość naszej estymacji.

Podsumowanie

Ostateczne wyniki testu T-studenta zamykają nasze rozważania na temat analizy danych w eksperymentach psychologicznych. Stworzony projekt umożliwia szybsze napisane raportu z badania. Dzięki niemu usprawnieniu ulegnie część odpowiedzialna za analizę danych, wizualizację danych, manipulację danymi, weryfikację hipotez statystycznych. Projekt pozwala nam równocześnie na lepszą estymację badanego zjawiska przy metodzie bootstrap.

Bibliografia

  • Kopczewski T. , Ćwiakowski P. , Materiały dydaktyczne z zajęć " Wizualizacja danych i Raportowanie danych w R" , Wydział Nauk Ekonomicznych UW, Semestr Letni 2020/2021.
  • Sienicka Z., Kuzio T., Szymczak W., Badanie wpływu znajomości konsekwencji na ocenę decyzji. Praca Zaliczeniowa w ramach Psychologii Eksperymentalnej, Wydział Filozofii UW.