---
# System prepended metadata

title: Portfolio 2

---

---
title: "Portfolio 2"
author: 
- Sofie Bøjgaard Thomsen (202206226) 
- Laura Sørine Voldgaard (202207128) 
- Ulyana Maslouskaya (202000242)  
- Daniella Varga (202204615)

output:
  pdf_document: default
  word_document: default
header-includes:
- \usepackage{fancyhdr}
- \usepackage{fvextra}
- \pagenumbering{gobble}
---

\pagestyle{fancy}
\fancyhead[HL]{Portfolio Assignment 2} 
\fancyhead[HR]{202206226, 202000242, 202207128, 202204615} 
\DefineVerbatimEnvironment{Highlighting}{Verbatim}{breaklines,commandchars=\\\{\}}

#### (DV - 202204615)
```{r}
# setup chunk
knitr::opts_chunk$set(echo = TRUE, include = TRUE, message = FALSE, warning = FALSE)
knitr::opts_knit$set(root.dir = '/work/CogSci_Methods01/Portfolio_assignment_02_SofieBT2/Logfiles')
```



```{r}
# loading packages
pacman::p_load(tidyverse)
library("moments")
library("dplyr")
library("readr")
library("pastecs") 
library("WRS2") 
library("stringr")

# clear environment.
# rm(list=ls()) # clears global workspace.

```

```{r}

# load files
files <- 
  list.files(path = "/work/CogSci_Methods01/portfolio_assignment_02-soerinev/logfiles", pattern = "*.csv", full.names = TRUE)
  
```

```{r}

# anonymize data - ONLY RUN THIS ONCE!!!
data_out <- list() 
  num_files <- length(files) 
  rand_ids <- sample(seq(1,num_files,1)) 
  cnt_f <- 0 
  for (f in files){
    cnt_f <- cnt_f + 1
    data_out[[f]] <- read_csv(file = f, col_names = TRUE)
    data_out[[f]]$ID <- paste(c("snew", rand_ids[cnt_f]), collapse = "")
    out_name <- paste(c('/work/CogSci_Methods01/portfolio_assignment_02-soerinev/logfiles', "/logfile_", unique(data_out[[f]]$ID[1]), ".csv"), collapse = "")
    write_csv(data_out[[f]], out_name, na = "NA")
    file.remove(f)
  }

```

```{r}

# compile all the logfiles into one csv file
files <- list.files(path =  getwd(), pattern = "*logfile_snew*", full.names = T) 
data <- map_dfr(files, read_csv)

```

```{r}
# removing punctuation from words in the data frame
data$word <- gsub(pattern = "[.,’“”]", replacement = "", as.character(data$word))
data$word = tolower(data$word)
names(data)[1] <- 'column_1'

# change the word length after removing this punctuation.
data$word_length <- nchar(data$word)

```

```{r}

# Prepping data for tests
data$condition <- as.factor(data$condition)
data$ID <- as.factor(data$ID)

```



### Detecting and removing outliers

```{r}

ggplot(data, aes(x = column_1, y = reading_time, colour = ID)) +
  geom_point() +
  ggtitle("Reading time for each word")
  
```


```{r}
boxplot_data <- ggplot(data, aes(x = ID, y = reading_time, fill = ID)) +
  geom_boxplot(outlier.shape = NA) +
  coord_cartesian(ylim = c(0.1, 1.5)) +
  stat_summary(fun = mean, geom = "point", shape = 23)
  
boxplot_stats <- ggplot_build(boxplot_data)$data[[1]]

lower <- print(round((boxplot_stats$lower - (1.5*(boxplot_stats$upper - boxplot_stats$lower))),4))

upper <- print(round((boxplot_stats$upper + (1.5*(boxplot_stats$upper - boxplot_stats$lower))),4))

df <- filter(data, reading_time <= max(upper), reading_time >= min(lower))

boxplot_data

```




# Correlational section
### Preparing data
```{r}
# New data frame where the rows with "haircut" and "tip" are eliminated.
df_filter <- filter(df, column_1 != 69)

#Adding transformations to the data frame
df_filter <- df_filter %>%  
  mutate(log_reading_time = log(reading_time), 
         sqrt_reading_time = sqrt(reading_time), 
         reading_time_1 = (1/reading_time))

#Creating a column with ordinal word number
ordinal_number <- df_filter$column_1 + 1
df_filter <- mutate(df_filter, ordinal_number)

```

### Testing for normality in reading time

```{r}
#Probability density histogram of reading time
  ggplot(df_filter, aes(x = reading_time)) +
  geom_histogram(aes(y = ..density..), binwidth = 0.1) +
  ggtitle("Probability Density of Reading Time") +
  stat_function(fun = dnorm, 
                args = list(mean = mean(df_filter$reading_time, na.rm = TRUE), 
                sd = sd(df_filter$reading_time, na.rm = TRUE)), 
                colour= "orange", size = 1) +
  theme_classic() + xlim(range(df_filter$reading_time))

#QQ-plot
qqnorm(df_filter$reading_time)
qqline(df_filter$reading_time)

#Shapiro-Wilk test
shapiro.test(df_filter$reading_time)
```


### Transforming the data
```{r}
df_transform <- select(df_filter, reading_time) %>%  
  mutate(log_reading_time = log(reading_time), 
         sqrt_reading_time = sqrt(reading_time), 
         reading_time_1 = (1/reading_time))

round(stat.desc(df_transform, norm = TRUE), digits = 2)[c("skew.2SE", "kurt.2SE", "normtest.p"), ]
```



### Transformation using logarithm
```{r}
ggplot(df_filter, aes(x = log_reading_time)) +
  geom_histogram(aes(y = ..density..), 
  binwidth = 0.1, 
  color= 'plum4', fill = 'lightblue2') +
  stat_function(fun = dnorm, 
                args = list( mean = mean(df_filter$log_reading_time, na.rm = TRUE), 
                sd = sd(df_filter$log_reading_time, na.rm = TRUE)),
                colour = "red", size = 1)
ggplot(df_filter, aes(sample = log_reading_time)) +
    stat_qq() +
    stat_qq_line() +
    labs(x = "Theoretical quantiles", y = "Sample quantiles")
```


### Transformation using squareroot 
```{r}
ggplot(df_filter, aes(x = sqrt_reading_time)) +
  geom_histogram(aes(y = ..density..), 
  binwidth = 0.1, 
  color= 'black', fill = 'pink') +
  stat_function(fun = dnorm, 
                args = list( mean = mean(df_filter$sqrt_reading_time, na.rm = TRUE), 
                sd = sd(df_filter$sqrt_reading_time, na.rm = TRUE)),
                colour = "blue", size = 1)
ggplot(df_filter, aes(sample = sqrt_reading_time)) +
    stat_qq() +
    stat_qq_line() +
    labs(x = "Theoretical quantiles", y = "Sample quantiles")
shapiro.test(df_filter$sqrt_reading_time)
```

### Transformation using reciprocity
```{r}
ggplot(df_filter, aes(x = reading_time_1)) +
  geom_histogram(aes(y = ..density..), 
  binwidth = 0.25, 
  color= 'plum4', fill = 'lightblue2') +
  stat_function(fun = dnorm, 
                args = list( mean = mean(df_filter$reading_time_1, na.rm = TRUE), 
                sd = sd(df_filter$reading_time_1, na.rm = TRUE)),
                colour = "darkblue", size = 1)
ggplot(df_filter, aes(sample = reading_time_1)) +
    stat_qq() +
    stat_qq_line() +
    labs(x = "Theoretical quantiles", y = "Sample quantiles")
```


### Word length vs reaction time
```{r}
#Creating a scatter plot with word length and the reading time
ggplot(df_filter, aes(x = word_length, y = reading_time)) +
  geom_point() +
  geom_smooth(method = lm, se=TRUE, colour = 'red') +
  labs(x = 'Word Length', y = 'Reading Time') +
  ggtitle('Word Length vs Reading Time') +
  theme_minimal()
```

### Testing correlational assumption
```{r}
cor.test(df_filter$word_length, df_filter$reading_time, method = 'kendall')
```


### Word frequency
```{r}
# Putting everything into one dataframe 
words_df <- read_table("word_frequency.txt", col_names = c("word", "frequency"))

df_freq <- left_join(df_filter, words_df, by = "word") 

#Creating data frame with mean reaction time of each word
mean_rt <- df_freq %>% 
          group_by(word) %>% 
          summarise(mean_rt = mean(reading_time), frequency) %>% 
          distinct()
```


### Plotting mean reading time in scatterplot
```{r}
#Creating scatter plot with assigned word frequency and reading time
ggplot(mean_rt, aes(x = log(frequency), y = mean_rt)) +
  geom_point() +
  geom_smooth(method = lm, se=TRUE, colour = 'orange') +
  labs(x = 'Word Frequency', y = 'Reading Time') +
  ggtitle('Word Frequency vs Reading Time') +
  theme_minimal()
#Note: here, we only plot reading time for unique words and not all the words in the text (i.e. we aggregated the data between subjects, but not within subjects)
```





### Testing correlational assumption
```{r}
cor.test(log(mean_rt$frequency), mean_rt$mean_rt, method = 'spearman')
```




### Additional: correlation between word frequency and word length

```{r}
ggplot(words_df, aes(x = log(words_df$frequency), y = nchar(words_df$word))) + 
  geom_point() + 
  geom_smooth(method = lm, se=TRUE, colour = 'red') + 
  labs(x = 'Word Frequency', y = 'Word Length') +
  ggtitle('Word Frequency vs Word Length') +
  theme_minimal()
  
cor.test(log(words_df$frequency), nchar(words_df$word), method = 'spearman')

```


### Ordinal word number vs reaction time
```{r}
#Creating scatter plot with ordinal word number and logged reaction time
ggplot(df_freq, aes(x = ordinal_number, y = reading_time)) +
  geom_point() +
  geom_smooth(method = lm, se=TRUE, colour = 'red') +
  labs(x = 'Ordinal Word Number', y = 'Logged Reading Time') +
  ggtitle('Ordinal Word Number vs Reading Time') +
  theme_minimal()
```

### Testing correlational assumption
```{r}
cor.test(df_freq$ordinal_number, df_freq$reading_time, method = 'spearman')
```


#### (UM - 202000242)
## Correlation section explanation

We conducted three analyses, correlating word reading times with word length, word frequency and ordinal words number. 
The assumptions of Pearson's correlation test is that both of the variables are normally distributed. Normality test on the reading time variable showed that its distribution differed from normality (W = 0.89, p-value < .000). We transformed the data using logarithm, square root and reciprocity, neither of which resulted in a normal distribution (p-value < .00). Based on this outcome, non-parametric tests were used.

Kendall's rank correlation was computed to assess the relationship between reading time and word length. There was a small positive correlation, r(9) = .05, p = .000, meaning that reading time increased with the increase of word length. The results of Spearman's correlation test for variables reading time and word frequency showed a significant negative correlation, r(82) = -0.08, p < 0.000. Thus, words with higher frequency required less reading time. 

For the analysis, word frequency variable was logarithmically transformed, which is justified by Zipf's law which states that frequencies of words in English language are inversely proportional to their rank. Given that the most frequent words in English language are the shortest, we conducted a correlation test between word length and word frequency variables. The results of Spearman's rank correlation confirmed our hypothesis that they are negatively correlated r(163) = -0.41, p = .000. Therefore, these variables are not independent and testing for one of them can predict the outcome of another. 

Finally, Spearman's correlation test was conducted for variables reading time and ordinal word number. The results showed a negative correlation, r(162) = -0.08, p < 0.000. This small but statistically significant correlation indicates that the further the participant was in the text, the less time it took for them to read each word.


#### (LSV - 202207128)
# Hypothesis section

Now we want to test if there is a significant difference in the mean reading times in the two conditions of our reading experiment. To do this we have formed two hypotheses:

Null-hypothesis H0: there is no difference in the mean reading time for the target word between the experimental and control condition.

Alternative hypothesis HA: mean reading time for the target word is higher for the experimental condition than for the control condition.


### Filtering out the target word (no. 69) and the following (no. 70) to be able to test on these:

```{r}
data_target <- filter(df, column_1 == 69)
data_following <- filter(df, column_1 == 70)
```

### Checking if the data for the target word is normally distributed to determine which t-test is appropriate. First we check the target word:

```{r}
# Making a probability density histogram of reading time on target word.
  ggplot(data_target, aes(x = reading_time)) +
  geom_histogram(aes(y = ..density..), binwidth = 0.1) +
  ggtitle("Probability Density of Reading Time for target word") +
  stat_function(fun = dnorm, 
                args = list(mean = mean(data_target$reading_time, na.rm = TRUE), 
                sd = sd(data_target$reading_time, na.rm = TRUE)), 
                colour= "orange", size = 1) +
  theme_classic() + xlim(range(data_target$reading_time))

#Making a QQ-plot of reading time on target word.
qqnorm(data_target$reading_time)
qqline(data_target$reading_time)

#Performing the Shapiro-Wilk test on reading time for target word.
shapiro.test(data_target$reading_time)
```
The Shapiro Wilk test p-value is > 0.05, which suggests that the data is normally distributed. Although, we have a very small data set and the plots don't visually indicate the same conclusion. So we will also perform Levene's test to make sure the data is normally distributed.

```{r}
# Performing Levene's test to make sure it is normally distributed.
car::leveneTest(data_target$reading_time, data_target$condition, center=mean)
```
The Levene's test, p-value > 0.05, so variances are equal and we can conclude that the data for the target word is normally distributed for further testing. 


### Now checking if the data for the following word is normally distributed to determine which t-test is appropriate:
```{r}
# Making a probability density histogram of reading time on following word.
  ggplot(data_following, aes(x = reading_time)) +
  geom_histogram(aes(y = ..density..), binwidth = 0.1) +
  ggtitle("Probability Density of Reading Time for following word") +
  stat_function(fun = dnorm, 
                args = list(mean = mean(data_following$reading_time, na.rm = TRUE), 
                sd = sd(data_following$reading_time, na.rm = TRUE)), 
                colour= "orange", size = 1) +
  theme_classic() + xlim(range(data_following$reading_time))

#Making a QQ-plot of reading time on following word.
qqnorm(data_following$reading_time)
qqline(data_following$reading_time)

#Performing the Shapiro-Wilks test on reading time for following word.
shapiro.test(data_following$reading_time)
```

For the Shapiro-Wilk test, the p-value is > 0.05, which suggests that the data is normally distributed. But once again we will perform the Levene's test to make sure. 

```{r}
# Performing Levene's test to make sure it is normally distributed.
car::leveneTest(data_following$reading_time, data_following$condition, center=mean)
```

Here, the p-value < 0.05, which indicates that variances of the data for the following word are not homogeneous, which violates the assumption of the t-test. We then choose to transform both the data for the target word and the following word to be able to perform a parametric t-test on it. Alternatively, we could have chosen to keep the data as original and perform a non-parametric t-test on it.

#### (SBT - 202206226)
### Transforming the data to see if we can make it normally distributed as a preperation for t-tests

```{r}
# transforming the target word reading time
trans_target <- data_target %>%  
  mutate(log_reading_time = log(reading_time), 
         sqrt_reading_time = sqrt(reading_time), 
         reading_time_1 = (1/reading_time))

# testing if the transformed data for the target word is normally distributed
round(stat.desc(cbind(Log = trans_target$log_reading_time, Sqrt = trans_target$sqrt_reading_time, reciprocity = trans_target$reading_time_1), basic = FALSE, norm = TRUE), digits = 3)

# transforming the following word reading time
trans_following <- data_following %>%  
  mutate(log_reading_time = log(reading_time), 
         sqrt_reading_time = sqrt(reading_time), 
         reading_time_1 = (1/reading_time))
 
# testing if the transformed data for the following word is normally distributed 
round(stat.desc(cbind(Log = trans_following$log_reading_time, Sqrt = trans_following$sqrt_reading_time, reciprocity = trans_following$reading_time_1), basic = FALSE, norm = TRUE), digits = 3)

# performing Levene's test on the following word to make sure that variance is equal
car::leveneTest(trans_following$log_reading_time, trans_following$condition, center=mean)

```

In all the transformations for both the target word data and the following word data, the Shapiro-Wilk test's p-value > 0.05, which indicates a normal distribution that we then can perform a t-test on. The results of the Levene's test for the following word indicate that the log-transformation made the variances equal.
Since we are comparing the mean reading time for the two conditions in a between-subjects design, we will use a independent samples t-test.

### Performing t-test

```{r}
# t-test for the transformed target word data
t.test(log_reading_time ~ condition, trans_target, var.equal = TRUE)

# t-test for the transformed following word data
t.test(log_reading_time ~ condition, trans_following, var.equal = TRUE )
```

## Hypothesis section explanation

For the following word, the mean reading time between the two conditions is not significantly different, t(10) = -1.54, p-value > .05. But, for the target word participants took a significantly longer time reading the semantically surprising word (M = -0.38) than the control word (M = -0.84), t(11) = -2.45, p-value < .05. 
*Note, that this mean is calculated from the logarithmic reading times.*
We have formulated a directional hypothesis where we only predict a higher mean reading time for the experimental condition, and not a lower one. And this generates the threshold for the one-tailed significance level at 0.05, which our p-value is below. Therefore these results allow us to reject the null hypothesis being "there is no difference in the mean reading time for the target word between the experimental and control condition." 

Note: 
For this assignment, we did not aggregate our data within participants. But if we had, we would have aggregated words of the same groups (i.e. length and frequency) within subjects and then words of the same group between subjects. This order of aggregation is to make sure it doesn't affect the variance. 
