Testing Simple Slopes in Multilevel Models with a Continuous Predictor

A friend asked me how to test simple slopes from interaction effects in multilevel models when you have a continuous predictor, since the only packages that seem to compute them are for one-level/fixed effects models. Here I’ll show you how to (1) test whether simple slopes are different from 0, and (2) test whether simple slopes are different from each other.

A big thing to note up front—the random intercepts or slopes in your model have no impact on how you compute the simple slopes for your fixed effects. You can use this method for normal regression too.

tl;dr

You can recode your predictors so that the slopes in your model are the simple slopes you want to test.

Motivating example

I want to test whether people who are more Agreeable empathize more. But I also think the effect of Agreeableness on empathy is moderated by how strongly the target of empathy expresses their emotions (are they crying or just a little bummed?) and how serious the target’s situation is (did their mother die or did they just spill their coffee?).

I run an experiment with a 2 (Expression: strong vs. weak) x 2 (Situation: serious vs. trivial) mixed design, with Situation as the within-subjects factor. Every research subject sees two targets—one in a serious situation, one in a trivial situation—but both targets express the same amount of emotion—strong or weak. I measure subjects’ empathy for the targets and their Agreeableness.

Here are my sample data.

# load the packages we need
library(dplyr)
library(car)
library(lme4)

set.seed(123) # set the seed so you can get the same results

SubjectID <- rep(1:140, each = 2) 
expression <- rep(c('strong', 'weak'), each = 140)
situation <- rep(c('serious', 'trivial'), times= 140)
agree <- rep(sample(1:7, 140, replace = TRUE), each = 2)
empathy <- sample(1:7, 280, replace = TRUE)

emp.data <- data.frame(SubjectID, expression, situation, agree, empathy)

Here’s a sample of what the data look like

emp.data[c(1:4, 141:144), ]

##     SubjectID expression situation agree empathy
## 1           1     strong   serious     3       3
## 2           1     strong   trivial     3       3
## 3           2     strong   serious     6       1
## 4           2     strong   trivial     6       2
## 141        71       weak   serious     6       6
## 142        71       weak   trivial     6       1
## 143        72       weak   serious     5       3
## 144        72       weak   trivial     5       2

Testing Main Effects and Interactions

First, I want to test the main effects and interaction effects. I could do this with repeated-measures ANOVA, but for this example I’ll do it with multilevel modeling, with a random intercept nested within Subjects.

To test this model, I need to mean-center Agreeableness and I need to compute contrasts for the two experimental conditions.

emp.data$agree.cent <- emp.data$agree - mean(emp.data$agree)
emp.data$expr.cont <- recode(emp.data$expression, as.factor.result = FALSE, "'strong' = 1; 'weak' = -1")
emp.data$sit.cont <- recode(emp.data$situation, as.factor.result = FALSE, "'serious' = 1; 'trivial' = -1")

Now I can run my multilevel model.

m <- lmer(empathy ~ expr.cont*sit.cont*agree.cent + (1|SubjectID), data = emp.data) # I'm using %>% from dplyr to make this example easier to read
summary(m) %>% # get the model summary
  coef() %>% # get the coefficients from the model summary
  round(digits = 2) # round the coefficients to two digits

##                               Estimate Std. Error t value
## (Intercept)                       3.89       0.12   33.19
## expr.cont                         0.09       0.12    0.77
## sit.cont                         -0.05       0.12   -0.42
## agree.cent                        0.01       0.06    0.16
## expr.cont:sit.cont               -0.09       0.12   -0.74
## expr.cont:agree.cent             -0.08       0.06   -1.34
## sit.cont:agree.cent               0.07       0.06    1.19
## expr.cont:sit.cont:agree.cent    -0.07       0.06   -1.15

In this model, the slope for expr.cont is the main effect of Expression, the slope for sit.cont is the main effect of Situation, the slope for expr.cont:agree.cent is the interaction of Expression and Agreeableness, and so on.

Testing Whether Simple Slopes are Different from 0

Perhaps I’d like to test whether there is an effect of Agreeableness on empathy in the trivial Situation/strong Expression condition. How can I do this?

We can find the full effect of Agreeableness by taking all the slopes from the model that include Agreeableness and multiplying them by the variables that contribute to those slopes. I’m going to use the letter B to represent these slopes instead of using the actual values from our model, so that the final lesson will make more sense.

B_Agree * Agreeableness + B_Agree*Exp * Agreeableness * Expression + B_Agree*Sit * Agreeableness * Situation + B_{Agree*Exp*Sit} * Agreeableness * Expression * Situation

When we factor out Agreeableness we get

Agreeableness * (B_Agree + B_Agree*Exp * Expression + B_Agree*Sit * Situation + B_{Agree*Exp*Sit} * Expression * Situation)

What’s the effect of Agreeableness in the strong/trivial condition? We can find this by plugging in the proper values for that condition. Currently, the Expression variable is coded so 1 = strong and -1 = weak, and the Situation variable is coded so 1 = serious and -1 = trivial. In the strong/trivial condition we have Expression = 1 and Situation = -1

When we plug these into our equation we get

Agreeableness * (B_Agree + B_Agree*Exp * 1 + B_Agree*Sit * -1 + B_{Agree*Exp*Sit} * 1 * -1)
= Agreeableness * (B_Agree + B_Agree*Exp – B_Agree*Sit – B_{Agree*Exp*Sit})

Great, we can compute the effect of Agreeableness in the strong/trivial condition. But is it significantly different from 0? We don’t know, because we don’t have significance tests for sums and differences of slopes. To get past this problem, we’re going to replace the current Expression and Situation variables with dummy codes that have 0’s for the strong and trivial conditions and 1’s for the weak and serious conditions. This will change the meaning of our slopes in a useful way.

emp.data$expr.dum <- recode(emp.data$expression, as.factor.result = FALSE, "'strong' = 0; 'weak' = 1")
emp.data$sit.dum <- recode(emp.data$situation, as.factor.result = FALSE, "'serious' = 1; 'trivial' = 0")

Now let’s run our model again with these dummy variables.

m.dummy <- lmer(empathy ~ expr.dum*sit.dum*agree.cent + (1|SubjectID), data = emp.data) summary(m.dummy) %>% 
  coef() %>% 
  round(digits = 2)

##                             Estimate Std. Error t value
## (Intercept)                     4.11       0.23   17.56
## expr.dum                       -0.35       0.33   -1.06
## sit.dum                        -0.27       0.33   -0.82
## agree.cent                     -0.07       0.12   -0.62
## expr.dum:sit.dum                0.35       0.47    0.74
## expr.dum:agree.cent             0.02       0.17    0.13
## sit.dum:agree.cent              0.00       0.16    0.03
## expr.dum:sit.dum:agree.cent     0.27       0.23    1.15

Now again, the full effect of Agreeableness is

Agreeableness * (B_Agree + B_Agree*Exp * Expression + B_Agree*Sit * Situation + B_{Agree*Exp*Sit} * Expression * Situation)

But when we plug in the new values for the strong Expression/trivial Situation condition, we get

Agreeableness * (B_Agree + B_Agree*Exp * 0 + B_Agree*Sit * 0 + B_{Agree*Exp*Sit} * 0 * 0) \
= Agreeableness * (B_Agree)

In this model, the slope for Agreeableness is no longer the main effect of Agreeableness; instead, it’s the effect of Agreeableness in the strong Expression/trivial Situation condition, and we can get the significance test of the simple slope.

Testing Whether Simple Slopes are Different from Each Other

Now, perhaps I’d like to know whether the effect of Agreeableness on empathy in the trivial Situation/strong Expression condition is different from the effect of Agreeableness on empathy in the serious Situation/strong Expression condition. How can I do this?

In our original model, our Situation and Expression variables represent three orthogonal contrasts:

	C1	C2	C3
Strong Expression/Trivial Situation	1	1	1
Strong Expression/Serious Situation	1	-1	-1
Weak Expression/Trivial Situation	-1	1	-1
Weak Expression/Serious Situation	-1	-1	1

C1 is the main effect of Expression, C2 is the main effect of Situation, and C3 is the Expression x Situation interaction. We know these are orthogonal contrasts because if we multiple any two columns row-wise, the sum of the product will equal 0.

Now let’s say we wanted to test the difference in empathy between the strong/trivial condition and the strong/serious condition. We could replace those three contrasts with a different set of orthogonal contrasts.

	C1	C2	C3
Strong Expression/Trivial Situation	1	0	1
Strong Expression/Serious Situation	-1	0	1
Weak Expression/Trivial Situation	0	1	-1
Weak Expression/Serious Situation	0	-1	-1

Now C1 tests whether there’s a difference between the strong/trivial and strong/serious conditions, C2 tests whether there’s a difference between the weak/trivial and weak/serious conditions, and C3 tests the main effect of Expression.

We can test whether there’s a difference in the effect of Agreeableness on empathy between the strong/trivial and strong/serious conditions (and between the weak/trivial and weak/serious conditions) by entering the sum of these contrasts in the model and multiplying by Agreeableness.

emp.data$C1 <- with(emp.data, ifelse(expression == 'strong' & situation == 'trivial', 1, 
                              ifelse(expression == 'strong' & situation == 'serious', -1,
                                     0)))
emp.data$C2 <- with(emp.data, ifelse(expression == 'weak' & situation == 'trivial', 1, 
                              ifelse(expression == 'weak' & situation == 'serious', -1,
                                     0)))
emp.data$C3 <- with(emp.data, ifelse(expression == 'strong', 1, -1))

m <- lmer(empathy ~ agree.cent*(C1 + C2 + C3) + (1 | SubjectID), data = emp.data) summary(m) %>% 
  coef() %>% 
  round(digits = 2)

##               Estimate Std. Error t value
## (Intercept)       3.89       0.12   33.19
## agree.cent        0.01       0.06    0.16
## C1                0.14       0.17    0.82
## C2               -0.04       0.17   -0.22
## C3                0.09       0.12    0.77
## agree.cent:C1     0.00       0.08   -0.03
## agree.cent:C2    -0.14       0.08   -1.63
## agree.cent:C3    -0.08       0.06   -1.34

The agree.cent:C1 interaction tests whether the effect of Agreeableness is different in the strong/trivial and strong/serious condition.

How to Test Moderation when X is Continuous and M is Categorical with 3+ Levels

A colleague asked how he could test a model with a continuous predictor, a continuous dependent variable, and a moderator with 4 levels. For this example, I’m using real data where I had subjects read a letter from someone about a misfortune. I asked subjects either to take the other person’s perspective (perspective taking condition), to remain objective (objective condition), or I didn’t ask them to do anything (control condition). I also measured dispositional perspective taking using the Interpersonal Reactivity Index (IRI).

Let’s test a model where dispositional perspective taking is predicting subjects’ sympathy, and the experimental condition is moderating that effect. And let’s start with the easy case where we only have two categories in the moderator—perspective taking and objective.

# Remember to always mean-center your continuous predictor when 
# you're testing interactions, which is what moderation does

ptdata$iri_pt_centered <- ptdata$iri_pt - mean(ptdata$iri_pt, na.rm = TRUE)

# I'm going to change the factor for experimental condition to 
# a contrast code, which doesn't need to be centered
# I'll use the recode() function from the car package

library(car)
ptdata$pt_vs_objective <- recode(ptdata$ptfactor, as.factor.result = FALSE,
 "'perspective taking' = 1; 'objective' = -1")

m1 <- lm(symp ~ iri_pt_centered * pt_vs_objective, data = ptdata)
summary(m1)

Call:
lm(formula = symp ~ iri_pt_centered * pt_vs_objective, data = ptdata2)

Residuals:
    Min      1Q Median     3Q    Max
-3.3842 -0.5246 0.3288 0.7496 1.5048

Coefficients:
                                 Estimate Std. Error t value Pr(>|t|)
(Intercept)                      3.912979   0.091323  42.848  < 2e-16 ***
iri_pt_centered                  0.196312   0.254364   0.772  0.44170
pt_vs_objective                  0.305775   0.091323   3.348  0.00107 **
iri_pt_centered:pt_vs_objective -0.008959   0.254364  -0.035  0.97196
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.027 on 125 degrees of freedom
(4 observations deleted due to missingness)
Multiple R-squared: 0.09038, Adjusted R-squared: 0.06855
F-statistic: 4.14 on 3 and 125 DF, p-value: 0.007778

Here the test of moderation is the t test of the interaction term. So I would say that subjects sympathized more in the perspective taking condition than in the objective condition, B = .31, SE = .09, t(125) = 3.35, p = .001. However, dispositional perspective taking was unrelated to sympathy, B = .20, SE = .25, t(125) = .77, p = .44, and there was no moderation, B = -.01, SE = .25, t(125) = .04, p = .97.

So what if we include that third category? We have to do two things.

First, instead of having just one contrast code for the factors, we’ll have two (because you need one less contrast code than the number of groups to represent the full study design), and we’ll interact both of them with dispositional perspective taking (but not with each other).

Second, we have to use R’s routine for hierarchical regression, which involves saving one model that excludes the interaction terms, and one model that includes the interaction terms. Then we’ll compare the two models by using the anova() function.

# Again, always mean-center your continuous predictor when 
# you're testing interactions
# I already did this, but I'll do it again just to hammer in the point

ptdata$iri_pt_centered <- ptdata$iri_pt - mean(ptdata$iri_pt, 
  na.rm = TRUE)

# Now I'm going to change the factor for experimental condition to 
# two orthogonal contrast codes, which don't need to be centered

library(car)

# this first contrast tests whether the perspective taking and 
# control conditions are different from each other, and ignores 
# the objective condition 
ptdata$pt_vs_control <- recode(ptdata$ptfactor, as.factor.result = FALSE, 
  "'perspective taking' = 1; 'control' = -1; 'objective' = 0")

# the second contrast tests whether the objective condition is 
# different from the other two combined
ptdata$objective_vs_others <- recode(ptdata$ptfactor, 
  as.factor.result = FALSE, 
  "'perspective taking' = -1; 'control' = -1; 'objective' = 2")

# now make a model with no interactions
no.interactions <- lm(symp ~ iri_pt_centered + pt_vs_objective + 
  objective_vs_others, data = ptdata)

# then make a model with the interactions
yes.interactions <- lm(symp ~ iri_pt_centered * (pt_vs_objective + 
  objective_vs_others), data = ptdata)

# then compare the two of them with anova()
anova(no.interactions, yes.interactions)

Analysis of Variance Table

Model 1: symp ~ iri_pt_centered + pt_v_control + objective_v_others
Model 2: symp ~ iri_pt_centered * (pt_v_control + objective_v_others)
  Res.Df    RSS Df Sum of Sq     F Pr(>F)
1    192 200.35                          
2    190 198.80  2    1.5507 0.741  0.478

The F test on line 2 is for the difference in the two models, and it’s equivalent to the omnibus interaction tests in ANOVA with only categorical predictors. So I would conclude that there was no interaction between dispositional perspective taking and experimental condition, F(2, 190) = .74, p = .48.

If you look at the summary for the model with interactions, you can also see the separate slopes testing whether dispositional perspective taking was moderated by whether subjects were in the perspective taking vs. control group (our first contrast) and whether dispositional perspective taking was moderated by whether subjects were in the objective vs. other two groups (our second contrast). An equivalent (and clearer) way to say that is: we can check whether the difference between the perspective taking and control groups was influenced by dispositional perspective taking (contrast 1) and whether the difference between the objective group and the other two groups combined was influenced by dispositional perspective taking (contrast 2).

 

A Better Custom Function for Contrasts in R

For all of it’s glory, R sucks when it comes to contrasts. In my mind, this is awful, because I’ve come to believe that nine times out of ten we should skip ANOVAs, which we don’t care about, and just go straight to the contrasts that do answer the questions we care about.

You can, of course, run contrasts in a regression model, but it would be nice to have a simple function that computes group contrasts without the need to create contrast codes for regression.

Here’s a custom function I wrote (with Rich Gonzalez):

# Runs a contrast based on same variances and separate variances

t.contrast <- function(dv, groups, contrast) {

  means <- by(dv, groups, mean)
vars <- by(dv, groups, var)
Ns <- by(dv, groups, length)
ihat <- contrast %*% means

#classic
df.classic <- sum(Ns)-length(Ns)
mse <- sum(vars*(Ns-1))/df.classic
se.classic <- sqrt(mse*(contrast^2 %*% (1/Ns)))
t.classic <- ihat/se.classic
p.classic <- 2*(1-pt(abs(t.classic), df.classic))

#welch
df.welch <- (contrast^2 %*% (vars/Ns))^2/(contrast^2 %*% (vars^2/(Ns^2*(Ns-1))))
se.welch <- sqrt(contrast^2 %*% (vars/Ns))
t.welch <- ihat/se.welch
p.welch <- 2*(1-pt(abs(t.welch), df.welch))

#compute confidence intervals
t.ci.classic <- qt(.025, df = df.classic)
t.ci.welch <- qt(.025, df = df.welch)
classic.lb <- ihat-t.ci.classic*se.classic
classic.ub <- ihat+t.ci.classic*se.classic
welch.lb <- ihat-t.ci.welch*se.welch
welch.ub <- ihat+t.ci.welch*se.welch
output <- data.frame(t = c(t.classic, t.welch),
df = c(df.classic, df.welch),
p = c(p.classic, p.welch),
lb.95CI = c(min(classic.lb, classic.ub), min(welch.lb, welch.ub)),
ub.95CI = c(max(classic.lb, classic.ub), max(welch.lb, welch.ub)))
rownames(output) <- c(“Equal variances”, “Separate variances”)
output <- round(output, digits = 3)
return(output)
}

 

The function computes contrasts for Student’s t test and Welch’s separate variances t test, with 95% confidence intervals. It takes three arguments: the dependent variable, the groups, and a vector of weights for the contrast that matches the order of the groups. There’s no data argument, so you need to specify the location of the data. It also doesn’t handle missing data by default, so you need to specify that you only want to use non-missing values when you input the data. But here’s an example of response times from a study I ran with 5 groups, where I compare group 4 to groups 1 and 2 together, ignoring groups 3 and 5. (And yes, I know that response times are never normally distributed, but these are the data that I had easily accessible at the time of writing)

# select data from data.frame jenn with non-missing response times

with(jenn[!is.na(jenn$rt), ], t.contrast(rt, condition, c(1, 1, 0, -2, 0)))

                       t       df     p  lb.95CI ub.95CI
Equal variances    0.159 37183.00 0.874 -100.645 118.405
Separate variances 0.098 34329.15 0.922 -168.070 185.830

Quickly Find Orthogonal Contrasts for Any Number of Conditions

I’m going to blaspheme as a social psychologist and say that I don’t like to use ANOVA to analyze my data. The omnibus F tests don’t answer the questions that I want to answer, which is whether the means of each pair of conditions differ from each other. Instead, I prefer to create sets of orthogonal contrasts that complete the full factorial design and then running regression. I can run a few different regression models to test the contrast for each pair of conditions. For example, if I have three conditions, A, B, and C, then I can test if groups A and B differ by using the following contrasts in a regression model.

Condition	Contrast 1	Contrast 2
A	-1	-1
B	1	-1
C	0	2

What I care about is the first contrast that compares groups A and B. I don’t care at all about the second contrast, but I just need it to complete the experimental design (otherwise my error term could increase, which could make it harder for me to find the effect I care about).

When using this strategy, you want to make sure that these contrasts are orthogonal—meaning that they are uncorrelated. To check whether any pair of contrasts are orthogonal, you can multiple the values for each group, and them sum those products. If they sum to zero, then the contrasts are orthogonal.

When I teach people how to use contrasts in regression, one of the things they always get stuck on is how to find the full set of orthogonal contrasts. The easiest way to do this is with Helmert contrasts. You can start out with a contrast that compares two groups and leaves out every other group (by giving them a weight of 0). Then, for every subsequent contrast, you give all the groups that were involved in the previous contrast the same weight, and then you compare them to one other group that was not involved in the previous contrast.

In the example above, the first contrast compares groups A and B, and it leaves out group C. The second contrast gives A and B the same weight, and then compares them to C. If we wanted to add a fourth group, D, we would need one more contrast, and we would give groups A, B, and C the same weight an then compare them to D, which would get a weight of 0 in all the other contrasts.

Condition	Contrast 1	Contrast 2	Contrast 3
A	-1	-1	-1
B	1	-1	-1
C	0	2	-1
D	0	0	3

All of these contrasts are orthogonal to each other. And we can keep on going to any number of groups, and all of the Helmert contrasts will be orthogonal.

Condition	Contrast 1	Contrast 2	Contrast 3	Contrast 4	Contrast 5	Contrast 6
A	-1	-1	-1	-1	-1	-1
B	1	-1	-1	-1	-1	-1
C	0	2	-1	-1	-1	-1
D	0	0	3	-1	-1	-1
E	0	0	0	4	-1	-1
F	0	0	0	0	5	-1
G	0	0	0	0	0	6

Again, the ONLY contrast we care about is the first one. The rest of them are just there to represent the full experimental design. To compare any two other groups, just substitute them for groups A and B and compute a new set of contrast weights. For example, here’s how we could compare groups D and F instead of groups A and B.

Condition	Contrast 1	Contrast 2	Contrast 3	Contrast 4	Contrast 5	Contrast 6
D	-1	-1	-1	-1	-1	-1
F	1	-1	-1	-1	-1	-1
C	0	2	-1	-1	-1	-1
A	0	0	3	-1	-1	-1
E	0	0	0	4	-1	-1
B	0	0	0	0	5	-1
G	0	0	0	0	0	6

Write APA-Style Manuscripts Directly in RStudio

As penance for my long absence, I’m going to tackle a big topic today. I’m going to show you how to write your APA-Style manuscript directly in RStudio using LaTeX and Sweave. LaTeX is a typesetting program that uses some syntax to automatically format your manuscript, unlike MS Word where you have to format everything manually. Sweave just runs LaTeX from R and allows you to use R commands directly in your manuscript. My goal isn’t to make you an expert on everything that LaTeX can do. Quite frankly, I don’t know it that well. Instead, my goal is to give you an accessible framework so that you can go write your next manuscript in RStudio.

Quick note: It’s been a while since I set all of this up for the first time on my computer, so I can’t remember if you need to download some LaTeX packages, like the apa6 and apacite packages. If something doesn’t work for you, then shoot me an e-mail and I’ll try to help you through the problem and update this tutorial.

Motivation

Why write your manuscript in RStudio? I’ll give you three reasons.

1. You can do everything in one program. Instead of switching back and forth between windows for SPSS and MS Word, or between R and MS Word, you can do everything in a couple of RStudio tabs.

2. Never forget to update your statistics again. Sometimes you collect more data, or you come up with a better strategy for data analysis, and then you need to update your manuscript. This usually requires running the new analyses in your stats program, and then writing the results in your MS Word document. Well, sometimes you run the analysis and then forget to update your Word document. Or you enter a typo and so your manuscript numbers differ from your stats program numbers. By writing everything in RStudio, you can write the analysis directly in your manuscript. That means if you tweak the analysis, your manuscript is automatically updated to reflect the changes.

3. Never fidget with APA-style formatting again. There are packages in LaTeX that will give you an APA-style document, including references, automatically. This lets you just focus on the writing.

So let’s get started.

Install LaTeX

Sweave is a part of R and RStudio, but you need to install LaTeX. You can find installation instructions at the LaTeX Project website.

Set up RStudio

Now that you’ve installed LaTeX, it’s time to make sure that RStudio is set up properly. Find the RStudio Preferences. Then go to the Sweave section. Where it says “Weave Rnw files using:”, make sure to select Sweave. Where it says “Typeset LaTeX into PDF using:”, make sure to select pdfLaTeX. Save the changes, and now you’re ready to write your manuscript in RStudio!

Create a New R Sweave File

In RStudio, go to File -> New File -> R Sweave to create a new R Sweave file. The file extension will be Rnw, so if I refer to your R Sweave file or your Rnw file, it means exactly the same thing.

You’ll see a new file that should have the following code:

\documentclass{article}

\begin{document}
\SweaveOpts{concordance=TRUE}

\end{document}


This is LaTeX code. In LaTeX syntax, all functions begin with a backslash, then they have the function name, and then the arguments go in curly braces like this: \functionName{arguments}. The \documentclass{} function specifies what kind of document you’re writing, which will determine how LaTeX formats the document. The section between the \documentclass{} function and the beginning of the document is called the “preamble”. You can add functions to the preamble that will change the format of the document, which we’ll do next.

Setting Up Your APA-Style Manuscript

Now we’re going to change the R Sweave default code to give you the basic code for an APA-Style manuscript. We’ll do this by accessing the apa6 and apacite packages in LaTeX.

Replace the code in your R Sweave file with the following code:

\documentclass[man,a4paper,noextraspace,apacite]{apa6}
\usepackage{apacite}
\title{}
\shorttitle{}
\author{}
\affiliation{}

\abstract{}
\keywords{}

\authornote{}

\begin{document}
\maketitle

\bibliography{}
\bibliographystyle{apacite}

\end{document}


Okay, there’s a lot going on here, so let’s break it down.

\documentclass[man,noextraspace,apacite]{apa6} This line specifies that we want to use the apa6 package to create an apa6 document. That’s why apa6 appears in the curly braces. All the stuff in the square brackets specifies some formatting for the apa6 document.

• man formats the document as a manuscript, whereas jou formats it like a published article (you can try switching back and forth after we compile our pdf later)
• noextraspace fixes a formatting problem that sometimes comes up with the apa6 documents, though I can’t actually recall the specific problem
• apacite specifies that we’ll use the apacite package to format our references section

\usepackage{apacite} This line accesses the apacite package, which will give us an APA-style References section.

\title{}, \shorttitle{}, \author{}, \affiliation{}, \abstract{}, \keywords{} These are exactly what they sound like. You can write the title of your paper, the running head (shorttitle), your name, your affiliation, your entire abstract, and your keywords (separated with commas) in the curly braces of each function. I’ll show you a complete example in a moment.

\authornote{} This is also what it sounds like, but I set it apart because you need to have double spaces between different parts of the author note to create line breaks. You’ll see this in the example below.

\maketitle This function has no arguments, it just takes all of the information in your preamble and turns it into a title page.

\bibliography{} This line specifies the name of the file with your references. So in the curly braces I would enter something like mybibliography.bib, or whatever I happened to call my references file. It should be stored in the same file folder as your Rnw file. It has to have the .bib extension because LaTeX uses BibTeX files to format the references. I’ll show you how to format these files using the apacite package later.

\bibliographystyle{apacite} This line tells LaTeX that you’re using the apacite package to format the references section.

Now that we’ve gone through all of the specific functions, let me show you what a filled-in document looks like.

\documentclass[man,a4paper,noextraspace,apacite]{apa6}
\usepackage{apacite}
\title{Criminals are Punished Less When They Harm Statistical Victims}
\shorttitle{Criminals and Statistical Victims}
\author{Joshua D. Wondra and Phoebe C. Ellsworth}
\affiliation{University of Michigan}

\abstract{There is a wealth of literature showing that people feel less compassion for statistical victims than for identified victims. But do people punish criminals less when they harm statistical victims than when they harm identified victims? We don't know if this study has been run already, but Josh thinks it's an interesting idea.}
\keywords{statistical victims, law and emotion, sentencing}

\authornote{Joshua D. Wondra, Department of Psychology, University of Michigan.

Phoebe C. Ellsworth, Department of Psychology, University of Michigan.

We are grateful to PRlab for their comments on an earlier version of this manuscript.

Correspondence concerning this article should be addressed to Josh Wondra, Department of Psychology, University of Michigan, 530 Church St., Ann Arbor, MI 48109-1043.

Contact: jdwondra@umich.edu}

\begin{document}
\maketitle

\bibliography{StatisticalVictimsReferences.bib}
\bibliographystyle{apacite}

\end{document}


Go ahead and fill in all of the arguments with information from a manuscript that you might write. Then look right above the syntax window in RStudio and click the Compile PDF button. If all goes well, then you should see the beginning of your APA-style document, which should look something like this.

By the way, your R Sweave/Rnw file needs to be given a name without spaces, because for some reason LaTeX hates spaces.

Writing the Manuscript, Except for the Data

Before we go into how you add your analyses to the document, here’s how you set up the general framework for an empirical paper.

\begin{document}
\maketitle

This is the first paragraph of my manuscript.

Create separate paragraphs with double spaces.

You can \textit{italicize} text with the textit function, if needed. For some special symbol, such as \%, you need to use a backslash before the symbol. You can even use Greek letters like \alpha, \beta, \Sigma, and \sigma.

\section{Method}

The section function creates level 1 headings.

\subsection{Procedure}
The subsection function creates level 2 headings. Any idea how to create level 3 headings?

\subsubsection{First Phase of the Study}
Subjects did tasks A, B, and C.

\subsubsection{Second Phase of the Study}
Subjects did tasks X, Y, and Z.

\section{Results}

And here are the results!

\section{Discussion}

Now, let's discuss the results...

\bibliography{StatisticalVictimsReferences.bib}
\bibliographystyle{apacite}

\end{document}

Try filling it out on your own and then you can Compile the pdf again. It should look something like this.

Adding the Analyses with Sweave

Now we’re ready to add data. You can insert chunks of R syntax using the following code:


<<ChunkTitle, echo=FALSE>>=
# R code goes here
@

The @ symbol ends the chunk. The echo=FALSE argument makes it so that every piece of R output doesn’t show up in your manuscript.

I suggest creating an initial chunk of R code at the beginning of your manuscript or at the beginning of your Method section to read in the data, set up any packages that you want to use, and do anything else that will apply throughout your manuscript. Then you can create other chunks that will run your analyses and insert the results into your manuscript. Here’s an example:

\section{Method}
<<InitializeData, echo=FALSE>>=
# I'm going to generate random data, but if you had a real data file you could read it in here

# Generate data
set.seed(1234) # I'm setting the seed to 1234 so that you can get the same results if you decide to follow along
punishment <- rnorm(40, mean=4, sd=1)
group <- factor(rep(c('identifiable victim','statistical victim'), each=20))
gender <- sample(c('female','male'), replace=TRUE, size=40)
myData <- data.frame(punishment, group, gender)

# Load packages
library(ggplot2)
@

\subsection{Overview}
Subjects read about a legal case in which the criminal stole money from one identifiable victim, or from many statistical victims. Then they indicated how many years they thought the criminal should spend in prison.

\subsection{Subjects}
<<SubjectDemographics, echo=FALSE>>=
# In this chunk, I'll generate subject information
totalN <- length(myData$gender)
femaleN <- sum(myData$gender=='female')
@
Subjects were \Sexpr{totalN} students (\Sexpr{femaleN} female) who participated for course credit.

See what I did? In the Sweave chunks I wrote my R code and created variables that I could insert into my manuscript using the \Sexpr{} function. The name of the function is short for “S expression”, because the programming language R is based on the programming language S. You can do the same thing for the results section.

\section{Results}
<<Punishment, echo=FALSE>>=
t.punishment <- t.test(punishment ~ group, data=myData)

## NOTE: Everything is rounded to 2-3 digits

# Save the mean and standard deviation for each group as variables
id.mean <- round(mean(myData$punishment[which(myData$group=='identifiable victim')]), digits=2)
id.sd <- round(sd(myData$punishment[which(myData$group=='identifiable victim')]), digits=2)
st.mean <- round(mean(myData$punishment[which(myData$group=='statistical victim')]), digits=2)
st.sd <- round(sd(myData$punishment[which(myData$group=='statistical victim')]), digits=2)

# Save stats from the t test as variables
t.value <- round(t.punishment$statistic, digits=2)
df <- round(t.punishment$parameter, digits=2)
p <- round(t.punishment$p.value, digits=3)
@

Subjects' punishment was no different when they read about identifiable victims (\textit{M} = \Sexpr{id.mean}, \textit{SD} = \Sexpr{id.sd}) than when they read about statistical victims (\textit{M} = \Sexpr{st.mean}, \textit{SD} = \Sexpr{st.sd}), \textit{t}(\Sexpr{df}) = \Sexpr{t.value}, \textit{p} = \Sexpr{p}.


Notice that I used the \textit{} function a lot to italicize the letters for the statistics. I think there’s a way to make this process more automatic when you report standard statistics, but I don’t know it offhand.

Try it out, compile the pdf, and you should get something like this.

Adding Figures

You can add a Sweave chunk to create a figure anywhere, surround it with a little LaTeX syntax, and you’ll be all set.

\begin{figure}
<<Fig1, echo=FALSE, fig=TRUE>>=
# Make sure to add the fig=TRUE part to the beginning of the chunk!
# We don't need to load ggplot2 because we did that in our initial chunk

# Set up the plot data
plot.means <- c(id.mean, st.mean)
plot.sds <- c(id.sd, st.sd)
plot.ses <- plot.sds/sqrt(20)
plot.groups <- factor(c('identifiable victim','statistical victim'))
plot.data <- data.frame(plot.means, plot.ses, plot.groups)

# Create barplot with standard error bars
ggplot(plot.data, aes(y=plot.means, x=plot.groups)) +
geom_bar(stat='identity') +
geom_errorbar(stat='identity', aes(ymin=plot.means-plot.ses, ymax=plot.means+plot.ses))
@

\textit{Figure 1.} Average punishment by group. Bars represent standard errors.
\end{figure}


This can go anywhere in the document, as long as its below all the code chunks that it refers to. So I couldn’t enter this before my initial chunk, because I load the ggplot2 package in the initial chunk. And I couldn’t enter it before the descriptive statistics, because I use the variables for the means and standard deviations that I created in the previous chunk.

A Note on Tables

Writing tables with LaTeX syntax is a bit more involved. I don’t use tables very often in my own manuscripts, so I won’t go into it here. I might revisit this later if someone asks, but for now I recommend that you search online for help formatting tables.

In-Text Citations and References

Now we’re on to the last big part of writing the APA-style manuscript, your references.

Start by creating a text file in TextEdit, NotePad, or my favorite, TextMate. Save it with the name you want to use for your References document, but make sure to add the extension .bib at the end. For example, I might save it as StatisticalVictimsReferences.bib.

Inside the .bib file, you’ll use the apacite syntax to create references. Here’s the basic format for a journal article.

@article{,
author = {},
title = {},
journal = {},
year = {},
volume = {},
pages = {x--x},
doi = {http://dx.doi.org/},
}


Aside from filling in the blanks in the curly braces, you need to create your own citekey that you’ll use for your in-text citations. Here’s an example.

@article{Small2007,
author = {Small, Deborah A. and Loewenstein, George and Slovic, Paul},
title = {Sympathy and callousness: The impact of deliberative thought on donations to identifiable and statistical victims},
journal = {Organizational Behavior and Human Decision Processes},
year = {2007},
volume = {102},
pages = {143--153},
doi = {http://dx.doi.org/10.1016/j.obhdp.2006.01.005},
}


My citekey is Small2007. And here is an example of the apacite entry for a book:

@book{Box1973,
author = {Box, G. E. P. and Tiao, G. C.},
title = {Bayesian inference in statistical analysis},
year = {1973},
publisher = {Addison-Wesley Publishing Company},
city = {Reading, MA},
}

You can find other reference formats for the apacite package in the package documentation.

Most of the errors I run into when writing documents are in my .bib file. I usually forget a comma, add too many commas, or forget the double-hyphen between page numbers.

To do in-text citations, you use the \cite{} function and feed it the citekey that you created.

There is a wealth of literature showing that people feel less compassion for statistical victims than for identified victims \cite{Small2007}.

And that’s it. LaTeX will create the References section and the in-text citation when you compile the pdf. You might need to compile twice for it to work, at least the first time. For other cases of in-text citations, see the apacite documentation to get the syntax.

The bad news: you need to type in all the information for each reference you use. The good news: once you create a .bib file, you can keep re-using it for every new manuscript you write, and just use the citekeys in future papers without any other typing.

Conclusion

So that’s it! You should now be able to write a manuscript in RStudio. Please send me an e-mail at jdwondra@umich.edu if questions come up, and I’m happy to help, and to revise this tutorial to make it better.

Assigning the Same Value to Multiple Objects in R

This was a simple, stupid problem to have, but I’m posting how to assign the same value to multiple objects in R because I couldn’t find a straightforward answer in my searches online.

The short answer is that if you do x <- y <- 5 then both x and y will be set equal to 5. But I’ll elaborate on the problem that I was trying to solve and my first failed attempt in case that information is useful to you.

Some Context

I have some data where subjects grouped emotions (sympathy, compassion, tenderness, empathy, and pity) together if they thought the emotion words meant the same thing. My dataset has variables designating which group each emotion was assigned to by a participant, but I wanted to create a matrix to display how many times the emotions were grouped together. So to see whether each participant grouped pity and sympathy together, I could do this:

pity_group==symp_group


For each participant who grouped them together, I would get TRUE as the output, and for each participant who grouped them separately, I would get FALSE. To get the total number of participants who grouped the two emotions together, I could do this:

sum(pity_group==symp_group)


Because R treats TRUE values equal to 1 and FALSE values equal to 0, the output is the number of TRUE values (i.e., the number of people who grouped the two emotions together). Next, I wanted to create a 5 x 5 matrix with the emotion names as the rows and the columns so that I could display all of the groupings together.

group_matrix <- matrix(rep(NA, 25), nrow=5) # this creates a 5 x 5 matrix with missing values in all 25 slots
colnames(group_matrix) <- c('pity','compassion','sympathy','empathy','tenderness')
rownames(group_matrix) <- c('pity','compassion','sympathy','empathy','tenderness')
diag(group_matrix) <- rep(196, 5) # this sets the values along the diagonal, where each emotion is grouped with itself, to 196, which is the number of subjects in the study


The matrix looks like this:

	pity	compassion	sympathy	empathy	tenderness
pity	196	NA	NA	NA	NA
compassion	NA	196	NA	NA	NA
sympathy	NA	NA	196	NA	NA
empathy	NA	NA	NA	196	NA
tenderness	NA	NA	NA	NA	196

What I needed to do is to make this a symmetric matrix so that the values that are reflected across the diagonal match. So, for example, the value in the first row, second column should match the value in the second row, first column.

First, I tried this:

group_matrix[2,1]; group_matrix[1,2] <- sum(pity_group==comp_group)


In my head, this syntax made sense because I was thinking about it as a single equation. On the left side, there were the two symmetric elements of the matrix. On the right side, there was the value that I wanted to assign to both of them. But I kept getting only one entry in the matrix, and one “NA” as output in the console in RStudio.

It turns out that this idea was silly. In R syntax, what semi-colons do is allow you to write multiple lines of syntax on a single line and execute them all at once. So I could write this…

2 + 2
1 + 5


… or I could write this…

2 + 2; 1 + 5

… and it would amount to the same thing.

What that means is that the code I was trying to use…

group_matrix[2,1]; group_matrix[1,2] <- sum(pity_group==comp_group)


… was first printing the value of group_matrix[2,1], which is why I got an “NA” in the console, and then it was executing this group_matrix[1,2] <- sum(pity_group==comp_group) as a separate line, which is why one of the elements in the grouping matrix was assigned the correct number.

What ultimately worked to correct the problem was to use the assignment operator <- twice…

group_matrix[2,1] <- group_matrix[1,2] <- sum(pity_group==comp_group)


… which assigns the grouping matrix element [2,1] the same value as the grouping matrix element [1,2], which is assigned the number of times that pity and compassion were grouped together. So my final matrix looked like this:

	pity	compassion	sympathy	empathy	tenderness
pity	196	5	31	13	3
compassion	5	196	39	34	48
sympathy	31	39	196	34	11
empathy	13	34	34	196	10
tenderness	3	48	11	10	196

Lesson 13: Using One-Way ANOVA to Test for Differences in the Means of More than Two Independent Groups

Overview

In this lesson, I introduce analysis of variance (ANOVA). We’re going to focus on analyzing the data when you have a single factor, which is called one-way ANOVA. This lesson picks up where we left off in Lesson 10.

The Problem

You predict that people enjoy watching movies more when they watch with other people than when they watch by themselves. To test this hypothesis, you have subjects watch a movie on a computer. Some subjects watch the movie in different rooms by themselves. Other subjects watch the movie in the same room as others who are watching the same movie. After watching, you ask how much they enjoyed the movie on a 7-point Likert-type scale (1 = not at all, 7 = very much).

However, with only these two conditions you can’t tell if differences in enjoyment are due to the shared experience or if they are due to the mere presence of other people. So you run a third condition where participants watch different movies in the same room. If subjects in the same room/same movie condition enjoy it more than those in the other two conditions, then you will know that the shared experience increases enjoyment. But if subjects in the same room/different movies condition enjoy the movie more than those in the different rooms condition, then you will know that the mere presence of others increases enjoyment.

Step 1: Simulate the Data

For a detailed example of simulating these kind of data, revisit Lesson 10. If you feel comfortable deciphering this kind of R syntax, then you can use the following code and proceed with the lesson.

# create a factor for the groups
condition <- factor(rep(c(1,2,3), each=70), levels=c(1,2,3), labels=c('same movie/same room','different rooms','different movie/same room'))
# create the dependent variables and restrict the range from 1 to 7
# set the seed so that you get the same simulated data that I post here
# in these simulated data, we make it so that the two same room conditions have the same mean enjoyment, which is higher than enjoyment in the different rooms condition
set.seed(83015)
enjoyment <- c(rnorm(n=70, mean=5, sd=1), rnorm(n=70, mean=4.6, sd=1), rnorm(n=70, mean=5, sd=1))
enjoyment <- round(enjoyment, digits=0)
enjoyment <- ifelse(enjoyment > 7, 7, ifelse(enjoyment < 1, 1, enjoyment))
movie.data <- data.frame(condition, enjoyment)

Step 2: Check Assumptions

ANOVA makes the same assumptions as the two sample t test—independence, normality, and equal variances. As noted in Lesson 10, you can do this using a boxplot. Let’s use the ggplot2 package for the boxplot.

library(ggplot2)
ggplot(movie.data, aes(y=enjoyment, x=condition)) + geom_boxplot()

# you can also use the qplot function (q for quick)
qplot(data=movie.data, y=enjoyment, x=condition, geom='boxplot')

These data look pretty ugly. In particular, the box for the different rooms condition has the median (dark line) on one edge and looks thinner than the other two groups. However, what we really care about is that the populations that the distributions came from are normal and have equal variances (we know that they are because we simulated the data).

Step 3: Run the ANOVA

As usual, we want to start out by getting descriptive statistics (and perhaps by plotting the data).

with(movie.data, by(enjoyment, condition, mean))

condition: same movie/same room
[1] 5.071429
---------------------------------------------------------------------- 
condition: different rooms
[1] 4.657143
---------------------------------------------------------------------- 
condition: different movie/same room
[1] 4.928571

with(movie.data, by(enjoyment, condition, sd))

condition: same movie/same room
[1] 0.9827708
---------------------------------------------------------------------- 
condition: different rooms
[1] 0.9763243
---------------------------------------------------------------------- 
condition: different movie/same room
[1] 1.094405

There are a couple of ways to run the one-way ANOVA. The most straightforward is to use the aov() function.

# save the ANOVA model as a new object
m.enjoyment <- aov(enjoyment ~ condition, data=movie.data)

# if you just print the model, you get the following output
m.enjoyment

Call:
   aov(formula = enjoyment ~ condition, data = movie.data)

Terms:
                condition Residuals
Sum of Squares     6.2000  215.0571
Deg. of Freedom         2       207

Residual standard error: 1.019276
Estimated effects may be unbalanced

# now try using the summary() function
summary(m.enjoyment)

             Df Sum Sq Mean Sq F value Pr(>F)  
condition     2    6.2   3.100   2.984 0.0528 .
Residuals   207  215.1   1.039                 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

There are a couple of other ways to run the one-way ANOVA that require you to use the lm() function first. Here they are:

# first, we save the model as an lm() object
m.enjoyment <- lm(enjoyment ~ condition, data=movie.data)

# now enter the model into the anova() function
anova(m.enjoyment)

Analysis of Variance Table

Response: enjoyment
           Df Sum Sq Mean Sq F value  Pr(>F)  
condition   2   6.20  3.1000  2.9839 0.05278 .
Residuals 207 215.06  1.0389                  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

# next, we'll enter the model into the Anova() function in the car package
# install the car package first if you don't have it
install.packages('car')
library(car)
Anova(m.enjoyment)

Anova Table (Type II tests)

Response: enjoyment
          Sum Sq  Df F value  Pr(>F)  
condition   6.20   2  2.9839 0.05278 .
Residuals 215.06 207                  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

With a one-way ANOVA, the three functions (aov(), anova(), and Anova()) all give the same results. But when we move on to factorial ANOVAs, there will be some differences. In the long run, you’re better off using the Anova() function from the car package.

So what can you conclude from these data? The ANOVA tests the null hypothesis that all of the group means are equal. So if the test had been significant, then you would have concluded that there was a difference in enjoyment among the groups. Where is that difference? You have no idea. ANOVA can’t answer that question without follow-up tests.

It’s very common to see researchers analyze their data by first running an ANOVA, and then running follow-up tests (contrasts) if the ANOVA is statistically significant. This is a good strategy if you’re running exploratory analyses and you don’t have any specific hypotheses. But it’s a bad strategy when you do have specific hypotheses. First off, if the ANOVA is significant, it can’t tell you where the mean differences are in the groups, so it can’t actually answer your research question. Second off, if the ANOVA is NOT significant, then that does NOT mean that there are no significant differences among your groups! We’ll come back to this point in the next lesson where we’ll demonstrate with the same data that we’re using here.

Step 4: Write it Up

Here is a sample of how I would write up the results of our one-way ANOVA.

Subjects were equally* entertained when they watched the same movie in the same room (M = 5.07, SD = .98), when they watched the same movie in different rooms (M = 4.66, SD = .98), and when they watched different movies in the same room (M = 4.93, SD = 1.09), F(2,207) = 2.98, p = .05.

* This is a stylistic choice to keep things simple, and “equally entertained” shouldn’t be taken too literally. As noted above, a non-significant ANOVA doesn’t really tell us that there are no differences among the group means. Alternatively, I could have written “There was no evidence that subjects felt entertained to a different extent…”

Lesson 12: Writing Functions to Make Analyses More Efficient

Overview

In this lesson I’ll show you how to write your own functions. We’ll start out simple and then build our way up to show some of the diverse ways that you can use functions in your own research. Additionally, writing your own functions should help you understand what’s going on when you use other functions.

The Basics

The basic format of any function is this:

myFunction <- function(arg1, arg2, ... argN) {
*steps in function*
return(output)
}


myFunction is the object that you’re using to save the new function so that you can call it in the future. arg1 through argN are the arguments, or inputs, of the function. The space inside the brackets is where you write the steps that your function should take. At the end of those steps you can use the return() function to determine what output you receive at the end of the function. To call the function, you would enter myFunction(arg1=value1, arg2=value2, ... argN=valueN).

Let’s demonstrate by writing a function to find the mean of a set of numbers.

myMean <- function(numbers, denominator) {
total <- sum(numbers) #we use the built-in sum() function to add our numbers together
meanValue <- total/denominator #the denominator is how many numbers we're trying to average
return(meanValue)
}


After you run the syntax to save the function myMean(), you can use it to find the average of a set of numbers, like this:

myMean(numbers=c(1,2,3,4), denominator=4)
[1] 2.5


By default, the arguments are processed in the order that you wrote them, so writing myMean(c(1,2,3,4), 4) will work just fine, but writing myMean(4, c(1,2,3,4)) will treat it as though you want 4 to be the value for the numbers argument and c(1,2,3,4) to be the value of the denominator (try running it that way and see what happens).

There are at least three reasons to write the names of each argument. First, it makes it so that the order of the arguments doesn’t matter, so writing myMean(denominator=4, numbers=c(1,2,3,4) would work just fine. Second, when you have the argument names with their values, it helps you to understand what you’re doing, especially when you use functions that have a lot of arguments. And third, it helps other researchers understand what you’re doing when you share your syntax. You should never assume that you’re the only person who will look at your data and analyses.

NOTE: When you write a new function you should make sure that you don’t call it the same name as something that already exists. For example, I used the name myMean instead of mean because R already has a function called mean. If I write a new function called mean then it won’t erase the existing function, but every time I call the mean function R will use the one that I created.

There are Multiple Routes to the Same End

One of the joys (and pains) of programming functions is that there are multiple ways to accomplish the same thing. For example, we could write our myMean() function in the following alternative ways.

myMean <- function(numbers) {
total <- sum(numbers)
meanValue <- total/length(numbers) #instead of specifying the value of the denominator in our arguments, we use the built-in function length() to find it for us
return(meanValue)
}


myMean <- function(numbers) {
meanValue <- sum(numbers)/length(numbers) #here we accomplish the summation and division in a single line
return(meanValue)
}


myMean <- function(numbers) {
return(sum(numbers)/length(numbers)) #here we don't specify any values in advance
}


The choice between the more complex coding and the simpler coding is a subjective one. My loose rules are to (1) try to automatize things as much as possible (so I would use the length() function for the denominator instead of writing the denominator as an argument), (2) try to reduce the number of lines of code as much as possible, (3) make sure that the lines of code are simple enough that when I look at the function in six months I’ll understand what’s going on. Rules 2 and 3 often come into conflict. For example, writing meanValue as an object and then returning it as an output probably makes it easier for me to understand what’s going on than writing everything in a single line as we did in the last example.

Sample Function: Effect Sizes and Building Up Functions

Let’s suppose that you want to compute a standardized effect size for the mean between two groups. The common measure that people use is Cohen’s d, which is the difference in group means divided by the pooled standard deviation. It’s okay if you’re not familiar with Cohen’s d, just try to follow along enough to understand the function. We can write the function as follows:

cohens.d <- function(m1, m2, var1, var2, n1, n2) {
mean.diff <- m1-m2

# here we break down the pooled standard deviation into more manageable parts by computing the numerator and the denominator separately
num.sd <- (n1-1)*var1 + (n2-1)*var2
denom.sd <- n1+n2-2
pooled.sd <- sqrt(num.sd/denom.sd)

d <- mean.diff/pooled.sd
return(d)
}

m1 and m2 are the group means, var1 and var2 are the group variances, and n1 and n2 are the group sample sizes. If we try it out by running cohens.d(m1=2, m2=3, var1=3.5, var2=4.2, n1=20, n2=20) then we get an effect size of d = -.51. This is a good, simple starting point, but it requires us to compute the means, variances, and sample sizes before can get what we want. We want our function to be more versatile (i.e., we want to be lazier in the long run) and so it’s time to build it up. The way we’ll do this is to compute the means, standard deviations, and sample sizes within the function itself so that all we have to do is input the data and the groups. Here is our new function:

cohens.d <- function(data, groups) {
# first we compute the means, variances, and sample sizes
# the by() function will apply the mean(), var(), and length() functions to our data for each group separately
means <- by(data, groups, mean)
vars <- by(data, groups, var)
ns <- by(data, groups, length)

# the square brackets indicate which element of an object you want
# in this case the first element is the first group's mean and the second element is the second group's mean
mean.diff <- means[1]-means[2]

num.sd <- (ns[1]-1)*vars[1] + (ns[2]-1)*vars[2]
denom.sd <- sum(ns)-2
pooled.sd <- sqrt(num.sd/denom.sd)

d <- mean.diff/pooled.sd
return(d)
}

Now let’s create some data to use with the new function. We need to create our dependent variable and our groups.

dv <- c(2,4,3,2,4,5,6,2,3,5,5,4,3,6,5,4,4,3,7,2)
groups <- rep(c(1,2), each=10)
cohens.d(data=dv,groups=groups)
-0.4786344


Great! So now we have a more flexible function. We can build up our function even further by adding a t test so that we get the inferential statistics along with the effect size. For our output, we’ll tell the function to return a list with both the outcome of the t test and the effect size.

cohens.d <- function(data, groups) {
means <- by(data, groups, mean)
vars <- by(data, groups, var)
ns <- by(data, groups, length)

mean.diff <- means[1]-means[2]

num.sd <- (ns[1]-1)*vars[1] + (ns[2]-1)*vars[2]
denom.sd <- sum(ns)-2
pooled.sd <- sqrt(num.sd/denom.sd)

d <- mean.diff/pooled.sd

# here we add the t test
# Cohen's d assumes that you have equal variances, so we set the var.equal argument to TRUE
t <- t.test(data ~ groups, var.equal=TRUE)

return(list(t=t, d=d))
}

Now try running the function with the same data that we used in the last example. You’ll see that you get both the results of the t test and the value of Cohen’s d in a single run of the same function.

Sample Function: Plotting Group Means in a Bar Graph with ggplot and Optional Arguments with Default Values

In the last lesson I showed you how to plot group means in a Bar Graph with ggplot. It required quite a few lines of code overall, so we’re going to simplify the process by putting all of it into a function. Let’s start with a function for two groups. As with our function for effect sizes, we want to make it so that we only have to enter the raw dependent data and the groups to get our plot.

barplots <- function(data, groups) {
# start by finding the means and standard errors and putting them in a data frame with group codes
means <- by(data,groups,mean)
ses <- by(data,groups,function(x){sd(x)/sqrt(length(x))})
# the unique() function will create one code per group
# we're also going to turn the group code into a factor because it will behave better with ggplot
group.codes <- factor(unique(groups))

# the means and ses that we created with the by() function are a special class of objects that won't cooperate when we try to put them in a data frame
# we can use the as.numeric() function to make those objects behave as though they are vectors of numbers
# this is called "coercing" one class of data to another class of data
plot.data <- data.frame(means=as.numeric(means), ses=as.numeric(ses), group.codes)

# create the plot
library(ggplot2) #makes sure ggplot is loaded
ggplot(plot.data, aes(x=group.codes, y=means)) +
geom_bar(stat='identity') +
geom_errorbar(stat='identity', aes(ymin=means-ses, ymax=means+ses))
}


Now if you run the function with the data from our effect size example, barplots(data=dv, groups=groups), then you will get the following plot.

Great! Now whenever you have groups with one factor you can run the function, enter the data in a single line, and voilà! You have your bar graph with standard error bars!

Now let’s build in the option to add a second factor to your plot. There are certainly a few ways to do it, but here’s an example that I figured out that works.

# we begin by changing the groups argument to factor1 and adding a factor2 argument with the default value NULL
barplots <- function(data, factor1, factor2=NULL) {

# the first two lines are the same
means <- by(data,factor1,mean)
ses <- by(data,factor1,function(x){sd(x)/sqrt(length(x))})
# but we're going to change our group codes by making a data frame with factor1 as a variable before we extract the unique group values
group.codes <- unique(data.frame(factor1=as.factor(factor1)))

# now we add an if() qualification to redo the data if there is a second factor
# the function will check if the statement "factor2 is NULL" is FALSE, in which case the code in curly brackets will run
if(is.null(factor2)==FALSE) {
# you can use the by() function with a list of multiple factors
means <- by(data,list(factor1,factor2),mean)
ses <- by(data,list(factor1,factor2),function(x){sd(x)/sqrt(length(x))})
# you can use the unique() function to find unique combinations of factors as long as they are combined in a data frame or in a matrix
group.codes <- unique(data.frame(factor1=as.factor(factor1),factor2=as.factor(factor2)))
}

# we create the data frame as normal
plot.data <- data.frame(means=as.numeric(means), ses=as.numeric(ses), group.codes)

library(ggplot2)

# save the plot into an object instead of just running it
plot.object <- ggplot(plot.data, aes(x=factor1, y=means)) +
geom_bar(stat='identity') +
geom_errorbar(stat='identity', aes(ymin=means-ses, ymax=means+ses))

# if there's a second factor, override the plot
if(is.null(factor2)==FALSE) {
# specify the second factor as the "fill" aesthetic
plot.object <- ggplot(plot.data, aes(x=factor1, fill=factor2, y=means)) +
# add the position='dodge' aesthetic to geom_bar and geom_errorbar (check out what happens when you leave this part out too)
geom_bar(stat='identity', position='dodge') +
geom_errorbar(stat='identity', position='dodge', aes(ymin=means-ses, ymax=means+ses))
}

# now print the plot
print(plot.object)
}


Give it a shot with the data we used in the previous function and it should produce the same graph because factor2 is NULL by default. Then, try this.

f1 <- rep(c(1,2), 10)
f2 <- rep(c(1,2), each=10)
barplot(dv,f1,f2)


It should give you the following plot.

Now let’s do one last thing and add optional arguments that will let you change the color of the bars. All you need to do is to change two lines. First, add one more argument to the function, so the first line will read barplots <- function(data, factor1, factor2=NULL, colors=c('blue','pink')) {. Second, change the last line of the function to print(plot.object + scale_fill_manual(values=colors)). Now if you make a two factor plot the color of the bars will be blue and pink by default, but you can change that with the colors argument.

You could continue building more arguments into this function until you find something that suits your general needs. And once it’s built, you can use it for all of your analyses and you never have to do it again.

A Quick Note on Debugging

As you build more complicating functions, you’re going to make mistakes (I did it a lot while writing the functions for this post). Debugging is going to deserve it’s own treatment with some examples, but my general recommendation is to try running each line of code in your function separately to see where things break down.

Lesson 11: Plotting Group Means in a Bar Graph

Overview

In this lesson, I’ll show you how to use the ggplot2 package to plot group means in a bar graph. We’ll also plot standard error bars.

Short Introduction to ggplot

The ggplot2 package allows you to do some complicated plotting, but we’ll start really simple. First, you’ll need to install the ggplot2 package and load it.

install.packages('ggplot2')
library(ggplot2)


The basic pieces of graphs with ggplot are the data, a geometric object (or “geom” for short) that defines what kind of plot you’ll get, and the aesthetics that customize the way the plot looks.

First, you need some data in a data frame. Let’s use the data from Lesson 10 where we had two groups—high power and low power—and our dependent variable was how sympathetic the subjects felt. Here is the R code to generate those data:

groups <- factor(rep(c('powerless','powerful'), each=60))

#generate dependent variables
set.seed(1111)
sympathetic <- c(rnorm(60,3,.5), rnorm(60,2.5,.5))
compassionate <- c(rnorm(60,2.6,.5), rnorm(60,2.3,.5))
empathetic <- c(rnorm(60,2.2,.5), rnorm(60,1.4,.5))

#restrict the maximum value to 4
sympathetic <- ifelse(sympathetic > 4, 4, sympathetic)
compassionate <- ifelse(compassionate > 4, 4, compassionate)
empathetic <- ifelse(empathetic > 4, 4, empathetic)

#round to the nearest integer
sympathetic <- round(sympathetic, digits=0)
compassionate <- round(compassionate, digits=0)
empathetic <- round(empathetic, digits=0)

#create an average of the three emotions
sympathy <- rowMeans(cbind(sympathetic, compassionate, empathetic))

#put the group and emotion variables in a data frame
power.data <- data.frame(groups, sympathetic, compassionate, empathetic, sympathy)


To begin plotting, we use the ggplot() function, which defines the data and some simple aesthetics.

ggplot(power.data, aes(y=sympathy, x=groups))


When you run the ggplot() function by itself, you’ll get an error and no plot. This is because we didn’t say what kind of geom we want. We can add the geom using the + character. Let’s start with a boxplot.

ggplot(power.data, aes(y=sympathy, x=groups)) + geom_boxplot()


For really basic plots like this example, we can use the qplot() function (q for “quick”) to do the same thing.

qplot(power.data, y=sympathy, x=groups, geom='boxplot)

You can make many different kinds of plots just by changing the y-value, x-value, and the geom. Check out this post from Noam Ross for a different introduction and some helpful references.

Plotting Group Means with ggplot

Plotting group means with ggplot takes a couple of extra steps. Let me show you what I mean by trying to plot a bar graph using the raw data.

#R will require you to add stat='identity' inside the geom_bar() function
ggplot(power.data, aes(y=sympathy, x=groups)) + geom_bar(stat='identity')


Those values are too high to be the group means, and they are. Without taking some extra steps, we can only plot the sum of all the values in the data, or we can remove the y-value and change stat=’identity’ to stat=’bin’ and plot the number of subjects in each condition. That’s not what we want. Instead, we’ll create a new data frame that has the group means and then plot them.

# find the group means
with(power.data, by(sympathy, groups, mean))
groups: powerful
[1] 2.077778
---------------------------------------------------------------
groups: powerless
[1] 2.622222

# save the group means and group labels into a new data frame
sympathy.means <- c(2.08, 2.62)
sympathy.groups <- factor(c('powerful','powerless'))
means.data <- data.frame(sympathy.means, sympathy.groups)

# plot the means using ggplot
ggplot(means.data, aes(y=sympathy.means, x=sympathy.groups)) + geom_bar(stat='identity')

Now if we want standard error bars, then we need to compute the standard error for each group. Recall that the standard error is equal to the standard deviation divided by the square root of the sample size. R doesn’t have a built-in function to compute the standard error, so we can write one.

standard.error <- function(x){sd(x)/sqrt(length(x))}


This standard.error() function takes one input, x, and then divides the standard deviation of x by the square root of its length. The length of an object is the number of elements in it, so that gives us the sample size. Let’s apply this to each group separately using the by() function.

with(power.data, by(sympathy, groups, standard.error))
groups: powerful
[1] 0.049076
---------------------------------------------------------------
groups: powerless
[1] 0.03998326


Now we can add the standard errors to our data frame with the means and plot our standard errors.

means.data$ses <- c(.05, .04)
# add the errorbar geom where the minimum is the mean - 1 standard error and the maximum is the mean + 1 standard error
ggplot(means.data, aes(y=sympathy, x=groups)) + geom_bar(stat='identity') + geom_errorbar(ymin=sympathy.means-ses, ymax=sympathy.means+ses)

Great! We have successfully plotted the means and standard errors! This plot looks ugly though. We need better labels, a bigger y-axis that represents the scale of our data and doesn’t cut off the top error bar, and maybe some nicer colors. Here is where we add layers to our plot to change the aesthetics. I’ll plot each layer on a new line so it’s easier to read.

ggplot(means.data, aes(y=sympathy.means, x=sympathy.groups)) +
geom_bar(stat='identity', fill='#02D4F1') + # the fill argument changes the bar color; can use hexadecimal codes for colors
geom_errorbar(ymin=sympathy.means-ses, ymax=sympathy.means+ses) +
scale_y_continuous(limits=c(0,4)) + # set the limits on the y-axis
labs(title='Figure 1. Average Sympathy by Power Condition', x='Condition', y='Average Sympathy') # adds plot title and axis labels


Lesson 10: Testing Differences in the Means of Two Independent Groups

Overview

In this lesson, I introduce the independent samples t test, which you use to test for a difference in the means of two independent samples. You will generate fake data to practice using the t test.

The Problem

You predict that people who feel powerless will feel more sympathy for someone else who is upset than people who feel powerful.

To test this prediction, you run an experiment with two conditions. In the powerless condition (n = 60), subjects write about a time when they were subordinate to someone else. In the powerful condition (n = 60), subjects write about a time when they had power over someone else.

After the power manipulation, subjects read a letter from someone else who says that he is upset because he recently lost his job. You ask subjects to report how much they feel sympathetic, compassionate, and empathetic on a five-point Likert-type scale (0=not at all, 4=extremely).

Your hypothesis is that subjects in the powerless condition will feel more sympathetic than those in the powerful condition on average. Your null hypothesis is that subjects in the two conditions will feel equally sympathetic on average.

Step 1: Generate the Fake Data

First, let’s create the grouping variables. We have 60 subjects per condition, so we need to create a variable that has 60 observations in the powerless condition and 60 observations in the powerful condition.

Remember that you can make different kinds of objects in R. When you’re coding groups, you have a choice to make: should you make the grouping variable a numeric object or a factor? There are some conditions where you will want your group variable to be numeric, but for this simple problem it’s probably easier to code it as a factor so that you can see the group labels. The good news is that when you’re using R it’s easy to convert a factor into a numeric object later if needed.

There are a couple of ways that you can create a factor, so just use whatever makes the most sense to you.

# Option 1: Create a string variable and then convert it into a factor
> groups <- c(rep('powerless',60), rep('powerful',60))
> groups <- factor(groups)


If you look at the help file for the rep() function, you’ll see that two arguments are rep(x, times), meaning that it will replicate whatever you put into the first argument (‘powerless’ or ‘powerful’) the number of times that you enter into the second argument (60 times). Alternatively, you can use rep(x, each) to replicate each element of the first argument 60 times, as in the following:

> groups <- rep(c('powerless','powerful'), each=60)
> groups <- factor(groups)

# Option 2: Create the factor directly
> groups <- factor(rep(c('powerless','powerful'), each=60))


Second, let’s create the dependent variables. We need to create responses to the three emotions: sympathetic, compassionate, empathetic. We’ll start by using the rnorm(n,mean,sd) function to randomly sample data from normal distributions. The first argument n sets the sample size. The second argument mean sets the mean of the population that we’ll sample from. The third argument sd sets the standard deviation of the population that we’ll sample from. Our first 60 observations will be for the powerless group and the second 60 observations will be for the powerful group.

Before we generate the dependent variable data, we’ll set the seed to a specific number. This number will affect the way that your computer generates random data. In particular, if you set the seed to the number that I specify here, you should get the same “random” results that I post here.

> set.seed(1111)

> sympathetic <- c(rnorm(60,3,.5), rnorm(60,2.5,.5)) # the powerless group has mean=3 and sd=.5; the powerful group has mean=2.5 and sd=.5
> compassionate <- c(rnorm(60,2.6,.5), rnorm(60,2.3,.5))
> empathetic <- c(rnorm(60,2.2,.5), rnorm(60,1.4,.5))


If you look at your results, you should notice a couple of problems. We said that subjects completed a 0-4 Likert-type scale, but some of the scores are greater than 4 and all of the observations are decimals rather than integers. We’ll have to fix this before we move on. Each time, we’ll make new variables for sympathetic, compassionate, and empathetic by applying some function to the existing variables.

First, make it so that if an observation is greater than 4, it is set equal to 4. We’ll use the ifelse(test, yes, no) function. The first argument is a testing condition that can be TRUE or FALSE. The second argument specifies what happens if the testing condition is TRUE. The third argument specifies what happens if the testing condition is FALSE.

> sympathetic <- ifelse(sympathetic > 4, 4, sympathetic) # if the value of an element of the variable sympathetic is greater than 4, then set it equal to 4; otherwise, set it equal to its original value
> compassionate <- ifelse(compassionate > 4, 4, compassionate)
> empathetic <- ifelse(empathetic > 4, 4, empathetic)


Next, we’ll use the round(x, digits) function to round the decimals to integers.

> sympathetic <- round(sympathetic, digits=0)
> compassionate <- round(compassionate, digits=0)
> empathetic <- round(empathetic, digits=0)


Finally, let’s combine the grouping variable and the dependent variables into a single data frame. We’ll also clean up our workspace by using the rm() function to get rid of the original variables.

> power.data <- data.frame(groups, sympathetic, compassionate, empathetic)
> rm(groups, sympathetic, compassionate, empathetic)


From now on if you want to access the variables, you’ll have to pull them from the power.data data frame.

> sympathetic # just entering the variable name returns "Error: object 'sympathetic' not found" because it is inside the power.data data frame
> power.data$sympathetic # this works because the dollar sign says "Look for sympathetic inside the object power.data"
> with(power.data, sympathetic) # this also works because the with() function says "Do everything inside the object specified by the first argument," which in this case is the data frame power.data


Now let’s say that you think the items sympathetic, compassionate, and empathetic are all measuring the same emotional reaction in slightly different ways. You’ll want to combine them by summing them or by finding the average of all three. It doesn’t really matter which, but to keep our overall sympathy measure on the same scale as the individual items (0-4), we’ll create a new variable in our data frame that has the mean of the three emotion items.

Your data are organized in the data frame as rows of participants and columns of variables. We want to combine the three emotion items across participants, so we want our new variable to show the average of the three items for each row (i.e., each participant). We can’t do this with the mean() function because it will give us the average of all three variables combined across all participants rather than the average for each participant separately. Instead, we’ll use the rowMeans() function.

> power.data$sympathy <- with(power.data, rowMeans(cbind(sympathetic, compassionate, empathetic))) # the cbind() function treats the variables as columns of data, which is required by the rowMeans() function; there is also an rbind() function that would treat them as rows of data


Step Two: Check Assumptions

The independent samples t test makes three assumptions

1. Independence. All observations are independent of each other, meaning that the probability of an observation taking on a specific value does not depend on the value of any other observation.
2. Normality. The data in each group are normally distributed.
3. Equal variances. The data from the two groups are equally variable. The complicated term that you’ll sometimes see for this is “homoscedasticity.”

In order for the results of the t test to be meaningful, we need to check these assumptions.

Checking independence in this case has more to do with your research design than a statistical check. How might your data violate the independence assumption? One clear example is if the same person participates in your study twice. If that person has a tendency to be more sympathetic in general, then they are likely to be extra sympathetic both times that they participate. Those two observations of the person’s sympathy are linked and independence is violated.

There are a few ways that you can check for normality and equal variances, but the easiest thing to do is to plot a boxplot of our composite sympathy score by group.

> boxplot(sympathy ~ group, data=power.data)




Recall that our normal distribution looks like this.



You can think of the boxplots as looking at the distribution of our data from above instead of from the side. Each box represents the middle 50% of the data for each group, which is called the interquartile range. The dark line is the median (the middle number for each group). The top and bottom whiskers reach out to the farthest points that are within 1.5 times the length of the interquartile range from the edges of the box. If you see any dots outside of the whiskers, then those are potential outliers (there are none in this case).

For the data in each group to be normally distributed, the boxplot for each group should look symmetric. This means that the medians (the dark lines) should be in the center of each box and the top and bottom whiskers within each group should be the same length. In our plot here the medians look okay, but the whiskers are a bit off. In the powerless group there is no top whisker because the box hits the maximum value in the data (which is 3). In the powerful group the top whisker is a bit longer than the bottom whisker, though the difference is not terrible. There might be a violation of the normality assumption in these data, but for the sake of this lesson we’ll pretend that things are okay.

For equal variances to hold, the interquartile ranges and whiskers should be the same length across the two groups. In our plot here, the interquartile ranges and the bottom whiskers look okay. However, the top whisker is longer in the powerful group than in the powerless group. There might be a violation of the equal variances assumption in these data. In practice, this will be okay. The classic t test assumes that variances are equal across groups, but there is also the Welch t test that does not assume equal variances. I’ve done some research showing that you tend to make better decisions if you just use the Welch t test at all times, so I’ll encourage you to use that t test exclusively.

Step Three: Run the t test

If our data are normally distributed and have equal variances across groups, then we can run our t test to compare the means of the two groups. First, we’ll look at the means and standard deviations of our groups using the by(x, group, function) function.

# Remember to use the with() function to access the variables that are inside the power.data data frame
> with(power.data, by(sympathy, groups, mean))
groups: powerless
[1] 2.622222
---------------------------------------------------
groups: powerful
[1] 2.077778

> with(power.data, by(sympathy, groups, sd))
groups: powerless
[1] 0.309709
---------------------------------------------------
groups: powerful
[1] 0.3801411


To run the classic t test in R, you use the t.test(y ~ x, var.equal=TRUE, data=myData) function.

> t.test(sympathy ~ groups, var.equal=TRUE, data=power.data)

Two Sample t-test

data: sympathy by groups
t = 8.6008, df = 118, p-value = 3.948e-14
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
0.4190897 0.6697992
sample estimates:
mean in group powerless mean in group powerful
2.622222 2.077778


In these data, the classic t test shows that there is a difference between the means of the two groups. The p value is less than the .05 cutoff and the 95% confidence interval does not contain 0 (more on confidence intervals in a future lesson). Now let’s run the Welch t test, which does not assume separate variances. All we need to do is change the value of the var.equal argument in to FALSE.

> t.test(sympathy ~ groups, var.equal=FALSE, data=power.data)

Welch Two Sample t-test

data: sympathy by groups
t = 8.6008, df = 113.37, p-value = 5.066e-14
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
0.4190366 0.6698522
sample estimates:
mean in group powerless mean in group powerful
2.622222 2.077778


What changed? The Welch t test adjusts the degrees of freedom to the extent that the variances of the two groups are unequal, and the degrees of freedom can have a decimal now. Additionally, the Welch t test uses a different standard error. As a consequence, the p value and confidence interval are slightly different, even though the t value did not change. In this case, the Welch t test also shows that there is a difference between the two groups. In practice, the two tests will almost always lead to the same decision when either the sample sizes or variances are equal, but when they lead to different decisions it is always safer to go with the Welch t test.

Step Four: Write it Up

Here is a sample of how I would write the results (you have my permission to use this format in your own writing if you would like). Using the Welch t test is not standard, even though it’s the better decision, so you need to specify that you used it in a paper.

To test whether there was a difference in sympathy between subjects in the powerless condition and those in the powerful condition, I ran a separate variances t test. Subjects in the powerless condition felt more sympathy (M = 2.62, SD = .31) than subjects in the powerful condition (M = 2.08, SD = .38, t(113.37) = 8.60, p < .001, 95% CI [.42, .67]).