---
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')