ViewSlide("figures/Rcomic.png")
#### Install packages (if required) ####
install.packages("psych")
install.packages("fBasics")
#install.packages("lawstat")
install.packages("car")
install.packages("lattice")
install.packages("compute.es")
#### Writing/Compiling a Function ####
# A simple example
simplefunc <- function(inputxt) {
# Print the input text to the screen
print(inputxt)
}
simplefunc("R is FUN!")
simplefunc() # Since we did not provide a default value for inputxt, this throws an error
# Another simple example
# Here we added default values for the input (though they don't make much sense!)
adder <- function(x=1, y=2) {
x + y
}
adder(5,7)
adder(y=9)
adder(x=9)
out <- adder(2,3)
out
# A fancy example that does something useful
lsp <- function(package, all.names = FALSE, pattern)
{
package <- deparse(substitute(package))
ls(
pos = paste("package", package, sep = ":"),
all.names = all.names,
pattern = pattern
)
}
lsp(psych)
lsp(psych, TRUE)
lsp(psych, pattern = "^sim")
lsp(psych, pattern = "^win")
#### Utility Function for the Gentle R Course ####
source("GentleR.R") # Some functions I made up for viewing slides, not important to learn!
# If you get the following message, make sure you have set the R working
# directory to the directory containing the file "GentleR.R"
# "Error in file(filename, "r", encoding = encoding) :
# cannot open the connection"
#### Explore Data ####
# This example comes from the following website:
# http://lib.stat.cmu.edu/DASL/Stories/FusionTime.html
# Results from an experiment in visual perception using
# random dot sterograms.
# Group VV was given verbal and visual information about the embedded picture
# Group NV was either not given information or just verbal information
# Time to fuse the image was recorded and compared between groups
ViewSlide("figures/exdata.png")
fusion <- read.csv("data/fusion.csv")
head(fusion)
str(fusion)
fusion <- cbind(ID=1:nrow(fusion),fusion) # Adding row numbers as IDs
head(fusion)
str(fusion)
#*# Attach the data to make time and group visible #*#
attach(fusion)
#*# We should remember to detach it when we are done #*#
# Explore data separated by group #
library(psych)
# NV Group
describe(ftime[group=="NV"])
# var n mean sd median trimmed mad min max range skew kurtosis se
# 1 1 43 8.56 8.09 6.9 7.19 5.63 1.7 47.2 45.5 **2.68** **9.66** 1.23
# VV Group
describe(ftime[group=="VV"])
# var n mean sd median trimmed mad min max range skew kurtosis se
# 1 1 35 5.55 4.8 3.6 4.8 3.26 1 19.7 18.7 1.42 1.16 0.81
? boxplot
# Boxplots are useful tools for assessing symmetry of distributions
# and comparison of variance across groups
boxplot(ftime ~ group, data=fusion)
# Adding a label to the y-axis and a main title
boxplot(ftime ~ group, data=fusion,
main="Fusion Times for Random Dot Stereograms",
ylab="Time (s)")
mtext("No Visual (NV) and Visual/Verbal Cue (VV)") # adds a sub-title
# Identify lets you pick out points on the boxplot
identify(fusion$group, fusion$ftime, fusion$ID)
fusion[1,] # First value is a serious outlier
# First a frequency histogram with the default settings
hist(ftime[group=="NV"])
? hist
# We can improve the quality of this plot by changing some settings
hist(ftime[group=="NV"],
breaks=seq(0,50,1), # sets break points at 1, 2, 3, ..., 48, 49, 50
xlim=c(0,50), # Limits on the X axis
ylim=c(0,8), # Limits on the Y axis
main="Group NV", # Changes the main title
xlab="Time(s)") # Changes the x-label
hist(ftime[group=="VV"], breaks=seq(0,50,1), xlim=c(0,50), ylim=c(0,8), main="Group VV", xlab="Time(s)")
# The lattice() package can be used to make very nice graphs!
library(lattice)
histogram(~ ftime|group, data=fusion, type="count", breaks=seq(0,50,1),
main="Fusion Time by Group", xlab="Time(s)")
# These samples clearly deviate from normality!!
# We can create a density plot by setting the histogram 'freq' parameter to FALSE
# And overlay a normal distribution set with mean/sd of our sample to compare using 'curve'
hist(ftime[group=="NV"], freq=FALSE, breaks=seq(0,50,2), xlim=c(0,50))
curve(dnorm(x, mean=mean(ftime[group=="NV"]), sd=sd(ftime[group=="NV"])), add=TRUE, lty=2)
hist(ftime[group=="VV"], freq=FALSE, breaks=seq(0,50,2), xlim=c(0,50))
curve(dnorm(x, mean=mean(ftime[group=="VV"]), sd=sd(ftime[group=="VV"])), add=TRUE, lty=2)
## Univariate Statistical Tests of Normality ##
ViewSlide("figures/normality.png")
# Null Hypothesis of normality
# Thus rejection of null == non-normality
# (large sample == high power == rejection == conclude non-normal despite only small deviation from normality!)
# (small sample == low power == no rejection == conclude normal)
# Kolmogorov-Smirnov Test
# Best for large samples
ks.test(ftime[group=="NV"], "pnorm")
ks.test(ftime[group=="VV"], "pnorm")
# Shaprio-Wilk Test
# For small to medium sized samples (up to N=1000)
shapiro.test(ftime[group=="NV"])
shapiro.test(ftime[group=="VV"])
# D’Agostino’s test
# Recommended test, and only available in R!
# Small to large samples
library(fBasics)
dagoTest(ftime[group=="NV"])
dagoTest(ftime[group=="VV"])
## Univariate Tests of Equality of Variance ##
# Note: both methods for detecting variance heterogeneity are affected by nonnormality
## Levene Test ##
ViewSlide("figures/variance.png")
# Levene’s Test: Levene developed a test of variance homogeneity
# that tests the null hypothesis that the group variances are equal
# Therefore, a significant Levene test
# (i.e., p ≤α, where α is usually set at .10) indicates
# that the variances are not equal
library(car)
leveneTest(ftime, group)
# Test Statistic = 2.0185, p-value = 0.1595
# Non-significant suggests difference in variances non-significant
## Variance Ratios ##
# Look at the ratio of the largest to smallest *variance*
# Ratios larger than 2:1 (for unequal ns) or 4:1 (for equal ns)
# indicate variance heterogeneity
var(ftime[group=="NV"]) # 65.37388, n=43
var(ftime[group=="VV"]) # 23.0567, n=35
# Compute the ratio:
var(ftime[group=="NV"]) / var(ftime[group=="VV"]) # 2.835353
#### Student's t-test for 2 independent groups ####
## Assumptions:
# 1. The two samples are independently and randomly drawn from the source population(s)
# 2. The scale of measurement for both samples has the properties of an equal interval scale.
# 3. The source population(s) can be reasonably supposed to have a normal distribution.
# 4. The source population(s) can be reasonably supposed equal variances.
# From looking at our data through boxplots and histograms we likely have a
# problem with the 3rd and possibly 4th assumptions
# But we will conduct a t-test just to see how it works out
ViewSlide("figures/studentt.png")
# One option is to pass the two samples as vectors:
t.test(ftime[group=="VV"],ftime[group=="NV"], var.equal=TRUE)
# Another option is to use the 'formula' notation which is used by linear models
t.test(ftime ~ group, data=fusion, var.equal=TRUE) # Note: Will work without attaching the data frame, because we specified the parameter 'data'
# t = 1.9395, df = 76, p-value = 0.05615
t.test(ftime ~ group, data=fusion, var.equal=TRUE,
alternative="greater") # Note: added option for one-tailed test
## Effect Size ##
# There is ongoing controversy about the appropriateness of null-hypothesis testing
# Many recommend that null-hypothesis tests ought always be presented with some
# measure of effect size
# (effect size guidelines are arbitrary, and better to use values from substantive theory!)
ttestout <- t.test(ftime ~ group, data=fusion, var.equal=TRUE)
str(ttestout)
tstat <- as.numeric(ttestout$statistic) # getting out the t-statistic
df <- as.numeric(ttestout$parameter) # degrees of freedom
mean_diff <- abs(as.numeric(ttestout$estimate[1] - ttestout$estimate[2])) # absolute difference between the means
# Cohen's d
# Cohen's d guidelines:
# <0.2 is a small effect
# 0.2-0.5 is a medium effect
# 0.8+ a large effect
ViewSlide("figures/cohend.png")
se <- mean_diff / tstat
sp <- se / sqrt(1/length(ftime[group=="NV"]) + 1/length(ftime[group=="VV"]))
cohen_d <- mean_diff / sp
cohen_d # 0.4415395
# Eta-Squared (point-biserial correlation)
# Cohen's guidelines for eta^2:
# .01-.05 is a small association
# .06-.14 is a medium association
# .15 + is a large association
ViewSlide("figures/eta2.png")
eta2 <- (tstat)^2 / ((tstat)^2 + df)
eta2 # ~5% of the variability in fusion time attributed to visual/verbal cues (*if* data had actually met test assumptions)
sqrt(eta2) # r
# Omega-Squared
# Some argue eta^2 is biased, and recommend a modified version omega^2
ViewSlide("figures/omega2.png")
omega2 <- ((tstat)^2 - 1) / ((tstat)^2 + df + 1)
omega2 # ~3.4% of the variability in fusion time attributed to visual/verbal cues (*if* data had actually met test assumptions)
# Summary of lots of effect sizes
library(compute.es)
tes(tstat,length(ftime[group=="NV"]),length(ftime[group=="VV"]))
# r^2 == eta^2
#### Two-independent groups Welch t-test ####
# Assumptions:
# 1. The two samples are independently and randomly drawn from the source population(s)
# 2. The scale of measurement for both samples has the properties of an equal interval scale.
# 3. The source population(s) can be reasonably supposed to have a normal distribution.
# The Welch's t-test is an adaptation of Student's t-test
# intended for use with two samples having possibly unequal variances
# Unlike in Student's t-test, the denominator is not based on a pooled variance estimate.
ViewSlide("figures/welcht1.png")
# The degrees of freedom ν associated with this variance estimate is approximated
# using the Welch-Satterthwaite equation
ViewSlide("figures/welcht2.png")
# When variances are equal the Welch t (t’) is only slightly less powerful
# than the Student t
# When variances are unequal the Welch t maintains Type I error rates
# at the nominal level (i.e., α)
# Recommendation: Always use Welch's t!
# The Welch's t-test is the R default, we get this test if we remove the var.equal=TRUE
# statement, or set var.equal=FALSE
# Welch's t-test is known to be ...
t.test(ftime ~ group, data=fusion)
# t = 2.0384, df = 70.039, p-value = 0.04529
# (Notice the df is not a whole number, as it is computed using the Welch-Satterthwaite equ'n)
ViewSlide("figures/ttest_r_vs_spss.png")
#!@$* Exercise! *$@!#
#!@$* Compute the effect sizes for the Welch t-test?
#!@$* cohen_d = 0.4640548 standard deviation (*if* data had actually met test assumptions)
#!@$* eta2 = 5.6% of the variability in fusion time attributed to visual/verbal cues (*if* data had actually met test assumptions)
#!@$* omega2 = 4.2% of the variability in fusion time attributed to visual/verbal cues (*if* data had actually met test assumptions)
#### Non-Normality & Parametric Tests ####
ViewSlide("figures/nonnorm1.png")
ViewSlide("figures/nonnorm2.png")
ViewSlide("figures/nonnorm3.png")
#### Transformations ####
# Sometimes transformations can be used to improve normality
# The idea behind a transformation is to make the distributions
# more normal (or make the variances of the distributions more similar)
# At a statistical level, transformations can help improve the behavior
# of the test statistics
# It is okay to “fish around” for an appropriate transformation that will
# normalize your variables
# Interpreting the results after transformations becomes more complicated
# The transformed means (and mean differences) do not necessarily
# correspond to the original observations
# For +ve skew (skewed to the right):
# sqrt(x), ln(x), log10(x), -1*(1/x)
# For -ve skew (skewed to the right):
# x2, x3, or exp(x)
# or reflect the data [(Xhighest + 1) – X] then use a transform. for +ve skew
# The car() package has a function symbox for trying the 'ladder of powers'
# Ladder of Powers:
# -1 == inverse
# -1/2 == inverse square root
# 0 == log
# 1/2 == square root
# 1 == No Transformation
# 2 == Square
library(car)
symbox(~ ftime, data=fusion)
symbox(~ ftime, data=fusion[group=="NV",])
symbox(~ ftime, data=fusion[group=="VV",])
symbox(~ ftime, data=fusion[group=="VV",], powers = c(-2, -1, -0.5, 0, 0.5, 1, 2))
# Both overall and by group the 'log' transformation seems to
# do a good job making the fusion times more symmetric
#*# Since we are adding a variable to the data frame, we detach and attach again later #*#
detach(fusion)
fusion <- cbind(fusion, logft=log(fusion$ftime))
attach(fusion)
str(fusion)
# Original Data #
histogram(~ ftime|group, data=fusion, type="count", breaks=seq(0,50,1),
main="Fusion Time by Group", xlab="Time(s)", layout=c(1,2))
# Transformed Data #
histogram(~ logft|group, data=fusion, type="count", breaks=15,
main="Fusion Time by Group", xlab="Log Time", layout=c(1,2))
# Original Data #
boxplot(ftime ~ group, data=fusion,
main="Fusion Times for Random Dot Stereograms",
ylab="Time (s)")
# Transformed Data #
boxplot(logft ~ group, data=fusion,
main="Fusion Times for Random Dot Stereograms",
ylab="Log Time")
# Check homogeneity of variance for transformed time
leveneTest(logft, group) # p-value = 0.9894 <- non-significant
var(logft[group=="NV"]) # 0.6621488
var(logft[group=="VV"]) # 0.6688039 <- very good
var(logft[group=="NV"]) / var(logft[group=="VV"]) # 0.9900491 <- very good
# Welch's t-test on transformed data
t.test(logft ~ group, data=fusion)
# t = 2.3178, df = 72.673, p-value = 0.02328
# Based on this result we can much more comfortably assert that
# a statistically significant different is present between
# means of log(fusion time), where the group given visual and verbal information
# about the stimulus had on average smaller log(fusion time) than those who did not
#### Non-Parametric Mann-Whitney-Wilcoxon Test ####
# The nonparametric analogue to the two independent-samples t test
# Tests whether one of two samples of independent observations tends to have
# larger values than the other
# Except for differences in the treatment of ties,
# the Wilcoxon Rank Sum test is equivalent to the t-test
# performed on the ranks of the data
ViewSlide("figures/nonpar1.png")
ViewSlide("figures/nonpar2.png")
wilcox.test(ftime ~ group, data=fusion)
# W = 973, p-value = 0.02706
# Note that except for differences in the treatment of ties,
# the Wilcoxon Rank Sum test is equivalent to the t-test performed
# on the ranks of the data
rankftime <- rank(ftime)
t.test(rankftime[group=="VV"],rankftime[group=="NV"], var.equal=TRUE)
## Effect Size ##
# Glass’ rank biserial correlation coefficient #
ViewSlide("figures/glass.png")
rg <- 2*(mean(rankftime[group=="NV"]) - mean(rankftime[group=="VV"]))/length(rankftime)
rg
#### Independent Sample t-test for Trimmed Means ####
# “trim” or “remove” observations from the tails of the distribution
# Especially effective for reducing the effect of extreme observations and
# provide a clearer picture of the performance of a group as a whole
# Symmetric trimming, or the equal trimming of observations from both tails
# A trimmed mean is computed by simply removing [X]% of the observations from each
# tail of the sample, and averaging the remaining observations
# Winsorized variance, is computed by replacing the “trimmed” observations
# in the tails with the highest untrimmed value (upper tail) or the lowest
# untrimmed value (lower tail)
# The degrees of freedom are h1+h2-2, where h represents the “effective” sample size
ViewSlide("figures/trimt.png")
source("TrimFunctions.R")
yuen(ftime[group=="NV"], ftime[group=="VV"], tr=0.2)
# t-statistic = 2.092864
# df = 45.7124
# p-value = 0.0419407
## Effect Size ##
trimt <- yuen(ftime[group=="NV"], ftime[group=="VV"], tr=0.2)
trimt$teststat
eta2 <- (trimt$teststat)^2 / ((trimt$teststat)^2 + df)
eta2 # ~5.5% of the variability in fusion time attributed to visual/verbal cues (*if* data had actually met test assumptions)
#*# detach the data set again #*#
detach(fusion)
#*# sub-objects are no longer visible #*#
#### Writing to a File ####
# Saving the modified version of fusion for future use
write.table(fusion, file="data/newfusion", row.names=FALSE)
write.csv(fusion, file="data/newfusion.csv", row.names=FALSE)
#### Explore Data ####
# This example comes from the following website:
# http://lib.stat.cmu.edu/DASL/Stories/WomenintheLaborForce.html
ViewSlide("figures/exdata2.png")
lf <- read.csv("data/labforce.csv")
lf
hist(lf$Y1972)
boxplot(lf$Y1972)
describe(lf$Y1972)
hist(lf$Y1968)
boxplot(lf$Y1968)
describe(lf$Y1968)
#### Paired Sample t-test ####
# 1. The observations are independently and randomly drawn from the source population(s)
# 2. The scale of measurement for both samples has the properties of an equal interval scale.
# 3. The source population(s) of the differences can be reasonably supposed to have a normal distribution.
ViewSlide("figures/pairedt.png")
lf <- cbind(lf, diff = lf$Y1972 - lf$Y1968)
hist(lf$diff)
boxplot(lf$diff)
describe(lf$diff)
# var n mean sd median trimmed mad min max range skew kurtosis se
# 1 1 19 0.03 0.06 0.01 0.03 0.03 -0.1 0.15 0.25 0.26 0.02 0.01
t.test(lf$Y1972, lf$Y1968, paired=TRUE)
# t = 2.4577, df = 18, p-value = 0.02435
# mean of the differences 0.03368421
# Statistically significant and positive, suggesting increase in mean labour
# force participation rate
#### Wilcoxon signed ranks test ####
# If the distribution of difference scores is not normally distributed, the
# Wilcoxon signed ranks test can be much more powerful than the
# paired samples t-test
wilcox.test(lf$Y1972, lf$Y1968, paired=TRUE)
# V = 104, p-value = 0.01324
# Again indicates a significant difference
#### Explore Data ####
# The data are from an experimental study conducted by Baumann and Jones,
# as reported by Moore and McCabe (1993) Students were randomly
# assigned to one of three experimental groups in which different
# methods of teaching were used to teach reading (Basal, DTRA and Strat)
# Where Basal is the standard method, and Strat & DRTA are novel
# Researchers conducted two pretests and three post-tests for reading comprehension
# We will focus on the 3rd post-test
# This example comes from
# Fox, J., Weisberg, S. (2011). An R Companion to Applied Regression 2nd Edition. Sage Publications.
data(Baumann)
head(Baumann)
str(Baumann)
#*# Remember to detach Baumann later! #*#
attach(Baumann)
#*# Remember to detach Baumann later! #*#
# Create a contingency table (optionally a sparse matrix)
# from cross-classifying factors
xtabs(~ group, data=Baumann)
describe(post.test.3[group=="Basal"])
# var n mean sd median trimmed mad min max range skew kurtosis se
# 1 1 22 41.05 5.64 41 40.89 5.93 32 54 22 0.16 -0.55 1.2
describe(post.test.3[group=="DRTA"])
# var n mean sd median trimmed mad min max range skew kurtosis se
# 1 1 22 46.73 7.39 48.5 47.5 7.41 30 57 27 -0.78 -0.27 1.58
describe(post.test.3[group=="Strat"])
# var n mean sd median trimmed mad min max range skew kurtosis se
# 1 1 22 44.27 5.77 45 44.67 5.19 33 53 20 -0.62 -0.69 1.23
boxplot(post.test.3 ~ group, data=Baumann, xlab="Group", ylab="Reading Score")
# Levene's Test for Homogeneity of Variance
leveneTest(post.test.3 ~ group, data=Baumann)
# Df F value Pr(>F)
# group 2 0.3429 0.711 <- non-significant, do not reject null of homogeneity
# Ratio of largest variance to smallest variance
var(post.test.3[group=="DRTA"]) / var(post.test.3[group=="Basal"])
# = 1.718803, also pretty small, < 4 which is 'rule of thumb' for equal N's
# Stacked histogram of the data
histogram(~ post.test.3|group, data=Baumann, type="count", breaks=10,
main="Frequency of Reading Scores by Group", xlab="Reading Score @ Post-test 3", layout=c(1,3))
#### One-Way ANOVA ####
# Fit an ANOVA model to test the hypothesis of a difference amongst the group means
# Recall that an ANOVA is actually a linear regression model fit by least squares
# Hence, we use the lm() function in R to fit the model
# We should check to be sure that the grouping variable is of type “factor”
class(Baumann$group)
# Quick aside:
# If you have a data set where you called groups arbitrarily 1, 2 and 3, but they are in fact nominal
# R may well treat these as numeric when you import them
ex <- c(1, 1, 3, 2, 3, 3, 3, 2, 1, 1)
ex <- as.factor(ex) # You can convert them to a factor variable this way!
ex
# METHOD 1 #
baum.mod <- lm(post.test.3 ~ group, data=Baumann)
summary(baum.mod)
# F-statistic: 4.481 on 2 and 63 DF, p-value: 0.01515
anova(baum.mod)
# Df Sum Sq Mean Sq F value Pr(>F)
# group 2 357.3 178.652 4.4811 0.01515 *
# Residuals 63 2511.7 39.868
# METHOD 2 #
baum.mod.aov<-aov(post.test.3 ~ group, data=Baumann)
summary(baum.mod)
# METHOD 3 #
oneway.test(post.test.3 ~ group, data=Baumann, var.equal=T)
## Follow-up Pairwise Comparisons ##
TukeyHSD(baum.mod.aov) # Note: this only works for object from aov
# diff lwr upr p adj
# DRTA-Basal 5.681818 1.112137 10.251499 0.0111135 ** Honestly significant difference!
# Strat-Basal 3.227273 -1.342408 7.796953 0.2149995
# Strat-DRTA -2.454545 -7.024226 2.115135 0.4064363
# We can also use the pairwise.t.test to do all pairwise comparisons
pairwise.t.test(post.test.3, group, p.adj="none") # No adjustment for multiple comparisons
pairwise.t.test(post.test.3, group, p.adj="bonferroni") # Bonferroni adjustment (conservative)
pairwise.t.test(post.test.3, group, p.adj="holm") # Holm adjustment (a little less conservative), sequential procedure
pairwise.t.test(post.test.3, group, p.adj="BY") # Benjamini & Yekuteli -- false discovery rate
#### One-Way ANOVA under Variance Inequality ####
# Welch’s Independent Groups ANOVA
# As with the t-test, there is a Welch's adjusts the ANOVA weakening the assumption
# of homogeneity of variance
# As with the t-test, Welch’s test is actually the R default
# We can get this test by removing the var.equal parameter, or by setting it to FALSE
# METHOD 3 #
oneway.test(post.test.3 ~ group, data=Baumann)
# F = 4.3102, num df = 2.00, denom df = 41.49, p-value = 0.01992
## Follow-up Pairwise Comparisons ##
# Setting pool.sd to FALSE ensures use of Welch tests
pairwise.t.test(post.test.3, group, p.adj="none", pool.sd=FALSE)
pairwise.t.test(post.test.3, group, p.adj="bonferroni", pool.sd=F)
pairwise.t.test(post.test.3, group, p.adj="holm", pool.sd=FALSE)
pairwise.t.test(post.test.3, group, p.adj="fdr", pool.sd=FALSE)
#### One-Way ANOVA under Non-Normality ####
# The Kruskal-Wallis is a nonparametric version of ANOVA
# Useful if each ‘cell’ of the ANOVA table has the same distribution shape
kruskal.test(post.test.3 ~ group)
## Follow-up Pairwise Comparisons ##
# Unfortunately not as convienient
# We need to do the pairwise comparisons ourselves
# then we can use the function p.adjust to impose some
# kind of multiplicity control, such as Bonferroni
p1 <- wilcox.test(post.test.3[group=="Basal"],post.test.3[group=="DRTA"])$'p.value'
p2 <- wilcox.test(post.test.3[group=="Basal"],post.test.3[group=="Strat"])$'p.value'
p3 <- wilcox.test(post.test.3[group=="DRTA"],post.test.3[group=="Strat"])$'p.value'
pvals <- c(p1,p2,p3)
p.adjust(pvals, method="bonferroni")
#### One-Way ANOVA for Trimmed Means (Non-Normality & Variance Heterogeneity) ####
# Rand Wilcox also has a function (t1way) for computing a Welch omnibus test
# on trimmed means
# This test is much more reliable than a standard one-way ANOVA
# when the normality and variance homogeneity assumptions are violated
# The data are assumed to be stored in $x$ in list mode, so we need
# to do a little shuffling around
data <- list(post.test.3[group=="Basal"], post.test.3[group=="DRTA"], post.test.3[group=="Strat"])
t1way(data)
# test-statistic 5.747549, df1 = 2, df2 = 25.58675, p-value = 0.00867752
#### Factorial ANOVA ####
# Adding another factor variable, 'gender'
gender<- as.factor(c("F","M","F","M","F","F","F","F","M","M","F","F","F","F","M","F","F","M","M","M","F","M","M","F","M","M","M","M","M","F","M","M","M","M","M","F","F","M","M","M","M","F","M","M","F","F","F","M","F","F","F","M","F","F","F","F","M","M","F","M","F","F","M","F","F","F"))
# First we can add gender as a second factor
mod.fact <- lm(post.test.3 ~ group + gender, data=Baumann)
anova(mod.fact)
# Next we can test the model using group, gender and the interaction between the two
# We can specify the interaction explicitly using the :
mod.fact <- lm(post.test.3 ~ group + gender + group:gender, data=Baumann)
anova(mod.fact)
# Or use * and R will automatically expand to include all combinations of variables
mod.fact <- lm(post.test.3 ~ group * gender, data=Baumann)
anova(mod.fact)
# NOTE: We can skip a step
anova(lm(post.test.3 ~ group * gender, data=Baumann))
# An interaction plot can illustrate nicely where the differences lie
interaction.plot(group, gender, post.test.3)
## Types of Sums of Squares ##
# For balanced ANOVA (ie. same number of observations in each 'cell')
# the Type I, Type II and Type III methods are equivalent
# However, in it is frequently not the case, particularly in Psychology
# and other social sciences data!
# Note that the frequencies are not 'balanced'
# that is to say, there are different numbers of observations in the
# different levels of the variable.
xtabs(~group + gender)
# When the data is 'unbalanced' the results of the three Types of SS will differ
# (Note: this assumes each cell has at least one observation! If this is
# not the case still other methods are required, ie. SAS Type IV)
# Type I SS
# Sequential sum of squares
# Control only for variables already included in the model
ViewSlide("figures/typeISS.png")
# TYPE I SS #
# The anova function returns Type I (or 'sequential') sums of squares
# This means that the order in which the variables are entered into
# the model will effect the results of the resulting statistical test
anova(lm(post.test.3 ~ gender * group, data=Baumann))
# TYPE II SS #
# Partial sums of squares, not controling for higher order terms
# In other words, not controling for the interaction, but controling
# each main effect for the other
# Type II is definitely powerful if there is NO interaction in the population
# If there is an interaction in the population, but it goes undetected,
# sometimes Type III is more powerful, but main effects are tricky in that situation anyway
# Important to have enough power to detect interactions if they are substantial
ViewSlide("figures/typeIISS.png")
# To get Type II tests we can use the Anova() function from the package car()
# This is a DIFFERENT function! Note the Capital A!
Anova(lm(post.test.3 ~ group * gender, data=Baumann))
Anova(lm(post.test.3 ~ gender * group, data=Baumann))
# Order of entry does not affect the statistical test anymore
# TYPE III SS #
# Partial sums of squares, controling for the interaction and controling
# each main effect for the other
# Type III might be more powerful if there is an interaction in the population
ViewSlide("figures/typeIiISS.png")
# Type III SS are very unpopular among R circles, and aren't so easy to acquire
# You need to makes sure that the factor contrasts are set to the right
# type --- sum to zero coding --- or the analysis will not work properly!
# R defaults to dummy coding in general
contrasts(gender)
contrasts(group)
options(contrasts = c("contr.sum","contr.poly")) # Changing global settings
contrasts(gender) # These should have 1's and -1's
contrasts(group) # These should have 1's and -1's
Anova(lm(post.test.3 ~ group * gender, data=Baumann), type="III")
## Significant Interactions ##
# The SS allocated to the highest order term (the interaction in this case)
# is the same as for all three methods Type I, II & III.
# So if you have a significant interaction, it will be the case regardless of the
# method used
# You should then follow up with appropriate interaction contrasts
#*# Back to dummy coding #*#
options(contrasts = c("contr.treatment","contr.poly")) # Changing global settings
#*# detach the data set again #*#
detach(Baumann)
#*# sub-objects are no longer visible #*#
#### Explore Data ####
data(Davis) # Davis is a data set in the package car(), the data function makes it visible in the workspace
head(Davis)
str(Davis)
attach(Davis)
# 200 subjects engaged in regular exercise
# sex --- Female (F) and Male (M)
# weight --- weight in kg
# height --- height in cm
# repwt --- reported weight in kg
# repht --- reported height in cm
# the pairs() function creates a scatterplot matrix
pairs(~sex+weight+height+repwt+repht, data=Davis) # There appears to be a strange point in weight by height
# Looking more closely at Weight on Height
xyplot(weight ~ height, data=Davis)
trellis.focus("panel", 1, 1) # xyplot is a part of the lattice() package
panel.identify() # Identifing the outlier, seems to be a problem with observation #12
trellis.unfocus()
# Another version of scatter-plot
scatterplot(weight ~ height, data=Davis, smooth=FALSE, id.n=2)
# Option id.n=2 will automatically identify the 2 most extreme points
#### Correlation Test ####
# Test for association between paired samples, using one of Pearson's
# product moment correlation coefficient (default), Kendall's tau or Spearman's rho.
cor.test(weight, height)
# t = 2.7179, df = 198, p-value = 0.007152
# Test for monotonic relationship
cor.test(weight, height, method="spearman")
# S = 308751.6, p-value < 2.2e-16
#### Simple Linear Regression ####
davis.mod <- lm(weight ~ height, data=Davis) # Returns a model object
summary(davis.mod) # the summary function returns standard output
# Note: R-squared = 0.03597 <<<--- Not very good!
# The stray point #12 turns out to be a typo
Davis[12,]
# Note weight is 166kg, but reported weight is 56kg!
# We would prefer to fit our model without it.
xyplot(weight ~ height, data=Davis, subset=-12)
scatterplot(weight ~ height, data=Davis, smooth=FALSE, subset=-12)
davis.mod <- lm(weight ~ height, data=Davis, subset=-12)
summary(davis.mod)
# Multiple R-squared: 0.594 <<--- Much better fit!!
#### Multiple Linear Regression ####
# Extending to multiple linear regression is really quite easy
# simply add more variables to the formula expression!!
# Additive Model:
# mod <- lm( Y ~ X1 + X2 + X3, data = dataset)
# Multiple linear regression (with interactions):
# mod <- lm( Y ~ X1 + X2 + X1:X2, data = dataset)
# or equivalently
# mod <- lm( Y ~ X1*X2, data = dataset)
davis.mod2 <- lm(weight ~ height + repwt, data=Davis, subset=-12)
summary(davis.mod2)
# Interpretation
# repwt is a highly significant predictor
# as such height does not contribute to the prediction once repwt has been taken into account
# If we had wanted an interaction term too
davis.mod3 <- lm(weight ~ height * repwt, data=Davis, subset=-12)
summary(davis.mod3)