• ggplot() specifies what data to use and what variables will be mapped to where
• inside ggplot(), aes(x = , y = , color = ) specify what variables correspond to what aspects of the plot in general
• layers of plots can be combined using the + at the end of lines
• use geom_line() and geom_point() to add lines and points
• sometimes you need to add a group element to aes() if your plot looks strange
• make sure you are plotting what you think you are by checking the numbers!
• facet_grid(~variable) and facet_wrap(~variable) can be helpful to quickly split up your plot

• the factor class allows us to have a different order from alphanumeric for categorical data
• we can change data to be a factor variable using mutate(), as_factor() (in the forcats package), or factor() functions and specifying the levels with the levels argument
• fct_reorder({variable_to_reorder}, {variable_to_order_by}) helps us reorder a variable by the values of another variable
• arranging, tabulating, and plotting the data will reflect the new order

We will cover how to use R to compute some of basic statistics and fit some basic statistical models.

• Correlation
• T-test
• Linear Regression / Logistic Regression

There are plenty of resources online for learning more about these methods.

Check out www.opencasestudies.org for deeper dives on some of the concepts covered here and the resource page for more resources.

The correlation coefficient is a summary statistic that measures the strength of a linear relationship between two numeric variables.

• The strength of the relationship - based on how well the points form a line
• The direction of the relationship - based on if the points progress upward or downward

source

See this case study for more information.

Function cor() computes correlation in R.

cor(x, y = NULL, use = c("everything", "complete.obs"),
    method = c("pearson", "kendall", "spearman"))

• provide two numeric vectors of the same length (arguments x, y), or
• provide a data.frame / tibble with numeric columns only
• by default, Pearson correlation coefficient is computed

Function cor.test() also computes correlation and tests for association.

cor.test(x, y = NULL, alternative(c("two.sided", "less", "greater")),
    method = c("pearson", "kendall", "spearman"))

• provide two numeric vectors of the same length (arguments x, y), or
• provide a data.frame / tibble with numeric columns only
• by default, Pearson correlation coefficient is computed
• alternative values:
  • two.sided means true correlation coefficient is not equal to zero (default)
  • greater means true correlation coefficient is > 0 (positive relationship)
  • less means true correlation coefficient is < 0 (negative relationship)

First, we compute correlation by providing two vectors.

# x and y must be numeric vectors
y1980 <- yearly_co2_emissions %>% pull(`1980`)
y1985 <- yearly_co2_emissions %>% pull(`1985`)

Like other functions, if there are NAs, you get NA as the result. But if you specify use only the complete observations, then it will give you correlation using the non-missing data.

cor(y1980, y1985, use = "complete.obs")

[1] 0.9936257

cor.test(y1980, y1985)

    Pearson's product-moment correlation

data:  y1980 and y1985
t = 114.59, df = 169, p-value < 0.00000000000000022
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.9913844 0.9952853
sample estimates:
      cor 
0.9936257 

The broom package helps make stats results look tidy

library(broom)
cor_result <- tidy(cor.test(y1980, y1985))
glimpse(cor_result)

Rows: 1
Columns: 8
$ estimate    <dbl> 0.9936257
$ statistic   <dbl> 114.5851
$ p.value     <dbl> 0.00000000000000000000000000000000000000000000000000000000…
$ parameter   <int> 169
$ conf.low    <dbl> 0.9913844
$ conf.high   <dbl> 0.9952853
$ method      <chr> "Pearson's product-moment correlation"
$ alternative <chr> "two.sided"

We can compute correlation for all pairs of columns of a data frame / matrix. This is often called, “computing a correlation matrix”.

Columns must be all numeric!

co2_subset <- yearly_co2_emissions %>%
  select(c(`1950`, `1980`, `1985`, `2010`))

head(co2_subset)

# A tibble: 6 × 4
  `1950` `1980` `1985` `2010`
   <dbl>  <dbl>  <dbl>  <dbl>
1   84.3   1760   3510   8460
2  297     5170   7880   4600
3 3790    66500  72800 119000
4   NA       NA     NA    517
5  187     5350   4700  29100
6   NA      143    249    524

We can compute correlation for all pairs of columns of a data frame / matrix. This is often called, “computing a correlation matrix”.

cor_mat <- cor(co2_subset, use = "complete.obs")
cor_mat

          1950      1980      1985      2010
1950 1.0000000 0.9228253 0.8818288 0.5415047
1980 0.9228253 1.0000000 0.9935477 0.7270839
1985 0.8818288 0.9935477 1.0000000 0.7827256
2010 0.5415047 0.7270839 0.7827256 1.0000000

corrplot package can make correlation matrix plots

library(corrplot)
corrplot(cor_mat)

source

The commonly used are:

• one-sample t-test – used to test mean of a variable in one group
• two-sample t-test – used to test difference in means of a variable between two groups (if the “two groups” are data of the same individuals collected at 2 time points, we say it is two-sample paired t-test)

The t.test() function in R is one to address the above.

t.test(x, y = NULL,
       alternative = c("two.sided", "less", "greater"),
       mu = 0, paired = FALSE, var.equal = FALSE,
       conf.level = 0.95, ...)

t.test(y1980, y1985)

    Welch Two Sample t-test

data:  y1980 and y1985
t = -0.090533, df = 341, p-value = 0.9279
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -95902.79  87462.97
sample estimates:
mean of x mean of y 
 109769.0  113988.9 

The broom package has a tidy() function that can organize results into a data frame so that they are easily manipulated (or nicely printed)

result <- t.test(y1980, y1985)
result_tidy <- tidy(result)
glimpse(result_tidy)

Rows: 1
Columns: 10
$ estimate    <dbl> -4219.909
$ estimate1   <dbl> 109769
$ estimate2   <dbl> 113988.9
$ statistic   <dbl> -0.09053303
$ p.value     <dbl> 0.9279168
$ parameter   <dbl> 340.999
$ conf.low    <dbl> -95902.79
$ conf.high   <dbl> 87462.97
$ method      <chr> "Welch Two Sample t-test"
$ alternative <chr> "two.sided"

• wilcox.test() – Wilcoxon signed rank test, Wilcoxon rank sum test
• shapiro.test() – Shapiro test
• ks.test() – Kolmogorov-Smirnov test
• var.test()– Fisher’s F-Test
• chisq.test() – Chi-squared test
• aov() – Analysis of Variance (ANOVA)

• Use cor() to calculate correlation between two vectors, cor.test() can give more information.
• corrplot() is nice for a quick visualization!
• t.test() one sample test to test the difference in mean of a single vector from zero (one input)
• t.test() two sample test to test the difference in means between two vectors (two inputs)
• tidy() in the broom package is useful for organizing and saving statistical test output
• Remember to adjust p-values with p.adjust() when doing multiple tests on data

🏠 Class Website

💻 Lab

Linear regression is a method to model the relationship between a response and one or more explanatory variables.

Most commonly used statistical tests are actually specialized regressions, including the two sample t-test, see here for more.

Here is some of the notation, so it is easier to understand the commands/results.

\[ y_i = \alpha + \beta x_{i} + \varepsilon_i \] where:

• \(y_i\) is the outcome for person i
• \(\alpha\) is the intercept
• \(\beta\) is the slope (also called a coefficient) - the mean change in y that we would expect for one unit change in x (“rise over run”)
• \(x_i\) is the predictor for person i
• \(\varepsilon_i\) is the residual variation for person i

To fit regression models in R, we use the function glm() (Generalized Linear Model).

You may also see lm() which is a more limited function that only allows for normally/Gaussian distributed error terms (aka typical linear regressions).

We typically provide two arguments:

• formula – model formula written using names of columns in our data
• data – our data frame

Model formula \[ y_i = \alpha + \beta x_{i} + \varepsilon_i \] In R translates to

y ~ x

Model formula \[ y_i = \alpha + \beta x_{i} + \varepsilon_i \] In R translates to

y ~ x

In practice, y and x are replaced with the names of columns from our data set.

For example, if we want to fit a regression model where outcome is income and predictor is years_of_education, our formula would be:

income ~ years_of_education

Model formula \[ y_i = \alpha + \beta_1 x_{i1} + \beta_2 x_{i2} + \beta_3 x_{i3} + \varepsilon_i \] In R translates to

y ~ x1 + x2 + x3

In practice, y and x1, x2, x3 are replaced with the names of columns from our data set.

For example, if we want to fit a regression model where outcome is income and predictors are years_of_education, age, and location then our formula would be:

income ~ years_of_education + age + location

We will use our the calenviroscreen dataset from the dasehr package to examine how traffic estimates predict diesel particulate emissions.

The summary() function returns a list that shows us some more detail

summary(fit)

Call:
glm(formula = DieselPM ~ TrafficPctl, data = calenviroscreen)

Coefficients:
              Estimate Std. Error t value             Pr(>|t|)    
(Intercept) 0.04245151 0.00529915   8.011   0.0000000000000013 ***
TrafficPctl 0.00363651 0.00009177  39.625 < 0.0000000000000002 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for gaussian family taken to be 0.05614936)

    Null deviance: 537.24  on 7999  degrees of freedom
Residual deviance: 449.08  on 7998  degrees of freedom
  (35 observations deleted due to missingness)
AIC: -330.9

Number of Fisher Scoring iterations: 2

The broom package can help us here too!

The estimate is the coefficient or slope.

for one change in the traffic percentile, we see 0.003637 more Diesel particulate emissions. The error for this estimate is relatively small at 0.00009. This relationship appears to be significant with a small p value < 2e-16.

tidy(fit) %>% glimpse()

Rows: 2
Columns: 5
$ term      <chr> "(Intercept)", "TrafficPctl"
$ estimate  <dbl> 0.042451513, 0.003636512
$ std.error <dbl> 0.00529915058, 0.00009177366
$ statistic <dbl> 8.011003, 39.624789
$ p.value   <dbl> 0.0000000000000012983396857238743156732797112774494962140342…

Can also use tidy and glimpse to see the output nicely.

fit2 %>%
  tidy() %>%
  glimpse()

Rows: 3
Columns: 5
$ term      <chr> "(Intercept)", "TrafficPctl", "Ozone"
$ estimate  <dbl> 0.23025068, 0.00355094, -3.77418894
$ std.error <dbl> 0.01347753532, 0.00009067243, 0.24967137612
$ statistic <dbl> 17.08403, 39.16229, -15.11663
$ p.value   <dbl> 0.0000000000000000000000000000000000000000000000000000000000…

Factors get special treatment in regression models - lowest level of the factor is the comparison group, and all other factors are relative to its values.

Let’s create a variable that tells us whether a census tract has a high, middle, or low percentage of the population below the poverty line.

calenviroscreen <- calenviroscreen %>% mutate(
  PovertyPctl_level = case_when(
    PovertyPctl > 0.75 ~ "high",
    PovertyPctl > 0.25 & PovertyPctl <= 0.75 ~ "middle",
    PovertyPctl <= 0.25 ~ "low",
    TRUE ~ NA
  )
)

Generalized linear models (GLMs) allow for fitting regressions for non-continuous/normal outcomes. Examples include: logistic regression, Poisson regression.

Add the family argument – a description of the error distribution and link function to be used in the model. These include:

• binomial(link = "logit") - outcome is binary
• poisson(link = "log") - outcome is count or rate
• others

Very important to use the right test!

See this case study for more information.

See ?family documentation for details of family functions.

See this case study for more information.

An odds ratio (OR) is a measure of association between an exposure and an outcome. The OR represents the odds that an outcome will occur given a particular exposure, compared to the odds of the outcome occurring in the absence of that exposure.

Check out this paper.

Some final notes:

• Researcher’s responsibility to understand the statistical method they use – underlying assumptions, correct interpretation of method results

• Researcher’s responsibility to understand the R software they use – meaning of function’s arguments and meaning of function’s output elements

• glm() fits regression models:
  • Use the formula = argument to specify the model (e.g., y ~ x or y ~ x1 + x2 using column names)
  • Use data = to indicate the dataset
  • Use family = to do a other regressions like logistic, Poisson and more
  • summary() gives useful statistics
• oddsratio() from the epitools package can calculate odds ratios (outside of logistic regression - which allows more than one explanatory variable)
• this is just the tip of the iceberg!

For more check out:

• this chapter on modeling in this tidyverse book
  
• this chart on when to do what test
• opencasestudies.org

Content for similar topics as this course can also be found on Leanpub.

🏠 Class Website

💻 Lab

Image by Gerd Altmann from Pixabay