---
title: "Viab_Exp2"
output: html_notebook
---
INDEX
-General preparation
-Creating a PROP Data.Frame
-GLM
-Data Analysis
-Pairwise Comparisons
-PLOTS
-TPS plots
-GAMs plots
-Difference plot
# General preparation
```{r}
setwd('C:/Users/tkut0001/Documents/R Scripts')
Viab<-read.table('Viab_Exp2.txt',head=T)
head(Viab)
```
```{r}
# Packages
require("gdata")
require("data.table")
require("ggplot2")
require("plyr")
require("lsmeans")
require("mgcv")
require("nlme")
#require("AICcmodavg")
require("fields")
require("multcompView")
require('cowplot')
require('car')
```
# Creatnig PROP data frame (with PROP as dependent var. - The proprtion of survived- to use in tps and gams)
```{r}
#I need one row per vial, PROP has to be the "nrow" where SURV==1
temp0<-Viab[which(Viab$SURV==1),]
head(temp0)
temp0<-aggregate(SURV~VIAL+PROT+CARB+CONC+RATIO+TEMP, temp0, length)
head(temp0)
#I have 36 diets at 2 temperatures (72treats) and 5 vials per treat, so I should have
5*72
#360 vials in total and temp0 has 1 row per vial
360-nrow(temp0)
#I am missing 2 vials, because they had 0 flies.
temp1<-Viab[which(Viab$SURV==0),] #Data.frame counting the dead flies only
head(temp1)
temp1<-aggregate(SURV~VIAL+PROT+CARB+CONC+RATIO+TEMP, temp1 , length) #number of dead flies in each vial.
temp1[20:30,]
#I want only those rows where all flies died (i.e. temp0$SURV==20)
rm(temp2)
temp2<-temp1[which(temp1$SURV==20),]
head(temp2) # 2 vials
temp2$SURV<-rep(0,nrow(temp2)) #Substitute the SURV (which is the number of rows where a flie died per vial (in this case all of them, so 20)
temp2
```
```{r}
# Join together temp0 -all vials wiht at least 1 fly- with temp2 -all vials with 0 flies-
Prop<-rbind(temp0,temp2)
head(Prop)
nrow(Prop)
#360 rows: 1 per vial
Prop$PROP<-Prop$SURV/20 #creating a new variable called PROP : proportion that survived out of the 20 eggs.
head(Prop)
library(xlsx)
#write.xlsx(Prop, "C:/Users/tkut0001/Documents/Prop_Exp2.xlsx")
```
```{r}
# just double-checking
nrow(Prop)
Prop[which(Prop$PROP==0),]
```
# GLM
## Data Analysis
```{r}
# For binomial I need a column for DIED
head(Viab)
Viab$DIED<-1:nrow(Viab)
Viab$DIED[which(Viab$SURV==1)]=0
Viab$DIED[which(Viab$SURV==0)]=1
#Viab<-cbind(Viab,DIED)
```
```{r}
#Data exploration
head(Viab)
is.data.frame(Viab)
str(Viab)
Viab$TEMP<-factor(Viab$TEMP)
Viab$VIAL<-factor(Viab$VIAL)
hist(Prop$PROP)
```
###Binomial family 25 degrees
```{r}
#Binomial family
#25oC degrees
Viab25.bin<-glm(cbind(SURV,DIED) ~ PROT + CARB + I(CARB^2) + I(PROT^2) + PROT:CARB, family = binomial, data = Viab, subset=TEMP=='25')
summary(Viab25.bin)
#(Dispersion parameter for quasibinomial family taken to be 1.015429. lets do binomial
par(mfrow=c(2,2))
plot(Viab25.bin)
```
###Binomial family 28 degrees
Dispersion parameter for quasibinomial family taken to be 1.03730
```{r}
Viab28.2.bin<-glm(cbind(SURV,DIED) ~ PROT + CARB + I(CARB^2) + I(PROT^2) + PROT:CARB, family = binomial, data = Viab, subset=TEMP=='28')
summary(Viab28.2.bin)
```
```{r}
Viab25.2.bin<-glm(cbind(SURV,DIED) ~ PROT + CARB + I(CARB^2) + I(PROT^2) + PROT:CARB, family = binomial, data = Viab, subset=TEMP=='25')
summary(Viab25.2.bin)
```
```{r}
y<-summary(Viab28.bin2)$coefficients
as.data.frame(y)
```
##Pairwise Comparisons
```{r}
#Comparing "wihtout temperature interactions"" to "with temperature interactions""
Viab.bin <- glm(cbind(SURV,DIED) ~ TEMP + PROT + CARB + I(CARB^2) + I(PROT^2) + PROT:CARB, family = binomial, data = Viab)
Viab.Temp.bin <- glm(cbind(SURV,DIED) ~ TEMP + PROT + CARB + I(CARB^2) + I(PROT^2) + PROT:CARB +
TEMP:PROT + TEMP:CARB + TEMP:PROT:CARB, family = binomial, data = Viab)
anova(Viab.bin, Viab.Temp.bin, test = "Chisq")
```
```{r}
Viab.bin <- glm(cbind(SURV,DIED) ~ TEMP + PROT + CARB + I(CARB^2) + I(PROT^2) + PROT:CARB, family = binomial, data = Viab)
ViabTempOnly.bin <- glm(cbind(SURV,DIED) ~ TEMP, family = binomial, data = Viab)
summary(ViabTempOnly.bin)
Anova(ViabTempOnly.bin, test.statistic = "F")
anova(Viab.bin, ViabTempOnly.bin, test = "Chisq")
anova(ViabTempOnly.bin, test = 'Chisq')
```
```{r}
x<-anova(Viab.bin, Viab.Temp.bin, test = "Chisq")
as.data.frame(x)
```
#PLOTS
##Preparation
```{r}
#Creating a numeric variable for "ratio"
Prop$RatioNo<-NULL
#levels(Prop$RATIO)
RatioNo<-c(1:nrow(Prop))
RatioNo[which(Prop$RATIO=='3:2')]=1.328
RatioNo[which(Prop$RATIO=='2:3')]=1/1.543
RatioNo[which(Prop$RATIO=='1:3')]=1/2.96
RatioNo[which(Prop$RATIO=='1:4')]=1/4
RatioNo[which(Prop$RATIO=='1:8')]=1/8
RatioNo[which(Prop$RATIO=='1:16')]=1/16
levels(factor(RatioNo))
Prop<-cbind(Prop,RatioNo)
head(Prop)
```
```{r}
# THIN PLATE SPLINE PLOT
##functions to generate plot space for heatmap/levelplot-----
genspace<-function(x){
with(x,{
full.space<-expand.grid(CARB=seq(min(CARB),max(CARB),length.out=100),
PROT=seq(min(PROT),max(PROT),length.out=100))
maxconc<-x[CONC==max(CONC),]
minconc<-x[CONC==min(CONC),]
xymax<-with(maxconc,coef(lm(CARB~PROT)))
xymin<-with(minconc,coef(lm(CARB~PROT)))
space.fit<-subset(full.space,(CARB>(xymin[2]*PROT+xymin[1]*0.99)&
CARB<(xymax[2]*PROT+xymax[1]*1.01))&
((PROT/CARB)>min(x$RatioNo*0.99)&
(PROT/CARB)