---
title: "One-way ANOVA and contrast analysis"
author: "Tera Letzring"
date: "November 2017"
output: html_document
---

```{r setup, include=FALSE}

knitr::opts_chunk$set(echo = TRUE)
library(pastecs)     #for stat.desc
library(multcomp)    #for Tukey posthoc test
library(compute.es)  #for effect sizes for pairs of groups
library(pwr)         #to compute power
library(gmodels)     #for contrast analysis
```

```{r read in data}

#set the working directory (wd) to the folder that contains the script file, and read in the data
x <- getwd()
setwd(x)

mydata = read.table("Expectation.csv", header=T, sep=",")  
attach(mydata)
head(mydata)
```

```{r calculate descriptive statistics and check homogeneity of variance assumption}

#Descriptive statistics by groups
tapply(rating, instruction, stat.desc)
tapply(rating, instruction, mean)
tapply(rating, instruction, sd)
table(instruction)

#Check for homogeneity of variance
tapply(rating, instruction, var)
bartlett.test(rating ~ instruction, data=mydata)
```

```{r anova, check assumptions, and post hocs}

#One-way ANOVA with equal variances
Anova.1 <- aov(rating ~ instruction, data=mydata)
summary(Anova.1)

#One-way ANOVA with unequal variances; uses the Welch correction
oneway.test(rating~instruction)

#Figures to check assumptions (you have to run the model before the plot will work)
plot(Anova.1)

#Bonferroni
pairwise.t.test(rating, instruction, paired=FALSE, p.adjust.method="bonferroni")

#unadjusted
pairwise.t.test(rating, instruction, paired=FALSE)
                
#Tukey
Tukey<-glht(Anova.1, linfct = mcp(instruction = "Tukey"))
summary(Tukey)   #output results
confint(Tukey)   #output confidence intervals
plot(confint(Tukey))   #create plots of the confidence intervals
```

```{r effect sizes}

#Calculate an overall effect size
summary.lm(Anova.1)   #eta-squared = Multiple R-squared

#calculate effect sizes for pairs
mes(-7.6875, -3.225806, 17.57094, 10.72912, 32, 31)  #good vs. scientific
```

```{r power}

#k=number of groups, n=number of observations per group, f=effect size similar to Cohen's d
#pwr.anova.test(k = NULL, n = NULL, f = NULL, sig.level = 0.05, power = NULL)
#to determine sample size

pwr.anova.test(k = 3, n = NULL, f = .25, sig.level = 0.05, power = .80)
```

```{r contrast analysis}

###planned contrasts 
fit.contrast(Anova.1, instruction, c(0, -1, 1), conf.int=.95, df=T)
```