diff --git a/ANOVA vs Regression.pdf b/ANOVA vs Regression.pdf
new file mode 100644
index 0000000..5e31f1d
Binary files /dev/null and b/ANOVA vs Regression.pdf differ
diff --git a/Peer Review Comments.Rmd b/Peer Review Comments.Rmd
new file mode 100644
index 0000000..8bedaab
--- /dev/null
+++ b/Peer Review Comments.Rmd	
@@ -0,0 +1,34 @@
+---
+title: "Peer Review"
+output: html_document
+---
+
+# Part 1:
+
+## ANOVA
+
+We like that this group used the following code:
+
+```{r}
+pchisq(teststat,df, lower=F)
+```
+
+Instead of doing 1-pchisq. Nifty.
+
+
+For lines 79, 81, and 83, include the results in hashtags, for more convienent comparisons. 
+
+Add an abline in your regression plot. 
+
+
+## Regression
+
+Looks great.
+
+# Part 2:
+
+Combine scripts for part 2 into one script. 
+
+In the ANOVA Monte carlo, delete the bottom lines to reduce reader confusion.
+
+At the end, if time, build a graph to visually display the p-values. 
\ No newline at end of file
diff --git a/Power Analysis.pdf b/Power Analysis.pdf
new file mode 100644
index 0000000..4ea994e
Binary files /dev/null and b/Power Analysis.pdf differ
diff --git a/anovaVSregression.Rmd b/anovaVSregression.Rmd
new file mode 100755
index 0000000..be0e1ee
--- /dev/null
+++ b/anovaVSregression.Rmd
@@ -0,0 +1,159 @@
+---
+title: "ANOVA vs Regression"
+output: html_document
+---
+
+# PART 1 REGRESSION VS ANOVA
+
+## An ANOVA-design and a regression-design experiment
+
+```{r include=FALSE}
+rm(list=ls())
+setwd("/Users/emilynonnamaker/Box/personalStuff/School/PhD/biocomputingFall2018/biocomputing_StatsGroupProject")
+library(ggplot2)
+library(arm)
+library(tidyr)
+library(reshape2)
+antibiotics <- read.csv("antibiotics.csv")
+sugar <- read.csv("sugar.csv")
+
+#playing around - check out the data
+boxplot(growth ~ trt, data=antibiotics)
+lm <- lm(growth ~ sugar, data=sugar)
+aov <- aov(growth ~ trt, data=antibiotics)
+
+# how a normal person would do this
+summary(aov)
+posthoc <- TukeyHSD(aov, "trt", conf.level=0.95) #a, b, ac, bd. 
+posthoc
+```
+
+## Now to do the actual assingment:
+
+```{r}
+# set up your data: 
+
+N <- length(antibiotics$growth)
+y <- antibiotics$growth
+x<- antibiotics$trt
+antibiotics$x1 <- ifelse(x=="ab1", 1, 0) #dummy variables
+antibiotics$x2 <- ifelse(x=="ab2", 1, 0) #dummy variables
+antibiotics$x3 <- ifelse(x=="ab3", 1, 0) #dummy variables
+
+#Build null model 
+nullmod<-function(p,x,y){
+  B0=p[1]
+  sigma=exp(p[2])
+  
+  expected=B0
+
+null=-sum(dnorm(x=y,mean=expected,sd=sigma, log = TRUE))
+return(null)
+}
+
+#Build full model
+fullmod<-function(p,x,y){
+B0=p[1]
+B1=p[2]
+B2=p[3]
+B3=p[4]
+sigma=exp(p[5])
+
+expected=B0+B1*x[,3]+B2*x[,4]+B3*x[,5]
+
+full=-sum(dnorm(x=y,mean=expected,sd=sigma,log = TRUE))
+return(full)
+}
+
+# Check fit
+nullguess <- c(1, 2)
+fullguess <- c(18, -12, -3, -4, 1)
+fitnull=optim(par=nullguess,fn=nullmod,y=antibiotics$growth) # converges
+fitfull=optim(par=fullguess,fn=fullmod,x=antibiotics,y=antibiotics$growth) # converges
+
+# Get t.statistic . 
+t.stat <- 2*(fitnull$value - fitfull$value)
+t.stat # 37.90134
+
+# get degrees of freedom
+df <- length(fitfull$par) - length(fitnull$par)
+
+# and p value
+1-pchisq(t.stat, df=df) # 2.965739e-08
+```
+
+Looks like there's a significant difference between the null and full anova model (t = 25.69, p < .001), suggesting that there are differences between treatments. 
+
+## Graphical visualization!
+
+Using a posthoc Tukey test (see hidden code) we can find significant differences between groups, and add these differences as letters to our visualization (see graph below)
+
+```{r}
+p <- ggplot(antibiotics, aes(x=trt, y=growth, fill=trt)) + 
+  geom_boxplot()
+p + geom_jitter(shape=16, position=position_jitter(0.2), alpha=0.3) + theme_classic() + labs(x = "Treatment", y= "Growth") + annotate("text", x = c(1, 2, 3, 4), y = c(23, 23, 23, 23), label = c("a", "b", "ac", "bd"))
+```
+
+
+## Regression
+
+Build the null model:
+ 
+```{r}
+lmnullmodreg<-function(p,x,y){
+  B0=p[1]
+  sigma=exp(p[2])
+  
+  expected=B0
+
+null=-sum(dnorm(x=y,mean=expected,sd=sigma, log = TRUE))
+return(null)
+}
+```
+
+Now build the extension:
+
+```{r}
+lmfullmodreg<-function(p,x,y){
+  B0=p[1]
+  B1=p[2]
+  sigma=exp(p[3])
+  
+  expected=B0+B1*x
+
+  null=-sum(dnorm(x=y,mean=expected,sd=sigma, log = TRUE))
+  return(null)
+}
+```
+
+Give it a guess:
+
+```{r}
+lmnullregguess <- c(1, 1)
+lmfullregguess <- c(1, 2, 3)
+lmfitregnull=optim(par=lmnullregguess,fn=lmnullmodreg,x=sugar$sugar,y=sugar$growth) # converges
+lmfitregfull=optim(par=lmfullregguess,fn=lmfullmodreg,x=sugar$sugar,y=sugar$growth) # converges
+```
+
+Now to compare:
+
+```{r}
+# Get t statistic
+t.statreg <- 2*(lmfitregnull$value - lmfitregfull$value)
+t.statreg # 39.92512
+
+# get degrees of freedom
+dfreg <- length(lmfitregfull$par) - length(lmfitregnull$par)
+
+#and p value
+1-pchisq(t.stat, df=dfreg) # 7.441429e-10, woot
+```
+
+Comparing between the null model of no effect of sugar concentration on growth to the model taking this affect into account, we can see that sugar concentration does have a significant affect on growth (t = 39.92512, p < .0001). Looking at the graph below, we can tell that sugar concentration has a positive affect, with higher concentrations of sugar increasing growth. 
+
+## Visualization!
+
+```{r}
+a <- ggplot(data=sugar,aes(x=sugar,y=growth))
+a + geom_point() + coord_cartesian() + labs( x = "sugar concentration", y = "growth") + theme_classic() + geom_smooth(method = "lm")
+```
diff --git a/statPowerAnalysis.Rmd b/statPowerAnalysis.Rmd
new file mode 100755
index 0000000..6cf52fe
--- /dev/null
+++ b/statPowerAnalysis.Rmd
@@ -0,0 +1,211 @@
+
+---
+title: "Power Analysis"
+output: html_notebook
+---
+
+# Part 2 Statistical Power Analysis
+
+```{r, include=FALSE}
+library(ggplot2)
+library(dplyr)
+library(reshape2)
+```
+
+### Step 1 - build the models
+
+```{r}
+
+# Null model
+nullmod<- function(p,x,y){
+  B0=p[1]
+  sigma=exp(p[2])
+  
+  expected=B0
+  
+  null=-sum(dnorm(x=y,mean=expected,sd=sigma, log = TRUE))
+  return(null)
+}
+
+# Regression model
+lin<-function(p,x,y){
+  B0=p[1]
+  B1=p[2]
+  sigma=exp(p[3])
+  
+  expected=B0+B1*x
+
+  lmmod=-sum(dnorm(x=y,mean=expected,sd=sigma, log = TRUE))
+  return(lmmod)
+}
+
+# T-test model
+ttest<-function(p,x,y){
+  B0=p[1]
+  B1=p[2]
+  sigma=exp(p[3])
+  x1=x
+  expected=B0+B1*x1
+
+  tmod=-sum(dnorm(x=y,mean=expected,sd=sigma, log = TRUE))
+  return(tmod)
+}
+
+# ANOVA - four groups
+anovamod<-function(p,x,y){
+B0=p[1]
+B1=p[2]
+B2=p[3]
+B3=p[4]
+sigma=exp(p[5])
+
+expected=B0+B1*x[,1]+B2*x[,2]+B3*x[,3]
+
+anova=-sum(dnorm(x=y,mean=expected,sd=sigma,log = TRUE))
+return(anova)
+}
+
+# Now an ANOVA with 8 groups
+eightanovamod<-function(p,x,y){
+  B0=p[1]
+  B1=p[2]
+  B2=p[3]
+  B3=p[4]
+  B4=p[5]
+  B5=p[6]
+  B6=p[7]
+  B7=p[8]
+  sigma=exp(p[9])
+  
+  expected=B0+B1*x[,1]+B2*x[,2]+B3*x[,3]+B4*x[,4]+B5*x[,5]+B6*x[,6]+B7*x[,7]
+  
+  eightanovamod=-sum(dnorm(x=y,mean=expected,sd=sigma, log = TRUE))
+  return(eightanovamod)
+}
+
+```
+
+### Step 2: Make your for loops!
+
+```{r}
+# Set up an empty data frame in which to store your p-values for each sigma and model.
+dfp <- setNames(data.frame(matrix(ncol = 5, nrow = 8)), c("Sigmas", "Regression", "T.test", "Anova",  "Eight.Anova"))
+
+# Set up your vector of sigmas to loop through
+sigmas <- c(1,2,4,6,8,12,16,24)
+
+# Write your nested for loop:
+for (i in 1:length(sigmas)){ # sigma runs
+  # set up your vectors for your p-values from each monte-carlo run. 
+  reg.p <- numeric(10) # changes each time you use a new sigma
+  t.p <- numeric(10)
+  anova.p <- numeric(10)
+  eightanova.p <- numeric(10)
+  for (j in 1:10){ # monte carlo runs
+    x <- runif(24, 0, 50) # build set of x's
+    e <- rnorm(24, 0, sd = sigmas[i]) # put in your error
+    y <- 10 + 0.4*x + e # write your linear equation
+    df <- cbind(x, y) # bind x and y into a data frame
+    df <- as.data.frame(df)
+    df <- df[order(df$x),] # order for later sorting. 
+    
+    # dummy coding:
+    
+    # For t-test
+    tdummy <- matrix(0, 24, 1)
+    tdummy[13:24,1] <- 1
+  
+    # For anova
+    dummy <- matrix(0, 24, 3)
+    dummy[7:12,1] <- 1
+    dummy[13:18,2] <-1
+    dummy[19:24,3] <-1
+    
+    #for eight-level anova
+    edummy <- matrix(0, 24, 7)
+    edummy[4:6,1]<-1
+    edummy[7:9,2]<-1
+    edummy[10:12,3]<-1
+    edummy[13:15,4]<-1
+    edummy[16:18,5]<-1
+    edummy[19:21,6]<-1
+    edummy[22:24,7]<-1
+    
+    # guesses
+    nullmodguess <- c(1,2) # 2
+    lmmodguess <- c(1, 2, 3) # 3 
+    tmodguess<-c(1,1,1,1,1) # 5
+    aovguess<- c(9,10,-6,8,2) # 5 
+    eightaovguess<-c(1, 12.5,15,17.5,20,22.5,25,27.5,30) # 9
+    
+    #optims
+    nullfitmod=optim(par=nullmodguess,fn=nullmod,x=df$x,y=df$y) # converges
+    lmfitmod=optim(par=lmmodguess,fn=lin,x=df$x,y=df$y) # converges 
+    tfitmod=optim(par=tmodguess,fn=ttest,x=tdummy,y=df$y) # converges
+    anovafitmod=optim(par=aovguess,fn=anovamod,x=dummy,y=df$y, control=list(maxit=1e5)) # converges
+    eightanovafitmod=optim(par=eightaovguess,fn=eightanovamod,x=edummy,y=df$y, control=list(maxit=1e5)) # converges
+   
+    # degrees of freedom 
+     degfreg <- length(lmfitmod$par) - length(nullfitmod$par) 
+     degft <- length(tfitmod$par) - length(nullfitmod$par) 
+     degfanova <- length(anovafitmod$par) - length(nullfitmod$par) 
+     degfeightanova <- length(eightanovafitmod$par) - length(nullfitmod$par) 
+     
+     # test stats
+    t.statreg <- 2*(nullfitmod$value - lmfitmod$value)
+    t.statt <- 2*(nullfitmod$value - tfitmod$value)
+    t.statanova <- 2*(nullfitmod$value - anovafitmod$value)
+    t.stateightanova <- 2*(nullfitmod$value - eightanovafitmod$value)
+    
+    # pvalues 
+    reg.p[j] <- 1-pchisq(t.statreg, df=degfreg)
+    t.p[j]<- 1-pchisq(t.statt, df=degft)
+    anova.p[j]<- 1-pchisq(t.statanova, df=degfanova)
+    eightanova.p[j]<- 1-pchisq(t.stateightanova, df=degfeightanova)
+  }
+  
+  # the mean p values for each of ten monte carlo runs, and each of 8 sigmas into a data frame.
+  dfp[,1]<- c(1, 2, 4, 6, 8, 12, 16, 24)
+  dfp[i,2] <- mean(reg.p)
+  dfp[i,3] <- mean(t.p)
+  dfp[i,4] <- mean(anova.p)
+  dfp[i,5] <- mean(eightanova.p)
+}
+```
+
+Now reshape your dataframe of p-values into something you can plot!
+
+```{r}
+mframe <- melt(dfp)
+mframe <- mframe[9:40,]
+allsigmas <- c(1,  2,  4,  6,  8, 12, 16, 24, 1,  2,  4,  6,  8,12, 16, 24, 1,  2,  4,  6,  8, 12, 16, 24, 1,  2,  4,  6,  8, 12, 16, 24)
+mframe$sigma <- allsigmas
+```
+
+### Step 3: Visualize!!
+
+```{r}
+mframe <- as.data.frame(mframe)
+p <- ggplot(mframe, aes(x=mframe$sigma, y=mframe$value, colour=variable)) + geom_point (size=2) + xlab("Sigma Values") + ylab("Mean p-values") 
+p
+     
+```
+
+### Step 4: Interpret
+
+The regression model performs best across the sigmas, as evidenced by it's low p-value scores (see figure above). We originally believed that the eight level anova (Eight.Anova) should perform close to the regression, when in fact it performed the worst out of all four models. The four level anova (Anova) performs in between the t-test (T.test) and the regression. Finally the t-test, as expected, performed worse than either the regression or the four level ANOVA. This all makes sense, as our data is linearly correlated, so a linear regression should fit best.
+
+Upon talking to Stuart, with the eight level ANOVA, we are starting to make finer and finer comparisons. When there's a linear relationship with noise and you split the x finer, the p-values should be high. With fewer groups, as in the four-way ANOVA, the groups shouldn't overlap as much, and the p-values should be lower.
+
+It's clearly best to try to fit linear data to a linear model, but when resources are constrained, it seems the best strategy is to fit as many non-overlapping groups as possible. In this case, you need at least four experimental units (four way anova), to detect a significant signal within the noise. 
+
+
+```{r, include=FALSE}
+# get values for the eight level anova guess:
+x <- c(6.25, 12.5, 18.72, 25, 31.25, 37.5, 43.75, 50)
+y.vec <- rep(NA, 8)
+for (i in 1:length(x)){
+   y[i] <- 10 + 0.4*x[i]
+   y.vec <- y
+}
+```