## Loading data

```{r, warning=FALSE}
data("HairEyeColor")
hairEyeColor = as.data.frame(HairEyeColor)

gombe = read.csv(file="./data/gombe_128.csv", header=TRUE)

horseKicks = read.table(file="./data/HorseKicks.txt", header=TRUE, sep="\t")

library(readxl)
folktales = read_xlsx(path="./data/rsos150645supp1.xlsx",
                      sheet=1, range="A2:JP52")
folktales = as.data.frame(folktales)
colnames(folktales)[1] = "society"
```


## Manipulating data

### Subsetting

Subsetting is one of the most common tasks when you're dealing with data. There are multiple ways to subset data, specifically when it is in a data frame.

#### Using references

The first you've already done: just select the rows or columns you want.

```{r}
horseKicks[,c(1:2)]
```

You can do so by name:

```{r}
horseKicks[c(1,5:10),c("C1","C2","C3")]
```

You can also exclude rows or columns you **don't** want by putting a `-` before the vector. Here I use a sequence to exclude odd years:

```{r}
horseKicks[-seq(1, nrow(horseKicks), 2),]
```

#### Using conditionals

Remember `which()`? You already used it for subsetting. If we just want the rows that match a certain condition, we put the statement in the row reference spot:

```{r}
gombe[which(gombe$innov > 6),]
```

Good news is that you can do this without the `which()` function and it still works:

```{r}
gombe[gombe$innov > 6,]
```

Why is this? Let's see:

```{r}
which(gombe$innov > 6)
gombe$innov > 6
```

The first returns the row indices that match that condition, and the second returns a logical vector for every row in the data frame. Remember that either can be used to reference a data frame.

You can get rows that must match more than one condition with the `&` operator:

```{r}
gombe[gombe$innov > 6 & gombe$sex == 1,]
```

Or return rows that match at least one of the conditions with the `|` (vertical pipe) "OR" operator:

```{r}
gombe[gombe$innov >= 6.5 | gombe$invt >= 6.5,]
```

You can also use (or make) another vector that contains the names of the columns you want, and keep only the values that match. This uses the useful operator `%in%`, which asks "Is X `%in%` Y?" and returns `TRUE` or `FALSE` for each element in X.

```{r}
varsToKeep = c("chimpcode","dom","soc")
head(gombe[,colnames(gombe) %in% varsToKeep])
```

#### Using `subset()`

Finally, sometimes the most straightforward way is to use the `subset()` function. You can pull up its help documentation and have a look at the syntax.

```{r}
subset(hairEyeColor, Hair == "Black" & Eye == "Blue")
```

Notice that I didn't specify the names of the arguments here (`x=`, `subset=`). If you enter your arguments in the same order they appear in the function definition, you don't have to specify which is which. For functions like subset, it's easy to tell what each argument is referring to (and `x=Hair ==...` would look really confusing). Just make sure you're entering them in the right order, or you may get unexpected results.


<big>Remember you need to **assign** your subset if you want to keep it! Otherwise your data frame is unchanged.</big>

You should always keep your original data frame and assign subsets to new data frames, so that you can always go back and see the original or subset it differently.


## Cleaning data

### Tidy data

**Three rules:**  

1. Each variable forms a column.  
2. Each observation forms a row.  
3. Each type of observational unit forms a table.  

Visit: https://opendata.fcgov.com/Neighborhood-Livability-and-Social-Health/Fort-Collins-Shelter-Service-Data/u8nn-nj59

On the top right, click the "Export" button and then the "CSV" button to download the data. Place it in your class `/data` directory.

Read it into R:

```{r}
focoShelter = read.csv(file="./data/Fort_Collins_Shelter_Service_Data.csv", header=TRUE)
```

Let's take a look at the data.

```{r}
head(focoShelter)
str(focoShelter)
colnames(focoShelter)
```

So... are these data "tidy"?

**No.** Why not?

We've got variables spread across multiple columns. One variable is whether it was a Rescue Mission shelter or a Catholic Charities shelter. Another is whether it is referring to men, women, families, or all. The other variables are split apart by each of these.

So how do we fix it? Let's consult the [Cheat Sheet](https://www.rstudio.com/wp-content/uploads/2015/02/data-wrangling-cheatsheet.pdf).

Then we need to load the `tidyr` package. If you haven't installed the `tidyverse` package (which includes `tidyr`), you can do so now, or just install `tidyr` alone.

```{r}
library(tidyr)
```

First, let's get rid of the last two columns since we don't need them.

```{r}
focoShelter = focoShelter[,-c(12:13)]
```

Right now our data are in what's called **wide** format. It's called this because all of your variables are columns, so if you have a lot of variables then your data frame will stretch out horizontally. Because of Tidy Data Rule #1, most "tidy" data will be wide.

Our goal is to get these data cleaned up and back into the wide format. But first we have to change it to **tall** format to split some of our variables apart, and then convert it back to wide.

To get our data into the **tall** format, we use the `gather()` function from the Cheat Sheet. Prior to the `tidyr` package being developed, this was called "melting" (by using the `melt()` function from `reshape2`).

```{r}
focoShelterTidy = gather(focoShelter, "variable", "value", 2:11)
```

We can use `unique()` to see all the unique values of the `variable` variable, to make sure all of our columns made it.

```{r}
unique(focoShelterTidy$variable)
```

They did!

But we can see from our variable names that not all of them are split by "Men," "Women," and "Families." Some of them represent total numbers.

So, for our next step, let's split our data frame into two: one that contains only the "Men," "Women," and "Families" variables; and another that has all the others.

One way we can do this is to use `grepl()`, the logical version of the `grep()` function, which searches for partial character matches. We use it to search for any values of `variable` that contain the text "Men," "Women," or "Families," and then we subset the data frame to the rows that contain any of them (with the `|` "OR" operator).

```{r}
focoShelterMWF = focoShelterTidy[grepl("Men", focoShelterTidy$variable, fixed=T) |
                               grepl("Women", focoShelterTidy$variable, fixed=T) |
                               grepl("Families", focoShelterTidy$variable, fixed=T),]
```

For our second subset that contains all the other rows, all we have to do is put a `!` before the `grepl()` functions to negate the logical vector that comes out of it (all "TRUE"s become "FALSE" and vice versa). This gives us the complementary set of rows to what we subset before.

```{r}
focoShelterNot = focoShelterTidy[!grepl("Men", focoShelterTidy$variable, fixed=T) &
                               !grepl("Women", focoShelterTidy$variable, fixed=T) &
                               !grepl("Families", focoShelterTidy$variable, fixed=T),]
```

Now let's make sure both of our subsets add up to the whole.

```{r}
nrow(focoShelterTidy)
nrow(focoShelterMWF) + nrow(focoShelterNot)
```

And we can take a look at our "MWF" data frame.

```{r}
head(focoShelterMWF)
str(focoShelterMWF)
colnames(focoShelterMWF)
```

Our next task is to split the "Sex" element off of each variable. We can do this easily by using the "..." part of each variable name just before the sex it's referencing. The trouble is that some of the variable names have another "..." earlier in the character string and some don't.

One way to deal with this would be to just rename the variables (you could use `recode()` or `gsub()`) to make them easier to split.

Here, I use a nasty bit of regex (["regular expressions"](https://en.wikipedia.org/wiki/Regular_expression)) to tell the function that I want it to separate based on the last occurrence of "..." for each variable name.

Some functions, like `grepl()` above, use regex but also allow you to just type in the match you want (with the `fixed=TRUE` argument). Unfortunately, `separate()` doesn't, but will often still work with plain text.

```{r}
focoShelterMWF = separate(focoShelterMWF, 2, into=c("variable","Sex"),
                          sep="\\.\\.\\.(?!.*\\.\\.\\.)")
```

After doing that, we want to merge our two subset data frames back together. But to do that, they both need to have the same structure: the same columns in the same order.

So, for the "Not" data frame, we need to add the "Sex" column and, when we do, we can assign a value of "Total" for all of those variables, since they all refer to the total number of people served regardless of sex or gender.

Then, we reorder the columns to match the "MWF" data frame.

```{r}
focoShelterNot$Sex = "Total"
focoShelterNot = focoShelterNot[,c(1,2,4,3)]
```

Let's check them out to make sure they look compatible.

```{r}
head(focoShelterMWF)
head(focoShelterNot)
```

They do. So now we can use `rbind()` to join them together row-wise (the second data frame underneath the first). If we needed to join extra columns, we could use `cbind()`, or there are a number of `dplyr` solutions for joining variables from different data sets on the Cheat Sheet.

```{r}
focoShelterTidy = rbind(focoShelterMWF, focoShelterNot)
str(focoShelterTidy)
```

Looks good. Our last step is to convert it all back to wide format again.

```{r}
focoShelterTidy = spread(focoShelterTidy, variable, value)
head(focoShelterTidy)
```

Our data are now (mostly) tidy! (Note: We actually still have one more variable to split off: the type of shelter being Rescue Mission or Catholic Charities, but we'll leave it for now.)

There are a couple of final details we could take care of, though. Our `Month` variable and `Sex` variable have gotten out of order compared to the original data frame. We can fix them by reassigning them as factors and specifying the order of the levels. For `Month`, we can use the order given in the original data frame. For `Sex` we can write it out manually.

Then, we sort the data frame by those two columns. One way, used below, is to use the `with()` and `order()` functions, but there are others.

```{r}
focoShelterTidy$Month = factor(focoShelterTidy$Month,
                               levels=as.character(focoShelter$Month))
focoShelterTidy$Sex = factor(focoShelterTidy$Sex,
                             levels=c("Men","Women","Families","Total"))
focoShelterTidy = focoShelterTidy[with(focoShelterTidy, order(Month, Sex)),]

head(focoShelterTidy)
str(focoShelterTidy)
```

There is more we could do to these data, but we will leave it for now. For one, we can see there are a lot of missing data values. This is just because none of the variables had all of those levels defined in the first place. It would probably make sense for us to sum the other values of `Sex` to the level `Total` for variables that don't have `Total` defined. That way we could use that level for analysis, and look at the other levels if we wanted finer breakdowns by `Men`, `Women`, and `Families` for the variables that have them.


## Problems with data

### Missing data

#### Finding missing data

As we have seen, some data sets can have a lot of missing data, and this can be a problem when we want to run analyses.

For another example, we can load the built-in `airquality` data set.

```{r}
data("airquality")
head(airquality)
```

We can see right away that `airquality` has some `NA` values.

A quick way to see how many missing values our data has is to use the `summary()` function, which lists the number of `NA`s in each column below some summary statistics.

```{r}
summary(airquality)
```

We can also do this for the `focoShelter` data that we just tidied.

```{r}
summary(focoShelterTidy)
```

There are other ways to find missing values in a data set. For specific columns, we can subset using the `is.na()` function, which returns `TRUE` if it finds an `NA` value.

```{r}
airquality[is.na(airquality$Ozone),]
```

If we want to see the rows that have ANY missing values across all columns, we can use `rowSums()` with `is.na()` to return rows where at least one `NA` value exists.

```{r}
airquality[rowSums(is.na(airquality)) > 0,]
```

With smaller data sets, it can also be easy to pick out `NA` values by scanning the data frame in a spreadsheet format, with `View()` (note the capital "V"):

```{r, eval=FALSE}
View(airquality)
```


#### Fixing missing data

Now that we've found our missing values, what do we do with them? Well, there usually isn't a good answer, and it depends on your data and what you plan to do with it.

**Removal**

If you need to do something about them because your analyses require you to, the most straightforward answer is to remove them. However, if you don't have many data points or they were especially expensive or difficult to collect, you may not want to throw them away.

If we want to remove values from just one variable, we can negate the statement that we used before to find the missing values to instead return the inverse where there are no missing values, using `!`.

(Note: I use `head()` below just to avoid flooding the screen with the full data set every time. It is not required for the task.)

```{r}
head(airquality[!is.na(airquality$Ozone),])
```

Likewise, we can return only rows that have zero `NA` values by slightly modifying our statement from before:

```{r}
head(airquality[rowSums(is.na(airquality)) == 0,])
```

There is also a function called `complete.cases()` that will return only the observations that have no missing values:

```{r, error=TRUE}
head(airquality[complete.cases(airquality)])
```

Oops. I forgot to include the `,` after my `complete.cases()` expression to indicate that I'm subsetting by rows and want to keep all the columns. If I fix that, it should work:

```{r}
head(airquality[complete.cases(airquality),])
```

Another way to return _only_ missing values is to negate this condition:

```{r}
airquality[!complete.cases(airquality),]
```

Let's go ahead and remove all rows with any missing values from `airquality` and see how many rows were removed.

```{r}
airqualityRem = airquality[complete.cases(airquality),]
nrow(airquality)
nrow(airqualityRem)
nrow(airquality) - nrow(airqualityRem)
(nrow(airquality) - nrow(airqualityRem)) / nrow(airquality)
```

We lost 42 rows, or 27.4% of the original data set. That's a pretty sizeable chunk. So, again, keep in mind what you're throwing away when you remove missing data casewise.

**Replacement**

Another option is replacement. Sometimes this is a good idea if you need to fix missing values by inserting the value that should have been there.

```{r}
airqualityFix = airquality
airqualityFix[is.na(airqualityFix$Ozone),][1,1] = 17
airqualityFix[is.na(airqualityFix$Solar.R),][1,2] = 142
head(airqualityFix)
```

Usually, though, you don't want to replace missing values with some other arbitrary value, because R won't treat them as missing values anymore.

```{r}
airqualityRep = airquality
airqualityRep[is.na(airqualityRep) == T] = 0
head(airqualityRep)
```

**Imputation**

Another more advanced method for dealing with missing data, usually when you really need to be able to use those values for analysis, is called imputation. There are a variety of imputation methods that calculate a value to substitute in for each missing value.

One of the most basic is to use the mean of the variable:

```{r}
airqualityImp = airquality
mean(airqualityImp$Ozone, na.rm=T)
airqualityImp$Ozone[is.na(airqualityImp$Ozone)] = mean(airqualityImp$Ozone, na.rm=T)
airqualityImp$Ozone
```

Whether this is a statistically valid solution is another question. My inclination is that it's usually a bad idea. The saying goes: "garbage in, garbage out." If you're running analyses using imputed data, you can't really trust the results that come out, because your data were fabricated in the first palce. This is especially questionable if at least a third of your data are imputed. Below that, you might be able to get away with it, but I still think it's usually a bad idea.

More sophisticated imputation methods are available in the `mice` package.


### Collinearity and confounding

Collinearity and confounding can be huge problems in your data that you should be aware of. Sometimes they're more of a theoretical concern in that they may not throw errors in your analyses, but will give you spurious results whether you're aware of them or not.

Collinearity is usually only an issue when you get to running an ANOVA or multiple regression and throw in a bunch of variables to see what best predicts some outcome variable. Collinear variables are a problem because they are highly correlated and represent the same information.

```{r}
cor(gombe$dom, gombe$dominance)
```

In the `gombe` data set, both of these variables represent dominance (`dom` is an individual attitudinal rating and `dominance` is a combination of several variables). If you were running analyses, you'd want to pick one and use it, not both.

Confounding is, in my opinion, one of the most important problems to understand about data, especially when you're analyzing social data.

The short definition of confounding is when you have variables that influence your predictors and your outcome, so that you think there is a relationship between your variables that is really caused by another factor.

Some good (hilarious) examples are available at: http://www.tylervigen.com/spurious-correlations

Not all of these are examples of confounding, but it's easy to imagine how some third factor might be causing the trends in both correlated variables.

A classic example is that homicide rates are highly correlated with ice cream sales. Why? Does ice cream make people kill each other? No, of course not! It's because both homicide rates and ice cream sales are higher in the summer months.

Confounding is common in social science studies when people are looking for significant relationships between social variables (like [genes that cause "success"](https://phys.org/news/2018-07-genetic-variants-chances-success-life.html)), and forget to think about things like ethnicity and socioeconomic status that are associated with many other demographic and life history variables.


### Outliers

```{r}
# windows()  # Command to open a new plot window on Windows
# quartz()  # Command to open a new plot window on Mac
# X11()  # Command to open a new plot window on Linux
```

```{r}
boxplot(airquality$Wind)
```

```{r}
boxplot.stats(airquality$Wind)
```

```{r}
names(boxplot.stats(airquality$Wind))
```

```{r}
boxplot.stats(airquality$Wind)$out
```

"Tukey's Fences"

Quartile plus or minus 1.5 times the interquartile range:

    Lower outlier: < Q1 - 1.5 (Q3 - Q1)  
    Upper outlier: > Q3 + 1.5 (Q3 - Q1)

```{r}
quantile(airquality$Wind)
IQR(airquality$Wind)
```

```{r}
airquality$Wind[airquality$Wind < (quantile(airquality$Wind)[[2]] - 1.5 * IQR(airquality$Wind)) |
                  airquality$Wind > (quantile(airquality$Wind)[[4]] + 1.5 * IQR(airquality$Wind))]
```

```{r}
plot(airquality$Wind, airquality$Temp)
```

```{r}
boxplot(Temp ~ Wind, data=airquality)
```

Identifying outliers when you're looking at a distribution of more than one variable is known as multivariate outlier detection. One package that handles this is `mvoutlier`.


## Exploratory data analysis

EDA is NOT "fishing," "data dredging," or "_p_-hacking." EDA is an important part of understanding your data.

### Descriptive statistics

```{r}
summary(gombe)
mean(gombe$extraversion)
median(gombe$extraversion)
```

Note: There is no built-in function to calculate the mode in R, but you can find one online or in a variety of packages.

```{r}
range(gombe$extraversion)
min(gombe$extraversion)
max(gombe$extraversion)
quantile(gombe$extraversion)
sd(gombe$extraversion)
var(gombe$extraversion)
```

We can verify that variance is the SD squared, and that SD is the square root of the variance.

```{r}
sd(gombe$extraversion)^2
sqrt(var(gombe$extraversion))
```


### Univariate plots

#### Histograms and density plots

Let's look at all the values from the `horseKicks` data together, regardless of column.

```{r, error=TRUE}
horseKicks
hist(horseKicks[,c(2:ncol(horseKicks))])
```

Weird. Okay, let's try to coerce it to numeric.

```{r, error=TRUE}
hist(as.numeric(horseKicks[,c(2:ncol(horseKicks))]))
```

For some strange reason, it's stored as a list. To get it out, we can use `unlist()`:

```{r}
unlist(horseKicks[,c(2:ncol(horseKicks))])
```

And we can use `unname()` to get just the values of the vector without the column names:

```{r}
horseKicksVals = unname(unlist(horseKicks[,c(2:ncol(horseKicks))]))
horseKicksVals
```

```{r}
hist(horseKicksVals)
hist(horseKicksVals, breaks=4)
```

```{r}
density(horseKicksVals)
```

```{r}
plot(density(horseKicksVals))
```


#### Box plots and violin plots

```{r}
boxplot(airquality$Wind)
```

A violin plot displays the same data, but shows the distribution of values on the horizontal (X) axis.

```{r, eval=FALSE}
install.packages("vioplot")
```

```{r}
library(vioplot)
```

```{r}
vioplot(airquality$Wind)
```


### Bivariate plots

#### Scatterplots

```{r}
plot(gombe$openness, gombe$extraversion)
```

We can add a regression line:

```{r}
plot(gombe$openness, gombe$extraversion)
abline(lm(gombe$extraversion ~ gombe$openness), col="red")
```

Or a smoothed LOESS/LOWESS line:

```{r}
plot(gombe$openness, gombe$extraversion)
lines(lowess(gombe$openness, gombe$extraversion), col="blue")
```

```{r}
pairs(~ dominance + extraversion + conscientiousness +
        agreeableness + neuroticism + openness,
      data=gombe)
```

#### Box plots and violin plots

```{r}
boxplot(Temp ~ Wind, data=airquality)
```

```{r}
vioplot(airquality$Wind[airquality$Month == 5],
        airquality$Wind[airquality$Month == 6],
        airquality$Wind[airquality$Month == 7])
```


### Association

```{r}
cor(gombe$openness, gombe$extraversion)
```

```{r}
cor.test(gombe$openness, gombe$extraversion)
```

```{r}
cor(gombe[,28:33])
gombeCorMat = cor(gombe[,28:33])
heatmap(gombeCorMat)
```

```{r, eval=FALSE}
install.packages("corrplot")
```

```{r}
library(corrplot)
corrplot.mixed(gombeCorMat, lower="ellipse", upper="number")
```


<br><br>


([pdf](./lecture04.pdf) / [Rmd](./lecture04.Rmd))

<div class="tocify-extend-page" data-unique="tocify-extend-page" style="height: 0;"></div>