Statistics

The most important thing about choosing statistics, is being able to justify the statistics that you choose.

There are many tools in the toolbox. There are many recommendations about what to do or not to do. It is a good idea to understand the statistics that you use, so that you can justify why you are using them for your analysis.

t-tests

one sample t-test

Set mu to test against a population value

x <- c(1,4,3,2,3,4,3,2,3,2,3,4)

#non-directional
t.test(x, mu = 2)

## 
##  One Sample t-test
## 
## data:  x
## t = 3.0794, df = 11, p-value = 0.01048
## alternative hypothesis: true mean is not equal to 2
## 95 percent confidence interval:
##  2.237714 3.428952
## sample estimates:
## mean of x 
##  2.833333

t.test(x, mu = 2, alternative = "two.sided") # the default

## 
##  One Sample t-test
## 
## data:  x
## t = 3.0794, df = 11, p-value = 0.01048
## alternative hypothesis: true mean is not equal to 2
## 95 percent confidence interval:
##  2.237714 3.428952
## sample estimates:
## mean of x 
##  2.833333

# directional
t.test(x, mu = 2, alternative = "greater")

## 
##  One Sample t-test
## 
## data:  x
## t = 3.0794, df = 11, p-value = 0.005241
## alternative hypothesis: true mean is greater than 2
## 95 percent confidence interval:
##  2.34734     Inf
## sample estimates:
## mean of x 
##  2.833333

t.test(x, mu = 2, alternative = "less")

## 
##  One Sample t-test
## 
## data:  x
## t = 3.0794, df = 11, p-value = 0.9948
## alternative hypothesis: true mean is less than 2
## 95 percent confidence interval:
##      -Inf 3.319326
## sample estimates:
## mean of x 
##  2.833333

paired-sample t-test

Set paired=TRUE

x <- c(1,4,3,2,3,4,3,2,3,2,3,4)
y <- c(3,2,5,4,3,2,3,2,1,2,2,3)

#non-directional
t.test(x, y, paired=TRUE)

## 
##  Paired t-test
## 
## data:  x and y
## t = 0.37796, df = 11, p-value = 0.7126
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.8038766  1.1372099
## sample estimates:
## mean of the differences 
##               0.1666667

t.test(x, y, paired=TRUE, alternative = "two.sided") # the default

## 
##  Paired t-test
## 
## data:  x and y
## t = 0.37796, df = 11, p-value = 0.7126
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.8038766  1.1372099
## sample estimates:
## mean of the differences 
##               0.1666667

# directional
t.test(x, y, paired=TRUE, alternative = "greater")

## 
##  Paired t-test
## 
## data:  x and y
## t = 0.37796, df = 11, p-value = 0.3563
## alternative hypothesis: true difference in means is greater than 0
## 95 percent confidence interval:
##  -0.6252441        Inf
## sample estimates:
## mean of the differences 
##               0.1666667

t.test(x, y, paired=TRUE, alternative = "less")

## 
##  Paired t-test
## 
## data:  x and y
## t = 0.37796, df = 11, p-value = 0.6437
## alternative hypothesis: true difference in means is less than 0
## 95 percent confidence interval:
##       -Inf 0.9585774
## sample estimates:
## mean of the differences 
##               0.1666667

Independent sample t-test

Note: Default is Welch test, set var.equal=TRUE for t-test

x <- c(1,4,3,2,3,4,3,2,3,2,3,4)
y <- c(3,2,5,4,3,2,3,2,1,2,2,3)

#non-directional
t.test(x, y, var.equal=TRUE)

## 
##  Two Sample t-test
## 
## data:  x and y
## t = 0.40519, df = 22, p-value = 0.6893
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.6863784  1.0197117
## sample estimates:
## mean of x mean of y 
##  2.833333  2.666667

t.test(x, y, var.equal=TRUE, alternative = "two.sided") # the default

## 
##  Two Sample t-test
## 
## data:  x and y
## t = 0.40519, df = 22, p-value = 0.6893
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.6863784  1.0197117
## sample estimates:
## mean of x mean of y 
##  2.833333  2.666667

# directional
t.test(x, y, var.equal=TRUE, alternative = "greater")

## 
##  Two Sample t-test
## 
## data:  x and y
## t = 0.40519, df = 22, p-value = 0.3446
## alternative hypothesis: true difference in means is greater than 0
## 95 percent confidence interval:
##  -0.5396454        Inf
## sample estimates:
## mean of x mean of y 
##  2.833333  2.666667

t.test(x, y, var.equal=TRUE, alternative = "less")

## 
##  Two Sample t-test
## 
## data:  x and y
## t = 0.40519, df = 22, p-value = 0.6554
## alternative hypothesis: true difference in means is less than 0
## 95 percent confidence interval:
##       -Inf 0.8729787
## sample estimates:
## mean of x mean of y 
##  2.833333  2.666667

printing

library(broom)

t_results <- t.test(x, y, var.equal=TRUE)
knitr::kable(tidy(t_results))

estimate1	estimate2	statistic	p.value	parameter	conf.low	conf.high	method	alternative
2.833333	2.666667	0.4051902	0.6892509	22	-0.6863784	1.019712	Two Sample t-test	two.sided

writing

The contents of R variables can be written into R Markdown documents.

t(22) =0.41, p = 0.689

# write your own function
report_t <- function(x){
  return(paste(c("t(",round(x$parameter, digits=2),") ",
        "= ",round(x$statistic, digits=2),", ",
        "p = ",round(x$p.value, digits=3)), collapse=""))
}
report_t(t_results)

## [1] "t(22) = 0.41, p = 0.689"

t_results

## 
##  Two Sample t-test
## 
## data:  x and y
## t = 0.40519, df = 22, p-value = 0.6893
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.6863784  1.0197117
## sample estimates:
## mean of x mean of y 
##  2.833333  2.666667

The results of my t-test were t(22) = 0.41, p = 0.689.

ANOVA

Requires data to be in long-format.

# example creation of 2x2 data 
Subjects <- rep(1:10,each=4)
Factor1 <- rep(rep(c("A","B"), each = 2), 10)
Factor2 <- rep(rep(c("1","2"), 2), 10)
DV <- rnorm(40,0,1)

all_data <- data.frame(Subjects = as.factor(Subjects),
                       DV,
                       Factor1,
                       Factor2)

between-subjects 1 factor

# run anova
aov_out <- aov(DV~Factor1, all_data)

# summary
summary(aov_out)

##             Df Sum Sq Mean Sq F value Pr(>F)
## Factor1      1   0.45  0.4467   0.386  0.538
## Residuals   38  44.01  1.1581

# kable xtable printing
library(xtable)
knitr::kable(xtable(summary(aov_out)))

	Df	Sum Sq	Mean Sq	F value	Pr(>F)
Factor1	1	0.44674	0.446740	0.385746	0.538252
Residuals	38	44.00854	1.158119	NA	NA

# print means
print(model.tables(aov_out,"means"), format="markdown")

## Tables of means
## Grand mean
##           
## 0.1031354 
## 
##  Factor1 
## Factor1
##        A        B 
## -0.00255  0.20882

between-subjects 2 factor

# run anova
aov_out <- aov(DV~Factor1*Factor2, all_data)

# summary
summary(aov_out)

##                 Df Sum Sq Mean Sq F value Pr(>F)
## Factor1          1   0.45  0.4467   0.383  0.540
## Factor2          1   0.13  0.1340   0.115  0.737
## Factor1:Factor2  1   1.86  1.8601   1.594  0.215
## Residuals       36  42.01  1.1671

# kable xtable printing
library(xtable)
knitr::kable(xtable(summary(aov_out)))

	Df	Sum Sq	Mean Sq	F value	Pr(>F)
Factor1	1	0.4467400	0.4467400	0.3827886	0.5400100
Factor2	1	0.1340183	0.1340183	0.1148335	0.7366757
Factor1:Factor2	1	1.8601091	1.8601091	1.5938323	0.2148966
Residuals	36	42.0144124	1.1670670	NA	NA

# print means
print(model.tables(aov_out,"means"), format="markdown")

## Tables of means
## Grand mean
##           
## 0.1031354 
## 
##  Factor1 
## Factor1
##        A        B 
## -0.00255  0.20882 
## 
##  Factor2 
## Factor2
##       1       2 
## 0.16102 0.04525 
## 
##  Factor1:Factor2 
##        Factor2
## Factor1 1       2      
##       A  0.2710 -0.2761
##       B  0.0511  0.3666

repeated-measures 1 factor

# run anova
aov_out <- aov(DV~Factor1 + Error(Subjects/Factor1), all_data)

# summary
summary(aov_out)

## 
## Error: Subjects
##           Df Sum Sq Mean Sq F value Pr(>F)
## Residuals  9  15.45   1.716               
## 
## Error: Subjects:Factor1
##           Df Sum Sq Mean Sq F value Pr(>F)
## Factor1    1  0.447  0.4467   0.242  0.635
## Residuals  9 16.640  1.8489               
## 
## Error: Within
##           Df Sum Sq Mean Sq F value Pr(>F)
## Residuals 20  11.92  0.5961

# kable xtable printing
library(xtable)
knitr::kable(xtable(summary(aov_out)))

	Df	Sum Sq	Mean Sq	F value	Pr(>F)
Residuals	9	15.44616	1.7162404	NA	NA
Factor1	1	0.44674	0.4467400	0.2416261	0.6347993
Residuals1	9	16.64000	1.8488893	NA	NA
Residuals2	20	11.92237	0.5961186	NA	NA

# print means
print(model.tables(aov_out,"means"), format="markdown")

## Tables of means
## Grand mean
##           
## 0.1031354 
## 
##  Factor1 
## Factor1
##        A        B 
## -0.00255  0.20882

repeated-measures 2 factor

# run anova
aov_out <- aov(DV~Factor1*Factor2 + Error(Subjects/(Factor1*Factor2)), all_data)

# summary
summary(aov_out)

## 
## Error: Subjects
##           Df Sum Sq Mean Sq F value Pr(>F)
## Residuals  9  15.45   1.716               
## 
## Error: Subjects:Factor1
##           Df Sum Sq Mean Sq F value Pr(>F)
## Factor1    1  0.447  0.4467   0.242  0.635
## Residuals  9 16.640  1.8489               
## 
## Error: Subjects:Factor2
##           Df Sum Sq Mean Sq F value Pr(>F)
## Factor2    1  0.134  0.1340   0.527  0.486
## Residuals  9  2.290  0.2544               
## 
## Error: Subjects:Factor1:Factor2
##                 Df Sum Sq Mean Sq F value Pr(>F)
## Factor1:Factor2  1  1.860  1.8601   2.192  0.173
## Residuals        9  7.639  0.8487

# kable xtable printing
library(xtable)
knitr::kable(xtable(summary(aov_out)))

	Df	Sum Sq	Mean Sq	F value	Pr(>F)
Residuals	9	15.4461637	1.7162404	NA	NA
Factor1	1	0.4467400	0.4467400	0.2416261	0.6347993
Residuals1	9	16.6400035	1.8488893	NA	NA
Factor2	1	0.1340183	0.1340183	0.5267850	0.4864057
Residuals2	9	2.2896723	0.2544080	NA	NA
Factor1:Factor2	1	1.8601091	1.8601091	2.1916373	0.1728958
Residuals	9	7.6385729	0.8487303	NA	NA

# print means
print(model.tables(aov_out,"means"), format="markdown")

## Tables of means
## Grand mean
##           
## 0.1031354 
## 
##  Factor1 
## Factor1
##        A        B 
## -0.00255  0.20882 
## 
##  Factor2 
## Factor2
##       1       2 
## 0.16102 0.04525 
## 
##  Factor1:Factor2 
##        Factor2
## Factor1 1       2      
##       A  0.2710 -0.2761
##       B  0.0511  0.3666

papaja

library(papaja)

apa_stuff <- apa_print.aov(aov_out)

The main effect for factor 1 was, \(F(1, 9) = 0.24\), \(\mathit{MSE} = 1.85\), \(p = .635\). The main effect for factor 2 was, \(F(1, 9) = 0.53\), \(\mathit{MSE} = 0.25\), \(p = .486\). The interaction was, \(F(1, 9) = 2.19\), \(\mathit{MSE} = 0.85\), \(p = .173\)

Linear Regression

lm(DV~Factor1, all_data)

## 
## Call:
## lm(formula = DV ~ Factor1, data = all_data)
## 
## Coefficients:
## (Intercept)     Factor1B  
##   -0.002546     0.211362

summary(lm(DV~Factor1, all_data))

## 
## Call:
## lm(formula = DV ~ Factor1, data = all_data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.56677 -0.52568 -0.01787  0.70304  2.19385 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.002546   0.240637  -0.011    0.992
## Factor1B     0.211362   0.340312   0.621    0.538
## 
## Residual standard error: 1.076 on 38 degrees of freedom
## Multiple R-squared:  0.01005,    Adjusted R-squared:  -0.016 
## F-statistic: 0.3857 on 1 and 38 DF,  p-value: 0.5383

lm(DV~Factor1+Factor2, all_data)

## 
## Call:
## lm(formula = DV ~ Factor1 + Factor2, data = all_data)
## 
## Coefficients:
## (Intercept)     Factor1B     Factor22  
##     0.05534      0.21136     -0.11577

summary(lm(DV~Factor1+Factor2, all_data))

## 
## Call:
## lm(formula = DV ~ Factor1 + Factor2, data = all_data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.50889 -0.53309 -0.00142  0.70919  2.17161 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept)  0.05534    0.29822   0.186    0.854
## Factor1B     0.21136    0.34435   0.614    0.543
## Factor22    -0.11577    0.34435  -0.336    0.739
## 
## Residual standard error: 1.089 on 37 degrees of freedom
## Multiple R-squared:  0.01306,    Adjusted R-squared:  -0.04028 
## F-statistic: 0.2449 on 2 and 37 DF,  p-value: 0.7841

lm(DV~Factor1*Factor2, all_data)

## 
## Call:
## lm(formula = DV ~ Factor1 * Factor2, data = all_data)
## 
## Coefficients:
##       (Intercept)           Factor1B           Factor22  
##            0.2710            -0.2199            -0.5471  
## Factor1B:Factor22  
##            0.8626

summary(lm(DV~Factor1*Factor2, all_data))

## 
## Call:
## lm(formula = DV ~ Factor1 * Factor2, data = all_data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.2932 -0.6123  0.1350  0.5750  1.9560 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)
## (Intercept)         0.2710     0.3416   0.793    0.433
## Factor1B           -0.2199     0.4831  -0.455    0.652
## Factor22           -0.5471     0.4831  -1.132    0.265
## Factor1B:Factor22   0.8626     0.6832   1.262    0.215
## 
## Residual standard error: 1.08 on 36 degrees of freedom
## Multiple R-squared:  0.05491,    Adjusted R-squared:  -0.02385 
## F-statistic: 0.6972 on 3 and 36 DF,  p-value: 0.5599

Randomization Test

A <- c(1,2,3,4,5,6,7,8,9,10)
B <- c(2,4,6,8,10,12,14,16,18,20)

all <- c(A,B)

mean_difference <- c()
for(i in 1:10000){
  shuffle <- sample(all)
  newA <- shuffle[1:10]
  newB <- shuffle[11:20]
  mean_difference[i] <- mean(newB)-mean(newA)
}

observed <- mean(B)-mean(A)
length(mean_difference[mean_difference >= observed])/10000

## [1] 0.0106

library(EnvStats)

twoSamplePermutationTestLocation(x=B, y=A, fcn = "mean", 
                                 alternative = "greater", 
                                 mu1.minus.mu2 = 0, 
                                 paired = FALSE, 
                                 exact = FALSE, 
                                 n.permutations = 10000, 
                                 seed = NULL, 
                                 tol = sqrt(.Machine$double.eps))

## 
## Results of Hypothesis Test
## --------------------------
## 
## Null Hypothesis:                 mu.x-mu.y = 0
## 
## Alternative Hypothesis:          True mu.x-mu.y is greater than 0
## 
## Test Name:                       Two-Sample Permutation Test
##                                  Based on Differences in Means
##                                  (Based on Sampling
##                                  Permutation Distribution
##                                  10000 Times)
## 
## Estimated Parameter(s):          mean of x = 11.0
##                                  mean of y =  5.5
## 
## Data:                            x = B
##                                  y = A
## 
## Sample Sizes:                    nx = 10
##                                  ny = 10
## 
## Test Statistic:                  mean.x - mean.y = 5.5
## 
## P-value:                         0.0122

Correlation

x <- rnorm(10,0,1)
y <- rnorm(10,0,1)
cor(x,y)

## [1] 0.01187039

ranks1 <- sample(1:10) 
ranks2 <- sample(1:10)
cor(ranks1,ranks2, method = "spearman")

## [1] -0.430303

Chi-square

M <- as.table(rbind(c(762, 327, 468), c(484, 239, 477)))
dimnames(M) <- list(gender = c("F", "M"),
                    party = c("Democrat","Independent", "Republican"))
(Xsq <- chisq.test(M))

## 
## Results of Hypothesis Test
## --------------------------
## 
## Alternative Hypothesis:          
## 
## Test Name:                       Pearson's Chi-squared test
## 
## Data:                            M
## 
## Test Statistic:                  X-squared = 30.07015
## 
## Test Statistic Parameter:        df = 2
## 
## P-value:                         2.953589e-07

binomial test

number_of_successes <- 100
trials <- 150
binom.test(x=number_of_successes, n=trials, p = 0.5,
           alternative = "two.sided",
           conf.level = 0.95)

## 
## Results of Hypothesis Test
## --------------------------
## 
## Null Hypothesis:                 probability of success = 0.5
## 
## Alternative Hypothesis:          True probability of success is not equal to 0.5
## 
## Test Name:                       Exact binomial test
## 
## Estimated Parameter(s):          probability of success = 0.6666667
## 
## Data:                            number_of_successes and trials
## 
## Test Statistic:                  number of successes = 100
## 
## Test Statistic Parameter:        number of trials = 150
## 
## P-value:                         5.447533e-05
## 
## 95% Confidence Interval:         LCL = 0.5851570
##                                  UCL = 0.7414436

Post-hoc tests

Many kinds of post-hoc tests can be done in R using base R functions, or functions from other packages.

At the same time, it is common for post-hoc comparison functions to be limited and specific in their functionality. So, if you can’t find an existing function to complete the analysis you want, you may have to modify, extend, or write your own function to do the job.

pairwise t-tests

Performs all possible t-tests between the levels of a grouping factor. Can implement different corrections for multiple comparisons.

df <- airquality #loads the airquality data from base R
df$Month <- as.factor(df$Month)

# default uses holm correction
pairwise.t.test(df$Ozone,df$Month)

## 
##  Pairwise comparisons using t tests with pooled SD 
## 
## data:  df$Ozone and df$Month 
## 
##   5       6       7       8      
## 6 1.00000 -       -       -      
## 7 0.00026 0.05113 -       -      
## 8 0.00019 0.04987 1.00000 -      
## 9 1.00000 1.00000 0.00488 0.00388
## 
## P value adjustment method: holm

# use bonferroni correction
pairwise.t.test(df$Ozone,df$Month, p.adj="bonf")

## 
##  Pairwise comparisons using t tests with pooled SD 
## 
## data:  df$Ozone and df$Month 
## 
##   5       6       7       8      
## 6 1.00000 -       -       -      
## 7 0.00029 0.10225 -       -      
## 8 0.00019 0.08312 1.00000 -      
## 9 1.00000 1.00000 0.00697 0.00485
## 
## P value adjustment method: bonferroni

tukey HSD test

notes: the warpbreaks data.frame should be available in base R, you do not need to load it to be able to use it.

Doesn’t work with repeated measures ANOVAs.

summary(fm1 <- aov(breaks ~ wool + tension, data = warpbreaks))

##             Df Sum Sq Mean Sq F value  Pr(>F)   
## wool         1    451   450.7   3.339 0.07361 . 
## tension      2   2034  1017.1   7.537 0.00138 **
## Residuals   50   6748   135.0                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

TukeyHSD(fm1, "tension", ordered = TRUE)

##   Tukey multiple comparisons of means
##     95% family-wise confidence level
##     factor levels have been ordered
## 
## Fit: aov(formula = breaks ~ wool + tension, data = warpbreaks)
## 
## $tension
##          diff        lwr      upr     p adj
## M-H  4.722222 -4.6311985 14.07564 0.4474210
## L-H 14.722222  5.3688015 24.07564 0.0011218
## L-M 10.000000  0.6465793 19.35342 0.0336262

Fishers LSD

Explanation of Fisher’s least significant difference test. https://www.utd.edu/~herve/abdi-LSD2010-pretty.pdf

library(agricolae)

data(sweetpotato)
model <- aov(yield~virus, data=sweetpotato)
out   <- LSD.test(model,"virus", p.adj="bonferroni")

Linear Contrasts

# Create data for a one-way BW subjects, 4 groups
A <- rnorm(n=10, mean=100, sd=25)
B <- rnorm(n=10, mean=120, sd=25)
C <- rnorm(n=10, mean=140, sd=25)
D <- rnorm(n=10, mean=80, sd=25)

DV <- c(A,B,C,D)
Conditions <- as.factor(rep(c("A","B","C","D"), each=10))
df <- data.frame(DV,Conditions)

# one-way ANOVA
aov_out <- aov(DV~Conditions, df)
summary(aov_out)

##             Df Sum Sq Mean Sq F value   Pr(>F)    
## Conditions   3  26220    8740    19.5 1.12e-07 ***
## Residuals   36  16137     448                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# look at the order of levels
levels(df$Conditions)

## [1] "A" "B" "C" "D"

# set up linear contrasts
c1 <- c(.5, -.5, -.5, .5) # AD vs. BC
c2 <- c(0, 0, 1, -1) # C vs. D
c3 <- c(0, -1, 1, 0) # B vs. C

# create a contrast matrix
mat <- cbind(c1,c2,c3)

# assign the contrasts to the group
contrasts(df$Conditions) <- mat

# run the ANOVA
aov_out <- aov(DV~Conditions, df)

# print the contrasts, add names for the contrasts
summary.aov(aov_out, split=list(Conditions=list("AD vs. BC"=1, "C vs. D" = 2, "B vs. C"=3)))

##                         Df Sum Sq Mean Sq F value   Pr(>F)    
## Conditions               3  26220    8740  19.498 1.12e-07 ***
##   Conditions: AD vs. BC  1  18961   18961  42.302 1.48e-07 ***
##   Conditions: C vs. D    1   5745    5745  12.816  0.00101 ** 
##   Conditions: B vs. C    1   1514    1514   3.377  0.07438 .  
## Residuals               36  16137     448                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1