Watching Driving Data Exploration and MANOVA
Author: S. Wu
# Check and load required R packages
# car: for Levene’s test, compute.es: for effect sizes, multcomp: for post hoc tests,
# ez: for ANOVA (ezANOVA()), pastecs: for descriptive statistics, WRS: for robust tests
pkg<-c("knitr", "ggplot2", "grid", "gridExtra", "plyr", "tidyr", "dplyr", "Hmisc", "reshape2", "compute.es", "car", "ez", "multcomp", "pastecs", "MASS", "mvoutlier", "mvnormtest", "heplots", "broom");
pkgCheck<-pkg %in% rownames(installed.packages());
for(i in 1:length(pkg)) {
if(pkgCheck[i]==FALSE) {
install.packages(pkg[i])
}
library(pkg[i],character.only = TRUE, quietly = TRUE)
};
rm(i,pkg,pkgCheck);
# install.packages("akima")
# install.packages("WRS", repos="http://R-Forge.R-project.org");
library(WRS)
options(scipen=999); #remove scientific notation in printingData Exploration
Load Data
data.flat<- read.csv('flat file.csv', stringsAsFactors = FALSE);
names(data.flat)## [1] "ss"
## [2] "condition"
## [3] "order"
## [4] "safe_agg"
## [5] "law_abiding_agg"
## [6] "speed_agg"
## [7] "competent_agg"
## [8] "satisfaction_agg"
## [9] "recommend_agg"
## [10] "buy_agg"
## [11] "style_agg"
## [12] "rideagain_agg"
## [13] "safe_pass"
## [14] "law_abiding_pass"
## [15] "speed_pass"
## [16] "competent_pass"
## [17] "satisfaction_pass"
## [18] "recommend_pass"
## [19] "buy_pass"
## [20] "style_pass"
## [21] "rideagain_pass"
## [22] "prefer"
## [23] "whichagg"
## [24] "age"
## [25] "gender"
## [26] "liscence"
## [27] "driveregularly"
## [28] "oftendrive"
## [29] "Have.you.ever.been.in.an.autonomous..self.driving..vehicle."
## [30] "usedcriuse"
## [31] "usedadvcriuse"
## [32] "usedvoice"
## [33] "automationcar"
## [34] "oftencriuse"
## [35] "oftenvoice"names(data.flat)<- gsub('law_abiding','lawAbiding',names(data.flat));
names(data.flat)[29]<- 'selfDrive.Veh.Experience';
names(data.flat)[3]<- 'firstDrive'
Ratings Overview by Drive/Order/Condition
# Long flatfile, different rows for passive/aggressive drives & variables
data.ratings<- data.flat[,-c(22:35)];
data.ratings<- data.ratings %>%
mutate(ss= factor(ss),
condition=factor(condition, labels=c('untimed','timed')),
firstDrive=factor(firstDrive, labels=c('pass','agg')),
rideagain_agg=abs(rideagain_agg-2), # transform to 0=no, 1=yes
rideagain_pass=abs(rideagain_pass-2)) %>%
gather(var_drive, rating, 4:21) %>%
separate(var_drive, c('var', 'drive')) %>%
mutate(order= factor(ifelse(drive==firstDrive, 1, 2), labels=c('first','second')),
drive= factor(drive, levels=c('pass', 'agg')))
# dplyr::select(-firstDrive);fnum=fnum+1;
# exclude "rideagain", which is a yes/no question
data.ratingsSub1<- subset(data.ratings, var!='rideagain');
# average across all ratings excluding 'rideagain'
data.ratingsAverage1 <- ddply(data.ratingsSub1, .(ss, condition, firstDrive, drive, order), summarise, rating=mean(rating));
ggplot(data.ratingsAverage1, aes(x=drive, y=rating, colour=order)) +
geom_boxplot() +
labs(x = "Drive", y = "Average Rating") +
facet_grid(. ~ condition) +
theme(panel.grid.minor = element_blank()) +
scale_colour_manual("Order", values=c("steelblue1", "steelblue")) +
scale_y_continuous(limits=c(0,10), breaks=seq(0,10,2.5)) +
scale_x_discrete(labels=c('passive', 'aggressive')) +
theme_bw();
Figure 1 Boxplots of the watching driving data I
fnum=fnum+1;
labels <- c(pass='Passive Drive First', agg='Aggressive Drive First')
ggplot(data.ratingsAverage1, aes(x=drive, y=rating, colour=condition)) +
geom_boxplot() +
labs(x = "Drive", y = "Average Rating") +
facet_grid(. ~ firstDrive, labeller=labeller(firstDrive=labels)) +
theme(panel.grid.minor = element_blank()) +
scale_colour_manual("Condition", values=c("steelblue1", "steelblue")) +
scale_y_continuous(limits=c(0,10), breaks=seq(0,10,2.5)) +
scale_x_discrete(labels=c('passive', 'aggressive')) +
theme_bw();
Figure 2 Boxplots of the watching driving data II
fnum=fnum+1;
# Definition of whiskers
f <- function(x) {
r <- quantile(x, probs = c(0.05, 0.25, 0.5, 0.75, 0.95))
names(r) <- c("ymin", "lower", "middle", "upper", "ymax")
r
}
# Add points outside of whiskers
o <- function(x) {
subset(x, x < quantile(x,probs=0.05) | quantile(x,probs=0.95) < x)
}
ggplot(data.ratingsSub1, aes(x=var, y=rating, colour = drive)) +
stat_summary(fun.data = f, geom = 'boxplot', position=position_dodge(width=1), size = 0.6) +
stat_summary(fun.y = o, geom="point", position=position_dodge(width=1)) +
labs(x = "", y = "Rating") +
scale_y_continuous(limits=c(1,10), breaks=seq(1,10,1)) +
scale_colour_manual("", values=c("blue", "red"),
labels=c("passive drive", "aggressive drive")) +
theme_bw() #+ 
# theme(axis.text.x = element_text(angle = 45, hjust = 1)) #legend.position="none", Figure 3 Boxplots of the watching driving data III
# graph
fnum=fnum+1
g1<- ggplot(data.ratingsSub1, aes(x=var, y=rating, colour=drive)) +
stat_summary(fun.data = "mean_cl_normal", size = 0.6) +
labs(x = "", y = "Rating") +
scale_y_continuous(limits=c(0,10), breaks=seq(0,10,2.5)) +
scale_colour_manual("", values=c("yellowgreen", "steelblue"),
labels=c("passive drive", "aggressive drive")) +
theme_bw() +
theme(legend.position="none", axis.text.x = element_text(angle = 45, hjust = 1))
g2<- ggplot(data.ratingsSub1, aes(x=var, y=rating, colour=drive)) +
stat_summary(fun.data = "mean_cl_normal", size = 0.6) +
facet_grid(condition ~ order) + #, scales="free_y"
labs(x = "Catagory", y = "Rating") +
scale_y_continuous(limits=c(0,10), breaks=seq(0,10,2.5)) +
scale_colour_manual("", values=c("yellowgreen", "steelblue"),
labels=c("passive drive", "aggressive drive")) +
theme_bw() +
theme(legend.position="bottom", axis.text.x = element_text(angle = 45, hjust = 1))
g3<- ggplot(data.ratingsAverage1, aes(x=drive, y=rating, colour=drive)) +
stat_summary(fun.data = "mean_cl_normal", size = 0.6) +
labs(x = "Drive", y = "Average Rating") +
scale_y_continuous(limits=c(0,10), breaks=seq(0,10,2.5)) +
scale_colour_manual("", values=c("yellowgreen", "steelblue"),
labels=c("passive drive", "aggressive drive")) +
theme_bw() +
theme(legend.position="none")
grid.arrange(grid.arrange(g3, g1, ncol=2), g2, nrow=2)
Figure 4 Ratings by Drive / Order / Condition
# graph - rideagain
fnum=fnum+1
# 'ridegagin' ratings only
data.ratingsSub2<- subset(data.ratings, var=='rideagain');
# f <- function(x) {
# r <- quantile(x, probs = c(0.10, 0.25, 0.5, 0.75, 0.90))
# names(r) <- c("ymin", "lower", "middle", "upper", "ymax")
# r
# }
g1<- ggplot(data.ratingsSub2, aes(x=drive, y=rating, colour=drive)) +
stat_summary(fun.data = "mean_cl_normal", size = 1) +
labs(x = "", y = "Would you ride in this car again?") +
scale_y_continuous(limits=c(-0.3,1.3), breaks=c(0, 0.5, 1),
labels=c('No', '', 'Yes')) +
scale_x_discrete(labels=c('passive', 'aggressive')) +
scale_colour_manual("", values=c("yellowgreen", "steelblue"),
labels=c("passive drive", "aggressive drive")) +
theme_bw() +
theme(legend.position="none")
g2<- ggplot(data.ratingsSub2, aes(x=drive, y=rating, colour=drive)) +
stat_summary(fun.data = "mean_cl_normal", size = 1) +
# stat_summary(fun.data = f, geom='boxplot', position=position_dodge(width=1.5)) +
facet_grid(condition ~ order) + #, scales="free_y"
labs(x = "", y = "") +
scale_y_continuous(limits=c(-0.3,1.3), breaks=c(0, 0.5, 1),
labels=c('No', '', 'Yes')) +
scale_x_discrete(labels=c('passive', 'aggressive')) +
scale_colour_manual("", values=c("yellowgreen", "steelblue"),
labels=c("passive drive", "aggressive drive")) +
theme_bw() +
theme(legend.position="none")
grid.arrange(g1, g2, ncol=2)
Figure 5 Would you ride in this car again?
MANOVA + Univariate
1. Descriptive Stats
# Wide flatfile, different rows for passive/aggressive drives, different columns for each variables
data.ratingsWide<- dcast(data.ratings, ss+condition+drive+firstDrive+order ~ var, value.var = 'rating');
options(contrasts=c("contr.helmert","contr.poly"))
#drive
stats<- NULL;
for (i in 1:9) {
stats[[i]]<- by(data.ratingsWide[[i+5]], list(data.ratingsWide$drive), stat.desc, basic = FALSE);
}
names(stats)<- names(data.ratingsWide[c(6:14)]); print(stats)## $buy
## : pass
## median mean SE.mean CI.mean.0.95 var
## 3.0000000 3.5312500 0.4065601 0.8291847 5.2893145
## std.dev coef.var
## 2.2998510 0.6512852
## --------------------------------------------------------
## : agg
## median mean SE.mean CI.mean.0.95 var
## 7.0000000 6.6875000 0.4697235 0.9580073 7.0604839
## std.dev coef.var
## 2.6571571 0.3973319
##
## $competent
## : pass
## median mean SE.mean CI.mean.0.95 var
## 2.0000000 2.3125000 0.3026679 0.6172953 2.9314516
## std.dev coef.var
## 1.7121482 0.7403884
## --------------------------------------------------------
## : agg
## median mean SE.mean CI.mean.0.95 var
## 5.5000000 5.0625000 0.4305321 0.8780759 5.9314516
## std.dev coef.var
## 2.4354572 0.4810780
##
## $lawAbiding
## : pass
## median mean SE.mean CI.mean.0.95 var
## 1.0000000 2.2500000 0.3620061 0.7383162 4.1935484
## std.dev coef.var
## 2.0478155 0.9101402
## --------------------------------------------------------
## : agg
## median mean SE.mean CI.mean.0.95 var
## 6.0000000 5.8125000 0.3899894 0.7953886 4.8669355
## std.dev coef.var
## 2.2061132 0.3795464
##
## $recommend
## : pass
## median mean SE.mean CI.mean.0.95 var
## 2.0000000 3.0312500 0.4567689 0.9315864 6.6764113
## std.dev coef.var
## 2.5838752 0.8524125
## --------------------------------------------------------
## : agg
## median mean SE.mean CI.mean.0.95 var
## 7.5000000 6.6562500 0.4848706 0.9889002 7.5231855
## std.dev coef.var
## 2.7428426 0.4120702
##
## $rideagain
## : pass
## median mean SE.mean CI.mean.0.95 var
## 1.00000000 0.87500000 0.05939887 0.12114479 0.11290323
## std.dev coef.var
## 0.33601075 0.38401229
## --------------------------------------------------------
## : agg
## median mean SE.mean CI.mean.0.95 var
## 0.00000000 0.31250000 0.08324929 0.16978804 0.22177419
## std.dev coef.var
## 0.47092907 1.50697304
##
## $safe
## : pass
## median mean SE.mean CI.mean.0.95 var
## 1.0000000 1.9062500 0.3154598 0.6433845 3.1844758
## std.dev coef.var
## 1.7845100 0.9361364
## --------------------------------------------------------
## : agg
## median mean SE.mean CI.mean.0.95 var
## 7.0000000 6.2812500 0.4217723 0.8602103 5.6925403
## std.dev coef.var
## 2.3859045 0.3798455
##
## $satisfaction
## : pass
## median mean SE.mean CI.mean.0.95 var
## 2.0000000 3.1250000 0.4439331 0.9054076 6.3064516
## std.dev coef.var
## 2.5112649 0.8036048
## --------------------------------------------------------
## : agg
## median mean SE.mean CI.mean.0.95 var
## 6.0000000 6.0937500 0.4172669 0.8510214 5.5715726
## std.dev coef.var
## 2.3604179 0.3873506
##
## $speed
## : pass
## median mean SE.mean CI.mean.0.95 var
## 4.0000000 3.9062500 0.2391154 0.4876790 1.8296371
## std.dev coef.var
## 1.3526408 0.3462760
## --------------------------------------------------------
## : agg
## median mean SE.mean CI.mean.0.95 var
## 7.0000000 6.8750000 0.3412205 0.6959239 3.7258065
## std.dev coef.var
## 1.9302348 0.2807614
##
## $style
## : pass
## median mean SE.mean CI.mean.0.95 var
## 2.0000000 2.7812500 0.4229657 0.8626442 5.7247984
## std.dev coef.var
## 2.3926551 0.8602805
## --------------------------------------------------------
## : agg
## median mean SE.mean CI.mean.0.95 var
## 7.5000000 7.2187500 0.4253423 0.8674914 5.7893145
## std.dev coef.var
## 2.4060994 0.3333125#condition
stats<- NULL;
for (i in 1:9) {
stats[[i]]<- by(data.ratingsWide[[i+5]], list(data.ratingsWide$condition), stat.desc, basic = FALSE);
}
names(stats)<- names(data.ratingsWide[c(6:14)]); print(stats)## $buy
## : untimed
## median mean SE.mean CI.mean.0.95 var
## 4.0000000 4.5000000 0.5236349 1.0679605 8.7741935
## std.dev coef.var
## 2.9621265 0.6582503
## --------------------------------------------------------
## : timed
## median mean SE.mean CI.mean.0.95 var
## 5.0000000 5.7187500 0.4984540 1.0166037 7.9506048
## std.dev coef.var
## 2.8196817 0.4930591
##
## $competent
## : untimed
## median mean SE.mean CI.mean.0.95 var
## 3.0000000 3.7500000 0.4445004 0.9065646 6.3225806
## std.dev coef.var
## 2.5144742 0.6705265
## --------------------------------------------------------
## : timed
## median mean SE.mean CI.mean.0.95 var
## 3.0000000 3.6250000 0.4484516 0.9146231 6.4354839
## std.dev coef.var
## 2.5368255 0.6998139
##
## $lawAbiding
## : untimed
## median mean SE.mean CI.mean.0.95 var
## 3.5000000 4.0312500 0.4825261 0.9841184 7.4506048
## std.dev coef.var
## 2.7295796 0.6771050
## --------------------------------------------------------
## : timed
## median mean SE.mean CI.mean.0.95 var
## 4.0000000 4.0312500 0.5049840 1.0299216 8.1602823
## std.dev coef.var
## 2.8566208 0.7086191
##
## $recommend
## : untimed
## median mean SE.mean CI.mean.0.95 var
## 4.0000000 4.6875000 0.5340002 1.0891007 9.1250000
## std.dev coef.var
## 3.0207615 0.6444291
## --------------------------------------------------------
## : timed
## median mean SE.mean CI.mean.0.95 var
## 3.5000000 5.0000000 0.6074139 1.2388288 11.8064516
## std.dev coef.var
## 3.4360517 0.6872103
##
## $rideagain
## : untimed
## median mean SE.mean CI.mean.0.95 var
## 1.00000000 0.62500000 0.08695104 0.17733782 0.24193548
## std.dev coef.var
## 0.49186938 0.78699100
## --------------------------------------------------------
## : timed
## median mean SE.mean CI.mean.0.95 var
## 1.00000000 0.56250000 0.08909831 0.18171719 0.25403226
## std.dev coef.var
## 0.50401613 0.89602867
##
## $safe
## : untimed
## median mean SE.mean CI.mean.0.95 var
## 2.5000000 4.0937500 0.5281575 1.0771843 8.9264113
## std.dev coef.var
## 2.9877100 0.7298223
## --------------------------------------------------------
## : timed
## median mean SE.mean CI.mean.0.95 var
## 3.0000000 4.0937500 0.5542348 1.1303692 9.8296371
## std.dev coef.var
## 3.1352252 0.7658565
##
## $satisfaction
## : untimed
## median mean SE.mean CI.mean.0.95 var
## 5.0000000 4.5312500 0.5019807 1.0237964 8.0635081
## std.dev coef.var
## 2.8396317 0.6266773
## --------------------------------------------------------
## : timed
## median mean SE.mean CI.mean.0.95 var
## 5.0000000 4.6875000 0.5108451 1.0418754 8.3508065
## std.dev coef.var
## 2.8897762 0.6164856
##
## $speed
## : untimed
## median mean SE.mean CI.mean.0.95 var
## 5.0000000 5.6875000 0.4222575 0.8611999 5.7056452
## std.dev coef.var
## 2.3886492 0.4199823
## --------------------------------------------------------
## : timed
## median mean SE.mean CI.mean.0.95 var
## 5.0000000 5.0937500 0.3630054 0.7403544 4.2167339
## std.dev coef.var
## 2.0534687 0.4031350
##
## $style
## : untimed
## median mean SE.mean CI.mean.0.95 var
## 5.0000000 5.1250000 0.6086573 1.2413647 11.8548387
## std.dev coef.var
## 3.4430856 0.6718216
## --------------------------------------------------------
## : timed
## median mean SE.mean CI.mean.0.95 var
## 5.0000000 4.8750000 0.5531253 1.1281065 9.7903226
## std.dev coef.var
## 3.1289491 0.6418357#dirve x condition
stats<- NULL;
for (i in 1:9) {
stats[[i]]<- by(data.ratingsWide[[i+5]], list(data.ratingsWide$drive, data.ratingsWide$condition), stat.desc, basic = FALSE);
}
names(stats)<- names(data.ratingsWide[c(6:14)]); print(stats)## $buy
## : pass
## : untimed
## median mean SE.mean CI.mean.0.95 var
## 2.0000000 2.5000000 0.3763863 0.8022485 2.2666667
## std.dev coef.var
## 1.5055453 0.6022181
## --------------------------------------------------------
## : agg
## : untimed
## median mean SE.mean CI.mean.0.95 var
## 7.0000000 6.5000000 0.6770032 1.4429982 7.3333333
## std.dev coef.var
## 2.7080128 0.4166174
## --------------------------------------------------------
## : pass
## : timed
## median mean SE.mean CI.mean.0.95 var
## 4.0000000 4.5625000 0.6322496 1.3476082 6.3958333
## std.dev coef.var
## 2.5289985 0.5543010
## --------------------------------------------------------
## : agg
## : timed
## median mean SE.mean CI.mean.0.95 var
## 7.5000000 6.8750000 0.6700435 1.4281640 7.1833333
## std.dev coef.var
## 2.6801741 0.3898435
##
## $competent
## : pass
## : untimed
## median mean SE.mean CI.mean.0.95 var
## 2.0000000 2.1250000 0.3966001 0.8453332 2.5166667
## std.dev coef.var
## 1.5864005 0.7465414
## --------------------------------------------------------
## : agg
## : untimed
## median mean SE.mean CI.mean.0.95 var
## 5.5000000 5.3750000 0.5543389 1.1815455 4.9166667
## std.dev coef.var
## 2.2173558 0.4125313
## --------------------------------------------------------
## : pass
## : timed
## median mean SE.mean CI.mean.0.95 var
## 2.0000000 2.5000000 0.4654747 0.9921358 3.4666667
## std.dev coef.var
## 1.8618987 0.7447595
## --------------------------------------------------------
## : agg
## : timed
## median mean SE.mean CI.mean.0.95 var
## 5.5000000 4.7500000 0.6677075 1.4231849 7.1333333
## std.dev coef.var
## 2.6708301 0.5622800
##
## $lawAbiding
## : pass
## : untimed
## median mean SE.mean CI.mean.0.95 var
## 1.0000000 1.9375000 0.3815402 0.8132337 2.3291667
## std.dev coef.var
## 1.5261608 0.7876959
## --------------------------------------------------------
## : agg
## : untimed
## median mean SE.mean CI.mean.0.95 var
## 6.5000000 6.1250000 0.4819665 1.0272872 3.7166667
## std.dev coef.var
## 1.9278658 0.3147536
## --------------------------------------------------------
## : pass
## : timed
## median mean SE.mean CI.mean.0.95 var
## 1.0000000 2.5625000 0.6189288 1.3192156 6.1291667
## std.dev coef.var
## 2.4757154 0.9661328
## --------------------------------------------------------
## : agg
## : timed
## median mean SE.mean CI.mean.0.95 var
## 6.0000000 5.5000000 0.6191392 1.3196639 6.1333333
## std.dev coef.var
## 2.4765567 0.4502830
##
## $recommend
## : pass
## : untimed
## median mean SE.mean CI.mean.0.95 var
## 2.0000000 2.3750000 0.4269563 0.9100358 2.9166667
## std.dev coef.var
## 1.7078251 0.7190843
## --------------------------------------------------------
## : agg
## : untimed
## median mean SE.mean CI.mean.0.95 var
## 8.0000000 7.0000000 0.5322906 1.1345507 4.5333333
## std.dev coef.var
## 2.1291626 0.3041661
## --------------------------------------------------------
## : pass
## : timed
## median mean SE.mean CI.mean.0.95 var
## 2.0000000 3.6875000 0.7890857 1.6818964 9.9625000
## std.dev coef.var
## 3.1563428 0.8559574
## --------------------------------------------------------
## : agg
## : timed
## median mean SE.mean CI.mean.0.95 var
## 7.0000000 6.3125000 0.8201562 1.7481216 10.7625000
## std.dev coef.var
## 3.2806249 0.5197030
##
## $rideagain
## : pass
## : untimed
## median mean SE.mean CI.mean.0.95 var
## 1.0000000 0.9375000 0.0625000 0.1332156 0.0625000
## std.dev coef.var
## 0.2500000 0.2666667
## --------------------------------------------------------
## : agg
## : untimed
## median mean SE.mean CI.mean.0.95 var
## 0.0000000 0.3125000 0.1196784 0.2550884 0.2291667
## std.dev coef.var
## 0.4787136 1.5318834
## --------------------------------------------------------
## : pass
## : timed
## median mean SE.mean CI.mean.0.95 var
## 1.0000000 0.8125000 0.1007782 0.2148037 0.1625000
## std.dev coef.var
## 0.4031129 0.4961389
## --------------------------------------------------------
## : agg
## : timed
## median mean SE.mean CI.mean.0.95 var
## 0.0000000 0.3125000 0.1196784 0.2550884 0.2291667
## std.dev coef.var
## 0.4787136 1.5318834
##
## $safe
## : pass
## : untimed
## median mean SE.mean CI.mean.0.95 var
## 1.0000000 1.3750000 0.1547848 0.3299160 0.3833333
## std.dev coef.var
## 0.6191392 0.4502830
## --------------------------------------------------------
## : agg
## : untimed
## median mean SE.mean CI.mean.0.95 var
## 7.0000000 6.8125000 0.3788002 0.8073935 2.2958333
## std.dev coef.var
## 1.5152008 0.2224148
## --------------------------------------------------------
## : pass
## : timed
## median mean SE.mean CI.mean.0.95 var
## 1.0000000 2.4375000 0.5913878 1.2605134 5.5958333
## std.dev coef.var
## 2.3655514 0.9704826
## --------------------------------------------------------
## : agg
## : timed
## median mean SE.mean CI.mean.0.95 var
## 6.0000000 5.7500000 0.7444237 1.5867016 8.8666667
## std.dev coef.var
## 2.9776949 0.5178600
##
## $satisfaction
## : pass
## : untimed
## median mean SE.mean CI.mean.0.95 var
## 2.0000000 2.3125000 0.4254287 0.9067798 2.8958333
## std.dev coef.var
## 1.7017148 0.7358767
## --------------------------------------------------------
## : agg
## : untimed
## median mean SE.mean CI.mean.0.95 var
## 7.0000000 6.7500000 0.4518481 0.9630913 3.2666667
## std.dev coef.var
## 1.8073922 0.2677618
## --------------------------------------------------------
## : pass
## : timed
## median mean SE.mean CI.mean.0.95 var
## 3.0000000 3.9375000 0.7386291 1.5743506 8.7291667
## std.dev coef.var
## 2.9545163 0.7503534
## --------------------------------------------------------
## : agg
## : timed
## median mean SE.mean CI.mean.0.95 var
## 5.5000000 5.4375000 0.6768108 1.4425882 7.3291667
## std.dev coef.var
## 2.7072434 0.4978838
##
## $speed
## : pass
## : untimed
## median mean SE.mean CI.mean.0.95 var
## 4.0000000 3.8125000 0.3318980 0.7074237 1.7625000
## std.dev coef.var
## 1.3275918 0.3482208
## --------------------------------------------------------
## : agg
## : untimed
## median mean SE.mean CI.mean.0.95 var
## 7.5000000 7.5625000 0.3975838 0.8474299 2.5291667
## std.dev coef.var
## 1.5903354 0.2102923
## --------------------------------------------------------
## : pass
## : timed
## median mean SE.mean CI.mean.0.95 var
## 4.0000000 4.0000000 0.3535534 0.7535812 2.0000000
## std.dev coef.var
## 1.4142136 0.3535534
## --------------------------------------------------------
## : agg
## : timed
## median mean SE.mean CI.mean.0.95 var
## 6.0000000 6.1875000 0.5100551 1.0871568 4.1625000
## std.dev coef.var
## 2.0402206 0.3297326
##
## $style
## : pass
## : untimed
## median mean SE.mean CI.mean.0.95 var
## 2.0000000 2.5000000 0.5773503 1.2305930 5.3333333
## std.dev coef.var
## 2.3094011 0.9237604
## --------------------------------------------------------
## : agg
## : untimed
## median mean SE.mean CI.mean.0.95 var
## 8.0000000 7.7500000 0.5283622 1.1261775 4.4666667
## std.dev coef.var
## 2.1134490 0.2727031
## --------------------------------------------------------
## : pass
## : timed
## median mean SE.mean CI.mean.0.95 var
## 2.0000000 3.0625000 0.6289459 1.3405664 6.3291667
## std.dev coef.var
## 2.5157835 0.8214803
## --------------------------------------------------------
## : agg
## : timed
## median mean SE.mean CI.mean.0.95 var
## 7.0000000 6.6875000 0.6564980 1.3992923 6.8958333
## std.dev coef.var
## 2.6259919 0.3926717# by(data.ratingsWide[, 6:14], data.ratingsWide$drive, cov);
2. Construct Contrasts
contrasts(data.ratingsWide$drive)<- contr.helmert(2) #=c(-1, 1), passive=-1, aggresive=1
contrasts(data.ratingsWide$condition)<- contr.helmert(2) #=c(-1, 1), untimed=-1, timed=1
# contrasts(data.ratingsWide$order)<- contr.helmert(2) #=c(-1, 1), first=-1, second=1
# contrasts(data.ratingsWide$firstDrive)<- contr.helmert(2) #=c(-1, 1), passFirst=-1, aggFirst=1
attr(data.ratingsWide$drive, 'contrasts'); attr(data.ratingsWide$condition, 'contrasts');
# attr(data.ratingsWide$order, 'contrasts'); attr(data.ratingsWide$firstDrive, 'contrasts');
3. MANOVA
DV: 9 ratings: safe, lawAbiding, speed, competent, satisfaction, recommend, buy, style, rideagain(0-=no/1=yes)
IV: drive(within subject), condition(between subject)
outcome<- as.matrix(data.ratingsWide[,6:14]);
# model2<-manova(outcome ~ condition*drive, data = data.ratingsWide);
model<-manova(outcome ~ condition*drive + Error(ss/drive)+condition, data = data.ratingsWide);
modSum<- summary(model, intercept = TRUE, test = "Wilks") #default test=Pillai's trace, or, test='Hotelling', test='Roy'
modSumMatrix<-rbind(modSum[[1]][[4]], modSum[[2]][[4]]);
p<- round(modSumMatrix[,6], 3); df<- round(modSumMatrix[,4:5], 3); f<- round(modSumMatrix[,3], 3);
lambda<- round(modSumMatrix[,2], 3);
etaSQ<- round(summary.lm(model[[3]])$r.squared, 3);
# p<- round(modSum[[1]][,6], 3); df<- round(modSum[[4]][,4:5], 3); f<- round(modSum[[4]][,3], 3);
# lambda<- round(modSum[[4]][,2], 3);
# etaSQ<- round(summary.lm(model[2:3])$r.squared, 3);
modSum;##
## Error: ss
## Df Wilks approx F num Df den Df Pr(>F)
## condition 1 0.71007 0.99811 9 22 0.4698
## Residuals 30
##
## Error: ss:drive
## Df Wilks approx F num Df den Df Pr(>F)
## drive 1 0.20519 9.4688 9 22 0.000009455 ***
## condition:drive 1 0.69633 1.0660 9 22 0.424
## Residuals 30
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Interaction of condition and drive was not significant, Wilk’s Lambda = 0.696, F(9, 22) = 1.066, p = 0.424.
Main effect of condition was not significant, Wilk’s Lambda = 0.71, F(9, 22) = 0.998, p = 0.47.
Main effect of drive was significant, Wilk’s Lambda = 0.205, F(9, 22) = 9.469, p < 0.001.
\(\eta^2\) = 0.522 (This is not correct.)
4. Univariate Tests
#Univariate ANOVA table
modUnivariate<- summary.aov(model[[3]]);
# modUnivariate<- summary.aov(model);
# a<- bind_rows(modUnivariate, .id='response')
modUnivariate;## Response buy :
## Df Sum Sq Mean Sq F value Pr(>F)
## drive 1 159.391 159.391 28.0093 0.0000102 ***
## condition:drive 1 11.391 11.391 2.0016 0.1674
## Residuals 30 170.719 5.691
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Response competent :
## Df Sum Sq Mean Sq F value Pr(>F)
## drive 1 121 121.000 26.1151 0.00001706 ***
## condition:drive 1 4 4.000 0.8633 0.3602
## Residuals 30 139 4.633
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Response lawAbiding :
## Df Sum Sq Mean Sq F value Pr(>F)
## drive 1 203.06 203.062 51.7631 0.00000005247 ***
## condition:drive 1 6.25 6.250 1.5932 0.2166
## Residuals 30 117.69 3.923
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Response recommend :
## Df Sum Sq Mean Sq F value Pr(>F)
## drive 1 210.25 210.250 26.9840 0.00001345 ***
## condition:drive 1 16.00 16.000 2.0535 0.1622
## Residuals 30 233.75 7.792
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Response rideagain :
## Df Sum Sq Mean Sq F value Pr(>F)
## drive 1 5.0625 5.0625 31.1538 0.000004513 ***
## condition:drive 1 0.0625 0.0625 0.3846 0.5398
## Residuals 30 4.8750 0.1625
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Response safe :
## Df Sum Sq Mean Sq F value Pr(>F)
## drive 1 306.250 306.250 66.2461 0.000000004375 ***
## condition:drive 1 18.062 18.062 3.9072 0.05734 .
## Residuals 30 138.688 4.623
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Response satisfaction :
## Df Sum Sq Mean Sq F value Pr(>F)
## drive 1 141.016 141.016 24.6002 0.00002609 ***
## condition:drive 1 34.516 34.516 6.0213 0.02016 *
## Residuals 30 171.969 5.732
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Response speed :
## Df Sum Sq Mean Sq F value Pr(>F)
## drive 1 141.016 141.016 51.7686 0.00000005242 ***
## condition:drive 1 9.766 9.766 3.5851 0.06798 .
## Residuals 30 81.719 2.724
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Response style :
## Df Sum Sq Mean Sq F value Pr(>F)
## drive 1 315.062 315.062 38.6777 0.0000007569 ***
## condition:drive 1 10.563 10.563 1.2967 0.2638
## Residuals 30 244.375 8.146
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1results<- lapply(modUnivariate, function(x) {
x<- round(x, 3);
dfError<- x['Residuals', 1];
x$results<- paste0('Effect of ', row.names(x), ' is ', ifelse(x[[5]]>0.05, "not significant", "significant"),
': F(', x[[1]], ', ', dfError, ') = ', x[[4]],
', p ', ifelse(x[[5]]==0, "< 0.001", paste("=", x[[5]])), '.');
# return(paste(x$results[1:nrow(x)-1], collapse=' '));
return(x$results[1:nrow(x)-1]);
});
results## $` Response buy`
## [1] "Effect of drive is significant: F(1, 30) = 28.009, p < 0.001."
## [2] "Effect of condition:drive is not significant: F(1, 30) = 2.002, p = 0.167."
##
## $` Response competent`
## [1] "Effect of drive is significant: F(1, 30) = 26.115, p < 0.001."
## [2] "Effect of condition:drive is not significant: F(1, 30) = 0.863, p = 0.36."
##
## $` Response lawAbiding`
## [1] "Effect of drive is significant: F(1, 30) = 51.763, p < 0.001."
## [2] "Effect of condition:drive is not significant: F(1, 30) = 1.593, p = 0.217."
##
## $` Response recommend`
## [1] "Effect of drive is significant: F(1, 30) = 26.984, p < 0.001."
## [2] "Effect of condition:drive is not significant: F(1, 30) = 2.053, p = 0.162."
##
## $` Response rideagain`
## [1] "Effect of drive is significant: F(1, 30) = 31.154, p < 0.001."
## [2] "Effect of condition:drive is not significant: F(1, 30) = 0.385, p = 0.54."
##
## $` Response safe`
## [1] "Effect of drive is significant: F(1, 30) = 66.246, p < 0.001."
## [2] "Effect of condition:drive is not significant: F(1, 30) = 3.907, p = 0.057."
##
## $` Response satisfaction`
## [1] "Effect of drive is significant: F(1, 30) = 24.6, p < 0.001."
## [2] "Effect of condition:drive is significant: F(1, 30) = 6.021, p = 0.02."
##
## $` Response speed`
## [1] "Effect of drive is significant: F(1, 30) = 51.769, p < 0.001."
## [2] "Effect of condition:drive is not significant: F(1, 30) = 3.585, p = 0.068."
##
## $` Response style`
## [1] "Effect of drive is significant: F(1, 30) = 38.678, p < 0.001."
## [2] "Effect of condition:drive is not significant: F(1, 30) = 1.297, p = 0.264."
4. Graphs
fnum=fnum+1;
ggplot(data.ratingsAverage, aes(x=condition, y=rating, colour=firstDrive)) +
stat_summary(fun.y = mean, geom="point") +
stat_summary(fun.y = mean, geom="line", aes(group=firstDrive)) +
stat_summary(fun.data = mean_cl_boot, geom="errorbar", width=0.2) +
labs(x = "Condition", y = "Average Rating") +
theme_bw() +
theme(panel.grid.major = element_blank()) +
scale_colour_manual("firstDrive", values=c("steelblue", "black")) +
scale_y_continuous(limits=c(1, 10), breaks=seq(1, 10, by=1));
Figure 6 Graph of interaction between condition and firstDrive
fnum=fnum+1;
ggplot(data.ratingsAverage, aes(x=drive, y=rating, colour=firstDrive)) +
stat_summary(fun.y = mean, geom="point") +
stat_summary(fun.y = mean, geom="line", aes(group=firstDrive)) +
stat_summary(fun.data = mean_cl_boot, geom="errorbar", width=0.2) +
labs(x = "Drive", y = "Average Rating") +
theme_bw() +
theme(panel.grid.major = element_blank()) +
scale_colour_manual("firstDrive", values=c("steelblue", "black")) +
scale_y_continuous(limits=c(1, 10), breaks=seq(1, 10, by=1));
Figure 7 Graph of interaction between drive and firstDrive
fnum=fnum+1;
ggplot(data.ratingsAverage, aes(x=condition, y=rating, colour=drive)) +
stat_summary(fun.y = mean, geom="point") +
stat_summary(fun.y = mean, geom="line", aes(group=drive)) +
stat_summary(fun.data = mean_cl_boot, geom="errorbar", width=0.2) +
labs(x = "Condition", y = "Average Rating") +
theme_bw() +
theme(panel.grid.major = element_blank()) +
scale_colour_manual("Drive", values=c("steelblue", "black")) +
scale_y_continuous(limits=c(1, 10), breaks=seq(1, 10, by=1));
Figure 8 Graph of interaction between condition and drive
fnum=fnum+1;
ggplot(data.ratingsAverage, aes(x=drive, y=rating, colour=condition)) +
stat_summary(fun.y = mean, geom="point") +
stat_summary(fun.y = mean, geom="line", aes(group=condition)) +
stat_summary(fun.data = mean_cl_boot, geom="errorbar", width=0.2) +
labs(x = "Drive", y = "Average Rating") +
theme_bw() +
theme(panel.grid.major = element_blank()) +
scale_colour_manual("Condition", values=c("steelblue", "black")) +
scale_y_continuous(limits=c(1, 10), breaks=seq(1, 10, by=1));
Figure 9 Graph of interaction between condition and drive
fnum=fnum+1;
ggplot(data.ratingsAverage, aes(x=firstDrive, y=rating, colour=drive)) +
stat_summary(fun.y = mean, geom="point") +
stat_summary(fun.y = mean, geom="line", aes(group=drive)) +
stat_summary(fun.data = mean_cl_boot, geom="errorbar", width=0.2) +
labs(x = "First Drive", y = "Average Rating") +
facet_grid(. ~ condition) +
theme_bw() +
theme(panel.grid.major = element_blank()) +
scale_colour_manual("Drive", values=c("steelblue", "black"), labels=(c('passive', 'aggressive'))) +
scale_x_discrete(labels=c('passive first', 'aggressive first')) +
scale_y_continuous(limits=c(1, 10), breaks=seq(1, 10, by=1));
Figure 10 Graphs showing the condition by drive interaction for passive-first and aggressive-first drivers. Lines represent passive drive (light blue) and aggressive drive (dark blue)
fnum=fnum+1;
ggplot(data.ratingsAverage, aes(x=firstDrive, y=rating)) +
stat_summary(fun.y = mean, geom = "bar", fill = "White", colour = "Black") +
stat_summary(fun.data = mean_cl_normal, geom = "pointrange") +
labs(x = "First Drive", y = "Average Rating") +
scale_x_discrete(labels=c('passive', 'aggressive')) +
theme_bw() +
theme(panel.grid.major = element_blank()) +
scale_y_continuous(limits=c(0, 10), breaks=seq(0, 10, by=1));
Figure 11 Error bar graph of the main effect of firstDrive
fnum=fnum+1;
ggplot(data.ratingsAverage, aes(x=drive, y=rating)) +
stat_summary(fun.y = mean, geom = "bar", fill = "White", colour = "Black") +
stat_summary(fun.data = mean_cl_normal, geom = "pointrange") +
labs(x = "Drive", y = "Average Rating") +
theme_bw() +
theme(panel.grid.major = element_blank()) +
scale_y_continuous(limits=c(0, 10), breaks=seq(0, 10, by=1));
Figure 12 Error bar graph of the main effect of drive
fnum=fnum+1;
ggplot(data.ratingsAverage, aes(x=condition, y=rating)) +
stat_summary(fun.y = mean, geom = "bar", fill = "White", colour = "Black") +
stat_summary(fun.data = mean_cl_normal, geom = "pointrange") +
labs(x = "Condition", y = "Average Rating") +
theme_bw() +
theme(panel.grid.major = element_blank()) +
scale_y_continuous(limits=c(0, 10), breaks=seq(0, 10, by=1));
Figure 13 Error bar graph of the main effect of condition
Repeated-Measures ANOVA
DV: average rating (include ‘rideagain’)
IV: drive(within subject), condition(between subject)
1. Descriptive Stats
# levels(data.ratingsAverage1$firstDrive)<- c('passFirst', 'aggFirst');
by(data.ratingsAverage$rating, list(data.ratingsAverage$drive), stat.desc, basic = FALSE)## : pass
## median mean SE.mean CI.mean.0.95 var
## 2.2222222 2.6354167 0.2453893 0.5004748 1.9269091
## std.dev coef.var
## 1.3881315 0.5267218
## --------------------------------------------------------
## : agg
## median mean SE.mean CI.mean.0.95 var
## 5.8888889 5.6666667 0.2832349 0.5776614 2.5671047
## std.dev coef.var
## 1.6022187 0.2827445by(data.ratingsAverage$rating, list(data.ratingsAverage$drive, data.ratingsAverage$condition), stat.desc, basic = FALSE)## : pass
## : untimed
## median mean SE.mean CI.mean.0.95 var
## 2.0000000 2.2083333 0.1641398 0.3498557 0.4310700
## std.dev coef.var
## 0.6565592 0.2973098
## --------------------------------------------------------
## : agg
## : untimed
## median mean SE.mean CI.mean.0.95 var
## 5.8888889 6.0208333 0.2596331 0.5533948 1.0785494
## std.dev coef.var
## 1.0385323 0.1724898
## --------------------------------------------------------
## : pass
## : timed
## median mean SE.mean CI.mean.0.95 var
## 2.5000000 3.0625000 0.4445566 0.9475499 3.1620885
## std.dev coef.var
## 1.7782262 0.5806453
## --------------------------------------------------------
## : agg
## : timed
## median mean SE.mean CI.mean.0.95 var
## 5.611111 5.312500 0.497444 1.060277 3.959208
## std.dev coef.var
## 1.989776 0.374546# by(data.ratingsAverage$rating, list(data.ratingsAverage$drive, data.ratingsAverage$order), stat.desc, basic = FALSE)
# by(data.ratingsAverage$rating, list(data.ratingsAverage$drive, data.ratingsAverage$condition, data.ratingsAverage$order), stat.desc, basic = FALSE)
2. Construct Contrasts
contrasts(data.ratingsAverage$drive)<- contr.helmert(2) #=c(-1, 1), passive=-1, aggresive=1
contrasts(data.ratingsAverage$condition)<- contr.helmert(2) #=c(-1, 1), untimed=-1, timed=1
contrasts(data.ratingsAverage$order)<- contr.helmert(2) #=c(-1, 1), first=-1, second=1
contrasts(data.ratingsAverage$firstDrive)<- contr.helmert(2) #=c(-1, 1), passFirst=-1, aggFirst=1
attr(data.ratingsAverage$drive, 'contrasts'); attr(data.ratingsAverage$condition, 'contrasts');
attr(data.ratingsAverage$order, 'contrasts'); attr(data.ratingsAverage$firstDrive, 'contrasts');
3. ANOVA
model1<-ezANOVA(data = data.ratingsAverage, dv = .(rating), wid = .(ss),
between = .(condition), within = .(drive),
type = 3, detailed = TRUE); #or type=2(default), type=1, specifying the Sums of Squares "type"
n<- nrow(model1[["ANOVA"]]);
df<- model1[["ANOVA"]][c(2:n), c('DFn', 'DFd')];
p<- round(model1[["ANOVA"]][c(2:n), 'p'], 3);
f<- round(model1[["ANOVA"]][c(2:n), 'F'], 3);
etaSQ<- round(model1[["ANOVA"]][c(2:n), 'ges'], 3); # generalized eta-squared
row.names(df)<- model1[["ANOVA"]][c(2:n), 'Effect'];
names(p)<- names(f)<- names(etaSQ)<- model1[["ANOVA"]][c(2:n), 'Effect'];
if(any(grep('Sphericity', names(model1)))) {
pMauchly<- round(model1[["Mauchly's Test for Sphericity"]]$p, 3);
if (any(pMauchly<=0.05)) {
msgMauchly<- "Where violations of sphericity were indicated by Mauchly's Test, degrees of freedom have been adjusted using the Greenhouse-Geisser correction.";
GGe<- round(model1[["Sphericity Corrections"]]$GGe, 3);
names(pMauchly)<- names(GGe)<- model1[["Sphericity Corrections"]]$Effect;
for (i in 1:length(GGe)) {
v<- names(GGe[i])
if (pMauchly[v]<=0.05) {
df[v,]<-df[v,]*GGe[v];
p[v]<- pMauchly[v];
}
}
} else { msgMauchly<- "" }
} else { msgMauchly<- "" }
model1;## $ANOVA
## Effect DFn DFd SSn SSd F
## 1 (Intercept) 1 30 1102.79340278 54.15239 610.9388870
## 2 condition 1 30 0.08506944 54.15239 0.0471278
## 3 drive 1 30 147.01562500 75.31134 58.5631407
## 4 condition:drive 1 30 9.76562500 75.31134 3.8901013
## p p<.05 ges
## 1 0.000000000000000000001675235 * 0.8949377280
## 2 0.829608332673397952383709253 0.0006566594
## 3 0.000000015519897797940396229 * 0.5317417735
## 4 0.057851210232398718458224351 0.0701405582Interaction of condition and drive was not significant, F(1, 30) = 3.89, p = 0.058, \(\eta^2\) = 0.07.
Main effect of condition was not significant, F(1, 30) = 0.047, p = 0.83, \(\eta^2\) = 0.001.
Main effect of drive was significant, F(1, 30) = 58.563, p < 0.001, \(\eta^2\) = 0.532.
4. Graphs
fnum=fnum+1;
ggplot(data.ratingsAverage, aes(x=condition, y=rating, colour=drive)) +
stat_summary(fun.y = mean, geom="point") +
stat_summary(fun.y = mean, geom="line", aes(group=drive)) +
stat_summary(fun.data = mean_cl_boot, geom="errorbar", width=0.2) +
labs(x = "Condition", y = "Average Rating") +
theme_bw() +
theme(panel.grid.major = element_blank()) +
scale_colour_manual("Drive", values=c("steelblue", "black")) +
scale_y_continuous(limits=c(1, 10), breaks=seq(1, 10, by=1));
Figure 6 Graph of interaction between condition and drive I
fnum=fnum+1;
ggplot(data.ratingsAverage, aes(x=drive, y=rating, colour=condition)) +
stat_summary(fun.y = mean, geom="point") +
stat_summary(fun.y = mean, geom="line", aes(group=condition)) +
stat_summary(fun.data = mean_cl_boot, geom="errorbar", width=0.2) +
labs(x = "Drive", y = "Average Rating") +
theme_bw() +
theme(panel.grid.major = element_blank()) +
scale_colour_manual("Condition", values=c("steelblue", "black")) +
scale_y_continuous(limits=c(1, 10), breaks=seq(1, 10, by=1));
Figure 7 Graph of interaction between condition and drive II
fnum=fnum+1;
ggplot(data.ratingsAverage, aes(x=drive, y=rating)) +
stat_summary(fun.y = mean, geom = "bar", fill = "White", colour = "Black") +
stat_summary(fun.data = mean_cl_normal, geom = "pointrange") +
labs(x = "Drive", y = "Average Rating") +
theme_bw() +
theme(panel.grid.major = element_blank()) +
scale_y_continuous(limits=c(0, 10), breaks=seq(0, 10, by=1));
Figure 8 Error bar graph of the main effect of drive
fnum=fnum+1;
ggplot(data.ratingsAverage, aes(x=condition, y=rating)) +
stat_summary(fun.y = mean, geom = "bar", fill = "White", colour = "Black") +
stat_summary(fun.data = mean_cl_normal, geom = "pointrange") +
labs(x = "Condition", y = "Average Rating") +
theme_bw() +
theme(panel.grid.major = element_blank()) +
scale_y_continuous(limits=c(0, 10), breaks=seq(0, 10, by=1));
Figure 9 Error bar graph of the main effect of condition