## Ecological analyses

Download the data from the annual Audubon Christmas Bird Count here: http://netapp.audubon.org/CBCObservation/Historical/ResultsByCount.aspx

Let's get all of the data available for Fort Collins: For "Start Year," select Count 1 in 1900; leave "End Year" at 2017; select "United States" and "Colorado" and flip through the pages until you find "Fort Collins" (was at the bottom of page 2 for me). Click the bubble, select CSV and Export.

Place the data in your `/data` folder.

If you open up the CSV and take a look, you'll see that the data are an absolute mess. Multiple tables are all provided one after another in the same spreadsheet. To get the data that we want, we need to skip some rows, read some rows, and skip some more rows. To make this easier, I opted to use `read_csv()` from the `readr` package, rather than the base `read.csv()` function, because it tries to figure out what the columns should be for messy data.

```{r, warning=FALSE}
library(readr)
```

```{r}
fcbird = as.data.frame(read_csv("./data/HistoricalResultsByCount [COFC-1901-2018].csv",
                    skip=208, n_max=18031))
```

```{r}
head(fcbird)
tail(fcbird)
```

Let's clean up our variables a bit.

```{r, warning=FALSE}
library(stringr)
```

```{r}
fcbird$SPEC_NAME = str_split_fixed(fcbird$COM_NAME, "\\r\\n", 2)[,2]
fcbird$SPEC_NAME = gsub("\\[|\\]", "", fcbird$SPEC_NAME)
fcbird$COM_NAME = str_split_fixed(fcbird$COM_NAME, "\r\n", 2)[,1]
fcbird$CountYear = as.integer(substr(fcbird$CountYear, 1, 4))

fcbird = fcbird[,c("COM_NAME","SPEC_NAME","CountYear","how_manyCW")]
```

```{r}
head(fcbird)
tail(fcbird)
```

Now, for the analyses we'll be doing (with the `vegan` package), we need our species as columns and our years (typically different sampling "sites") as rows. So we need to _spread_ the rows of our `COM_NAME` variable across columns, using `how_manyCW` as its values.

If we refer back to the [Data Wrangling Cheat Sheet](https://www.rstudio.com/wp-content/uploads/2015/02/data-wrangling-cheatsheet.pdf), we see that we need to use the `spread()` function from `tidyr`.

```{r, warning=FALSE}
library(tidyr)
```

```{r}
fcbirdW = spread(fcbird[,-2], "COM_NAME", "how_manyCW")
```

```{r}
fcbirdW[1:5,1:10]
```

Looks good. However, we can see there are a lot of missing values, and `vegan` can't deal with any missing values.

`complete.cases()` would remove every row, so we want to find a way to include as much of our data as we can while eliminating missing values. This is an optimization problem that R doesn't have a base function for, so I Googled it and came to [this thread](https://stackoverflow.com/questions/53613882/how-to-get-the-largest-possible-column-sequence-with-the-least-possible-row-nas) on StackOverflow.

From one of the responses, I copied the code below to find the "best" subset of the data:

```{r}
l1 = combn(2:length(fcbirdW[,-1]), 2, function(x) fcbirdW[,-1][x[1]:x[2]], simplify = FALSE)
# If you also need "combinations" of only single columns, then uncomment the next line
# l1 = c(d[-1], l1)
l2 = sapply(l1, function(x) sum(complete.cases(x)))

score = sapply(1:length(l1), function(i) NCOL(l1[[i]]) * l2[i])
best_score = which.max(score)
best = l1[[best_score]]
```
Source: [dww on StackOverflow, 12/4/18](https://stackoverflow.com/a/53616117/6408008)

And then I want to take the complete cases of those variables, so that we do have complete data with no missing values.

```{r}
rownames(best) = fcbirdW$CountYear
best = best[complete.cases(best),]
# best = apply(best, as.numeric)
best = data.frame(lapply(best, function(x) as.numeric(as.character(x))),
                  check.names=F, row.names=rownames(best))

head(best)
str(best)
```

Great. Now we can load up `vegan`.

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

```{r, warning=FALSE}
library(vegan)
```

<br>

### Diversity

```{r, eval=FALSE}
?diversity
```

```{r}
diversity(best, index="shannon")
diversity(best, index="simpson")
```

```{r}
bestDiv = diversity(best, index="simpson")

plot(as.numeric(names(bestDiv)), bestDiv)
abline(lm(bestDiv ~ as.numeric(names(bestDiv))), col="red")

cor.test(as.numeric(names(bestDiv)), bestDiv)
```

<br>

### Evenness

```{r}
diversity(best, index="shannon") / log(specnumber(best))
```

<br>

### Richness

```{r, eval=FALSE}
?rarefy
```

```{r}
rarefy(best, sample=10)
head(rarefy(best, sample=c(5, 15)))
```

```{r}
rarecurve(best)
```

We can also transpose our matrix with `t()` to look at it by species:

```{r}
rarecurve(t(best))
```

Species accumulation curves

```{r, eval=FALSE}
?specaccum
```

<br>

The `diverse` package has a number of different measures for "complex systems" research, which includes social sciences. A thorough description is available in [a published paper](https://journal.r-project.org/archive/2016/RJ-2016-033/RJ-2016-033.pdf) on the package.

<br><br>

## Ordinal data

For these exercises, let's get some Human Dimensions data for once! Go to the US Forest Service page for the 2004 visitor preference and usage data set for the Bob Marshall Wilderness Complex in Montana: https://www.fs.usda.gov/rds/archive/Product/RDS-2017-0016

At the bottom, click "Download data publication," which gives you a ZIP archive. Open it up, go into the "Data" folder and pull out both CSVs for your `/data` directory. You can hang on to the other files in the archive as well, for the metadata.

For now, let's load in the onsite data:

```{r}
bm = read.csv("./data/BMWC2004_onsitedata.csv", header=T, na.strings="88",
              stringsAsFactors=F)
```

```{r}
head(bm)
```

### Likert data

```{r}
summary(bm[,36:45])
```

<big>You can't take the mean of an ordinal variable!</big>

Why not? Remember: we can't assume the intervals between ordinal levels are equal. For example, for this survey, the participants' perceptions of the distance between "Not Important" and "Somewhat Important" may be different from the distance between "Somewhat Important" and "Very Important." So, a mean doesn't make sense.

But you can take the median.

We can also see that one of the Likert variables, `recent_f` has a 0 in it.

```{r}
bm$recent_f
```

Since there's only one 0, and the scale goes from 1 to 3 on the survey, this is likely a coding mistake. Might as well fix it and recode it as a missing value.

```{r}
bm$recent_f[bm$recent_f == 0] = NA
```

A good way of visualizing ordinal variables is through the use of a histogram.

```{r}
hist(bm$familiar)
```

A specialized package named `likert` has some additional options.

<br>

### Hypothesis testing

How do we test whether Likert responses are different by group?

#### **Permutation tests**

Permutation tests shuffle around the data to see how often the observed result occurs, to generate a _p_-value. They don't require the assumptions that normal parametric tests do, and can work regardless of the expected distribution, which is why they're good for ordinal data. This is a one-way test, but others are described in the [online textbook by Mangiafico](http://rcompanion.org/handbook/).

Functions for permutation tests in R are in the `coin` package:

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

```{r, warning=FALSE}
library(coin)
```

Let's examine whether the importance of familiarity was different for Montana residents versus visitors from elsewhere. So, we recode the `st` (state) variable to represent this. We also go ahead and declare our Likert variable as ordered, to make sure R treats it as ordinal.

```{r}
bmLik = bm
bmLik$st = factor(ifelse(bmLik$st != "MT", "Not MT", "MT"))
bmLik$familiar = ordered(bmLik$familiar)
```

We can take a look at the contingency table.

```{r}
table(bmLik$st, bmLik$familiar)
```

Difficult to tell, since the row sums are different.

```{r, eval=FALSE}
?independence_test
```

```{r}
independence_test(familiar ~ st, data=bmLik)
```

Not significant. Interesting, because we might have expected Montanans to care more about familiarity with the natural area.

As mentioned above, two-way tests, regression, etc. are available on the Mangiafico page.

<br>

### Polychoric correlations

For ordinal data, we can't use regular Pearson correlations (or Spearman, etc.). Instead, we need to calculate polychoric correlations, which assume that each variable is actually normally distributed, but represented ordinally in the data.

I like to use the `lavCor()` function in the `lavaan` package, because it can take a mix of variable types (numeric, ordinal) and calculate appropriate correlations for each. Other options are available, such as the `tetrachor()` function in the `psych` package:

```{r, eval=FALSE}
?psych::tetrachor
```

But for now we'll stick with `lavaan`, which we will also use later for structural equation modeling.

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

```{r, warning=FALSE}
library(lavaan)
```

Change all of the numeric Likert variables to ordered:

```{r}
bmLik[,36:45] = lapply(bmLik[,36:45], function(x) ordered(x))
str(bmLik)
```

And calculate our correlation matrix:

```{r, eval=FALSE}
?lavCor
```

```{r}
bmLikCor = lavCor(bmLik[,36:45])
bmLikCor
```

We can also plot it to have a look. I like `corrplot` for its different visualization options.

```{r, warning=FALSE}
library(corrplot)
```

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

<br>

### Treating ordinal data as continuous

If you have fewer than 5 levels (like we did here), **don't do it**. Your data are unlikely to meet the assumptions of the tests you want to run. You often won't get errors or warnings for doing so, and R will spit out a result, but it's statistically incorrect and your results won't mean anything.

If you have at least 5 levels and good sample size, you're usually okay. Data with 6 or 7 levels are essentially indistinguishable from continuous data. So, when you're designing a survey, go for 6 or 7.

See:

Rhemtulla, M., et al. (2012). When can categorical variables be treated as continuous? A comparison of robust continuous and categorical SEM estimation methods under suboptimal conditions. Psychological Methods, 17(3), 354. doi: [10.1037/a0029315](https://doi.org/10.1037/a0029315)

<br><br>

## Saving data

Before we go, let's save our cleaned bird count data for next time.

We can save it in CSV format, similar to the way we read CSVs in:

```{r, eval=FALSE}
?write.csv
```

```{r}
write.csv(fcbirdW, "./data/fcbirdW.csv")
write.csv(best, "./data/fcbirdbest.csv")
```

If we have a substantially larger data frame (or other object), and we know we'll only need to work with it in R or share it with others using R, we can save any R object as a compressed RDS file to save space:

```{r, eval=FALSE}
?saveRDS
```

```{r}
saveRDS(fcbirdW, "./data/fcbirdW.RDS")
saveRDS(best, "./data/fcbirdbest.RDS")
```

We could then reload it later with `readRDS()`.

Saving as RDS is also useful if you're working with and processing large files (geospatial raster layers, for example), and want to save the result to load later instead of having to do the processing steps every time.

<br><br>


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

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