--- title: "Gene Expression Analysis" author: "Alexander Shingleton" date: "11/11/2020" output: pdf_document editor_options: chunk_output_type: console --- ```{r setup, include=FALSE} knitr::opts_chunk$set(echo = TRUE,tidy.opts=list(width.cutoff=60),tidy=TRUE) ``` # Analysis of the Nutritional Geometry of Gene Expression in Males and Females Preamble ```{r,message=FALSE, warning=FALSE} setwd("~/Expression Analysis") suppressPackageStartupMessages({ require("gdata") require("car") require("data.table") require("ggplot2") require("plyr") require("lsmeans") require("mgcv") require("nlme") require("AICcmodavg") require("fields") require("lme4") require ("piecewiseSEM") require ("gridExtra") require("lmSupport") require("WebPower") require("pwr") require("formatR") }) source("R/functions.R") ``` Import the data and add additional columns specifying nutritional details. ```{r,message=FALSE, warning=FALSE} data <- read.csv("gene expression data log.csv") diets <- read.csv("diets gene expression.csv") all.data<-within(data,{ foodF<-factor(food) inter=interaction(ratio,foodF) logprot=log(prot) logcarb=log(carb) logratio=log(ratio_num) logfood=log(food) }) all.data$cal<-data$food*0.004 m.data<-subset(all.data, sex=="M") f.data<-subset(all.data, sex=="F") ``` ## Multivariate Analysis of Gene Expression Using Lee's model of nutritional geometry. ```{r,message=FALSE, warning=FALSE} gexp<-as.matrix(cbind(all.data[,c(9:16)])) all.manova0<-lm(cbind(X4eBP,InR,dILP2,dILP3,dILP5,dILP8,Ash2L,CG3071)~sex*(prot:carb+poly(carb,2)+poly(prot,2)), data=all.data) all.manova1<-lm(cbind(X4eBP,InR,dILP2,dILP3,dILP5,dILP8,Ash2L,CG3071)~sex+(prot:carb+poly(carb,2)+poly(prot,2)), data=all.data) anova(all.manova0,all.manova1) Man.all.manova0 <- Manova(all.manova0, multivariate=TRUE, type="III") Man.all.manova0 ``` It looks like there is a significant sex:gene expression interaction, but inspecting the factors, the prot:carb and sex:carb^2 terms are not significant. Try a simpler model with a linear carbohydrate term: ```{r,message=FALSE, warning=FALSE} all.manova0<-lm(cbind(X4eBP,InR,dILP2,dILP3,dILP5,dILP8,Ash2L,CG3071)~sex*(prot:carb+poly(carb,2)+poly(prot,2)), data=all.data) all.manova1<-lm(cbind(X4eBP,InR,dILP2,dILP3,dILP5,dILP8,Ash2L,CG3071)~sex*(poly(carb,1)+poly(prot,2)), data=all.data) all.manova2<-lm(cbind(X4eBP,InR,dILP2,dILP3,dILP5,dILP8,Ash2L,CG3071)~sex+(poly(carb,1)+poly(prot,2)), data=all.data) anova(all.manova2,all.manova1,all.manova0) Man.all.manova1 <- Manova(all.manova1, multivariate=TRUE, type="III") Man.all.manova1 ``` It looks like there is a significant quadratic protein and linear carbohydrate interaction with sex ## Gene Expression and Body Size We can also explore how body gene expression of each gene change changes with body size. Because we used different diets for the body size and analysis and gene expression analysis, we will use fitted values to test the relationship. First looking at gene expression and fitted values for overall body size. We can also get the fitted PCA1 for the diets used in the gene expression analysis. First load the body size data and fit the PCA values predicted by the body size model to the expression data diets. Load the body size data ```{r,message=FALSE, warning=FALSE} diets.gene.expression <- read.csv("~/Documents/Data/Nutritional Geometry of Gene Expression/Final Analysis/Expression Analysis/diets gene expression.csv") organsize <- read.csv("~/Documents/Data/Nutritional Geometry of Gene Expression/Final Analysis/Body Size/organsize.csv") ``` Now extract the females data, fit the complete model and add the predicted PCA values to the expression data. ```{r,message=FALSE, warning=FALSE} femsize<-(subset(organsize,sex=="F")) femsize$genital<-NULL femsize<-na.omit(femsize) size.model.lee.replicates<-lmer(PCA1~carb*prot+poly(carb,2)+poly(prot,2)+(1|repfull),data=femsize) GE_fem_size<-cbind(diets.gene.expression,(predict(size.model.lee.replicates,diets.gene.expression, re.form=NA))) names(GE_fem_size)[5]<-"lee" GE_fem_size[,c(3,4)]<-NULL f.data<-merge(GE_fem_size, f.data, by=c("carb","prot")) ``` Now do the same with the thin plate spline ```{r,message=FALSE, warning=FALSE} tps.f<-with(femsize,Tps(cbind(prot,carb),PCA1)) gen.tps.f<-within(diets.gene.expression,{ tps<-predict(tps.f,data.frame(prot,carb)) }) gen.tps.f[,c(3,4)]<-NULL f.data<-merge(gen.tps.f, f.data, by=c("carb","prot")) ``` Now males: ```{r,message=FALSE, warning=FALSE} malsize<-(subset(organsize,sex=="M")) malsize<-na.omit(malsize) size.model.lee.replicates<-lmer(PCA1~poly(prot,2)+(1|repfull),data=malsize) GE_mal_size<-cbind(diets.gene.expression,(predict(size.model.lee.replicates,diets.gene.expression, re.form=NA))) names(GE_mal_size)[5]<-"lee" GE_mal_size[,c(3,4)]<-NULL m.data<-merge(GE_mal_size, m.data, by=c("carb","prot")) tps.m<-with(malsize,Tps(cbind(prot,carb),PCA1)) gen.tps.m<-within(diets.gene.expression,{ tps<-predict(tps.m,data.frame(prot,carb)) }) gen.tps.m[,c(3,4)]<-NULL m.data<-merge(gen.tps.m, m.data, by=c("carb","prot")) ``` The expression of which gene correlates most strongly with body size? ```{r,message=FALSE, warning=FALSE} fem_lm_lee<-lm(lee~X4eBP+ InR + dILP2 + dILP3 + dILP5 + dILP8 + Ash2L + CG3071, data=f.data) Anova(fem_lm_lee) summary(fem_lm_lee) modelEffectSizes(fem_lm_lee) fem_lm_tps<-lm(tps~X4eBP+ InR + dILP2 + dILP3 + dILP5 + dILP8 + Ash2L + CG3071, data=f.data) Anova(fem_lm_tps) modelEffectSizes(fem_lm_tps) mal_lm_lee<-lm(lee~X4eBP+ InR + dILP2 + dILP3 + dILP5 + dILP8 + Ash2L + CG3071, data=m.data) Anova(mal_lm_lee) summary(mal_lm_lee) modelEffectSizes(mal_lm_lee) mal_lm_tps<-lm(tps~X4eBP+ InR + dILP2 + dILP3 + dILP5 + dILP8 + Ash2L + CG3071, data=m.data) Anova(mal_lm_tps) modelEffectSizes(mal_lm_tps) ``` Now lets examine the nutritional geometry of each gene's expression individually. ## Analysis of InR Test whether there is a difference between the sexes ```{r,message=FALSE, warning=FALSE} InRModel1<-lm(InR~sex+(carb:prot+poly(carb,2)+poly(prot,2)), data=all.data) InRModel2<-lm(InR~sex*(carb:prot+poly(carb,2)+poly(prot,2)), data=all.data) anova(InRModel1, InRModel2) ``` No difference. Test the effect of C and P on InR expression in males and females ```{r,message=FALSE, warning=FALSE} InR.m<-lm(InR~(carb:prot+poly(carb,2)+poly(prot,2)), data=m.data) InR.f<-lm(InR~(carb:prot+poly(carb,2)+poly(prot,2)), data=f.data) summary(InR.m) summary(InR.f) ``` There is a quadratic protein effect in males and females Reanalyse the data using the simplest model (the carb effect in males disappears when the carb*protein interaction is dropped, suggesting that it is not well supported) ```{r,message=FALSE, warning=FALSE} #Non-orthogoganol coef summary(lm(InR~(poly(prot,2,raw=TRUE)), data=m.data))[4] summary(lm(InR~(poly(prot,2,raw=TRUE)), data=f.data))[4] #Orthoganol test InR.m<-lm(InR~(poly(prot,2)), data=m.data) InR.f<-lm(InR~(poly(prot,2)), data=f.data) summary(InR.m)[4] summary(InR.f)[4] ``` Re-test whether there is a difference between the sexes ```{r,message=FALSE, warning=FALSE} InRModel0<-lm(InR~(poly(prot,2)), data=all.data) InRModel1<-lm(InR~sex+(poly(prot,2)), data=all.data) InRModel2<-lm(InR~sex*(poly(prot,2)), data=all.data) anova(InRModel0,InRModel1, InRModel2) summary(InRModel1) ``` There is still no effect of sex on InR expression althought it is slightly lower in males than females. Plot the model and the TPS ```{r, echo=FALSE, warning=FALSE} # Generate plot space based on female data. This will be used for all subsequent analyses. The geneexpression space specifies the coverage of the surface. The points give the dietary values used. genexpression.space<-genspace(f.data)[c(-3,-4,-5)] points<-genpoints(f.data) #generate the predictions lee.f<-predict(InR.f,newdata=genexpression.space,se.fit=F) gen.lee.f<-genexpression.space gen.lee.f$lee<-lee.f lee.m<-predict(InR.m,newdata=genexpression.space,se.fit=F) gen.lee.m<-genexpression.space gen.lee.m$lee<-lee.m #generate the TPS data tps.f<-with(f.data,Tps(cbind(prot,carb),InR)) gen.tps.f<-within(genexpression.space,{ tps<-predict(tps.f,data.frame(prot,carb)) }) tps.m<-with(m.data,Tps(cbind(prot,carb),InR)) gen.tps.m<-within(genexpression.space,{ tps<-predict(tps.m,data.frame(prot,carb)) }) #set z-axis attributes to span male and female expression for both model and tps m.min.lee<-min(gen.lee.m$lee) m.max.lee<-max(gen.lee.m$lee) f.min.lee<-min(gen.lee.f$lee) f.max.lee<-max(gen.lee.f$lee) m.min.tps<-min(gen.tps.m$tps) m.max.tps<-max(gen.tps.m$tps) f.min.tps<-min(gen.tps.f$tps) f.max.tps<-max(gen.tps.f$tps) z.range<-c(min(m.min.lee,m.min.tps,f.min.lee,f.min.tps),max(m.max.lee,m.max.tps,f.max.lee,f.max.tps)) #Model Plots gen.plot.lee.f<-ggplot(data=gen.lee.f,aes(x=prot,y=carb,z=lee,fill=lee))+ geom_tile(alpha=.9)+theme_bw()+theme(panel.grid=element_blank(), axis.title.x=element_text(size=16), axis.title.y=element_text(size=16), axis.text.x=element_text(size=14), axis.text.y=element_text(size=14), legend.text=element_text(size=14), legend.title=element_text(size=16))+ scale_x_continuous(limits=c(0,max(diets$prot*1.1)),expand=c(0,0))+ scale_y_continuous(expand=c(0,0),limits=c(0,max(diets$carb*1.1)))+ geom_contour(color="black",bins=10)+ xlab("protein (g/L)")+ylab("carbohydrate (g/L)")+ scale_fill_gradientn(colours=c("blue","yellow","red"), name="InR", limits=z.range)+ geom_abline(data=ratios,aes(slope=1/ratios, intercept=0),linetype=2,size=rel(0.8))+ geom_point(data=diets,aes(x=prot, y=carb, z=1),size=3,fill="black")+ guides(fill=guide_colorbar(order=1)) + ggtitle("female.lee") gen.plot.lee.m<-ggplot(data=gen.lee.m,aes(x=prot,y=carb,z=lee,fill=lee))+ geom_tile(alpha=.9)+theme_bw()+theme(panel.grid=element_blank(), axis.title.x=element_text(size=16), axis.title.y=element_text(size=16), axis.text.x=element_text(size=14), axis.text.y=element_text(size=14), legend.text=element_text(size=14), legend.title=element_text(size=16))+ scale_x_continuous(limits=c(0,max(diets$prot*1.1)),expand=c(0,0))+ scale_y_continuous(expand=c(0,0),limits=c(0,max(diets$carb*1.1)))+ geom_contour(color="black",bins=10)+ xlab("protein (g/L)")+ylab("carbohydrate (g/L)")+ scale_fill_gradientn(colours=c("blue","yellow","red"), name="InR", limits=z.range)+ geom_abline(data=ratios,aes(slope=1/ratios, intercept=0),linetype=2,size=rel(0.8))+ geom_point(data=diets,aes(x=prot, y=carb, z=1),size=3,fill="black")+ guides(fill=guide_colorbar(order=1)) + ggtitle("male.lee") gen.plot.tps.f<-ggplot(data=gen.tps.f,aes(x=prot,y=carb,z=tps,fill=tps))+ geom_tile(alpha=.9)+theme_bw()+theme(panel.grid=element_blank(), axis.title.x=element_text(size=16), axis.title.y=element_text(size=16), axis.text.x=element_text(size=14), axis.text.y=element_text(size=14), legend.text=element_text(size=14), legend.title=element_text(size=16))+ scale_x_continuous(limits=c(0,max(diets$prot*1.1)),expand=c(0,0))+ scale_y_continuous(expand=c(0,0),limits=c(0,max(diets$carb*1.1)))+ geom_contour(color="black",bins=10)+ xlab("protein (g/L)")+ylab("carbohydrate (g/L)")+ scale_fill_gradientn(colours=c("blue","yellow","red"), name="InR", limits=z.range)+ geom_abline(data=ratios,aes(slope=1/ratios, intercept=0),linetype=2,size=rel(0.8))+ geom_point(data=diets,aes(x=prot, y=carb, z=1),size=3,fill="black")+ guides(fill=guide_colorbar(order=1)) + ggtitle("female.tps") gen.plot.tps.m<-ggplot(data=gen.tps.m,aes(x=prot,y=carb,z=tps,fill=tps))+ geom_tile(alpha=.9)+theme_bw()+theme(panel.grid=element_blank(), axis.title.x=element_text(size=16), axis.title.y=element_text(size=16), axis.text.x=element_text(size=14), axis.text.y=element_text(size=14), legend.text=element_text(size=14), legend.title=element_text(size=16))+ scale_x_continuous(limits=c(0,max(diets$prot*1.1)),expand=c(0,0))+ scale_y_continuous(expand=c(0,0),limits=c(0,max(diets$carb*1.1)))+ geom_contour(color="black",bins=10)+ xlab("protein (g/L)")+ylab("carbohydrate (g/L)")+ scale_fill_gradientn(colours=c("blue","yellow","red"), name="InR", limits=z.range)+ geom_abline(data=ratios,aes(slope=1/ratios, intercept=0),linetype=2,size=rel(0.8))+ geom_point(data=diets,aes(x=prot, y=carb, z=1),size=3,fill="black")+ guides(fill=guide_colorbar(order=1)) + ggtitle("male.tps") # Print plots gen.plot.lee.f gen.plot.lee.m gen.plot.tps.f gen.plot.tps.m ``` ##Analysis of 4eBP Test whether there is a difference between the sexes ```{r,message=FALSE, warning=FALSE} X4eBPModel1<-lm(X4eBP~sex+(carb:prot+poly(carb,2)+poly(prot,2)), data=all.data) X4eBPModel2<-lm(X4eBP~sex*(carb:prot+poly(carb,2)+poly(prot,2)), data=all.data) anova(X4eBPModel1, X4eBPModel2) summary(X4eBPModel2) ``` There is a sex*diet interaction, but none of the indiviual parameters in the full model are significant apart from prot^2. Test the effect of C and P on 4EBP expression in males and females ```{r,message=FALSE, warning=FALSE} X4eBP.m<-lm(X4eBP~(carb:prot+poly(carb,2)+poly(prot,2)), data=m.data) X4eBP.f<-lm(X4eBP~(carb:prot+poly(carb,2)+poly(prot,2)), data=f.data) summary(X4eBP.m) summary(X4eBP.f) ``` There is no effect of diet on x4eBP expression except prot in females. Try a simpler model. ```{r,message=FALSE, warning=FALSE} X4eBP.m<-lm(X4eBP~prot, data=m.data) X4eBP.f<-lm(X4eBP~prot, data=f.data) summary(X4eBP.m) summary(X4eBP.f) ``` There is an effect in females but not in males. Compare sex effect again. ```{r,message=FALSE, warning=FALSE} X4eBPModel1<-lm(X4eBP~sex+prot, data=all.data) X4eBPModel2<-lm(X4eBP~sex*prot, data=all.data) anova(X4eBPModel1, X4eBPModel2) summary(X4eBPModel2) ``` There is a significant difference between sexes in the response to protein. Males appear to have lower expression in general, but there is an interaction with protein, such that there is no effect on expression in males. Plot the model and the TPS ```{r echo=FALSE, warning=FALSE} # Generate plot space based on female data. This will be used for all subsequent analyses. The geneexpression space specifies the coverage of the surface. The points give the dietary values used. genexpression.space<-genspace(f.data)[c(-3,-4,-5)] points<-genpoints(f.data) #generate the predictions #lee.f<-predict(X4eBP.f,newdata=genexpression.space,se.fit=F) lee.f<-predict(X4eBP.f,newdata=genexpression.space,se.fit=F) gen.lee.f<-genexpression.space gen.lee.f$lee<-lee.f #lee.m<-predict(X4eBP.m,newdata=genexpression.space,se.fit=F) gen.lee.m<-genexpression.space gen.lee.m$lee<-mean(m.data$X4eBP, na.rm=TRUE) #generate the TPS data tps.f<-with(f.data,Tps(cbind(prot,carb),X4eBP)) gen.tps.f<-within(genexpression.space,{ tps<-predict(tps.f,data.frame(prot,carb)) }) tps.m<-with(m.data,Tps(cbind(prot,carb),X4eBP)) gen.tps.m<-within(genexpression.space,{ tps<-predict(tps.m,data.frame(prot,carb)) }) #set z-axis attributes to span male and female expression for both model and tps m.min.lee<-min(gen.lee.m$lee) m.max.lee<-max(gen.lee.m$lee) f.min.lee<-min(gen.lee.f$lee) f.max.lee<-max(gen.lee.f$lee) m.min.tps<-min(gen.tps.m$tps) m.max.tps<-max(gen.tps.m$tps) f.min.tps<-min(gen.tps.f$tps) f.max.tps<-max(gen.tps.f$tps) z.range<-c(min(m.min.lee,m.min.tps,f.min.lee,f.min.tps),max(m.max.lee,m.max.tps,f.max.lee,f.max.tps)) #Model Plots gen.plot.lee.f<-ggplot(data=gen.lee.f,aes(x=prot,y=carb,z=lee,fill=lee))+ geom_tile(alpha=.9)+theme_bw()+theme(panel.grid=element_blank(), axis.title.x=element_text(size=16), axis.title.y=element_text(size=16), axis.text.x=element_text(size=14), axis.text.y=element_text(size=14), legend.text=element_text(size=14), legend.title=element_text(size=16))+ scale_x_continuous(limits=c(0,max(diets$prot*1.1)),expand=c(0,0))+ scale_y_continuous(expand=c(0,0),limits=c(0,max(diets$carb*1.1)))+ geom_contour(color="black",bins=10)+ xlab("protein (g/L)")+ylab("carbohydrate (g/L)")+ scale_fill_gradientn(colours=c("blue","yellow","red"), name="4EBP", limits=z.range)+ geom_abline(data=ratios,aes(slope=1/ratios, intercept=0),linetype=2,size=rel(0.8))+ geom_point(data=diets,aes(x=prot, y=carb, z=1),size=3,fill="black")+ guides(fill=guide_colorbar(order=1)) + ggtitle("female.lee") gen.plot.lee.m<-ggplot(data=gen.lee.m,aes(x=prot,y=carb,z=lee,fill=lee))+ geom_tile(alpha=.9)+theme_bw()+theme(panel.grid=element_blank(), axis.title.x=element_text(size=16), axis.title.y=element_text(size=16), axis.text.x=element_text(size=14), axis.text.y=element_text(size=14), legend.text=element_text(size=14), legend.title=element_text(size=16))+ scale_x_continuous(limits=c(0,max(diets$prot*1.1)),expand=c(0,0))+ scale_y_continuous(expand=c(0,0),limits=c(0,max(diets$carb*1.1)))+ geom_contour(color="black",bins=10)+ xlab("protein (g/L)")+ylab("carbohydrate (g/L)")+ scale_fill_gradientn(colours=c("blue","yellow","red"), name="4EBP", limits=z.range)+ geom_abline(data=ratios,aes(slope=1/ratios, intercept=0),linetype=2,size=rel(0.8))+ geom_point(data=diets,aes(x=prot, y=carb, z=1),size=3,fill="black")+ guides(fill=guide_colorbar(order=1)) + ggtitle("male.lee") gen.plot.tps.f<-ggplot(data=gen.tps.f,aes(x=prot,y=carb,z=tps,fill=tps))+ geom_tile(alpha=.9)+theme_bw()+theme(panel.grid=element_blank(), axis.title.x=element_text(size=16), axis.title.y=element_text(size=16), axis.text.x=element_text(size=14), axis.text.y=element_text(size=14), legend.text=element_text(size=14), legend.title=element_text(size=16))+ scale_x_continuous(limits=c(0,max(diets$prot*1.1)),expand=c(0,0))+ scale_y_continuous(expand=c(0,0),limits=c(0,max(diets$carb*1.1)))+ geom_contour(color="black",bins=10)+ xlab("protein (g/L)")+ylab("carbohydrate (g/L)")+ scale_fill_gradientn(colours=c("blue","yellow","red"), name="4EBP", limits=z.range)+ geom_abline(data=ratios,aes(slope=1/ratios, intercept=0),linetype=2,size=rel(0.8))+ geom_point(data=diets,aes(x=prot, y=carb, z=1),size=3,fill="black")+ guides(fill=guide_colorbar(order=1)) + ggtitle("female.tps") gen.plot.tps.m<-ggplot(data=gen.tps.m,aes(x=prot,y=carb,z=tps,fill=tps))+ geom_tile(alpha=.9)+theme_bw()+theme(panel.grid=element_blank(), axis.title.x=element_text(size=16), axis.title.y=element_text(size=16), axis.text.x=element_text(size=14), axis.text.y=element_text(size=14), legend.text=element_text(size=14), legend.title=element_text(size=16))+ scale_x_continuous(limits=c(0,max(diets$prot*1.1)),expand=c(0,0))+ scale_y_continuous(expand=c(0,0),limits=c(0,max(diets$carb*1.1)))+ geom_contour(color="black",bins=10)+ xlab("protein (g/L)")+ylab("carbohydrate (g/L)")+ scale_fill_gradientn(colours=c("blue","yellow","red"), name="4EBP", limits=z.range)+ geom_abline(data=ratios,aes(slope=1/ratios, intercept=0),linetype=2,size=rel(0.8))+ geom_point(data=diets,aes(x=prot, y=carb, z=1),size=3,fill="black")+ guides(fill=guide_colorbar(order=1)) + ggtitle("male.tps") # Print plots gen.plot.lee.f gen.plot.lee.m gen.plot.tps.f gen.plot.tps.m ``` ##Analysis of Ash2L Test whether there is a difference between the sexes ```{r,message=FALSE, warning=FALSE} Ash2LModel1<-lm(Ash2L~sex+(carb:prot+poly(carb,2)+poly(prot,2)), data=all.data) Ash2LModel2<-lm(Ash2L~sex*(carb:prot+poly(carb,2)+poly(prot,2)), data=all.data) anova(Ash2LModel1, Ash2LModel2) ``` There is a significant sex by diet interaction Test the effect of C and P on Ash2L expression in males and females ```{r,message=FALSE, warning=FALSE} Ash2L.m<-lm(Ash2L~(carb:prot+poly(carb,2)+poly(prot,2)), data=m.data) Ash2L.f<-lm(Ash2L~(carb:prot+poly(carb,2)+poly(prot,2)), data=f.data) summary(Ash2L.m) summary(Ash2L.f) ``` There is an effect of protein in females but not in males. The quadratic is marginally significant in females. Re-test the effect of P on expression in males and females using the simpler models ```{r,message=FALSE, warning=FALSE} #Non-Orthoganol Coef summary(lm(Ash2L~(poly(carb,1)+poly(prot,2, raw=TRUE)), data=f.data))[4] #Orthoganol Test Ash2L.m<-lm(Ash2L~(poly(carb,1)), data=m.data) Ash2L.f<-lm(Ash2L~(poly(carb,1)+poly(prot,2)), data=f.data) summary(Ash2L.m) summary(Ash2L.f) ``` There is a carb effect in both males and female and a quadratic protein effect in females Re-test the effect of sex using the simplest model ```{r,message=FALSE, warning=FALSE} Ash2LModel1<-lm(Ash2L~sex*(poly(carb,1)+poly(prot,2)), data=all.data) Ash2LModel1_5<-lm(Ash2L~sex*(carb), data=all.data) Ash2LModel2<-lm(Ash2L~sex+(carb), data=all.data) anova(Ash2LModel2,Ash2LModel1_5,Ash2LModel1) summary(Ash2LModel1_5) ``` There is a significiant difference in response between males and females, but it lies in the response to carbohydrates Plot the model and the TPS ```{r, echo=FALSE, warning=FALSE} #Generate plot space based on female data. This will be used for all subsequent analyses. The geneexpression space specifies the coverage of the surface. The points give the dietary values used. genexpression.space<-genspace(f.data)[c(-3,-4,-5)] points<-genpoints(f.data) #generate the predictions lee.f<-predict(Ash2L.f,newdata=genexpression.space,se.fit=F) gen.lee.f<-genexpression.space gen.lee.f$lee<-lee.f lee.m<-predict(Ash2L.m,newdata=genexpression.space,se.fit=F) gen.lee.m<-genexpression.space gen.lee.m$lee<-lee.m #generate the TPS data tps.f<-with(f.data,Tps(cbind(prot,carb),Ash2L)) gen.tps.f<-within(genexpression.space,{ tps<-predict(tps.f,data.frame(prot,carb)) }) tps.m<-with(m.data,Tps(cbind(prot,carb),Ash2L)) gen.tps.m<-within(genexpression.space,{ tps<-predict(tps.m,data.frame(prot,carb)) }) #set z-axis attributes to span male and female expression for both model and tps m.min.lee<-min(gen.lee.m$lee) m.max.lee<-max(gen.lee.m$lee) f.min.lee<-min(gen.lee.f$lee) f.max.lee<-max(gen.lee.f$lee) m.min.tps<-min(gen.tps.m$tps) m.max.tps<-max(gen.tps.m$tps) f.min.tps<-min(gen.tps.f$tps) f.max.tps<-max(gen.tps.f$tps) z.range<-c(min(m.min.lee,m.min.tps,f.min.lee,f.min.tps),max(m.max.lee,m.max.tps,f.max.lee,f.max.tps)) #Model Plots gen.plot.lee.f<-ggplot(data=gen.lee.f,aes(x=prot,y=carb,z=lee.f,fill=lee))+ geom_tile(alpha=.9)+theme_bw()+theme(panel.grid=element_blank(), axis.title.x=element_text(size=16), axis.title.y=element_text(size=16), axis.text.x=element_text(size=14), axis.text.y=element_text(size=14), legend.text=element_text(size=14), legend.title=element_text(size=16))+ scale_x_continuous(limits=c(0,max(diets$prot*1.1)),expand=c(0,0))+ scale_y_continuous(expand=c(0,0),limits=c(0,max(diets$carb*1.1)))+ geom_contour(color="black",bins=10)+ xlab("protein (g/L)")+ylab("carbohydrate (g/L)")+ scale_fill_gradientn(colours=c("blue","yellow","red"), name="Ash2L", limits=z.range)+ geom_abline(data=ratios,aes(slope=1/ratios, intercept=0),linetype=2,size=rel(0.8))+ geom_point(data=diets,aes(x=prot, y=carb, z=1),size=3,fill="black")+ guides(fill=guide_colorbar(order=1)) + ggtitle("female.lee") gen.plot.lee.m<-ggplot(data=gen.lee.m,aes(x=prot,y=carb,z=lee,fill=lee))+ geom_tile(alpha=.9)+theme_bw()+theme(panel.grid=element_blank(), axis.title.x=element_text(size=16), axis.title.y=element_text(size=16), axis.text.x=element_text(size=14), axis.text.y=element_text(size=14), legend.text=element_text(size=14), legend.title=element_text(size=16))+ scale_x_continuous(limits=c(0,max(diets$prot*1.1)),expand=c(0,0))+ scale_y_continuous(expand=c(0,0),limits=c(0,max(diets$carb*1.1)))+ geom_contour(color="black",bins=10)+ xlab("protein (g/L)")+ylab("carbohydrate (g/L)")+ scale_fill_gradientn(colours=c("blue","yellow","red"), name="Ash2L", limits=z.range)+ geom_abline(data=ratios,aes(slope=1/ratios, intercept=0),linetype=2,size=rel(0.8))+ geom_point(data=diets,aes(x=prot, y=carb, z=1),size=3,fill="black")+ guides(fill=guide_colorbar(order=1)) + ggtitle("male.lee") gen.plot.tps.f<-ggplot(data=gen.tps.f,aes(x=prot,y=carb,z=tps,fill=tps))+ geom_tile(alpha=.9)+theme_bw()+theme(panel.grid=element_blank(), axis.title.x=element_text(size=16), axis.title.y=element_text(size=16), axis.text.x=element_text(size=14), axis.text.y=element_text(size=14), legend.text=element_text(size=14), legend.title=element_text(size=16))+ scale_x_continuous(limits=c(0,max(diets$prot*1.1)),expand=c(0,0))+ scale_y_continuous(expand=c(0,0),limits=c(0,max(diets$carb*1.1)))+ geom_contour(color="black",bins=10)+ xlab("protein (g/L)")+ylab("carbohydrate (g/L)")+ scale_fill_gradientn(colours=c("blue","yellow","red"), name="Ash2L", limits=z.range)+ geom_abline(data=ratios,aes(slope=1/ratios, intercept=0),linetype=2,size=rel(0.8))+ geom_point(data=diets,aes(x=prot, y=carb, z=1),size=3,fill="black")+ guides(fill=guide_colorbar(order=1)) + ggtitle("female.tps") gen.plot.tps.m<-ggplot(data=gen.tps.m,aes(x=prot,y=carb,z=tps,fill=tps))+ geom_tile(alpha=.9)+theme_bw()+theme(panel.grid=element_blank(), axis.title.x=element_text(size=16), axis.title.y=element_text(size=16), axis.text.x=element_text(size=14), axis.text.y=element_text(size=14), legend.text=element_text(size=14), legend.title=element_text(size=16))+ scale_x_continuous(limits=c(0,max(diets$prot*1.1)),expand=c(0,0))+ scale_y_continuous(expand=c(0,0),limits=c(0,max(diets$carb*1.1)))+ geom_contour(color="black",bins=10)+ xlab("protein (g/L)")+ylab("carbohydrate (g/L)")+ scale_fill_gradientn(colours=c("blue","yellow","red"), name="Ash2L", limits=z.range)+ geom_abline(data=ratios,aes(slope=1/ratios, intercept=0),linetype=2,size=rel(0.8))+ geom_point(data=diets,aes(x=prot, y=carb, z=1),size=3,fill="black")+ guides(fill=guide_colorbar(order=1)) + ggtitle("male.tps") # Print plots gen.plot.lee.f gen.plot.lee.m gen.plot.tps.f gen.plot.tps.m ``` ## Analysis of CG3071 Test whether there is a difference between the sexes ```{r,message=FALSE, warning=FALSE} CG3071Model1<-lm(CG3071~sex+(carb:prot+poly(carb,2)+poly(prot,2)), data=all.data) CG3071Model2<-lm(CG3071~sex*(carb:prot+poly(carb,2)+poly(prot,2)), data=all.data) anova(CG3071Model1, CG3071Model2) summary(CG3071Model2) ``` There is a sex*diet interaction. It's complex, but males have (generally) a lower expression level, although you can really interpret this because of the higher order interactions. Test the effect of C and P on CG3071 expression in males and females ```{r,message=FALSE, warning=FALSE} #Non-orthoganol Coef summary(lm(CG3071~(carb:prot+poly(carb,2, raw=TRUE)+poly(prot,2, raw=TRUE)), data=f.data))[4] #Orthoganol Test CG3071.m<-lm(CG3071~(carb:prot+poly(carb,2)+poly(prot,2)), data=m.data) CG3071.f<-lm(CG3071~(carb:prot+poly(carb,2)+poly(prot,2)), data=f.data) summary(CG3071.m) summary(CG3071.f) ``` There is an effect of protein and carbohydrates in females but not in males Plot the model and the TPS ```{r, echo=FALSE, warning=FALSE} # Generate plot space based on female data. This will be used for all subsequent analyses. The geneexpression space specifies the coverage of the surface. The points give the dietary values used. genexpression.space<-genspace(f.data)[c(-3,-4,-5)] points<-genpoints(f.data) #generate the predictions lee.f<-predict(CG3071.f,newdata=genexpression.space,se.fit=F) gen.lee.f<-genexpression.space gen.lee.f$lee<-lee.f gen.lee.m<-genexpression.space gen.lee.m$lee<-mean(m.data$CG3071, na.rm=TRUE) #generate the TPS data tps.f<-with(f.data,Tps(cbind(prot,carb),CG3071)) gen.tps.f<-within(genexpression.space,{ tps<-predict(tps.f,data.frame(prot,carb)) }) tps.m<-with(m.data,Tps(cbind(prot,carb),CG3071)) gen.tps.m<-within(genexpression.space,{ tps<-predict(tps.m,data.frame(prot,carb)) }) #set z-axis attributes to span male and female expression for both model and tps m.min.lee<-min(gen.lee.m$lee) m.max.lee<-max(gen.lee.m$lee) f.min.lee<-min(gen.lee.f$lee) f.max.lee<-max(gen.lee.f$lee) m.min.tps<-min(gen.tps.m$tps) m.max.tps<-max(gen.tps.m$tps) f.min.tps<-min(gen.tps.f$tps) f.max.tps<-max(gen.tps.f$tps) z.range<-c(min(m.min.lee,m.min.tps,f.min.lee,f.min.tps),max(m.max.lee,m.max.tps,f.max.lee,f.max.tps)) #Model Plots gen.plot.lee.f<-ggplot(data=gen.lee.f,aes(x=prot,y=carb,z=lee.f,fill=lee))+ geom_tile(alpha=.9)+theme_bw()+theme(panel.grid=element_blank(), axis.title.x=element_text(size=16), axis.title.y=element_text(size=16), axis.text.x=element_text(size=14), axis.text.y=element_text(size=14), legend.text=element_text(size=14), legend.title=element_text(size=16))+ scale_x_continuous(limits=c(0,max(diets$prot*1.1)),expand=c(0,0))+ scale_y_continuous(expand=c(0,0),limits=c(0,max(diets$carb*1.1)))+ geom_contour(color="black",bins=10)+ xlab("protein (g/L)")+ylab("carbohydrate (g/L)")+ scale_fill_gradientn(colours=c("blue","yellow","red"), name="CG3071", limits=z.range)+ geom_abline(data=ratios,aes(slope=1/ratios, intercept=0),linetype=2,size=rel(0.8))+ geom_point(data=diets,aes(x=prot, y=carb, z=1),size=3,fill="black")+ guides(fill=guide_colorbar(order=1)) + ggtitle("female.lee") gen.plot.lee.m<-ggplot(data=gen.lee.m,aes(x=prot,y=carb,z=lee,fill=lee))+ geom_tile(alpha=.9)+theme_bw()+theme(panel.grid=element_blank(), axis.title.x=element_text(size=16), axis.title.y=element_text(size=16), axis.text.x=element_text(size=14), axis.text.y=element_text(size=14), legend.text=element_text(size=14), legend.title=element_text(size=16))+ scale_x_continuous(limits=c(0,max(diets$prot*1.1)),expand=c(0,0))+ scale_y_continuous(expand=c(0,0),limits=c(0,max(diets$carb*1.1)))+ geom_contour(color="black",bins=10)+ xlab("protein (g/L)")+ylab("carbohydrate (g/L)")+ scale_fill_gradientn(colours=c("blue","yellow","red"), name="CG3071", limits=z.range)+ geom_abline(data=ratios,aes(slope=1/ratios, intercept=0),linetype=2,size=rel(0.8))+ geom_point(data=diets,aes(x=prot, y=carb, z=1),size=3,fill="black")+ guides(fill=guide_colorbar(order=1)) + ggtitle("male.lee") gen.plot.tps.f<-ggplot(data=gen.tps.f,aes(x=prot,y=carb,z=tps,fill=tps))+ geom_tile(alpha=.9)+theme_bw()+theme(panel.grid=element_blank(), axis.title.x=element_text(size=16), axis.title.y=element_text(size=16), axis.text.x=element_text(size=14), axis.text.y=element_text(size=14), legend.text=element_text(size=14), legend.title=element_text(size=16))+ scale_x_continuous(limits=c(0,max(diets$prot*1.1)),expand=c(0,0))+ scale_y_continuous(expand=c(0,0),limits=c(0,max(diets$carb*1.1)))+ geom_contour(color="black",bins=10)+ xlab("protein (g/L)")+ylab("carbohydrate (g/L)")+ scale_fill_gradientn(colours=c("blue","yellow","red"), name="CG3071", limits=z.range)+ geom_abline(data=ratios,aes(slope=1/ratios, intercept=0),linetype=2,size=rel(0.8))+ geom_point(data=diets,aes(x=prot, y=carb, z=1),size=3,fill="black")+ guides(fill=guide_colorbar(order=1)) + ggtitle("female.tps") gen.plot.tps.m<-ggplot(data=gen.tps.m,aes(x=prot,y=carb,z=tps,fill=tps))+ geom_tile(alpha=.9)+theme_bw()+theme(panel.grid=element_blank(), axis.title.x=element_text(size=16), axis.title.y=element_text(size=16), axis.text.x=element_text(size=14), axis.text.y=element_text(size=14), legend.text=element_text(size=14), legend.title=element_text(size=16))+ scale_x_continuous(limits=c(0,max(diets$prot*1.1)),expand=c(0,0))+ scale_y_continuous(expand=c(0,0),limits=c(0,max(diets$carb*1.1)))+ geom_contour(color="black",bins=10)+ xlab("protein (g/L)")+ylab("carbohydrate (g/L)")+ scale_fill_gradientn(colours=c("blue","yellow","red"), name="CG3071", limits=z.range)+ geom_abline(data=ratios,aes(slope=1/ratios, intercept=0),linetype=2,size=rel(0.8))+ geom_point(data=diets,aes(x=prot, y=carb, z=1),size=3,fill="black")+ guides(fill=guide_colorbar(order=1)) + ggtitle("male.tps") # Print plots gen.plot.lee.f gen.plot.lee.m gen.plot.tps.f gen.plot.tps.m ``` ## Analysis of dILP2 Test whether there is a difference between the sexes ```{r,message=FALSE, warning=FALSE} dILP2Model1<-lm(dILP2~sex+(carb:prot+poly(carb,2)+poly(prot,1)), data=all.data) dILP2Model2<-lm(dILP2~sex*(carb:prot+poly(carb,2)+poly(prot,1)), data=all.data) anova(dILP2Model1, dILP2Model2) summary(dILP2Model1) ``` There is no sex*diet interaction, and males and females have the same expression level Test the effect of C and P on dILP2 expression in males and females. ```{r,message=FALSE, warning=FALSE} dILP2.m<-lm(dILP2~(carb:prot+poly(carb,2)+poly(prot,2)),data=m.data) dILP2.f<-lm(dILP2~(carb:prot+poly(carb,2)+poly(prot,2)), data=f.data) summary(dILP2.m) summary(dILP2.f) ``` There is no effect in males but there is an effect in females, although the quadratic effect in females is not significant. Retest with simpler ```{r,message=FALSE, warning=FALSE} #Non-orthoganol Ceof summary(lm(dILP2~(carb:prot+poly(carb,2, raw=TRUE)+poly(prot,1)), data=f.data))[4] #Orthoganol Test dILP2.m<-lm(dILP2~(carb:prot+poly(carb,2)+poly(prot,1)),data=m.data) dILP2.f<-lm(dILP2~(carb:prot+poly(carb,2)+poly(prot,1)), data=f.data) summary(dILP2.m) summary(dILP2.f) ``` Re-testing sex interaction with the simpler model ```{r,message=FALSE, warning=FALSE} dILP2Model0<-lm(dILP2~(carb:prot+poly(carb,2)+poly(prot,1)), data=all.data) dILP2Model1<-lm(dILP2~sex+(carb:prot+poly(carb,2)+poly(prot,1)), data=all.data) dILP2Model2<-lm(dILP2~sex*(carb:prot+poly(carb,2)+poly(prot,1)), data=all.data) anova(dILP2Model0, dILP2Model1, dILP2Model2) summary(anova(dILP2Model1, dILP2Model2)) summary(dILP2Model1) ``` We can test for the power to detect a sex-by-diet interaction ```{r,message=FALSE, warning=FALSE} RsqFull<-summary(dILP2Model2)$r.squared RsqRed<-summary(dILP2Model1)$r.squared fsq<- (RsqFull-RsqRed)/(1-RsqFull) predFull<-length(dILP2Model2$coefficients) predRed<-length(dILP2Model1$coefficients) wp.regression(n = nobs(dILP2Model2), p1= predFull, p2=predRed, f2 = fsq) ``` Plot the model and the TPS ```{r, echo=FALSE, warning=FALSE} # Generate plot space based on female data. This will be used for all subsequent analyses. The geneexpression space specifies the coverage of the surface. The points give the dietary values used. genexpression.space<-genspace(f.data)[c(-3,-4,-5)] points<-genpoints(f.data) #generate the predictions lee.f<-predict(dILP2.f,newdata=genexpression.space,se.fit=F) gen.lee.f<-genexpression.space gen.lee.f$lee<-lee.f gen.lee.m<-genexpression.space gen.lee.m$lee<-mean(m.data$dILP2, na.rm=TRUE) #generate the TPS data tps.f<-with(f.data,Tps(cbind(prot,carb),dILP2)) gen.tps.f<-within(genexpression.space,{ tps<-predict(tps.f,data.frame(prot,carb)) }) tps.m<-with(m.data,Tps(cbind(prot,carb),dILP2)) gen.tps.m<-within(genexpression.space,{ tps<-predict(tps.m,data.frame(prot,carb)) }) #set z-axis attributes to span male and female expression for both model and tps m.min.lee<-min(gen.lee.m$lee) m.max.lee<-max(gen.lee.m$lee) f.min.lee<-min(gen.lee.f$lee) f.max.lee<-max(gen.lee.f$lee) m.min.tps<-min(gen.tps.m$tps) m.max.tps<-max(gen.tps.m$tps) f.min.tps<-min(gen.tps.f$tps) f.max.tps<-max(gen.tps.f$tps) z.range<-c(min(m.min.lee,m.min.tps,f.min.lee,f.min.tps),max(m.max.lee,m.max.tps,f.max.lee,f.max.tps)) #Model Plots gen.plot.lee.f<-ggplot(data=gen.lee.f,aes(x=prot,y=carb,z=lee.f,fill=lee))+ geom_tile(alpha=.9)+theme_bw()+theme(panel.grid=element_blank(), axis.title.x=element_text(size=16), axis.title.y=element_text(size=16), axis.text.x=element_text(size=14), axis.text.y=element_text(size=14), legend.text=element_text(size=14), legend.title=element_text(size=16))+ scale_x_continuous(limits=c(0,max(diets$prot*1.1)),expand=c(0,0))+ scale_y_continuous(expand=c(0,0),limits=c(0,max(diets$carb*1.1)))+ geom_contour(color="black",bins=10)+ xlab("protein (g/L)")+ylab("carbohydrate (g/L)")+ scale_fill_gradientn(colours=c("blue","yellow","red"), name="dILP2", limits=z.range)+ geom_abline(data=ratios,aes(slope=1/ratios, intercept=0),linetype=2,size=rel(0.8))+ geom_point(data=diets,aes(x=prot, y=carb, z=1),size=3,fill="black")+ guides(fill=guide_colorbar(order=1)) + ggtitle("female.lee") gen.plot.lee.m<-ggplot(data=gen.lee.m,aes(x=prot,y=carb,z=lee,fill=lee))+ geom_tile(alpha=.9)+theme_bw()+theme(panel.grid=element_blank(), axis.title.x=element_text(size=16), axis.title.y=element_text(size=16), axis.text.x=element_text(size=14), axis.text.y=element_text(size=14), legend.text=element_text(size=14), legend.title=element_text(size=16))+ scale_x_continuous(limits=c(0,max(diets$prot*1.1)),expand=c(0,0))+ scale_y_continuous(expand=c(0,0),limits=c(0,max(diets$carb*1.1)))+ geom_contour(color="black",bins=10)+ xlab("protein (g/L)")+ylab("carbohydrate (g/L)")+ scale_fill_gradientn(colours=c("blue","yellow","red"), name="dILP2", limits=z.range)+ geom_abline(data=ratios,aes(slope=1/ratios, intercept=0),linetype=2,size=rel(0.8))+ geom_point(data=diets,aes(x=prot, y=carb, z=1),size=3,fill="black")+ guides(fill=guide_colorbar(order=1)) + ggtitle("male.lee") gen.plot.tps.f<-ggplot(data=gen.tps.f,aes(x=prot,y=carb,z=tps,fill=tps))+ geom_tile(alpha=.9)+theme_bw()+theme(panel.grid=element_blank(), axis.title.x=element_text(size=16), axis.title.y=element_text(size=16), axis.text.x=element_text(size=14), axis.text.y=element_text(size=14), legend.text=element_text(size=14), legend.title=element_text(size=16))+ scale_x_continuous(limits=c(0,max(diets$prot*1.1)),expand=c(0,0))+ scale_y_continuous(expand=c(0,0),limits=c(0,max(diets$carb*1.1)))+ geom_contour(color="black",bins=10)+ xlab("protein (g/L)")+ylab("carbohydrate (g/L)")+ scale_fill_gradientn(colours=c("blue","yellow","red"), name="dILP2", limits=z.range)+ geom_abline(data=ratios,aes(slope=1/ratios, intercept=0),linetype=2,size=rel(0.8))+ geom_point(data=diets,aes(x=prot, y=carb, z=1),size=3,fill="black")+ guides(fill=guide_colorbar(order=1)) + ggtitle("female.tps") gen.plot.tps.m<-ggplot(data=gen.tps.m,aes(x=prot,y=carb,z=tps,fill=tps))+ geom_tile(alpha=.9)+theme_bw()+theme(panel.grid=element_blank(), axis.title.x=element_text(size=16), axis.title.y=element_text(size=16), axis.text.x=element_text(size=14), axis.text.y=element_text(size=14), legend.text=element_text(size=14), legend.title=element_text(size=16))+ scale_x_continuous(limits=c(0,max(diets$prot*1.1)),expand=c(0,0))+ scale_y_continuous(expand=c(0,0),limits=c(0,max(diets$carb*1.1)))+ geom_contour(color="black",bins=10)+ xlab("protein (g/L)")+ylab("carbohydrate (g/L)")+ scale_fill_gradientn(colours=c("blue","yellow","red"), name="dILP2", limits=z.range)+ geom_abline(data=ratios,aes(slope=1/ratios, intercept=0),linetype=2,size=rel(0.8))+ geom_point(data=diets,aes(x=prot, y=carb, z=1),size=3,fill="black")+ guides(fill=guide_colorbar(order=1)) + ggtitle("male.tps") # Print plots gen.plot.lee.f gen.plot.lee.m gen.plot.tps.f gen.plot.tps.m ``` ## Analysis of dILP3 Test whether there is a difference between the sexes ```{r,message=FALSE, warning=FALSE} dILP3Model1<-lm(dILP3~sex+(carb:prot+poly(carb,2)+poly(prot,2)), data=all.data) dILP3Model2<-lm(dILP3~sex*(carb:prot+poly(carb,2)+poly(prot,2)), data=all.data) anova(dILP3Model1, dILP3Model2) summary(dILP3Model1) ``` There is no sex*diet interaction, and expression level is the same in both sexes Test the effect of C and P on dILP3 expression in males and females ```{r,message=FALSE, warning=FALSE} dILP3.m<-lm(dILP3~(carb:prot+poly(carb,2)+poly(prot,2)), data=m.data) dILP3.f<-lm(dILP3~(carb:prot+poly(carb,2)+poly(prot,2)), data=f.data) summary(dILP3.m) summary(dILP3.f) ``` Simplify the models: ```{r,message=FALSE, warning=FALSE} dILP3.m<-lm(dILP3~prot, data=m.data) dILP3.f<-lm(dILP3~prot, data=f.data) summary(dILP3.m) summary(dILP3.f) ``` Re-test using the simplest model ```{r,message=FALSE, warning=FALSE} dILP3Model0<-lm(dILP3~(prot), data=all.data) dILP3Model1<-lm(dILP3~sex+(prot), data=all.data) dILP3Model2<-lm(dILP3~sex*(prot), data=all.data) anova(dILP3Model0,dILP3Model1, dILP3Model2) summary(dILP3Model2) ``` Still no difference between males and females There is a marginal effect of protein in females but not in females We can do a power analysis on detecting a sex:protein interaction ```{r,message=FALSE, warning=FALSE} RsqFull<-summary(dILP3Model2)$r.squared RsqRed<-summary(dILP3Model1)$r.squared fsq<- (RsqFull-RsqRed)/(1-RsqFull) predFull<-length(dILP3Model2$coefficients) predRed<-length(dILP3Model1$coefficients) wp.regression(n = nobs(dILP2Model2), p1= predFull, p2=predRed, f2 = fsq) ``` We can do a power analysis on the linear effect of protein on expression in males ```{r,message=FALSE, warning=FALSE} anova(dILP3.m) summary(dILP3.m) pwr.f2.test(u = 1 , v= summary(dILP3.m)$df[2], f2 = (summary(dILP3.m)$r.squared/(1-summary(dILP3.m)$r.squared)), sig.level = 0.05) ``` Plot the model and the TPS ```{r, echo=FALSE, warning=FALSE} # Generate plot space based on female data. This will be used for all subsequent analyses. The geneexpression space specifies the coverage of the surface. The points give the dietary values used. genexpression.space<-genspace(f.data)[c(-3,-4,-5)] points<-genpoints(f.data) #generate the predictions lee.f<-predict(dILP3.f,newdata=genexpression.space,se.fit=F) gen.lee.f<-genexpression.space gen.lee.f$lee<-lee.f gen.lee.m<-genexpression.space gen.lee.m$lee<-mean(m.data$dILP3, na.rm=TRUE) #generate the TPS data tps.f<-with(f.data,Tps(cbind(prot,carb),dILP3)) gen.tps.f<-within(genexpression.space,{ tps<-predict(tps.f,data.frame(prot,carb)) }) tps.m<-with(m.data,Tps(cbind(prot,carb),dILP3)) gen.tps.m<-within(genexpression.space,{ tps<-predict(tps.m,data.frame(prot,carb)) }) #set z-axis attributes to span male and female expression for both model and tps m.min.lee<-min(gen.lee.m$lee) m.max.lee<-max(gen.lee.m$lee) f.min.lee<-min(gen.lee.f$lee) f.max.lee<-max(gen.lee.f$lee) m.min.tps<-min(gen.tps.m$tps) m.max.tps<-max(gen.tps.m$tps) f.min.tps<-min(gen.tps.f$tps) f.max.tps<-max(gen.tps.f$tps) z.range<-c(min(m.min.lee,m.min.tps,f.min.lee,f.min.tps),max(m.max.lee,m.max.tps,f.max.lee,f.max.tps)) #Model Plots gen.plot.lee.f<-ggplot(data=gen.lee.f,aes(x=prot,y=carb,z=lee.f,fill=lee))+ geom_tile(alpha=.9)+theme_bw()+theme(panel.grid=element_blank(), axis.title.x=element_text(size=16), axis.title.y=element_text(size=16), axis.text.x=element_text(size=14), axis.text.y=element_text(size=14), legend.text=element_text(size=14), legend.title=element_text(size=16))+ scale_x_continuous(limits=c(0,max(diets$prot*1.1)),expand=c(0,0))+ scale_y_continuous(expand=c(0,0),limits=c(0,max(diets$carb*1.1)))+ geom_contour(color="black",bins=10)+ xlab("protein (g/L)")+ylab("carbohydrate (g/L)")+ scale_fill_gradientn(colours=c("blue","yellow","red"), name="dILP3", limits=z.range)+ geom_abline(data=ratios,aes(slope=1/ratios, intercept=0),linetype=2,size=rel(0.8))+ geom_point(data=diets,aes(x=prot, y=carb, z=1),size=3,fill="black")+ guides(fill=guide_colorbar(order=1)) + ggtitle("female.lee") gen.plot.lee.m<-ggplot(data=gen.lee.m,aes(x=prot,y=carb,z=lee,fill=lee))+ geom_tile(alpha=.9)+theme_bw()+theme(panel.grid=element_blank(), axis.title.x=element_text(size=16), axis.title.y=element_text(size=16), axis.text.x=element_text(size=14), axis.text.y=element_text(size=14), legend.text=element_text(size=14), legend.title=element_text(size=16))+ scale_x_continuous(limits=c(0,max(diets$prot*1.1)),expand=c(0,0))+ scale_y_continuous(expand=c(0,0),limits=c(0,max(diets$carb*1.1)))+ geom_contour(color="black",bins=10)+ xlab("protein (g/L)")+ylab("carbohydrate (g/L)")+ scale_fill_gradientn(colours=c("blue","yellow","red"), name="dILP3", limits=z.range)+ geom_abline(data=ratios,aes(slope=1/ratios, intercept=0),linetype=2,size=rel(0.8))+ geom_point(data=diets,aes(x=prot, y=carb, z=1),size=3,fill="black")+ guides(fill=guide_colorbar(order=1)) + ggtitle("male.lee") gen.plot.tps.f<-ggplot(data=gen.tps.f,aes(x=prot,y=carb,z=tps,fill=tps))+ geom_tile(alpha=.9)+theme_bw()+theme(panel.grid=element_blank(), axis.title.x=element_text(size=16), axis.title.y=element_text(size=16), axis.text.x=element_text(size=14), axis.text.y=element_text(size=14), legend.text=element_text(size=14), legend.title=element_text(size=16))+ scale_x_continuous(limits=c(0,max(diets$prot*1.1)),expand=c(0,0))+ scale_y_continuous(expand=c(0,0),limits=c(0,max(diets$carb*1.1)))+ geom_contour(color="black",bins=10)+ xlab("protein (g/L)")+ylab("carbohydrate (g/L)")+ scale_fill_gradientn(colours=c("blue","yellow","red"), name="dILP3", limits=z.range)+ geom_abline(data=ratios,aes(slope=1/ratios, intercept=0),linetype=2,size=rel(0.8))+ geom_point(data=diets,aes(x=prot, y=carb, z=1),size=3,fill="black")+ guides(fill=guide_colorbar(order=1)) + ggtitle("female.tps") gen.plot.tps.m<-ggplot(data=gen.tps.m,aes(x=prot,y=carb,z=tps,fill=tps))+ geom_tile(alpha=.9)+theme_bw()+theme(panel.grid=element_blank(), axis.title.x=element_text(size=16), axis.title.y=element_text(size=16), axis.text.x=element_text(size=14), axis.text.y=element_text(size=14), legend.text=element_text(size=14), legend.title=element_text(size=16))+ scale_x_continuous(limits=c(0,max(diets$prot*1.1)),expand=c(0,0))+ scale_y_continuous(expand=c(0,0),limits=c(0,max(diets$carb*1.1)))+ geom_contour(color="black",bins=10)+ xlab("protein (g/L)")+ylab("carbohydrate (g/L)")+ scale_fill_gradientn(colours=c("blue","yellow","red"), name="dILP3", limits=z.range)+ geom_abline(data=ratios,aes(slope=1/ratios, intercept=0),linetype=2,size=rel(0.8))+ geom_point(data=diets,aes(x=prot, y=carb, z=1),size=3,fill="black")+ guides(fill=guide_colorbar(order=1)) + ggtitle("male.tps") # Print plots gen.plot.lee.f gen.plot.lee.m gen.plot.tps.f gen.plot.tps.m ``` ## Analysis of dILP5 Test whether there is a difference between the sexes ```{r,message=FALSE, warning=FALSE} dILP5Model1<-lm(dILP5~sex+(carb:prot+poly(carb,2)+poly(prot,2)), data=all.data) dILP5Model2<-lm(dILP5~sex*(carb:prot+poly(carb,2)+poly(prot,2)), data=all.data) anova(dILP5Model1, dILP5Model2) ``` There is no sex*protein interaction Test the effect of C and P on dILP5 expression in males and females ```{r,message=FALSE, warning=FALSE} dILP5.m<-lm(dILP5~(carb:prot+poly(carb,2)+poly(prot,2)), data=m.data) dILP5.f<-lm(dILP5~(carb:prot+poly(carb,2)+poly(prot,2)), data=f.data) summary(dILP5.m) summary(dILP5.f) ``` There is an effect of protein in females but not in males, but only when you include higher order interactions. Simplify: ```{r,message=FALSE, warning=FALSE} #Non-orthoganol Coef summary(lm(dILP5~(poly(carb,1)+poly(prot,2, raw=TRUE)), data=f.data))[4] #Orthoganol Test dILP5.m<-lm(dILP5~(poly(carb,1)+poly(prot,2)), data=m.data) dILP5.f<-lm(dILP5~(poly(carb,1)+poly(prot,2)), data=f.data) summary(dILP5.m) summary(dILP5.f) ``` There is an effect of protein in females but not in males, but only when you include carbohydrate. Re-test using the simplest model ```{r,message=FALSE, warning=FALSE} dILP5Model0<-lm(dILP5~(poly(carb,1)+poly(prot,2)), data=all.data) dILP5Model1<-lm(dILP5~sex+(poly(carb,1)+poly(prot,2)), data=all.data) dILP5Model2<-lm(dILP5~sex*(poly(carb,1)+poly(prot,2)), data=all.data) anova(dILP5Model0, dILP5Model1, dILP5Model2) summary(dILP5Model1) ``` Still no significant interaction, although males have higher expression We can do a power analysis on the sex interaction: ```{r,message=FALSE, warning=FALSE} RsqFull<-summary(dILP5Model2)$r.squared RsqRed<-summary(dILP5Model1)$r.squared fsq<- (RsqFull-RsqRed)/(1-RsqFull) predFull<-length(dILP5Model2$coefficients) predRed<-length(dILP5Model1$coefficients) wp.regression(n = nobs(dILP2Model2), p1= predFull, p2=predRed, f2 = fsq) ``` Plot the model and the TPS ```{r, echo=FALSE, warning=FALSE} # Generate plot space based on female data. This will be used for all subsequent analyses. The geneexpression space specifies the coverage of the surface. The points give the dietary values used. genexpression.space<-genspace(f.data)[c(-3,-4,-5)] points<-genpoints(f.data) #generate the predictions lee.f<-predict(dILP5.f,newdata=genexpression.space,se.fit=F) gen.lee.f<-genexpression.space gen.lee.f$lee<-lee.f gen.lee.m<-genexpression.space gen.lee.m$lee<-mean(m.data$dILP5,na.rm="TRUE") #generate the TPS data tps.f<-with(f.data,Tps(cbind(prot,carb),dILP5)) gen.tps.f<-within(genexpression.space,{ tps<-predict(tps.f,data.frame(prot,carb)) }) tps.m<-with(m.data,Tps(cbind(prot,carb),dILP5)) gen.tps.m<-within(genexpression.space,{ tps<-predict(tps.m,data.frame(prot,carb)) }) #set z-axis attributes to span male and female expression for both model and tps m.min.lee<-min(gen.lee.m$lee) m.max.lee<-max(gen.lee.m$lee) f.min.lee<-min(gen.lee.f$lee) f.max.lee<-max(gen.lee.f$lee) m.min.tps<-min(gen.tps.m$tps) m.max.tps<-max(gen.tps.m$tps) f.min.tps<-min(gen.tps.f$tps) f.max.tps<-max(gen.tps.f$tps) z.range<-c(min(m.min.lee,m.min.tps,f.min.lee,f.min.tps),max(m.max.lee,m.max.tps,f.max.lee,f.max.tps)) #Model Plots gen.plot.lee.f<-ggplot(data=gen.lee.f,aes(x=prot,y=carb,z=lee.f,fill=lee))+ geom_tile(alpha=.9)+theme_bw()+theme(panel.grid=element_blank(), axis.title.x=element_text(size=16), axis.title.y=element_text(size=16), axis.text.x=element_text(size=14), axis.text.y=element_text(size=14), legend.text=element_text(size=14), legend.title=element_text(size=16))+ scale_x_continuous(limits=c(0,max(diets$prot*1.1)),expand=c(0,0))+ scale_y_continuous(expand=c(0,0),limits=c(0,max(diets$carb*1.1)))+ geom_contour(color="black",bins=10)+ xlab("protein (g/L)")+ylab("carbohydrate (g/L)")+ scale_fill_gradientn(colours=c("blue","yellow","red"), name="dILP5", limits=z.range)+ geom_abline(data=ratios,aes(slope=1/ratios, intercept=0),linetype=2,size=rel(0.8))+ geom_point(data=diets,aes(x=prot, y=carb, z=1),size=3,fill="black")+ guides(fill=guide_colorbar(order=1)) + ggtitle("female.lee") gen.plot.lee.m<-ggplot(data=gen.lee.m,aes(x=prot,y=carb,z=lee,fill=lee))+ geom_tile(alpha=.9)+theme_bw()+theme(panel.grid=element_blank(), axis.title.x=element_text(size=16), axis.title.y=element_text(size=16), axis.text.x=element_text(size=14), axis.text.y=element_text(size=14), legend.text=element_text(size=14), legend.title=element_text(size=16))+ scale_x_continuous(limits=c(0,max(diets$prot*1.1)),expand=c(0,0))+ scale_y_continuous(expand=c(0,0),limits=c(0,max(diets$carb*1.1)))+ geom_contour(color="black",bins=10)+ xlab("protein (g/L)")+ylab("carbohydrate (g/L)")+ scale_fill_gradientn(colours=c("blue","yellow","red"), name="dILP5", limits=z.range)+ geom_abline(data=ratios,aes(slope=1/ratios, intercept=0),linetype=2,size=rel(0.8))+ geom_point(data=diets,aes(x=prot, y=carb, z=1),size=3,fill="black")+ guides(fill=guide_colorbar(order=1)) + ggtitle("male.lee") gen.plot.tps.f<-ggplot(data=gen.tps.f,aes(x=prot,y=carb,z=tps,fill=tps))+ geom_tile(alpha=.9)+theme_bw()+theme(panel.grid=element_blank(), axis.title.x=element_text(size=16), axis.title.y=element_text(size=16), axis.text.x=element_text(size=14), axis.text.y=element_text(size=14), legend.text=element_text(size=14), legend.title=element_text(size=16))+ scale_x_continuous(limits=c(0,max(diets$prot*1.1)),expand=c(0,0))+ scale_y_continuous(expand=c(0,0),limits=c(0,max(diets$carb*1.1)))+ geom_contour(color="black",bins=10)+ xlab("protein (g/L)")+ylab("carbohydrate (g/L)")+ scale_fill_gradientn(colours=c("blue","yellow","red"), name="dILP5", limits=z.range)+ geom_abline(data=ratios,aes(slope=1/ratios, intercept=0),linetype=2,size=rel(0.8))+ geom_point(data=diets,aes(x=prot, y=carb, z=1),size=3,fill="black")+ guides(fill=guide_colorbar(order=1)) + ggtitle("female.tps") gen.plot.tps.m<-ggplot(data=gen.tps.m,aes(x=prot,y=carb,z=tps,fill=tps))+ geom_tile(alpha=.9)+theme_bw()+theme(panel.grid=element_blank(), axis.title.x=element_text(size=16), axis.title.y=element_text(size=16), axis.text.x=element_text(size=14), axis.text.y=element_text(size=14), legend.text=element_text(size=14), legend.title=element_text(size=16))+ scale_x_continuous(limits=c(0,max(diets$prot*1.1)),expand=c(0,0))+ scale_y_continuous(expand=c(0,0),limits=c(0,max(diets$carb*1.1)))+ geom_contour(color="black",bins=10)+ xlab("protein (g/L)")+ylab("carbohydrate (g/L)")+ scale_fill_gradientn(colours=c("blue","yellow","red"), name="dILP5", limits=z.range)+ geom_abline(data=ratios,aes(slope=1/ratios, intercept=0),linetype=2,size=rel(0.8))+ geom_point(data=diets,aes(x=prot, y=carb, z=1),size=3,fill="black")+ guides(fill=guide_colorbar(order=1)) + ggtitle("male.tps") # Print plots gen.plot.lee.f gen.plot.lee.m gen.plot.tps.f gen.plot.tps.m ``` ## Analysis of dILP8 Test whether there is a difference between the sexes ```{r,message=FALSE, warning=FALSE} dILP8Model1<-lm(dILP8~sex+(carb:prot+poly(carb,2)+poly(prot,2)), data=all.data) dILP8Model2<-lm(dILP8~sex*(carb:prot+poly(carb,2)+poly(prot,2)), data=all.data) anova(dILP8Model1, dILP8Model2) summary(dILP8Model1) ``` There is no sex*diet interaction but male expression is higher. Test the effect of C and P on dILP8 expression in males and females ```{r,message=FALSE, warning=FALSE} dILP8.m<-lm(dILP8~(carb:prot+poly(carb,2)+poly(prot,2)), data=m.data) dILP8.f<-lm(dILP8~(carb:prot+poly(carb,2)+poly(prot,2)), data=f.data) summary(dILP8.m) summary(dILP8.f) ``` Protein effects dILP8 in females but not in males. Test the simpler model on dILP8 expression in males and females ```{r,message=FALSE, warning=FALSE} #Non-orthoganol Coef summary(lm(dILP8~(poly(prot,2, raw=TRUE)), data=f.data))[4] #Orthoganol Test dILP8.m<-lm(dILP8~(poly(prot,2)), data=m.data) dILP8.f<-lm(dILP8~(poly(prot,2)), data=f.data) summary(dILP8.m) summary(dILP8.f) ``` Still no effect in males Test using the simpler model ```{r,message=FALSE, warning=FALSE} dILP8Model1<-lm(dILP8~sex+(poly(prot,2)), data=all.data) dILP8Model2<-lm(dILP8~sex*(poly(prot,2)), data=all.data) anova(dILP8Model1, dILP8Model2) summary(dILP8Model2) ``` There is a sex*diet interaction. The negative linear effect of protein on dILP8 expression is less negative in males. Do we have the power to detect a female protein effect-size in males? ```{r,message=FALSE, warning=FALSE} anova(dILP8.m) summary(dILP8.m) pwr.f2.test(u = 1 , v= summary(dILP8.m)$df[2], f2 = (summary(dILP8.f)$r.squared/(1-summary(dILP8.f)$r.squared)), sig.level = 0.05) ``` Plot the model and the TPS ```{r, echo=FALSE, warning=FALSE} # Generate plot space based on female data. This will be used for all subsequent analyses. The geneexpression space specifies the coverage of the surface. The points give the dietary values used. genexpression.space<-genspace(f.data)[c(-3,-4,-5)] points<-genpoints(f.data) #generate the predictions lee.f<-predict(dILP8.f,newdata=genexpression.space,se.fit=F) gen.lee.f<-genexpression.space gen.lee.f$lee<-lee.f gen.lee.m<-genexpression.space gen.lee.m$lee<-mean(m.data$dILP8,na.rm="TRUE") #generate the TPS data tps.f<-with(f.data,Tps(cbind(prot,carb),dILP8)) gen.tps.f<-within(genexpression.space,{ tps<-predict(tps.f,data.frame(prot,carb)) }) tps.m<-with(m.data,Tps(cbind(prot,carb),dILP8)) gen.tps.m<-within(genexpression.space,{ tps<-predict(tps.m,data.frame(prot,carb)) }) #set z-axis attributes to span male and female expression for both model and tps m.min.lee<-min(gen.lee.m$lee) m.max.lee<-max(gen.lee.m$lee) f.min.lee<-min(gen.lee.f$lee) f.max.lee<-max(gen.lee.f$lee) m.min.tps<-min(gen.tps.m$tps) m.max.tps<-max(gen.tps.m$tps) f.min.tps<-min(gen.tps.f$tps) f.max.tps<-max(gen.tps.f$tps) z.range<-c(min(m.min.lee,m.min.tps,f.min.lee,f.min.tps),max(m.max.lee,m.max.tps,f.max.lee,f.max.tps)) #Model Plots gen.plot.lee.f<-ggplot(data=gen.lee.f,aes(x=prot,y=carb,z=lee,fill=lee))+ geom_tile(alpha=.9)+theme_bw()+theme(panel.grid=element_blank(), axis.title.x=element_text(size=16), axis.title.y=element_text(size=16), axis.text.x=element_text(size=14), axis.text.y=element_text(size=14), legend.text=element_text(size=14), legend.title=element_text(size=16))+ scale_x_continuous(limits=c(0,max(diets$prot*1.1)),expand=c(0,0))+ scale_y_continuous(expand=c(0,0),limits=c(0,max(diets$carb*1.1)))+ geom_contour(color="black",bins=10)+ xlab("protein (g/L)")+ylab("carbohydrate (g/L)")+ scale_fill_gradientn(colours=c("blue","yellow","red"), name="dILP8", limits=z.range)+ geom_abline(data=ratios,aes(slope=1/ratios, intercept=0),linetype=2,size=rel(0.8))+ geom_point(data=diets,aes(x=prot, y=carb, z=1),size=3,fill="black")+ guides(fill=guide_colorbar(order=1)) + ggtitle("female.lee") gen.plot.lee.m<-ggplot(data=gen.lee.m,aes(x=prot,y=carb,z=lee,fill=lee))+ geom_tile(alpha=.9)+theme_bw()+theme(panel.grid=element_blank(), axis.title.x=element_text(size=16), axis.title.y=element_text(size=16), axis.text.x=element_text(size=14), axis.text.y=element_text(size=14), legend.text=element_text(size=14), legend.title=element_text(size=16))+ scale_x_continuous(limits=c(0,max(diets$prot*1.1)),expand=c(0,0))+ scale_y_continuous(expand=c(0,0),limits=c(0,max(diets$carb*1.1)))+ geom_contour(color="black",bins=10)+ xlab("protein (g/L)")+ylab("carbohydrate (g/L)")+ scale_fill_gradientn(colours=c("blue","yellow","red"), name="dILP8", limits=z.range)+ geom_abline(data=ratios,aes(slope=1/ratios, intercept=0),linetype=2,size=rel(0.8))+ geom_point(data=diets,aes(x=prot, y=carb, z=1),size=3,fill="black")+ guides(fill=guide_colorbar(order=1)) + ggtitle("male.lee") gen.plot.tps.f<-ggplot(data=gen.tps.f,aes(x=prot,y=carb,z=tps,fill=tps))+ geom_tile(alpha=.9)+theme_bw()+theme(panel.grid=element_blank(), axis.title.x=element_text(size=16), axis.title.y=element_text(size=16), axis.text.x=element_text(size=14), axis.text.y=element_text(size=14), legend.text=element_text(size=14), legend.title=element_text(size=16))+ scale_x_continuous(limits=c(0,max(diets$prot*1.1)),expand=c(0,0))+ scale_y_continuous(expand=c(0,0),limits=c(0,max(diets$carb*1.1)))+ geom_contour(color="black",bins=10)+ xlab("protein (g/L)")+ylab("carbohydrate (g/L)")+ scale_fill_gradientn(colours=c("blue","yellow","red"), name="dILP8", limits=z.range)+ geom_abline(data=ratios,aes(slope=1/ratios, intercept=0),linetype=2,size=rel(0.8))+ geom_point(data=diets,aes(x=prot, y=carb, z=1),size=3,fill="black")+ guides(fill=guide_colorbar(order=1)) + ggtitle("female.tps") gen.plot.tps.m<-ggplot(data=gen.tps.m,aes(x=prot,y=carb,z=tps,fill=tps))+ geom_tile(alpha=.9)+theme_bw()+theme(panel.grid=element_blank(), axis.title.x=element_text(size=16), axis.title.y=element_text(size=16), axis.text.x=element_text(size=14), axis.text.y=element_text(size=14), legend.text=element_text(size=14), legend.title=element_text(size=16))+ scale_x_continuous(limits=c(0,max(diets$prot*1.1)),expand=c(0,0))+ scale_y_continuous(expand=c(0,0),limits=c(0,max(diets$carb*1.1)))+ geom_contour(color="black",bins=10)+ xlab("protein (g/L)")+ylab("carbohydrate (g/L)")+ scale_fill_gradientn(colours=c("blue","yellow","red"), name="dILP8", limits=z.range)+ geom_abline(data=ratios,aes(slope=1/ratios, intercept=0),linetype=2,size=rel(0.8))+ geom_point(data=diets,aes(x=prot, y=carb, z=1),size=3,fill="black")+ guides(fill=guide_colorbar(order=1)) + ggtitle("male.tps") # Print plots gen.plot.lee.f gen.plot.lee.m gen.plot.tps.f gen.plot.tps.m ``` ## Sample variablity It is useful to see how much variablity there is in gene expression among samples within diets. We can inspect the EMS for each gene after fitting a one-way ANOVA of gene expression against ```{r, echo=FALSE, warning=FALSE} var.m<-lm(InR~Sample, data=m.data) anova(var.m) var.f<-lm(InR~Sample, data=f.data) anova(var.f) var.test(var.m,var.f) var.m<-lm(X4eBP~Sample, data=m.data) anova(var.m) var.f<-lm(X4eBP~Sample, data=f.data) anova(var.f) var.test(var.m,var.f) var.m<-lm(CG3071~Sample, data=m.data) anova(var.m) var.f<-lm(CG3071~Sample, data=f.data) anova(var.f) var.test(var.m,var.f) var.m<-lm(Ash2L~Sample, data=m.data) anova(var.m) var.f<-lm(Ash2L~Sample, data=f.data) anova(var.f) var.test(var.m,var.f) var.m<-lm(dILP2~Sample, data=m.data) anova(var.m) var.f<-lm(dILP2~Sample, data=f.data) anova(var.f) var.test(var.m,var.f) var.m<-lm(dILP3~Sample, data=m.data) anova(var.m) var.f<-lm(dILP3~Sample, data=f.data) anova(var.f) var.test(var.m,var.f) var.m<-lm(dILP5~Sample, data=m.data) anova(var.m) var.f<-lm(dILP5~Sample, data=f.data) anova(var.f) var.test(var.m,var.f) var.m<-lm(dILP8~Sample, data=m.data) anova(var.m) var.f<-lm(dILP8~Sample, data=f.data) anova(var.f) var.test(var.m,var.f) ```