Extras

R Markdown for stats homework

In stats classes (especially graduate stats classes), you often have to submit written answers with your code and output from R. This can be a bit tedious to copy and paste into a Word document, and also might not be that easy to read.

What if there was a way for you to actually complete your whole assignment in R, without any copy and pasting, and have R format your answers and make them look pretty? Good news is…there is! It’s called R Markdown.

Here is a brief 5-minute tutorial on how you can use R markdown to do your homework more efficiently.

If you would like to learn more about R markdowns, we discuss them in Workshop 1 above, and R Studio has made a convenient R Markdown cheat sheet, which you can also open straight from inside RStudio.

The function sheet

Across the workshops we build six functions of our own. They are scattered through the chapters where each one was needed, so this section collects them in one place–and, more usefully, takes each one apart so you can see how it works.

Download all six functions as an R script

Save that file next to your own script and load it with:

source("r4psych-functions.R")

Between them these functions use dplyr, magrittr and psych, so load those first.

None of this is magic. Every one of these is built out of functions the workshops already taught you–the only new idea is wrapping them up so you can reuse them.

How a function is put together

Before the six, here is the shape they all share:

my_function <- function(argument1, argument2 = default_value) {
  
  # the body: do something with the arguments
  result <- argument1 + argument2
  
  result   # the last value is what gets handed back
}

my_function <- gives the function a name, exactly like naming any other object.
function(...) lists what the function needs. argument1 is required; argument2 = default_value has a default, so you can leave it out.
{ ... } is the body–the code that runs when you call it.
The last thing evaluated is what comes back out. You can write return(result) if you prefer being explicit.

Objects created inside a function live and die there. That is the point: result above does not turn up in your environment afterwards, so functions cannot quietly clobber your data.

`var.center()`

From Workshop 2. Mean-centres a variable.

Used in the book: Workshop 2, Section 2.3.11.2, where we centre the quality-of-life items before combining them.

var.center <- function(x) {
  scale(x, scale = FALSE)
}

Piece	What it does
`function(x)`	Takes one argument: the variable you want centred
`scale(x, scale = FALSE)`	`scale()` normally does two things–subtracts the mean and divides by the standard deviation. Setting `scale = FALSE` turns off the dividing, so you get mean-centering only

Centring subtracts the mean from every value, so the new mean is exactly 0. It changes where zero sits, not the shape of the distribution. That matters most in moderation (Workshop 8), where it makes the individual coefficients interpretable at the average.

`gen_comp()`

From Workshop 2. Builds a composite score–the mean of several items–for each person.

Used in the book: Workshop 2, Section 2.3.11.2, where we turn the four car-quality items into a single index.

gen_comp <- function(data, comp, vector){
   comp <- enquo(comp)
   data %>% 
       rowwise() %>% 
       mutate(!!quo_name(comp) := mean(c(!!!vector), na.rm = TRUE)) %>% 
       ungroup()
}

Piece	What it does
`enquo(comp)`	Captures the name you typed without evaluating it. Without this, R would look for an object with that name and fail
`rowwise()`	Makes the next step work across each row rather than down a column–we want each person’s average, not the average of everyone
`!!quo_name(comp) :=`	Uses the captured name as the new column’s name. `:=` is needed because the left side is a variable, not something typed literally
`c(!!!vector)`	`!!!` unpacks a list of item names into `c()`, so `mean()` sees the items themselves
`na.rm = TRUE`	Averages the items a person did answer, rather than returning `NA`
`ungroup()`	Undoes `rowwise()`, so later steps behave normally

The enquo/!!/:= machinery is called tidy evaluation. It is what lets dplyr functions take bare column names instead of quoted strings. You do not need to master it to use the function–but this is what it is for.

`MO_Detection()`

From Workshop 5. Finds multivariate outliers using Mahalanobis distance, and returns the data without them.

Used in the book: Workshop 5, Section 5.3.5.1, where we apply it to the belonging data and then vary alpha in the practice problem.

MO_Detection = function (CompleteDataset, AnalyzedDataset, alpha = 0.001) {
  Means = colMeans(AnalyzedDataset, na.rm = T)
  Covariance = cov(AnalyzedDataset, use = "pairwise.complete.obs")
  Distances = mahalanobis(AnalyzedDataset, Means, Covariance)
  cutoff = qchisq(1-alpha, ncol(AnalyzedDataset))
  remain = Distances < cutoff
  ...
}

Piece	What it does
`CompleteDataset`, `AnalyzedDataset`	Two arguments: the whole data set, and just the numeric columns to measure distance on. It returns the whole thing minus the outliers
`alpha = 0.001`	How extreme a case must be to count. The default is the convention
`colMeans()`, `cov()`	The centre of the data, and how the variables covary–the two ingredients Mahalanobis distance needs
`mahalanobis()`	For each row, how far it sits from the centre, accounting for the correlations between variables. Two points equally far from the mean are not equally unusual if the variables are correlated
`qchisq(1-alpha, ncol(...))`	The cutoff. Mahalanobis distances follow a chi-squared distribution with degrees of freedom equal to the number of variables
`remain = Distances < cutoff`	A TRUE/FALSE vector: TRUE means “keep this row”
`factor(remain, levels = c(FALSE, TRUE), ...)`	Labels those for the plot. Stating `levels` explicitly matters–without it the function breaks whenever every case falls on the same side

`alpha_table()`

From Workshop 5. Turns psych::alpha()’s very full output into a readable table.

Used in the book: Workshop 5, Section 5.4.5, for the efficacy scale–and again in that section’s practice problem, for the well-being items.

alpha_table <- function(data, digits = 3) {
  a <- psych::alpha(data)
  data.frame(
    Item = rownames(a$item.stats),
    obs  = a$item.stats$n,
    ...
  )
}

Piece	What it does
`a <- psych::alpha(data)`	Does the actual work. Everything else is just fishing the useful pieces back out
`a$item.stats`	Per-item statistics: `n`, `r.cor` (item-test), `r.drop` (item-rest)
`a$alpha.drop`	What the scale would look like without each item: `average_r` and `raw_alpha`
`round(..., digits)`	Trims the decimals. Making `digits` an argument means you can ask for more precision without editing the function
`data.frame(...)`	Assembles the pieces into a plain table, which `knitr::kable()` can then render in any format

The pattern here is worth noticing: R functions often return far more than they print. Use str() or names() on a model object to see everything hiding inside it.

`cor_table()`

From Workshop 7. An APA-style correlation table–descriptives, lower triangle, significance stars.

Used in the book: Workshop 7, Section 7.2.1.1, where we build it and then use it with and without the describe option.

cor_table <- function(data, vars, describe = TRUE, digits = 2) {
  d  <- data[, vars, drop = FALSE]
  d  <- d[complete.cases(d), , drop = FALSE]
  ct <- psych::corr.test(d)
  stars <- ifelse(ct$p < .001, "***",
           ifelse(ct$p < .01,  "**",
           ifelse(ct$p < .05,  "*", "")))
  ...
}

Piece	What it does
`data[, vars, drop = FALSE]`	Keeps only the requested columns. `drop = FALSE` stops R helpfully turning a single column into a plain vector
`complete.cases(d)`	Drops rows with missing values, so every correlation uses the same people
`psych::corr.test(d)`	Correlations and p-values in one object–`ct$r` and `ct$p`
nested `ifelse()`	Turns p-values into the stars you see in journals. Read it as a chain: if p < .001 use `***`, otherwise check .01, otherwise .05, otherwise nothing
`paste0(formatC(...), stars)`	Glues the rounded number and its stars into one cell
`upper.tri(cells, diag = TRUE) <- ""`	Blanks the upper triangle and the diagonal–the same information twice, plus a column of 1.00s
`if (describe)`	The optional mean and SD columns. This is how one argument gives you two different tables

`partial_cor_table()`

From Workshop 7. The same table, but each correlation controls for a set of covariates.

Used in the book: Workshop 7, Section 7.2.1.2, controlling the outcome correlations for age and gender.

partial_cor_table <- function(data, vars, control.vars, digits = 2) {
  d  <- data[, c(vars, control.vars), drop = FALSE]
  d  <- d[complete.cases(d), , drop = FALSE]
  pr <- psych::partial.r(d, vars, control.vars)
  n  <- nrow(d) - length(control.vars)
  p  <- psych::corr.p(pr, n = n)$p
  ...
}

Piece	What it does
`c(vars, control.vars)`	Keeps both sets of columns–you cannot control for a variable you have not kept
`psych::partial.r(d, vars, control.vars)`	The partial correlations: the relationship between each pair after removing what the control variables explain
`nrow(d) - length(control.vars)`	Each control variable costs a degree of freedom, so we adjust the sample size before testing significance
`psych::corr.p(pr, n = n)`	Gets p-values for a correlation matrix you already have

Everything after this point is identical to cor_table()–stars, lower triangle, assembly. That is deliberate: when two jobs are nearly the same, write the second by copying the first and changing only what differs.