--- title: "GrowthRoomIBD" output: html_document: default pdf_document: default --- ```{r setup, include=FALSE} knitr::opts_chunk$set(echo = TRUE) setwd("~/linktofilehere") ioexpRSB <- read.delim("~/Dropbox (University of Michigan)/InbreedingAnalysis-wMVE2020/GrowthRoom2019/ioexpRSB.txt") install.packages("gmodels", repos="http://cran.us.r-project.org") install.packages("dplyr", repos="http://cran.us.r-project.org") install.packages("ggplot2", repos="http://cran.us.r-project.org") install.packages("tidyr", repos="http://cran.us.r-project.org") install.packages("scales", repos="http://cran.us.r-project.org") install.packages("ggthemes", repos="http://cran.us.r-project.org") #install.packages("ztable", repos="http://cran.us.r-project.org") install.packages("lmerTest", repos="http://cran.us.r-project.org") install.packages("lme4", repos="http://cran.us.r-project.org") install.packages("wesanderson", repos="http://cran.us.r-project.org") install.packages("vctrs", repos="http://cran.us.r-project.org") install.packages("rlang", repos="http://cran.us.r-project.org") library(rlang) library(lme4) library(lmerTest) library(dplyr) # for data manipulation library(tidyr) # for reshaping data library(ggplot2) # plotting data library(scales) # for scale_y_continuous(label = percent) library(ggthemes) # for scale_fill_few('medium') #library(ztable) install.packages("viridis", repos="http://cran.us.r-project.org") install.packages("hrbrthemes", repos="http://cran.us.r-project.org") library(viridis) library(hrbrthemes) library(wesanderson) ioexpRSB$PLANTEXIST <- NULL ioexpRSB$CorrBiomass <- ioexpRSB$BIOMASS/ioexpRSB$AVGBIOMASSCONTROL ioexpRSB$HeightChange <- ioexpRSB$HEIGHTAFTER1/ioexpRSB$HEIGHTBEFORE ioexpRSB$LeafChange <- ioexpRSB$LEAFCOUNTAFTER1/ioexpRSB$LEAFCOUNTBEFORE herb1 <- subset(ioexpRSB, treatment == 1) herb2 <- subset(ioexpRSB, treatment == 2) ioexpRSB$treatment <- as.factor(ioexpRSB$treatment) ioexpRSB$MAT <- as.factor(ioexpRSB$MAT) herbonly <- subset(ioexpRSB, treatment !=3) #This orders the Rstatus alphabetically, do not use this command bc it removes the Rstatus, need to figure this out #herbonly$Rstatus <- factor(herbonly$Rstatus, levels = (levels(herbonly$Rstatus))) #One C entry had an outlier value bc control biomass super low, remove from analysis herbonly1 <- subset(herbonly, CorrBiomass < 3) herbonly2 <- subset(herbonly, LeafChange < 1.2) #The below is for alive/almost alive work herbonly$NewDead[herbonly$DEAD=="alive"] <- "Alive" herbonly$NewDead[herbonly$DEAD=="almost"] <- "Almost dead or dead" herbonly$NewDead[herbonly$DEAD=="dead"] <- "Almost dead or dead" test <- herbonly[-c(522),] #Removes pesky row with empty NA in NewDead highcounts <- subset(test, treatment == '1') lowcounts <- subset(test, treatment == '2') ``` ```{r stats Alive/almost dead} test$RecodeDead[test$NewDead=="Alive"] <- "1" test$RecodeDead[test$NewDead=="Almost dead or dead"] <- "0" test$RecodeDead <- as.numeric(test$RecodeDead) #Model not converging when mat line included, also interactions making it tricky for model fit #Test overall for treatment effect, no interactions significant, so remove lmOverall <- glm(RecodeDead ~ Rstatus + breeding + treatment, data = test, family = 'binomial') lmOverallmTrt <- glm(RecodeDead ~ Rstatus + breeding, data = test, family = 'binomial') lmOverallmRstatus <- glm(RecodeDead ~ breeding + treatment, data = test, family = 'binomial') lmOverallmbreeding <- glm(RecodeDead ~ Rstatus + treatment, data = test, family = 'binomial') anova(lmOverall, lmOverallmTrt, test='Chisq') #Test R vs C, R vs S, and S vs C in both environments RvC <- subset(test, Rstatus != 'S') RvS <- subset(test, Rstatus != 'C') SvC <- subset(test, Rstatus != 'R') chisq.test(RvC$Rstatus, RvC$RecodeDead) chisq.test(RvS$Rstatus, RvS$RecodeDead) chisq.test(SvC$Rstatus, SvC$RecodeDead) #Test I/O in each seln line across both treatments R_IO <- subset(test, Rstatus == 'R') S_IO <- subset(test, Rstatus == 'S') C_IO <- subset(test, Rstatus == 'C') lmROverall <- glm(RecodeDead ~ breeding + treatment, data = R_IO, family = 'binomial') lmRmbreeding <- glm(RecodeDead ~ treatment, data = R_IO, family = 'binomial') anova(lmROverall, lmRmbreeding, test='Chisq') lmSOverall <- glm(RecodeDead ~ breeding + treatment, data = S_IO, family = 'binomial') lmSmbreeding <- glm(RecodeDead ~ treatment, data = S_IO, family = 'binomial') anova(lmSOverall, lmSmbreeding, test='Chisq') lmCOverall <- glm(RecodeDead ~ breeding + treatment, data = C_IO, family = 'binomial') lmCmbreeding <- glm(RecodeDead ~ treatment, data = C_IO, family = 'binomial') anova(lmCOverall, lmCmbreeding, test='Chisq') #Very sig effect of treatment, so test within treatment level for Rstatus and breeding? test1.5 <- subset(test, treatment == '2') test3 <- subset(test, treatment == '1') lm1.5 <- glm(RecodeDead ~ Rstatus + breeding , data = test1.5, family = 'binomial') lm1.5mRstatus <- glm(RecodeDead ~ breeding , data = test1.5, family = 'binomial') anova(lm1.5, lm1.5mRstatus, test='Chisq') lm1.5mbreeding <- glm(RecodeDead ~ Rstatus , data = test1.5, family = 'binomial') anova(lm1.5, lm1.5mbreeding, test='Chisq') lm3 <- glm(RecodeDead ~ Rstatus + breeding, data = test3, family = 'binomial') lm3mRstatus <- glm(RecodeDead ~ breeding, data = test3, family = 'binomial') anova(lm3, lm3mRstatus, test='Chisq') lm3mbreed <- glm(RecodeDead ~ Rstatus, data = test3, family = 'binomial') anova(lm3, lm3mbreed, test='Chisq') #1.5 RvC <- subset(test1.5, Rstatus != 'S') RvS <- subset(test1.5, Rstatus != 'C') SvC <- subset(test1.5, Rstatus != 'R') chisq.test(RvC$Rstatus, RvC$RecodeDead) chisq.test(RvS$Rstatus, RvS$RecodeDead) chisq.test(SvC$Rstatus, SvC$RecodeDead) #3.0 RvC <- subset(test3, Rstatus != 'S') RvS <- subset(test3, Rstatus != 'C') SvC <- subset(test3, Rstatus != 'R') chisq.test(RvC$Rstatus, RvC$RecodeDead) chisq.test(RvS$Rstatus, RvS$RecodeDead) chisq.test(SvC$Rstatus, SvC$RecodeDead) test1.5bio <- subset(herbonly1, treatment =='2') test3bio <- subset(herbonly1, treatment =='1') test1.5bio$Rstatus <- factor(test1.5bio$Rstatus) test3bio$Rstatus <- factor(test3bio$Rstatus) test1.5bio$breeding <- factor(test1.5bio$breeding) test3bio$breeding <- factor(test3bio$breeding) ``` ```{r Biomass} #Biomass test over both treatments library(lme4) library(lmerTest) library(lsmeans) herbonly1$Rstatus <- as.factor(herbonly1$Rstatus) herbonly1$breeding <- as.factor(herbonly1$breeding) herbonly1$treatment <- as.factor(herbonly1$treatment) #Full model lm.fitBiomass <- lmer(log1p(CorrBiomass) ~ (1|MAT:Rstatus) + Rstatus + breeding + treatment + Rstatus*treatment + treatment*breeding + Rstatus*breeding + Rstatus*breeding*treatment, data = herbonly1) stepHam<-step(lm.fitBiomass) #Reduced model #Very significant effect of treatment in full model. lm.fitBiomassRed <- lmer(log1p(CorrBiomass) ~ Rstatus + breeding + treatment + (1|MAT:Rstatus) + breeding:treatment + Rstatus:breeding, data = herbonly1) lsmeansLT(lm.fitBiomassRed, "treatment") lm.fitBiomassRed <- lmer(CorrBiomass ~ Rstatus + breeding + treatment + (1|MAT:Rstatus) + breeding:treatment + Rstatus:breeding + Rstatus:treatment + Rstatus:breeding:treatment, data = herbonly1) lsmeansLT(lm.fitBiomassRed, "treatment") lsmeansLT(lm.fitBiomassRed, "breeding") lsmeansLT(lm.fitBiomassRed, "Rstatus") lsmeansLT(lm.fitBiomassRed, "Rstatus:breeding") lsmeansLT(lm.fitBiomassRed, "breeding:treatment") lsmeansLT(lm.fitBiomassRed, "Rstatus:breeding:treatment") lsmeansLT(lm.fitBiomassRed, "Rstatus:treatment") anova(lm.fitBiomassRed) rand(lm.fitBiomassRed) lm.fitBiomassRed <- lmer(CorrBiomass ~ Rstatus + breeding + treatment + (1|MAT:Rstatus) + breeding:treatment + Rstatus:breeding + Rstatus:breeding:treatment, data = herbonly1) anova(lm.fitBiomassRed) rand(lm.fitBiomassRed) #TukeyHSD(aov(CorrBiomass ~ Rstatus, data = herbonly1)) #In 1.5 only lm.fitBiomass1.5 <- lmer(log1p(CorrBiomass) ~ (1|MAT:Rstatus) + Rstatus + breeding + Rstatus*breeding, data = test1.5bio) anova(lm.fitBiomass1.5) rand(lm.fitBiomass1.5) #lsmeansLT(lm.fitBiomass1.5, "breeding") lsmeans(lm.fitBiomass1.5, pairwise ~ breeding) #In 3.0 only lm.fitBiomass3 <- lmer(log1p(CorrBiomass) ~ (1|MAT:Rstatus) + Rstatus + breeding + Rstatus*breeding, data = test3bio) anova(lm.fitBiomass3) rand(lm.fitBiomass3) #lsmeansLT(lm.fitBiomass3, "Rstatus") lsmeans(lm.fitBiomass3, pairwise ~ Rstatus) par(mfrow = c(2, 2)) plot(lm.fitBiomass) #assumptions of anova mostly met after log 1p transforming #Within R only, sig breeding effect test3bioR <- subset(test3bio, Rstatus =='R') lm.fitBiomass3 <- lmer(log1p(CorrBiomass) ~ (1|MAT:Rstatus) + breeding, data = test3bioR) anova(lm.fitBiomass3) rand(lm.fitBiomass3) lsmeans(lm.fitBiomass3, pairwise ~ breeding) #Multiple comparisons of Rstatus TukeyHSD(aov(CorrBiomass ~ Rstatus, data = test3bio)) #Multiple comparisons of breeding TukeyHSD(aov(CorrBiomass ~ breeding, data = test1.5bio)) ``` ```{r looking at data} #library('dplyr') # for data manipulation #library('tidyr') # for reshaping data #Percent alive and almost dead or dead by Rstatus : over both treatments herbonly_perc <- test %>% count(Rstatus, NewDead) %>% group_by(Rstatus,add=TRUE, .drop=TRUE) %>% mutate(perc = prop.table(n)*100) %>% # mutate count(n) into perc select(-n) %>% # remove the count... spread(Rstatus, perc) herbonly_perc <- as.data.frame(herbonly_perc) herbonly_perc #Percent alive and almost dead or dead by Rstatus : over both treatments herbonly_percTrt <- test %>% count(treatment, NewDead) %>% group_by(treatment,add=TRUE, .drop=TRUE) %>% mutate(perc = prop.table(n)*100) %>% # mutate count(n) into perc select(-n) %>% # remove the count... spread(treatment, perc) herbonly_percTrt <- as.data.frame(herbonly_percTrt) herbonly_percTrt #Percent alive overall by treatment herbonly_perc0 <- test %>% count(treatment,breeding,NewDead) %>% group_by(treatment, breeding, add=TRUE, .drop=TRUE) %>% mutate(perc = prop.table(n)*100) %>% # mutate count(n) into perc select(-n) %>% # remove the count... spread(treatment, perc) herbonly_perc0 #Percent alive overall by IO herbonly_percbreed <- test %>% count(breeding,NewDead) %>% group_by(breeding, add=TRUE, .drop=TRUE) %>% mutate(perc = prop.table(n)*100) %>% # mutate count(n) into perc select(-n) %>% # remove the count... spread(breeding, perc) herbonly_percbreed #Percent alive by treatment and Rstatus herbonly_perc1 <- test %>% count(Rstatus, treatment,NewDead) %>% group_by(Rstatus, treatment, add=TRUE, .drop=TRUE) %>% mutate(perc = prop.table(n)*100) %>% # mutate count(n) into perc select(-n) %>% # remove the count... spread(treatment, perc) herbonly_perc1 #Percent alive almost dead & dead at 3.0 herbonly3_perc <- highcounts %>% count(Rstatus, NewDead) %>% group_by(Rstatus,add=TRUE, .drop=TRUE) %>% mutate(perc = prop.table(n)*100) %>% # mutate count(n) into perc select(-n) %>% # remove the count... spread(Rstatus, perc) herbonly3_perc <- as.data.frame(herbonly3_perc) herbonly3_perc #Percent alive almost dead & dead at 1.5 herbonly1.5_perc <- lowcounts %>% count(Rstatus, NewDead) %>% group_by(Rstatus,add=TRUE, .drop=TRUE) %>% mutate(perc = prop.table(n)*100) %>% # mutate count(n) into perc select(-n) %>% # remove the count... spread(Rstatus, perc) herbonly1.5_perc <- as.data.frame(herbonly1.5_perc) herbonly1.5_perc #Percent alive and almost dead or dead by Rstatus and breeding values herbonly_percIO <- test %>% count(Rstatus, NewDead, breeding) %>% group_by(Rstatus,breeding, add=TRUE, .drop=TRUE) %>% mutate(perc = prop.table(n)*100) %>% # mutate count(n) into perc select(-n) %>% # remove the count... spread(Rstatus, perc) herbonly_percIO <- as.data.frame(herbonly_percIO) herbonly_percIO #Percent alive and almost dead or dead by Rstatus and breeding values, at 3.0 herbonly3_percIO <- highcounts %>% count(Rstatus, NewDead, breeding) %>% group_by(Rstatus,breeding, add=TRUE, .drop=TRUE) %>% mutate(perc = prop.table(n)*100) %>% # mutate count(n) into perc select(-n) %>% # remove the count... spread(Rstatus, perc) herbonly3_percIO <- as.data.frame(herbonly3_percIO) herbonly3_percIO #Percent alive and almost dead or dead by Rstatus and breeding values, at 1.5 herbonly1.5_percIO <- lowcounts %>% count(Rstatus, NewDead, breeding) %>% group_by(Rstatus,breeding, add=TRUE, .drop=TRUE) %>% mutate(perc = prop.table(n)*100) %>% # mutate count(n) into perc select(-n) %>% # remove the count... spread(Rstatus, perc) herbonly1.5_percIO <- as.data.frame(herbonly1.5_percIO) herbonly1.5_percIO ioexpRSB %>% group_by(Rstatus, breeding) %>% summarize(mean_CorrBiomass = mean(CorrBiomass, na.rm = TRUE), mean_HeightChange = mean(HeightChange, na.rm = TRUE), mean_LeafChange = mean(LeafChange, na.rm = TRUE), mean_NumGreen = mean(NumGreen, na.rm = TRUE), mean_NumDam50 = mean(NumdamageGT50, na.rm = TRUE)) herb1 %>% group_by(Rstatus, breeding) %>% summarize(mean_CorrBiomass = mean(CorrBiomass, na.rm = TRUE), mean_HeightChange = mean(HeightChange, na.rm = TRUE), mean_LeafChange = mean(LeafChange, na.rm = TRUE), mean_NumGreen = mean(NumGreen, na.rm = TRUE), mean_NumDam50 = mean(NumdamageGT50, na.rm = TRUE)) herb2 %>% group_by(Rstatus, breeding) %>% summarize(mean_CorrBiomass = mean(CorrBiomass, na.rm = TRUE), mean_HeightChange = mean(HeightChange, na.rm = TRUE), mean_LeafChange = mean(LeafChange, na.rm = TRUE), mean_NumGreen = mean(NumGreen, na.rm = TRUE), mean_NumDam50 = mean(NumdamageGT50, na.rm = TRUE)) ``` ```{r plot biomass IO} to_plot <- herbonly1 %>% group_by(treatment, Rstatus, breeding,MAT) %>% mutate (CorrBiomass = mean(CorrBiomass, na.rm=T)) %>% group_by(treatment, Rstatus, breeding,MAT) %>% summarize (avg=mean(CorrBiomass, na.rm=T), med = median(CorrBiomass, na.rm=T), sd=sd(CorrBiomass, na.rm=T), se = sd(CorrBiomass, na.rm=T)/sqrt(n())) write.table(to_plot, file="CorrBiomassMat.txt", sep = '\t')