---
title: 'Lecture 7: Getting Data'
author: "Statistical Computing, 36-350"
date: "Monday October 5, 2015"
---

In previous episodes
===

- Seen functions to load data in passing
- Learned about string manipulation and regexp

Outline
===

- Getting data into and out of the system when it's already 
  in R format
- Import and export when the data is already very structured 
  and machine-readable
- Dealing with less structured data
- Intro to regression modeling

Reading data from R
===

- You can load and save R objects
    + R has its own format for this, which is shared across operating 
      systems
    + It's an open, documented format if you really want to pry into it
- `save(thing, file="filename")` saves the object `thing` in a file called 
  `filename` (conventional extension: `RData`)
- `save(thing1, thing2, thing3, file="filename")` saves 3 objects
- `load("name")` loads the object or objects stored in the file called 
  `filename`, with their old variable names

---

```{r}
gmp = read.table("http://www.stat.cmu.edu/~ryantibs/statcomp-F15/lectures/gmp.dat")
gmp$pop = round(gmp$gmp/gmp$pcgmp)
save(gmp,file="gmp.RData")
rm(gmp)
exists("gmp")
load(file="gmp.RData")
colnames(gmp)
```

---

- We can load or save more than one object at once; this is how RStudio will load 
your whole workspace when you're starting, and offer to save it when you're done
- Many packages come with saved data objects; there's the convenience function 
`data()` to load them

```{r}
data(cats,package="MASS")
summary(cats)
```

Non-R data tables
===

- Tables full of data, just not in the R file format
- Main function: `read.table()`
    + Presumes space-separated fields, one line per row
    + Main argument is the file name or URL
    + Returns a data frame
    + Lots of options for things like field separator, 
    column names, forcing or guessing column types, skipping 
    lines at the start of the file...
- `read.csv()` is a short-cut to set the options for reading 
  comma-separated value (CSV) files

Example: prostate cancer data 
===

Data set from 97 men who have prostate cancer

```{r}
pros = read.table("http://statweb.stanford.edu/~tibs/ElemStatLearn/datasets/prostate.data")
nrow(pros)
head(pros)
```

Writing data tables
===

- Counterpart functions `write.table()`, `write.csv()` write a data frame 
  into a file
- Drawback: takes a lot more disk space than what you get from `load` or `save`
- Advantage: can communicate with other programs, and you can even edit it 
  manually
  
Example: prostate cancer data
===

Let's create a new fake variable, and print out a new table

```{r}
pros2 = pros                                      # Make a fresh copy 
pros2$bodywt = round(runif(nrow(pros2),150,300))  # Add fake body weights
head(pros2)
#setwd("path_to_dir")                             # Set working dir (optional)
write.table(pros2, file="prostate2.txt", quote=F) # Now go look at it!
```

Less friendly data formats
===

- We learned last week how to mine text files for data, using regular 
  expressions (an extremely useful skill!)
- The `foreign` package on CRAN has tools for reading data files from lots 
  of non-R statistical software
- Spreadsheets are special

Spreadsheets considered harmful
===

- Spreadsheets look like they should be data frames
- Real spreadsheets are full of ugly irregularities
    + Values or formulas?
    + Headers, footers, side-comments, notes
    + Columns change meaning half-way down
- Ought-to-be-notorious source of errors in both industry and academia

Spreadsheets, if you have to
===

- Save the spreadsheet as a CSV; `read.csv()`
- Save the spreadsheet as a CSV; edit in a text editor; `read.csv()`
- Use `read_excel()` from the `readxl` package
- You may still need to do a lot of frustrating hacking to get things
  to work

Let's do a bit of statistical modeling
===

Once we have data loaded into R, we can do fun and interesting
statistics with it. Now we'll discuss linear modeling, perhaps the most 
fundamental tool in applied statistics. You'll learn *a lot more* in future 
classes (especially 36-401). Here's an intro: we posit a relationship
$$
Y = \beta_0 + \beta_1 X_1 + \ldots + \beta_p X_p +  \mathrm{noise}
$$
between a *response* variable $Y$ and *predictor* variables
$X_1, \ldots X_p$.
Goal is to estimate the coefficients $\beta_0, \beta_1, \ldots \beta_p$ 

For example, we might think that the log PSA (prostate-specific antigen)
of the men can be linearly predicted from the other variables in the data 
frame (Stamey et al. 1989)

- `lpsa`: log PSA score
- `lcavol`: log cancer volume
- `lweight`: log prostate weight
- `age`: age of patient
- `lbph`: log of the amount of benign prostatic hyperplasia
- `svi`: seminal vesicle invasion
- `lcp`: log of capsular penetration
- `gleason`: Gleason score 
- ` pgg45`: percent of Gleason scores 4 or 5 

---

Let's just look at the data a bit before pursuing the 
linear model

```{r}
colnames(pros) # The first 8 columns are the predictor variables
#par(ask=T)    # Have R ask us before it makes each plot, uncomment for demo
for (j in 1:8) {
  hist(pros[,j], xlab=colnames(pros)[j],
       main=paste("Histogram of",colnames(pros[j])),
       col="lightblue",breaks=20)
}
par(ask=F)     # Back to normal
```

What did we learn? A bunch of things! 

- `lcp`, the log amount of capsular penetration, is very skewed;
  a bunch of men with little (or none?), then a big spread
- `svi`, the presence of seminal vesicle invasion, is binary
- `gleason`, seems to take only values larger than 6; so how does it
  relate to `pgg45`, the percentage of Gleason scores 4 or 5?
  
(At this point, we're confused perhaps and go to our doctor friend;
we find out that `pgg45` is the percentage of Gleason scores they
received *before* their final current Gleason score, stored in `gleason`, 
that were either 4 or 5; higher Gleason score is worse, so `pgg45` tells
us something about their histories ...)

---

Now let's look at some correlations between variables

```{r}
pros.cor = cor(pros[,1:8])
round(pros.cor,3)
```

Some strong-ish correlations here. Let's find the biggest one

```{r}
pros.cor[upper.tri(pros.cor,diag=T)] = 0 # Why the lower tri part?
pros.cor.sorted = sort(abs(pros.cor),decreasing=T)
pros.cor.sorted[1]
vars.big.cor = arrayInd(which(abs(pros.cor)==pros.cor.sorted[1]), 
                        dim(pros.cor))
colnames(pros)[vars.big.cor] 
```

This is not surprising, given what we know about `pgg45` and `gleason`;
essentially what this is telling us is the following: if their Gleason 
score is high now,then they likely had a bad history of Gleason scores

---

Let's look at the second biggest correlation, since the first one
was kind of trivial

```{r}
pros.cor.sorted[2]
vars.big.cor = arrayInd(which(abs(pros.cor)==pros.cor.sorted[2]), 
                        dim(pros.cor))
colnames(pros)[vars.big.cor] 
```

<small>(Exercise: write yourself a function to extract the 
$k$th biggest absolute correlation)</small>

This is more interesting! Seems like the log amount of capsular
penetration `lcp`, and the log cancer volume `lcavol` are pretty 
related

---

Now let's plot these two variables

```{r}
plot(pros[,vars.big.cor])
```

Seems like there were a bunch of men whose capsular penetration was
really small (we knew this already) and their cancer volumes are pretty
much evenly distributed; but then, as the amount of capsular penetration 
increases, it becomes highly related to the volume of the cancer

(Should go back to our doctor friend! Ask if this make sense, ask why, etc.
He tells us that when the recorded capsular penetration value is that small, it 
basically means that it couldn't be accurately measured, which is why we see 
this spread in the cancer volume)

---

So now we want to form a linear model of the log PSA score on these 
$p=8$ predictor variables. How?  Use the `lm()` function:

```{r}
pros.lm = lm(lpsa ~ lcavol + lweight + age + 
               lbph + svi + lcp + gleason + pgg45, 
             data=pros)
summary(pros.lm)
```

Can see the estimated coefficients, and their standard errors.
More next time ...

---

Summary
===

- Loading and saving R objects is very easy
- Reading and writing data tables is pretty easy
- Dealing with other formats are harder, but not impossible
- Once we have data loaded, we can do a lot of cool things with it