Novum R-ganum: ANOVA and between-subjects contrasts in R and SPSS

For my few posts I'm trying to come up with procedures that are essentially rituals in Psychology. By ritual I mean psychologists do them all the time without really thinking about them; they're the most common Google searches that end with, "...in R?" (e.g., How do I run ANOVAs in R?). Between-subjects ANOVA seems like a good choice.

Warning: This post is long because it turns out even basic ANOVA is complicated. Because I created the html in R Markdown, I tend to just list SPSS syntax and print the output (I assume you're more familiar with those steps) and then I go into more detail about the procedures when I reveal the R code.

I'll use data from Klenman, Stern, and Trope (2016) who posted their data on the Open Science Framework here. Specifically, I'll use data from Study 3, which you can download here.

Relevant information about the data (from Klenman, Stern, and Trope, 2016):

Metaphorical-conflict manipulation: "...participants assigned to the compatible condition (n = 136) used their left hand (“Q” key) to categorize a photograph as Clinton and their right hand (“P” key) to categorize a photograph as Bush. Participants assigned to the incompatible condition (n = 136) had the opposite assignment of candidate to response hand. Participants assigned to the control condition (n = 138) used the “T” key to categorize a photograph as Clinton and the “V” key to categorize a photograph as Bush. We chose these two keys because they are vertically aligned, and so using them to categorize liberal and conservative politicians would not prime any metaphorical association between spatial location and political beliefs."

Perceptions of candidates’ beliefs: "To measure perceptions of the former presidents’ ideologies, we had participants separately indicate the social, economic, and general ideologies of Clinton and Bush using scales ranging from 1 to 9 (1 = extremely liberal, 5 = moderate, 9 = extremely conservative). We created a single ideology score for each former president by combining the three ratings (α = .88 and α = .89 for Clinton and Bush, respectively). We then calculated our dependent variable of the perceived difference between the former presidents’ beliefs by subtracting Clinton’s score from Bush’s score. Higher numbers indicate that Clinton’s and Bush’s attitudes were perceived as more different."

Here's the syntax for loading the data in SPSS.

GET
FILE='C:\Users\nickmm\Desktop\nmmichalak-analysis-examples-dfd10b4aaf3b\Kleiman et al. (Study '+
'3).sav'.
DATASET NAME DataSet2 WINDOW=FRONT.

DATASET ACTIVATE DataSet1.

USE ALL. FILTER BY filter_$.

EXECUTE.

Commands used:

GET
DATA SET
USE ALL
FILTER

Here's the code for loading the data in R1.

study3_filepath <- "~/Desktop/analysis-examples/Kleiman et al. (Study 3).csv"

study3 <- read.csv(file = study3_filepath, header = TRUE)

study3_subset <- study3[study3$filter_. == "Selected", c("Pnum","Condition","filter_.","Bush_Clinton_difference")]

study3_subset[sample(x = 1:nrow(study3_subset),
                     size = 10,
                     replace = FALSE), ]

##     Pnum Condition filter_. Bush_Clinton_difference
## 238  238   Control Selected                3.000000
## 430  430 Congruent Selected                6.000000
## 403  403 Congruent Selected                3.333333
## 233  233   Control Selected                3.333333
## 369  369  Conflict Selected                3.333333
## 87    87   Control Selected                5.000000
## 52    52   Control Selected                3.333333
## 450  450  Conflict Selected                3.000000
## 306  306 Congruent Selected                8.000000
## 149  149  Conflict Selected                3.333333

Functions used:

read.csv() takes a .csv file from a folder on your computer or, generally, it takes comma-separated data and turns it into a data frame (i.e., R's version of an Excel spreadsheet).
sample() takes a random sample from some data you give it.
nrow() tells you how many rows are in a matrix or data frame.
I also use basic subsetting code. You can find an easy introductory tutorial for subsetting at Quick-R, but, essentially, the argument is structured like [ "rows you want", "columns you want" ] (I added spaces for clarity, but they don't matter).

Here's the syntax for numerical and graphical tests of assumptions in SPSS.

EXAMINE VARIABLES=Bush_Clinton_difference BY Condition
/PLOT HISTOGRAM NPPLOT BOXPLOT SPREADLEVEL
/STATISTICS NONE.

Tests of NormalityTests of Normality, table, 2 levels of column headers and 1 levels of row headers, table with 7 columns and 5 rows
	Kolmogorov-Smirnov^a			Shapiro-Wilk
	Statistic	df	Sig.	Statistic	df	Sig.
Bush_Clinton_difference	.093	410	.000	.956	410	.000
a. Lilliefors Significance Correction

Plot a Spread vs. Level Plot to determine how to transform non-normal data.

Transform the data and plot the new variable.

COMPUTE Bush_Clinton_differenceT = abs(Bush_Clinton_difference)**1.026. if (Bush_Clinton_difference >= 0) Bush_Clinton_differenceT = Bush_Clinton_difference**1.026.

GRAPH
/HISTOGRAM = Bush_Clinton_differenceT.

Commands used:

EXAMINE
COMPUTE
IF
ABS
GRAPH

Here's how you get the same tests and plots in R.

First, check the normality assumption of Bush_Clinton_difference.

if(!("ggplot2" %in% rownames(installed.packages())))
install.packages("ggplot2")
library(ggplot2)

ggplot(data = study3_subset, aes(x = Bush_Clinton_difference)) + geom_histogram()

The difference between average perceptions of Bush's and Clinton's ideologies is, obviously, a difference score: t-tests (thus ANOVAs) assume difference scores are normally distributed. Technically, this variable violates the normality assumption: you can see that in the left-skew of the histogram and when you test this assumption with the Shapiro-Wilk Normality Test. Below I run that test, but note that a test like this has its own asssumptions ripe for violating. In addition, because power to detect an effect increases with sample size, the test will almost certainly produce a significant p-value for large N studies. Don't do this test. It's not helpful.

with(study3_subset,
  shapiro.test(Bush_Clinton_difference)
)

## 
##  Shapiro-Wilk normality test
## 
## data:  Bush_Clinton_difference
## W = 0.95603, p-value = 1.031e-09

What you can do is run what is called a Spread-Level Plot, which, for,

"linear models, plots log(abs(studentized residuals) vs. log(fitted values); fits a line to the plot; and calculates a spread-stabilizing transformation from the slope of the line."

That's a lot of jargon. But what you can get from all that is pretty cool because it turns out that the slope of this model is related to the power (i.e., "to the power of...") that will make your data normal: power transformation = 1 - slope.

if(!("car" %in% rownames(installed.packages())))
install.packages("car")
library(car)

trnsfrm <- with(study3_subset,
     spreadLevelPlot(Bush_Clinton_difference ~ Condition)
)

trnsfrm

##           LowerHinge    Median UpperHinge Hinge-Spread
## Conflict    8.666667  9.666667   11.00000     2.333333
## Congruent   9.000000 10.000000   11.00000     2.000000
## Control     9.000000 10.333333   11.33333     2.333333
## 
## Suggested power transformation:  1.025685

study3_subset$Bush_Clinton_differenceT <- with(study3_subset,
                                               ifelse(Bush_Clinton_difference < 0, abs(Bush_Clinton_difference)^trnsfrm$PowerTransformation,
                                                      Bush_Clinton_difference^trnsfrm$PowerTransformation))

ggplot(data = study3_subset, aes(x = Bush_Clinton_differenceT)) + geom_histogram()

Compare this histogram to the untransformed histogram. Pretty friggn' neat. Also, notice I made a new column in the dataset that has this transformed variable, Bush_Clinton_differenceT (T = transformed). I can use this to test whether the transformation makes a substantive difference in how we interpret the results (it does maybe a little²).

Functions used

install.packages() downloads the name of the package you give it (R has so many packages and people develop more every day) and puts it in your library
rownames() outputs the row names of the data frame you give it
%in% matches the names of elements in one vector with the names of elements in another
installed.packages() gives you a list of all installed packages
library() activates a package in your library. Think of it as taking a package of the shelf.
ggplot() is the first element of a ggplot2 plot. If you think of ggplot2 plots as sentences with a grammatical structure, then ggplot() is the subject...kinda.
geom_histogram() adds a histogram to your ggplot
with() lets you call variable names from a dataset without having to include [name of dataframe]$ every time you want to call a variable.
shapiro.test() conducts a Shapiro-Wilk Normality Test on data you give it.
spreadLevelPlot() produces a spread-level plot of the data you give it
! means "not equal to" in R
^ is the power operater in R (i.e., "to the power of...")
$ is the extract function in R. It lets you look inside R objects. If an R list "mylist" has two objects in it named "x" and "y", then to access object "y", you'd run mylist$y

Packages used:

car John Fox and colleagues wrote this package of functions to accompany their book Fox, J., & Weisberg, S. (2010). An R companion to applied regression. Thousand Oaks, CA: Sage. This package is a go to for so much that is regression.
ggplot2 is how you plot, and I don't mean just in R. It's how you plot. Period.

Next, check the homogeneity of variance assumption.

ggplot(data = study3_subset, aes(x = Condition, y = Bush_Clinton_difference, color = Condition)) + geom_boxplot() + scale_color_brewer(palette = "Set1")

There are a couple points that fall outside the "whiskers" of the boxplots, but there seem to be about the same number of these for each condition. It's fine. Also, don't run the Levene's Test for the same reason I suggested you don't run Shapiro-Wilk: you'd be testing assumptions with tests that have their own assumptions, and big sample sizes can basically guarantee significance.

Here's how you run between-subjects ANOVAs in SPSS (only one factor for now).

Between-subjects ANOVA using Type I Sum of Squares

GLM Bush_Clinton_difference BY Condition
/METHOD=SSTYPE(1).

Tests of Between-Subjects EffectsTests of Between-Subjects Effects, table, Dependent Variable, Bush_Clinton_difference, 1 layers, 1 levels of column headers and 1 levels of row headers, table with 6 columns and 10 rows

Bush_Clinton_difference

Source

Type I Sum of Squares

Mean Square

Sig.

Corrected Model

33.077^a

16.539

3.950

.020

Intercept

6445.834

1539.506

.000

Condition

33.077

16.539

3.950

.020

Error

1704.089

407

4.187

Total

8183.000

410

Corrected Total

1737.166

409

a. R Squared = .019 (Adjusted R Squared = .014)

Between-subjects ANOVA using Type II Sum of Squares

GLM Bush_Clinton_difference BY Condition
/METHOD=SSTYPE(2).

Bush_Clinton_difference

Source

Type II Sum of Squares

Mean Square

Sig.

Corrected Model

33.077^a

16.539

3.950

.020

Intercept

6445.834

1539.506

.000

Condition

33.077

16.539

3.950

.020

Error

1704.089

407

4.187

Total

8183.000

410

Corrected Total

1737.166

409

a. R Squared = .019 (Adjusted R Squared = .014)

Between-subjects ANOVA using Type II Sum of Squares

GLM Bush_Clinton_difference BY Condition
/METHOD=SSTYPE(3).

Bush_Clinton_difference

Source

Type III Sum of Squares

Mean Square

Sig.

Corrected Model

33.077^a

16.539

3.950

.020

Intercept

6441.393

1538.445

.000

Condition

33.077

16.539

3.950

.020

Error

1704.089

407

4.187

Total

8183.000

410

Corrected Total

1737.166

409

a. R Squared = .019 (Adjusted R Squared = .014)

Commands used:

Here's how you run the same ANOVAs in R.

Between-subjects ANOVA using Type I Sum of Squares

# Type I sum of squares
summary(
  aov(formula = Bush_Clinton_difference ~ Condition,
      data = study3_subset,
      contrasts = list(Condition = contr.sum(3)))
)

##              Df Sum Sq Mean Sq F value Pr(>F)  
## Condition     2   33.1  16.539    3.95   0.02 *
## Residuals   407 1704.1   4.187                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Below I introuduce the functions for running models with different types of sums of squares procedures. These models don't have factorial designs, so you don't have to worry about which type to use. I'm merely presenting these because 1) the car package is baller and 2) you should see that the results are the same across types of sums of squares (as they were in the SPSS output).

For more information regarding types of sums of squares, see Chapter 7 of Maxwell & Delaney (2004)³. Or Google it.

Between-subjects ANOVA using Type II Sum of Squares⁴

# Type II sum of squares
Anova(mod = lm(formula = Bush_Clinton_difference ~ Condition,
      data = study3_subset,
      contrasts = list(Condition = contr.sum(3))),
  type = "II"
)

## Anova Table (Type II tests)
## 
## Response: Bush_Clinton_difference
##            Sum Sq  Df F value Pr(>F)  
## Condition   33.08   2    3.95   0.02 *
## Residuals 1704.09 407                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Between-subjects ANOVA using Type III Sum of Squares

# Type III sum of squares
Anova(mod = lm(formula = Bush_Clinton_difference ~ Condition,
      data = study3_subset,
      contrasts = list(Condition = contr.sum(3))),
  type = "III" 
)

## Anova Table (Type III tests)
## 
## Response: Bush_Clinton_difference
##             Sum Sq  Df F value Pr(>F)    
## (Intercept) 6441.4   1 1538.45 <2e-16 ***
## Condition     33.1   2    3.95   0.02 *  
## Residuals   1704.1 407                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Functions used:

Anova() has more features than the base R statistics function aov() (e.g., Type II SS)
aov() is nice for quick and dirty ANOVAs. It comes with R and works well with base R functions like TukeyHSD.
contrasts() lets you get contrasts from factors or apply contrasts to factors. You can give a factor contrasts independent of a model or you can set contrasts as an argument in a model function like lm().
list() makes a list for you. Lists are nice because anything and everything can go in them (e.g., dataframes, other lists, character strings, etc.).
contr.sum() is a coding system that compares the mean of the dependent variable for a given level to the overall mean of the dependent variable. Give it the numer of levels in your factor.
geom_boxplot() adds a boxplot to your ggplot.
scale_color_brewer() allows you to set pre-made or custom color palletes for your ggplot.

Here's how you run contrasts for between-subjects factors in SPSS.

Custom contrasts and 95% confidence intervals

GLM Bush_Clinton_difference BY Condition
/METHOD=SSTYPE(3)
/CONTRAST(Condition) = SPECIAL(-1 1 0
0 1 -1
-1 0 1).

Contrast Results (K Matrix)Contrast Results (K Matrix), table, 2 levels of column headers and 3 levels of row headers, table with 4 columns and 24 rows
Condition Special Contrast			Dependent Variable
			Bush_Clinton_difference
L1	Contrast Estimate		.656
	Hypothesized Value		0
	Difference (Estimate - Hypothesized)		.656
	Std. Error		.247
	Sig.		.008
	95% Confidence Interval for Difference	Lower Bound	.170
		Upper Bound	1.142
L2	Contrast Estimate		.127
	Hypothesized Value		0
	Difference (Estimate - Hypothesized)		.127
	Std. Error		.247
	Sig.		.609
	95% Confidence Interval for Difference	Lower Bound	-.359
		Upper Bound	.613
L3	Contrast Estimate		.529
	Hypothesized Value		0
	Difference (Estimate - Hypothesized)		.529
	Std. Error		.248
	Sig.		.033
	95% Confidence Interval for Difference	Lower Bound	.042
		Upper Bound	1.017

Here's how you run the same contrasts for between-subjects factors in R.

Custom contrasts and 95% confidence intervals

I wanted to test a linear contrast in the paper (they don't call it that), so below I start off by rearranging the order of the factor levels so that the control condition is in the middle and the conflict and congruent conditions are on the ends (i.e., -1, 0, 1 makes sense when the control condition gets the 0 weight). When I test contrasts, R applies weights in the order they're assigned in the dataset.

if(!("multcomp" %in% rownames(installed.packages())))
install.packages("multcomp")
library(multcomp)

study3_subset$Condition <- factor(x = study3_subset$Condition, levels = levels(study3_subset$Condition)[c(1, 3, 2)])

levels(study3_subset$Condition)

## [1] "Conflict"  "Control"   "Congruent"

study3_multcomp <- rbind("CntrlVCnflct" = c(-1, 1, 0),
                          "CntrlVCngrnt" = c(0, 1, -1),
                          "Linear" = c(-1, 0, 1))

summary(
    glht(model = lm(formula = Bush_Clinton_difference ~ Condition,
      data = study3_subset),
      linfct = mcp(Condition = study3_multcomp),
      alternative = "two.sided"),
    test = adjusted("none"))

## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Multiple Comparisons of Means: User-defined Contrasts
## 
## 
## Fit: lm(formula = Bush_Clinton_difference ~ Condition, data = study3_subset)
## 
## Linear Hypotheses:
##                   Estimate Std. Error t value Pr(>|t|)   
## CntrlVCnflct == 0   0.6560     0.2472   2.653  0.00828 **
## CntrlVCngrnt == 0   0.1266     0.2472   0.512  0.60889   
## Linear == 0         0.5294     0.2481   2.134  0.03348 * 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- none method)

confint(glht(model = lm(formula = Bush_Clinton_difference ~ Condition,
      data = study3_subset),
      linfct = mcp(Condition = study3_multcomp),
      alternative = "two.sided"),
    test = adjusted("none")
    )

## 
##   Simultaneous Confidence Intervals
## 
## Multiple Comparisons of Means: User-defined Contrasts
## 
## 
## Fit: lm(formula = Bush_Clinton_difference ~ Condition, data = study3_subset)
## 
## Quantile = 2.3522
## 95% family-wise confidence level
##  
## 
## Linear Hypotheses:
##                   Estimate lwr      upr     
## CntrlVCnflct == 0  0.65601  0.07444  1.23758
## CntrlVCngrnt == 0  0.12660 -0.45497  0.70816
## Linear == 0        0.52941 -0.05427  1.11310

Above I tested contrasts in probably the most intuitive way: the null hypothesis is that the means of each condition are equal. You might also define the null hypothesis such that all the differences from the grand mean are equal to zero. Below I use the lm() function to demonstrate contrasts as deviations from the grand mean. In the beginning I print the contrast weights using contr.sum() (these are equivalent to SPSS's Deviation contrasts).

contr.sum(3)

##   [,1] [,2]
## 1    1    0
## 2    0    1
## 3   -1   -1

summary(
  lm(formula = Bush_Clinton_difference ~ Condition,
     data = study3_subset,
     contrasts = list(Condition = contr.sum(3)))
)

## 
## Call:
## lm(formula = Bush_Clinton_difference ~ Condition, data = study3_subset, 
##     contrasts = list(Condition = contr.sum(3)))
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -9.235 -1.098  0.098  1.109  4.431 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   3.9638     0.1011  39.223  < 2e-16 ***
## Condition1   -0.3951     0.1431  -2.761  0.00601 ** 
## Condition2    0.2609     0.1426   1.830  0.06801 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.046 on 407 degrees of freedom
## Multiple R-squared:  0.01904,    Adjusted R-squared:  0.01422 
## F-statistic:  3.95 on 2 and 407 DF,  p-value: 0.02

confint(
  lm(formula = Bush_Clinton_difference ~ Condition,
     data = study3_subset,
     contrasts = list(Condition = contr.sum(3)))
)

##                   2.5 %     97.5 %
## (Intercept)  3.76510890  4.1624273
## Condition1  -0.67642796 -0.1138534
## Condition2  -0.01939425  0.5411334

Here's how you run the same Deviation contrasts for between-subjects factors in SPSS (version 21 .htm output, sorry).

GLM Bush_Clinton_difference BY Condition

/METHOD=SSTYPE(3)

/CONTRAST(Condition) = Deviation.

Contrast Results (K Matrix)

Condition Deviation Contrast^a Dependent Variable

Bush_Clinton_difference

Level 1 vs. Mean Contrast Estimate -.395

Hypothesized Value 0

Difference (Estimate - Hypothesized) -.395

Std. Error .143

Sig. .006

95% Confidence Interval for Difference Lower Bound -.676

Upper Bound -.114

Level 2 vs. Mean Contrast Estimate .261

Hypothesized Value 0

Difference (Estimate - Hypothesized) .261

Std. Error .143

Sig. .068

95% Confidence Interval for Difference Lower Bound -.019

Upper Bound .541

a. Omitted category = 3

Some quick notes on running omnibus tests before specified contrasts (what I just did): why? If you designed an experiment having no clue which comparisons you were interested in, then, yea, you should adjust your p-values to account for your blind statistical analysis⁵. If, however, you designed an experiment with some questions in mind, you should just test those contrasts and dust your hands. ANOVAs for ANOVA's sake just hide statistical differences that may be present in your data⁶.

Functions used:

factor() turns data you give it into a factor (R's version of SPSS's nominal variable type).
c() combines objects you give it, a lot like Excel's CONCATENATE function but with more applications.
levels() can be used to set or print the levels of a factor object in R.
rbind() stands for row bind: it binds vectors together by rows.
glht() stands for general linear hypothesis. You can use this function for contrasts in model objects (e.g., lm(), aov(), etc).
confint() gives you confidence intervals for model paramaters. You can use this for CIs for, say, a regression model.

Packages used:

multcomp is a handy package for conducting contrasts.

Here's how you compute effect size measures in R (past this point, I don't know how to do it in SPSS).

Effect size: Eta-squared

I think many psychologists believe that measures of effect size are married to certain tests. For example, you might hear that you use Cohen's d for t-tests, eta-squared or partial eta-squared for F tests and odds ratios for chi-squared tests. In reality, there exist formulas to convert them all into whichever measure you want. I don't think it matters that much which you use, but because I'm demonstrating ANOVA, below I compute eta-squared.

if(!("DescTools" %in% rownames(installed.packages())))
install.packages("DescTools")
library(DescTools)

## 
## Attaching package: 'DescTools'

## The following object is masked from 'package:car':
## 
##     Recode

EtaSq(x = aov(formula = Bush_Clinton_difference ~ Condition,
      data = study3_subset,
      contrasts = list(Condition = contr.sum(3))),
      type = 2)

##               eta.sq eta.sq.part
## Condition 0.01904082  0.01904082

Functions used:

EtaSq() Calculates eta-squared, partial eta-squared and generalized eta-squared for ANOVAs

Packages used:

DescTools stands for descriptive tools because it's full of functions for descriptive and exploratory procedures

Confession: There's obviously a way to compute eta-squared in SPSS or I wouldn't have shown you the R code. When you run the GLM function, you specify you want eta-squared in the PRINT subcommand, like this:

GLM Bush_Clinton_difference BY Condition /METHOD=SSTYPE(3)

/PRINT=ETAQ.

I don't have the output for this, but please just take my word for it that the SPSS results for Partial Eta Squared = 0.019041 (which matches the R output above). SPSS also gives Partial Eta Squared for the intercept and the condition estimate (thanks, SPSS).

Effect size: Cohen's d

Because, "Cohen's d is for t-tests."

if(!("compute.es" %in% rownames(installed.packages())))
install.packages("compute.es")
library(compute.es)

table(study3_subset$Condition)

## 
##  Conflict   Control Congruent 
##       136       138       136

levels(study3_subset$Condition)

## [1] "Conflict"  "Control"   "Congruent"

study3_multcomp

##              [,1] [,2] [,3]
## CntrlVCnflct   -1    1    0
## CntrlVCngrnt    0    1   -1
## Linear         -1    0    1

tes(t = as.vector(
  summary(
    glht(model = lm(formula = Bush_Clinton_difference ~ Condition,
      data = study3_subset),
      linfct = mcp(Condition = study3_multcomp),
      alternative = "two.sided"),
    test = adjusted("none"))$test$tstat),
  n.1 = as.vector(
    table(study3_subset$Condition))[c(2, 2, 1)],
  n.2 = as.vector(
    table(study3_subset$Condition))[c(1, 3, 3)],
  verbose = TRUE)

## Mean Differences ES: 
##  
##  d [ 95 %CI] = 0.32 0.06 0.26 [ 0.08 -0.18 0.02 , 0.56 0.3 0.5 ] 
##   var(d) = 0.01 0.01 0.01 
##   p-value(d) = 0.01 0.61 0.03 
##   U3(d) = 62.57 52.47 60.21 % 
##   CLES(d) = 58.97 51.74 57.26 % 
##   Cliff's Delta = 0.18 0.03 0.15 
##  
##  g [ 95 %CI] = 0.32 0.06 0.26 [ 0.08 -0.18 0.02 , 0.56 0.3 0.5 ] 
##   var(g) = 0.01 0.01 0.01 
##   p-value(g) = 0.01 0.61 0.03 
##   U3(g) = 62.54 52.46 60.18 % 
##   CLES(g) = 58.94 51.74 57.24 % 
##  
##  Correlation ES: 
##  
##  r [ 95 %CI] = 0.16 0.03 0.13 [ 0.04 -0.09 0.01 , 0.27 0.15 0.24 ] 
##   var(r) = 0 0 0 
##   p-value(r) = 0.01 0.61 0.03 
##  
##  z [ 95 %CI] = 0.16 0.03 0.13 [ 0.04 -0.09 0.01 , 0.28 0.15 0.25 ] 
##   var(z) = 0 0 0 
##   p-value(z) = 0.01 0.61 0.03 
##  
##  Odds Ratio ES: 
##  
##  OR [ 95 %CI] = 1.79 1.12 1.6 [ 1.16 0.73 1.04 , 2.76 1.72 2.47 ] 
##   p-value(OR) = 0.01 0.61 0.03 
##  
##  Log OR [ 95 %CI] = 0.58 0.11 0.47 [ 0.15 -0.32 0.03 , 1.02 0.54 0.9 ] 
##   var(lOR) = 0.05 0.05 0.05 
##   p-value(Log OR) = 0.01 0.61 0.03 
##  
##  Other: 
##  
##  NNT = 9.88 56.28 12.5 
##  Total N = 274 274 272

Functions used:

table() is a handy function for creating contingency tables.
tes() takes t-scores and condition sample sizes and gives you all of the effect sizes you'd ever want, with 95% CIs
as.vector() forces something to be a vector, which in R is basically a line of things. I use it a lot to turn labeled numbers into just the numbers.

Packages used:

compute.es is one of the best packages for computing effect sizes I've ever found, in R and on the web. For example, you can give it p-values and sample sizes and it'll give you all the effect sizes with 95% CIs.

Entire courses are taught on ANOVA, so the best I thought I could do with this post is demonstrate all the procedures you'd need to run a simple one between-subjects factor ANOVA in R and get all the numbers you'd need for a publication. While I was writing this, I realized how many functions I used to get to that point. It might seem like a lot, but remembering all these functions isn't so different from learning the buttons SPSS. I promise you, the coding will pay off when you need to run more complex procedures or when you want to run some creative analysis or customize a plot. Notice, by the way, a small but obvious limitation we've found in SPSS already: computing effect sizes. I noticed other less obvious limitations while I was searching the web for SPSS syntax error messages: few accessible web pages discuss them. The R community clearly dwarfs the SPSS community.

But that's it for now. I hope to cover more complex ANOVA designs in future posts. Thanks to the researchers who posted their data and materials on OSF. Also, thank you for reading.

Happy R,

Nick

Footnotes

For the life of me, I could not get the foreign package's read.spss() or the Hmisc package's get.spss() to knit in R Markdown without an error even though I was able to run it on the R console. After about 30 minutes, I just wrote the .sav file into a .csv I that's the object I use throughout the post.
I ran the contrasts reported in the paper and only the difference in the transformed variable between congruent and conflict conditions is reduced by an appreciable extent: t (407) = 1.70, p = .09. Given I know close to nothing about this research area and there are other studies reported in their paper whose results are consistent with their hypotheses, I make no claim about this effect.
Maxwell, S. E., & Delaney, H. D. (2004). Designing experiments and analyzing data: A model comparison perspective. New York, NY: Psychology Press
This SS procedure is equivalent to SPSS's Type III SS. Yes, Type II is Type III. Up is down and black is white. There is no obvious consistency for naming SS types so I suggest learning how the methods differ rather than how they're named.
To be sure, you can test interesting questions without having asked them ahead of time. In a future post I'll demonstrate post-hoc comparisons like Tukey, Scheffe, and Bonferroni.
Omnibus F = average of orthogonal contrasts = average of pairwise comparisons.

Thursday, June 2, 2016

ANOVA and between-subjects contrasts in R and SPSS

I'll use data from Klenman, Stern, and Trope (2016) who posted their data on the Open Science Framework here. Specifically, I'll use data from Study 3, which you can download here.

Relevant information about the data (from Klenman, Stern, and Trope, 2016):

Here's the syntax for loading the data in SPSS.

GET FILE='C:\Users\nickmm\Desktop\nmmichalak-analysis-examples-dfd10b4aaf3b\Kleiman et al. (Study '+ '3).sav'.DATASET NAME DataSet2 WINDOW=FRONT.DATASET ACTIVATE DataSet1.

Commands used:

Here's the code for loading the data in R1.

Functions used:

Here's the syntax for numerical and graphical tests of assumptions in SPSS.

Plot a Spread vs. Level Plot to determine how to transform non-normal data.

Transform the data and plot the new variable.

Commands used:

Here's how you get the same tests and plots in R.

First, check the normality assumption of Bush_Clinton_difference.

What you can do is run what is called a Spread-Level Plot, which, for,

"linear models, plots log(abs(studentized residuals) vs. log(fitted values); fits a line to the plot; and calculates a spread-stabilizing transformation from the slope of the line."

That's a lot of jargon. But what you can get from all that is pretty cool because it turns out that the slope of this model is related to the power (i.e., "to the power of...") that will make your data normal: power transformation = 1 - slope.

Functions used

Packages used:

Next, check the homogeneity of variance assumption.

Here's how you run between-subjects ANOVAs in SPSS (only one factor for now).

Between-subjects ANOVA using Type I Sum of Squares

Between-subjects ANOVA using Type II Sum of Squares

GLM Bush_Clinton_difference BY Condition /METHOD=SSTYPE(2).

Between-subjects ANOVA using Type II Sum of Squares

Commands used:

Here's how you run the same ANOVAs in R.

Between-subjects ANOVA using Type I Sum of Squares

For more information regarding types of sums of squares, see Chapter 7 of Maxwell & Delaney (2004)3. Or Google it.

Between-subjects ANOVA using Type II Sum of Squares4

Between-subjects ANOVA using Type III Sum of Squares

Functions used:

Here's how you run contrasts for between-subjects factors in SPSS.

Custom contrasts and 95% confidence intervals

Here's how you run the same contrasts for between-subjects factors in R.

Custom contrasts and 95% confidence intervals

Here's how you run the same Deviation contrasts for between-subjects factors in SPSS (version 21 .htm output, sorry).

Functions used:

Packages used:

Here's how you compute effect size measures in R (past this point, I don't know how to do it in SPSS).

Effect size: Eta-squared

Functions used:

Packages used:

Confession: There's obviously a way to compute eta-squared in SPSS or I wouldn't have shown you the R code. When you run the GLM function, you specify you want eta-squared in the PRINT subcommand, like this:

I don't have the output for this, but please just take my word for it that the SPSS results for Partial Eta Squared = 0.019041 (which matches the R output above). SPSS also gives Partial Eta Squared for the intercept and the condition estimate (thanks, SPSS).

Effect size: Cohen's d

Because, "Cohen's d is for t-tests."

Functions used:

Packages used:

But that's it for now. I hope to cover more complex ANOVA designs in future posts. Thanks to the researchers who posted their data and materials on OSF. Also, thank you for reading.

Happy R,

Nick

Footnotes

No comments:

Post a Comment

GET
FILE='C:\Users\nickmm\Desktop\nmmichalak-analysis-examples-dfd10b4aaf3b\Kleiman et al. (Study '+
'3).sav'.
DATASET NAME DataSet2 WINDOW=FRONT.

DATASET ACTIVATE DataSet1.

GLM Bush_Clinton_difference BY Condition
/METHOD=SSTYPE(2).

For more information regarding types of sums of squares, see Chapter 7 of Maxwell & Delaney (2004)³. Or Google it.

Between-subjects ANOVA using Type II Sum of Squares⁴