%load_ext rpy2.ipython

pwd

!cat /Users/sr320/Desktop/1190951/Lab\ experiments/ANOVA_Threshold.r

cd /Users/sr320/Desktop/1190951/Lab\ experiments

%%R
threshold.data <-read.csv('Threshold_data.csv', header = TRUE)
head(threshold.data)

%%R
th.anova<-lm(NormalizedExp~Treatment + Size, data=threshold.data)
summary(th.anova)
anova(th.anova)

pylab inline

%%R 
library(Hmisc) #for errorbars
#first calculate mean, sd, and count
th.mean<-tapply(data1$NormalizedExp,data1$Treatment, mean)
th.sd<-tapply(data1$NormalizedExp,data1$Treatment, sd)
th.count<-tapply(data1$NormalizedExp,data1$Treatment, length)
#calculate SE
th.se<-th.sd/(sqrt(th.count))
#Bind mean and SE for plotting
threshold.plot<-cbind(th.mean,th.se)
threshold.plot
#Data was exported and relative expression and SE were claculated

##Upload plot data
th.plot <-read.csv('Mean_threshold_plot.csv', header = TRUE)
head(th.plot )

yplus<-th.plot$Relative+th.plot$RelativeSE
yminus<-th.plot$Relative

barx<-barplot(th.plot$Relative, las=1, ylab="Relative HIF-1a mRNA Expression",
  ylim=c(0,3),cex.lab=1.1, cex.axis=1.11,col="gray",  axis.lty=1,
  names.arg=c("Control","Moderate","Hypoxia"), 
  cex=1.1)
errbar(barx,th.plot$Relative, yplus,yminus, type='n', add=T)

!cat /Users/sr320/Desktop/1190951/Lab\ experiments/Mean_threshold_plot.csv

%%R
workdir <-"/Users/sr320/Desktop/1190951/Field\ data"
setwd(workdir)

%%R
field.data <-read.csv('All_field_qPCR.csv', header = TRUE)
head(field.data)


%%R
library(lme4)
#Year model
yr.mod<- glm(Elevated ~ factor(Year) + Site + Size + TrawlDepth, family=binomial(link = "logit"), data=field.data)
summary(yr.mod)
anova(yr.mod, test="Chisq")