--- title: "R Notebook" output: html_notebook --- This is an [R Markdown](http://rmarkdown.rstudio.com) Notebook. When you execute code within the notebook, the results appear beneath the code. Try executing this chunk by clicking the *Run* button within the chunk or by placing your cursor inside it and pressing *Cmd+Shift+Enter*. ```{r} #library("gdata") library("tidyverse") library("car") library("ggplot2") library("plyr") library('dplyr') library("lsmeans") #require("mgcv") #require("nlme") #require("AICcmodavg") library("multcompView") library("tidyr") library("stringr") library(lme4) ``` ```{r} #setwd("/Users/acha0034/Dropbox/Sgro Lab/EXPERIMENTS/Avi/Feeding assays") DTFEEDmacro <- read_csv("~/Dropbox/Sgro Lab/PEOPLE/Avi/MANUSCRIPT CHAPTER2/Submission_JIP/Resubmission/Macronutrient_Intake2.csv", col_names = TRUE) DTFEEDmacro ``` ```{r} means <- DTFEEDING %>% group_by(Population,treatment, ratio ) %>% summarize(avgweightnormalisedfoodintake = mean(weightnormalisedfoodintake), sdweightnormalisedfoodintake = sd(weightnormalisedfoodintake)) means ``` ```{r} library(tidyverse) require(dunn.test) library (dunn.test) library(rstatix) DTFEEDmacro <- DTFEEDmacro %>% mutate(logtotal_intake=log10(total_intake+1), ratioF = as.factor(ratio)) ratios<-data.frame(logratios=c(log(1/5),log(1/3),log(2/3)), ratios=c(1/5,1/3,2/3)) intake_array_plot <- ggplot(data=DTFEEDmacro, aes(x= prot_intake, y= carb_intake, colour= as.factor(treatment), fill= as.factor (treatment), shape = as.factor(treatment)))+ theme_bw()+ theme(panel.grid=element_blank(),axis.title.x=element_text(size=12), axis.title.y=element_text(size=12), axis.text.y=element_text(size=12), axis.text.x=element_text(size=12), panel.background = element_rect(colour = "black"), legend.background = element_rect(), legend.key = element_rect(colour = "white"), legend.text=element_text(size=12))+ theme(legend.title=element_blank())+ theme(axis.title.x = element_text(vjust=-0.5), axis.title.y = element_text(vjust=1.5))+ coord_fixed(ratio=1)+ #scale_shape(solid=FALSE)+ xlab("Protein eaten (mg for 12 larvae)")+ ylab("Carbohydrate eaten (mg for 12 larvae)")+ scale_y_continuous(expand=c(0,0),limits=c(0,32), breaks=c(0,4,8,12,16,20,24,28,32))+ scale_x_continuous(expand=c(0,0),limits=c(0,15), breaks=c(0,4,8,12,15))+ #scale_colour_manual(values=c("brown4", "chocolate3"))+ geom_point(size=2, alpha=0.6)+ facet_grid(~Population) + geom_abline(data = ratios, aes(intercept = 0, slope=1/ratios),linetype=2,size=rel(0.25), colour="grey30")+ geom_blank() intake_array_plot intake_array_plotF <-intake_array_plot + theme(panel.spacing.x=unit(1.5, "lines")) pdf("./intake_array_plot.pdf", useDingbats=FALSE) grid.arrange(intake_array_plotF) dev.off() colors () ``` ```{r} ###Dispersion statistics ###test for differences in dispersion DTFEEDmacro2 <- read_csv("~/Dropbox/Sgro Lab/PEOPLE/Avi/MANUSCRIPT CHAPTER2/Submission_JIP/Resubmission/Macronutrient_Intake3.1.csv", col_names = TRUE) DTFEEDmacro2 ##grouping by pop and macro DTFEEDmacro2 <- DTFEEDmacro2 %>% mutate(Group = paste(Population,macro,sep="_")) fligner.Eaten<-fligner.test(intake ~ Group, data = DTFEEDmacro2) fligner.Eaten ##grouping by pop and diet type (treatment) DTFEEDmacro2.2 <- DTFEEDmacro2 %>% mutate(Group = paste(Population,treatment,sep="_")) fligner.Eaten2<-fligner.test(intake ~ Group, data = DTFEEDmacro2.2) fligner.Eaten2 ##GROUPING BY DIET TYPE AND MACRO DTFEEDmacro2.3 <- DTFEEDmacro2 %>% mutate(Group = paste(macro,treatment,sep="_")) fligner.Eaten3<-fligner.test(intake ~ Group, data = DTFEEDmacro2.3) fligner.Eaten3 #Fligner-Killeen test of homogeneity of variances, so if significant then that means variances across the groups arent the same? #data: Eaten by Group #Fligner-Killeen:med chi-squared = 127.0691, df = 3, p-value < 2.2e-16 #the function below substracts each data point from the median of each group to generate an estimate of variation within the data. It then uses three different tests to assess whether the dispersion differs between groups. You could choose either food type, macronutrient type, or population + macronutrient as your group. variability.trans = function(data) { ## find median for group of data ## subtract median; take absolute value a = sort(abs((data-median(data))/median(data))) ## if odd sample size, remove exactly one zero #if (length(a)%%2) #a[a != 0 | duplicated(a)] #else a } library(lme4) ## set up data frame with transformed data for anova V2 = lapply(split(DTFEEDmacro2[["intake"]], DTFEEDmacro2[["Group"]]), variability.trans) V1 = rep(seq(length(V2)), lapply(V2, length)) levdata = data.frame(V1 = factor(V1), V2 = unsplit(V2, V1)) ## perform one-way anova on transformed data cat("Overall ANOVA for Dispersion\n") model <- glm(V2 ~ V1, data = levdata, family ="quasipoisson") Anova (model) compare <- emmeans(model, ~ V1) multcomp::cld(compare) ## perform pairwise t-tests on transformed data cat("\nPairwise Dispersion Tests\n") pairwise.t.test(levdata$V2, levdata$V1, p.adjust = "holm") ## perform Kruskal-Wallis test on the transformed data cat("Overall Kruskal-Wallis Test for Dispersion\n") kruskal.test(V2 ~ V1, levdata) ## perform pairwise Wilcox tests on the transformed data cat("\nPairwise Dunn's Tests for Dispersion\n") dunn.test(levdata$V2, levdata$V1, method = "holm") #Parametric tests more powerful than non-parametric ##diagnostic checks hist(levdata$V2) qqPlot(levdata$V2) ##if you call V2 you'll see the order of the groups, which are relabelled in V1 and levdata. below are the groupings for macro and pop # 1 = Townsville_carb # 2 = Townsville_Prot # 3 = ballina_Carb # 4 = ballina_Prot #5 = melbourne_carb #6 = melbourne_prot # townsville_carb and ballina_carb group similar, ballina_prot/melb_prot/townsville_prot group similar; melb_carb similar to both these groups...protein leveraging across all populations more variation in carb intake #below are groupings for macro and diet type/treatment #1 = carb_cc #2 = carb_pc #3 = prot_cc #4 = prot_pc #more variation in the protein varies diet than carb varies? #carb_cc, carb_pc and prot_cc all different from each other. prot_pc similar to both prot_cc and carb_cc ``` ```{r} library(gridExtra) Fig1a <- ggplot(data=DTFEEDmacro , aes(x= weightnormalised_prot_intake,y=weightnormalised_carb_intake , colour= as.factor(treatment)))+ xlab("Protein Intake per 12L3s")+ ylab("Carb Intake per 12L3s ")+ scale_y_continuous(limits = c(0,30)) + scale_x_continuous(limits = c(0,30)) + theme_bw()+ theme(panel.grid=element_blank(),axis.title.x=element_text(size=16), axis.title.y=element_text(size=16), axis.text.y=element_text(size=16), axis.text.x=element_text(size=12), legend.background = element_rect(), legend.key = element_rect(colour = "white"), legend.text=element_text(face="italic", size=12))+ scale_colour_hue(name="Treatment", breaks=c("cc", "pc"), labels=c("Protein_varies", "Carb_varies")) + theme(legend.title=element_blank())+expand_limits(y=0)+ expand_limits(x=0.1) + #theme(panel.grid=element_blank(),axis.title.x=element_text(size=20), axis.title.y=element_text(size=20), #axis.text.y=element_text(size=16), axis.text.x=element_text(size=16), panel.background = element_rect(colour = "black"), legend.background = element_rect(), legend.key = element_rect(colour = "white"), legend.text=element_text(face="italic", size=14))+ #theme(axis.title.x = element_text(vjust=-0.5), axis.title.y = element_text(vjust=1.5))+ #theme(legend.title=element_blank())+ #scale_y_continuous(limits= c(0.00,0.13), breaks = seq (0.00, 0.13, 0.02))+ #geom_smooth( aes (fill = as.factor (ratio)))+ guides (fill=FALSE) + geom_point(alpha = 0.5)+ geom_line(aes(linetype= as.factor(ratio))) + facet_grid(~ Population)+ geom_blank() Fig1a pdf("./Macronutrient_Intake1.pdf", width=7, height=5, useDingbats=FALSE) grid.arrange(Fig1a, ncol=1) dev.off() #png("./DTFEEDING.png", width=10, height=8, useDingbats=FALSE) #grid.arrange(Fig1a, ncol=1) #dev.off() #expand_limits(y=0) ``` ``` Add a new chunk by clicking the *Insert Chunk* button on the toolbar or by pressing *Cmd+Option+I*. When you save the notebook, an HTML file containing the code and output will be saved alongside it (click the *Preview* button or press *Cmd+Shift+K* to preview the HTML file). The preview shows you a rendered HTML copy of the contents of the editor. Consequently, unlike *Knit*, *Preview* does not run any R code chunks. Instead, the output of the chunk when it was last run in the editor is displayed.