```{r, echo=F}
library(readr)
library(dplyr)
library(ggplot2)
library(babynames)
library(ez)
library(tidyr)
```

## Project for extension of lesson 3: Comparison of baby names in Austria and the US

What about baby names in Austria?

Below, we read in data from Austria, make it comparable to the US data set, and then run a t-test to compare the popularity of baby names in the two countries.

## Reading in and looking at data

### Data from Austria

-   I downloaded data about given names in Austria from 1984 to 2022 from [data.bv.at](https://www.data.gv.at/katalog/de/dataset/stat_vornamen#resources). Download the first file from there and put it into your data folder.

Continue working in your cleaning script from lesson 3, in the \# READ DATA \#### section.

-   To read the Austrian data, try the function `read_csv()`. Add after the lines that read in the data from the USA.

```{r}
data_aut <- readr::read_csv('data/OGDEXT_VORNAMEN_1.csv')
head(data_aut)
```

You can see that the data was not read in correctly. Why?

-   Cells in CSV-files are most commonly separated by a comma, with the point used for the decimal point.
-   In some European countries (including Austria) cells are instead separated by a semicolon, and a comma used as the decimal point.
-   Here, you can see that these data from Austria are indeed separated by a semicolon.
-   If you check the help section for `?read_csv`, you see that you can use the function `read_csv2()` to read in CSV-files with the field separator `;`.

```{r}
data_aut <- readr::read_csv2('data/OGDEXT_VORNAMEN_1.csv')
head(data_aut)
```

### Data from the US

Here, I load the baby names data set, a default package included in the library "babynames". To run the code below, you first need to install the package (via the console), and then load it by adding `library(babynames)` in your section for loading packages. This dataset, as opposed to the smaller one we used in lesson 3, includes all given names in the US, not just the ones from our course participants. Assign the dataset to `data_usa`.

```{r, read in US babynames data}
# load the dataset
data(babynames)
# assign it to data_usa
data_usa = babynames
```

## Get an overview of the dataset

To get an overview of this new data set, you could use any of the following functions:

```{r}
dplyr::glimpse(data_aut)
```

-   `summary()` shows us that the earliest year in the data is 1984, the latest one 2022.
    -   It also shows us that the column "Geschlecht" (sex) is not yet coded as a categorical variable.

```{r}
summary(data_aut)
```

-   The function `View()` opens the data in an extra window.

```{r, eval = F}
View(data_aut)
```

-   The function `str()` shows you the structure of the data set.

```{r}
str(data_aut) 
```

## Cleaning and organizing the data

### Renaming columns

Let's use more useful column names for coding, and choose the same names as in the US baby names data. Check the names of the US data with `names()`:

```{r}
names(data_usa)
```

We can use the `dplyr::rename()`function to change column names.

-   You simply assign each old variable name (in quotation marks) to the new variable name.
-   Add the following in your \# CLEAN DATA \#### section:

```{r}
data_aut_clean <- data_aut %>%
  rename(year = 'C-JAHR-0', 
         sex = 'C-GESCHLECHT-0',
         living_district = 'C-WOHNBEZIRK-0',
         name = 'F-VORNAME_NORMALISIERT',
         n = 'F-ANZAHL_LGEB') 
head(data_aut_clean)
```

### Selecting rows

If we compare the US and the Austrian data set, we can see that they cover different time spans.

```{r}
summary(data_usa)
summary(data_aut_clean)
```

Let's keep only years we have in both datasets. First, filtering the Austrian dataset:

```{r}
data_aut_clean <- data_aut %>%
  dplyr::rename(year = 'C-JAHR-0', 
                sex = 'C-GESCHLECHT-0',
                living_district = 'C-WOHNBEZIRK-0',
                name = 'F-VORNAME_NORMALISIERT',
                n = 'F-ANZAHL_LGEB') %>%
  
  #keep only years up to 2017
  filter(year <= 2017)

# Check the new data frame
summary(data_aut_clean)
```

-   Next, filter out all years before 1984 from the USA data set:

```{r}
data_usa_clean = data_usa %>% 
  #keep only years from 1984 on
  filter(year > 1983)

summary(data_usa_clean)
```

### Changing the format of columns

To run tests comparing both countries, we need to join both data frames into one. There are several steps to prepare our data sets for this.

-   First, columns in both data frames need the same format.
-   The US data had the following formats (see second column):

```{r}
dplyr::glimpse(data_usa_clean)
```

So, in the Austrian data, we need sex to be a character vector with levels "M" and "F" instead of "1" and "2".

-   We will use the `ifelse()` function, within a `mutate()` call to transform the variable `sex`.
-   If you check the help for `?ifelse`, the Usage section shows you how to use it: `ifelse(test, yes, no)`.
-   The arguments for the `ifelse()` function start with a condition (rows where the column sex is equal to "1"), then specifies what happens if the condition is met (use the label "M"), and then lists what should happen otherwise/else (use the label "F").
-   Now, add this as a last line to where you created `data_aut_clean`:

```{r}
data_aut_clean <- data_aut %>% 
  dplyr::rename(year = 'C-JAHR-0', 
         sex = 'C-GESCHLECHT-0',
         living_district = 'C-WOHNBEZIRK-0',
         name = 'F-VORNAME_NORMALISIERT',
         n = 'F-ANZAHL_LGEB') %>% 
  #keep only years up to 2017
  filter(year <= 2017) %>% 
  
  # Recode the levels of sex from 1 and 2, to M and F
  mutate(sex = ifelse(sex == "1", "M", "F"))

# check if there is now F instead of 1, and M instead of M: 
glimpse(data_aut_clean)
```

### Create a new variable

Next, we need to tell statistics from Austria and the US apart, once they are in a joint data frame. Therefore, let's add a variable that identifies the country.

-   We use `rep()` to repeat "Austria" n times, and set `nrow()` to the number of rows of the Austrian data.

```{r}
data_aut_clean <- data_aut %>% 
  dplyr::rename(year = 'C-JAHR-0', 
         sex = 'C-GESCHLECHT-0',
         living_district = 'C-WOHNBEZIRK-0',
         name = 'F-VORNAME_NORMALISIERT',
         n = 'F-ANZAHL_LGEB') %>% 
  #keep only years up to 2017
  filter(year <= 2017) %>% 
  # Recode the levels of sex from 1 and 2, to M and F
  mutate(sex = ifelse(sex == "1", "M", "F")) %>%
  
  #create a country variable for Austria
  mutate(country = rep("Austria", n=nrow(data_aut)))
```

Now, let's repeat this step for the USA data. Just add a line to where you started to clean the US data:

```{r}
data_usa_clean = data_usa %>%
  #keep only years from 1984 on
  filter(year > 1983) %>% 

  #Create a country variable for the US
  mutate(country = rep("USA", n=nrow(babynames)))
```

### Grouping data and summarising statistics

The US data did not separate names per living district.

-   To make our Austrian data similar, let's group names across districts, so that we get the total number of names per year and sex.
-   Using `group_by()` first group the data by year, sex and name, because we want the total n for each of these combined categories.
    -   We also add the column country, because we want to keep it, although it does not create any additional grouping.
-   We then `summarise()` the data by calculating the `sum()` of the number of names (n).
    -   `summarise()` is similar to `mutate()`, but creates one line of output per group. Here, it collapses across districts, reducing the row number for each name.
    -   `mutate()` creates as many lines of output as input.
-   Finally, we `ungroup()` the data again.

```{r}
data_aut_clean <- data_aut %>% 
  rename(year = 'C-JAHR-0', 
         sex = 'C-GESCHLECHT-0',
         living_district = 'C-WOHNBEZIRK-0',
         name = 'F-VORNAME_NORMALISIERT',
         n = 'F-ANZAHL_LGEB') %>% 
  #keep only years up to 2017
  filter(year <= 2017) %>% 
  # Recode the levels of sex from 1 and 2, to M and F
  mutate(sex = ifelse(sex == "1", "M", "F")) %>%
  #create a country variable for Austria
  mutate(country = rep("Austria", n=nrow(data_aut))) %>%
  
  #group by year, sex and name and calculate the new sum for each category
  group_by(year, sex, name, country) %>% 
  summarise(n = sum(n)) %>%
  ungroup()
```

Finally, we also want to calculate the proportion variable that exists in the US data set (propotion of babies getting each name).

-   N, the total number of babies, is obviously different because the US is larger.
-   Using the proportion allows us to meaningfully compare the popularity of names between countries.
-   We first want to calculate the sum of all babies born in each year (`n_total`), and then divide n of each name by the years' total.
-   To calculate the sum per year, we again use `group_by()`, grouping only by year.
    -   We again add country, also this does not create a group, to keep country variable.

```{r}
data_aut_clean <- data_aut %>% 
  rename(year = 'C-JAHR-0', 
         sex = 'C-GESCHLECHT-0',
         living_district = 'C-WOHNBEZIRK-0',
         name = 'F-VORNAME_NORMALISIERT',
         n = 'F-ANZAHL_LGEB') %>% 
  #keep only years up to 2017
  filter(year <= 2017) %>% 
  # Recode the levels of sex from 1 and 2, to M and F
  mutate(sex = ifelse(sex == "1", "M", "F")) %>%
  #create a country variable for Austria
  mutate(country = rep("Austria", n=nrow(data_aut))) %>%
  #group by year, sex and name and calculate the new sum for each category
  group_by(year, sex, name, country) %>% 
  summarise(n = sum(n)) %>%  # see lesson 3 page for documentation
  ungroup()%>%
  
  #calculate total number of birth per year, and then proportion for each name out of this total
  group_by(year, country) %>% 
  mutate(n_total = sum(n), 
         prop = n/n_total) %>% 
  ungroup()

head(data_aut_clean)
```

### Filter out one name

Once this is done, you can again filter one name. I recommend taking a name that occurs in both men and women (like "Kim"), so the results of the ANOVA below make more sense.

Add another filter() call to each of the data sets for this. After this, the cleaning and data organizing is finished.

Austrian data set complete data cleaning and organizing:

```{r}
data_aut_clean <- data_aut %>% 
  rename(year = 'C-JAHR-0', 
         sex = 'C-GESCHLECHT-0',
         living_district = 'C-WOHNBEZIRK-0',
         name = 'F-VORNAME_NORMALISIERT',
         n = 'F-ANZAHL_LGEB') %>% 
  #keep only years up to 2017
  filter(year <= 2017) %>% 
  # Re-code the levels of sex from 1 and 2, to M and F
  mutate(sex = ifelse(sex == "1", "M", "F")) %>%
  #create a country variable for Austria
  mutate(country = rep("Austria", n=nrow(data_aut))) %>%
  #group by year, sex and name and calculate the new sum for each category
  group_by(year, sex, name, country) %>% 
  summarise(n = sum(n)) %>%  # see lesson 3 page for documentation
  ungroup()%>%
  #calculate total number of birth per year, and then proportion for each name out of this total
  group_by(year, country) %>% 
  mutate(n_total = sum(n), 
         prop = n/n_total) %>% 
  ungroup() %>% 
  
  #filter out your name
  filter(name == "Kim")
```

USA data set complete data cleaning and organizing:

```{r}
data_usa_clean = data_usa %>%
  #keep only years from 1984 on
  filter(year > 1983) %>% 
  #Create a country variable for the US
  mutate(country = rep("USA", n=nrow(data_usa))) %>% 
      
  #filter out your name
  filter(name == "Kim")
```

Have a look at both organized data sets:

```{r}
head(data_usa_clean)
head(data_aut_clean)
```

## Combine two data frames

Now we can combine the 2 data frames. We'll learn about dplyr verbs that work on two data frames. Until now, all verbs (filter, mutate, summarise, rename, recode) manipulated only one data frame.

-   We will use the function `full_join`: it
    -   R looks for column they have in common
    -   then looks for rows that have the same matching values
    -   combines them into a new data frame, including all the other columns from both data frames
-   `full_join()` is one of 4 Mutating joins. It keeps all observations from both data frames.
-   Check out inner_join, right_join, and left_join with the Examples provided in the help on `?full_join` when you have a moment. It takes a bit of time to understand them.

```{r}
# join the two data frames, keeping all columns and entries from both data frames
data_combined <- data_usa_clean %>% 
  full_join(data_aut_clean) 
head(data_combined)
```

-   To see data points from both countries, you could sort values by year using `arrange()`.
-   In addition, lets transform the columns sex and country into categorical variables (factors).

```{r}
data_combined <- data_usa_clean %>%
  # join the two data frames
  full_join(data_aut_clean) %>% 
  # sort by year
  arrange(year)%>%
  # transform sex into a factor
  mutate(sex = factor(sex)) %>% 
  # transform country into a factor with 2 levels:
  mutate(country = factor(country))
head(data_combined)
summary(data_combined)
```

You might want to drop the column n_total, which only has values for Austria. You can do this with `select()`:

-   It selects the columns you list.
-   if you add a minus, it deletes the columns you list

```{r}
data_combined <- data_usa_clean %>%
  # join the two data frames
  full_join(data_aut_clean) %>% 
  # sort by year
  arrange(year)%>%
  # transform sex into a factor
  mutate(sex = factor(sex)) %>% 
  # transform country into a factor with 2 levels:
  mutate(country = factor(country)) %>% 
  
  # delete the column n_total
  select(-n_total) #to delete multiple columns: -c("column_1", "column_2")

data_combined
```

## Figures script

In the figures script, add the following figure:

```{r pet-plot, fig.cap="Figure 1. Scores by pet type and country."}
# Boxplot of name propotion by sex and country
ggplot(data_combined, aes(country, prop, fill = sex)) +
  geom_boxplot() +
  scale_fill_manual(values = c("orange", "dodgerblue")) +
  labs(x = "Country", y = "Proportion") +
  theme(text = element_text(size = 20, family = "Times"))+
  theme_bw()+
  ggtitle("Popularity of the name \"Kim\" by country and sex")
```

-   In both countries, there are more women than men with the name Kim (main effect of sex).
-   There is also a main effect of country: in the USA, the name "Kim" is less popular than in Austria.
-   In addition, there is an interaction: The difference in women named "Kim" between the two countries is larger than the difference for men.

## Statistics script

In the statistics script, rerun the line with `source()` to read in the new data from the cleaning script.

### T-test by country

Is your name more popular in Austria or the USA? Although this data is not normally distributed, I want to show you how to do a t-test:

```{r}
# t-test between countries
t.test(prop ~ country, data = data_combined)
```

-   You can see that "Kim" is significantly more popular in Austria compared to the US.

```{r}
# t-test between men and women
t.test(prop ~ sex, data = data_combined)
```

-   You can see that "Kim" is much more frequent in women than men.

### Multi-factorial ANOVA: Difference in name popularity by sex and country

If you want to investigate if your name is more popular for men or women, and if this difference varies across countries, you could run an ANOVA (analysis of variance) with two factors (sex and country).

I recommend using the package "ez" for this. Load it at the top of your statistics script.

```{r, eval = F}
# Install in the console
install.packages("ez")

# In your script: 
library(ez)
```

After your other analyses (the linear regression and t-test), you can now add the ANOVA in your statistics script, using the function `ezANOVA()`. It's arguments:

-   the data frame

-   dv: the dependent variable (proportion)

-   wid: the columns that identify unique rows in your data frame. Here these are 2 columns: each name exists for each year. The function ezANOVA requires that you enter multiple columns as a list: `.("column 1", "column 2")`.

    -   A list in R is a very flexible object type that can store all kinds of objects of different dimensions. The first object could be a string, the second a data frame, the third a vector.

-   between: the predictors. ezANOVA requires that you enter multiple predictors as a list: .(x, y, z)

-   Type 3 defines which sum of squares are calculated for the ANOVA. This is detailed statistics knowledge, and only slightly changes your result. In any case, SPSS by default gives you type 3, so we use it here.

```{r}
country_sex.aov = ezANOVA(data = data_combined,
        dv = prop, 
        wid = .("year", "name"),
        # predictors
        between = .("country", "sex"),
        type = 3)
country_sex.aov
```

### Make a table with statistics and export it to a csv-file

```{r}
# Create a table with number of babyes with the name Kim per country and sex, with mean and sd: 
summary_table <- data_combined %>%
  group_by(country, sex) %>%
  #summarise: calculate one number/statistic per group created above
  summarise(
    n = n(), 
    # calculate mean and standard deviation
    mean_proportion = mean(prop), 
    sd_proportion = sd(prop), 
    # transform to percent for readibility
    mean_percent = mean_proportion*100, 
    sd_percent = sd_proportion*100
  ) %>% 
  ungroup()
# print
summary_table

# To write results to an excel table, in order to include tables into papers, for example: 
readr::write_csv(summary_table, 'write_up/descriptives_table.csv')
```

## Transforming data from long to wide format

In this section, you'll learn about a few functions that you will regularly need - although we did not need them for the analyses above.

Sometimes you need your data in a wide format, that is, multiple columns for the dependent variable, as here separate columns for Austria and the US and this for each sex.

-   SPSS and other statistics software need data in this format by default.
-   In R, you can change the format of the data frame using the function `pivot_wider()` from the tidyverse package `tidyr`.
-   Install the package, and load the library at the start of your statistics script.
-   Arguments of pivot_wider():
    -   the data frame
    -   id_cols: the column or columns that identify a single unique row in your data frame. This could be you participant ID. Here, it is each name per year.
    -   names_from: the variables for which you want to create new columns (or column/variable names)
    -   values_from: the value/variable you want to fill these columns with, here proportion (log transformed)
-   Write the following code in your section for organizing data in the statistics script.

```{r}
## ORGANIZE DATA ####

# create data frame in wide format
data_combined_wide = tidyr::pivot_wider(data_combined, 
                                        id_cols = c("year", "name"), 
                                        names_from = c("country", "sex"), 
                                        values_from = "prop")
head(data_combined_wide)
```

-   This will create an extra column per country and gender (4 columns in total)
-   The function that transforms data frames in the reverse direction, is called `pivot_longer()`. Try to back transform your wide data frame to a long one as an exercise!

You can see that the wide data frame has empty values (NA = not available) for men in the USA. our pivot_wider call added a row for each year where the name Kim occurred in any country or sex, this is why cells with no data for Kim now have an NA.You could also check how many NAs are in the column with the line of code below.

```{r}
# are there any missing values? How many?
sum(is.na(data_combined_wide$Austria_M))
```

We actually know that in this case, NA simply means there were no babies with that name. So let's replace the NAs with 0 using `ifelse()`. It's arguments:

-   The condition for `if` is that Austria_M is an `NA`, you can check this with the function `is.na()` which will give TRUE for missing values (NAs).
-   if this condition is met (TRUE), then replace with a 0
-   else, keep the value of Austria_M.

```{r}
data_combined_wide = data_combined_wide %>% 
    # replace NAs with zero
  mutate(Austria_M = ifelse(is.na(Austria_M), 0, Austria_M))
data_combined_wide
```

## Feedback welcome!

If you find any errors or things that are not explained in enough detail in this lesson, please let me know: metzler\[\@\]csh.ac.at.
