################################################# #### Garzke, J; Connor, SJ; Sommer, U; O'Connor, MI #### Mixed effects models #### Code by Mary O'Connor and Jessica Garzke ################################################# ### load libraries library(lme4) library(MuMIn) library(plyr) library(broom) library(dplyr) library(nlme) ### set working directory and load data data <- read.csv("./temporal_data") dim(data) head(data) tail(data) data <- data[-(241:255),] ### weekly average of daily average temps tdata <- read.csv("./avgtemps") head(tdata) temp.data <- ddply(tdata, .(Week, Tank), summarise, mean(Temperature)) head(temp.data) names(temp.data) <- c('Week', 'Tank', 'wklyTemp') data2 <- merge(data, temp.data, by.x = c("week", "Tank"), by.y = c("Week", "Tank")) data <- data2 o.data <- read.csv("./oxygen_temp_temporal.csv") o.data <- o.data[,-(4:14)] o.data <- o.data[,-2] head(o.data) dim(o.data) data4 <- merge(data, o.data, by.x = c("week", "Tank"), by.y = c("week", "Tank")) data <- data4 ### following Barneche et al 2014 ### for each tank, we will want: ### size corrected biomass sdata <- read.csv("./zooplankton size for correction no eggs.csv") names(sdata) <- c('treatment', 'tank', 'week', 'species', 'stage', 'size') sdata$drywt <- ifelse(sdata$species == 'Daphnia', exp(1.468 + 2.829*log(sdata$size)), exp(1.821 + 0.654*log(sdata$size))) sdata$C.ug <- sdata$drywt*0.5 #add in notonectids sdata$noto.wt.corr <- ifelse(sdata$treatment == 'PZN', 2*((10.8/4)^.75), 0) #adding in a mass-corrected biomass estimate for notonectids mass.data <- ddply(sdata, .(week, tank, treatment), summarise, sum((C.ug^0.75))) # estimating a mass corrected biomass for adult zooplankton in all tanks. mass.data2 <- ddply(sdata, .(week, tank, treatment), summarise, sum((C.ug^0.75)+noto.wt.corr)) # estimating a mass corrected biomass for adult zooplankton incl Noto in all tanks. mass.data3 <- ddply(sdata, .(week, tank, treatment), summarise, sum((C.ug))) mass.data$M.corrz <- mass.data2$..1 # corrected total zooplankton mass.data$zpC.uncorr <- mass.data3$..1 # uncorrected total (no notonectids) names(mass.data) <- c('week', 'tank', 'treatment', 'M.corrT', 'M.corrz', 'M.uncorrz') ### load libraries, define variables and add columns k <- 8.617342*10^-5 # eV/K ## create test temp data$invT <- 1/((data$wklyTemp + 273)*k) data$invT <- as.numeric(as.character(data$invT)) data$PP.biomass <- (data$chla*55) #chla (ug/L)* 55 C in PP / 1 chla = ugPPC/L data$ZP.carbon1 <- ifelse(data$trophic.level=='P',0, data$M.uncorrz) # raw carbon data$HeteroB <- ifelse(data$trophic.level=='PZN', data$M.uncorrz + 10.8/2, data$ZP.carbon1) # add notonectid weight data$total.carbon <- data$PP.biomass + data$HeteroB data$HA <- data$HeteroB/data$PP.biomass ### estimating NPP and ER from the raw data: ## correct of water-air oxygen flux ### adjusting for effects of temperature on physical exchange C.star <- function(T) exp(7.7117 - 1.31403*log(T+45.93)) - 0.035 # the -0.035 is an approximate adjustment for elevation C.star(T) # yields the oxygen concentration expected at a given temperature (T in C). #calculate NPP and ER (hourly), in terms of umol O2 / l / hr, following Yvon-Durocher et al 2010. # O2 has molar mass of 32g/mol. so 1 umol = 32 ug. so take ug/32 data$NPP2 <- (((data$dusk - data$dawn1) - (C.star(data$temp.dusk) - C.star(data$temp.dawn1)))*1000)/(32) # oxygen produced umol / L /day, net all respiration. raw data is mg/L. Subtract o2 water/atm flux due to change in temperature data$ER2 <- -(24/data$hours2)*(((data$dawn2 - data$dusk) - (C.star(data$temp.dawn2) - C.star(data$temp.dusk)))*1000)/(32) # amount of oxygen consumed per day via respiaration. negative to get the change in oxygen umol / L /day; oxygen used in the dark and daylight. MeanER can be greater than meanNPP, because NPP reflects ER already. data$GPP <- data$NPP2+(data$ER2/24)*data$hours1 # daily oxygen production (NPP2) + estimated daytime community respiration (daily R / 24 * hours daylight) data$NEM <- data$ER2/data$GPP # following Yvon Durochers 2010. NEM > 1 means the system is respiring more than it's fixing per day. This does not need to be logged. data$NPP.mass <- data$NPP2 / (data$PP.biomass) # NPP on ummol 02/L/day/ugCPP data$ER.mass <- data$ER2/(data$M.corrz + data$PP.biomass) # ER on ummol 02/L/day/ugTPP divited by mass corrected zooplankton biomass (including notonectids) data <- data[data$week >= '4',] #### data prep complete #### #### ANALYSIS #### ## Does NPP vary with temperature? ## figures on invT #data1 <- data[(data$NPP2 >= 0.3),] # three negative values and one very small value now, not sure what to do about them. #hist(data[(data$NPP2 >= 0.5),]$NPP2) hist(data$NPP2) hist(log(data$NPP2)) plot(log(data$NPP2)~data$Tank, pch = 19, col = data$trophic.level) plot(log(data$NPP2)~data$week, pch = 19, col = data$trophic.level) plot(log(data$NPP2)~I(data$invT-mean(data$invT)), pch = 19, col = data$trophic.level, ylim = c(0,6)) abline(3.00, -1.14, lwd = 2, col = 1) abline((3.00+0.18), (-1.14-0.32), lwd = 2, col = 2) abline((3.00+0.23), (-1.14+0.06), lwd = 2, col = 3) ## analysis ## determine need for random effects in the full model: modNPP4a <- lmer(log(NPP2) ~ 1 + I(invT-mean(invT))*trophic.level + (I(invT-mean(invT))|week), data=data, REML = FALSE, na.action=na.omit) modNPP4b<-lmer(log(NPP2)~1+I(invT-mean(invT))*trophic.level + (1|week), data=data, REML = FALSE, na.action=na.omit) anova(modNPP4a, modNPP4b) modNPP0 <- lmer(log(NPP2) ~ 1 + (1|week), data=data, REML = FALSE, na.action=na.omit) modNPP1 <- lmer(log(NPP2) ~ 1 + I(invT-mean(invT)) + (1|week), data=data, REML = FALSE, na.action=na.omit) modNPP2 <- lmer(log(NPP2) ~ 1 + I(invT-mean(invT)) + trophic.level + (1|week), data=data, REML = FALSE, na.action=na.omit) modNPP4 <- lmer(log(NPP2) ~ 1 + I(invT-mean(invT))*trophic.level + (1|week), data=data, REML = FALSE, na.action=na.omit) model.sel(modNPP0, modNPP1, modNPP2, modNPP4) anova(modNPP4, modNPP1) anova(modNPP1, modNPP2) anova(modNPP2, modNPP0) # for model fitting: modNPP1r <- lmer(log(NPP2) ~ 1 + I(invT-mean(invT)) + (1|week), data=data, REML = TRUE, na.action=na.omit) modNPP2r <- lmer(log(NPP2) ~ 1 + I(invT-mean(invT)) + trophic.level + (1|week), data=data, REML = TRUE, na.action=na.omit) modNPP4r <- lmer(log(NPP2) ~ 1 + I(invT-mean(invT))*trophic.level + (1|week), data=data, REML = TRUE, na.action=na.omit) summary(modNPP4r) confint(modNPP4r) m.avg <- model.avg(modNPP2r, modNPP4r, modNPP1r) summary(m.avg) confint(m.avg) #Test for temporal autocorrelation with nlme modNPP <- lme(log(NPP2)~1+I(invT-mean(invT))*trophic.level, random= ~1|week, correlation=corAR1(), data=data, na.action=na.omit) ## Does mass-specific NPP vary with temperature? ## figures ### NPP in umol/L/day per ugC #data1 <- data[(data$NPP2 >= 0.3),] hist(data$NPP.mass) hist(log(data$NPP.mass)) # hist(log(data$NPP.mass+.001)) # no reason to add a number, there are no -inf values. plot(log(data$NPP.mass)~data$Tank, pch = 19, col = data$trophic.level) plot(log(data$NPP.mass)~data$week, pch = 19, col = data$trophic.level) plot(log(data$NPP.mass)~I(data$invT-mean(data1$invT)), pch = 19, col = data$trophic.level) abline(-2.83, -1.89, lwd = 2, col = 1) abline((-2.83+1.36), (-1.89-1.41), lwd = 2, col = 2) abline((-2.83+0.46), (-1.89-0.05), lwd = 2, col = 3) ## analysis ## determine need for random effects in the full model: modNPPma <- lmer(log(NPP.mass) ~ 1 + I(invT-mean(invT))*trophic.level + (I(invT-mean(invT))|week), data=data, REML = FALSE, na.action=na.omit) modNPPmb <- lmer(log(NPP.mass) ~ 1 + I(invT-mean(invT))*trophic.level + (1|week), data=data, REML = FALSE, na.action=na.omit) anova(modNPPma, modNPPmb) modNPPm0 <- lmer(log(NPP.mass) ~ 1 + (I(invT-mean(invT))|week), data=data, REML = FALSE, na.action=na.omit) modNPPm1 <- lmer(log(NPP.mass) ~ 1 + I(invT-mean(invT)) + (I(invT-mean(invT))|week), data=data, REML = FALSE, na.action=na.omit) modNPPm2 <- lmer(log(NPP.mass) ~ 1 + I(invT-mean(invT)) + trophic.level + (I(invT-mean(invT))|week), data=data, REML = FALSE, na.action=na.omit) modNPPm4 <- lmer(log(NPP.mass) ~ 1 + I(invT-mean(invT))*trophic.level + (I(invT-mean(invT))|week), data=data, REML = FALSE, na.action=na.omit) model.sel(modNPPm0, modNPPm1, modNPPm2, modNPPm4) anova(modNPPm4, modNPPm2) anova(modNPPm2, modNPPm1) anova(modNPPm1, modNPPm2) # for model fitting: modNPPm4r <- lmer(log(NPP.mass) ~ 1 + I(invT-mean(invT))*trophic.level + (I(invT-mean(invT))|week), data=data, REML = TRUE, na.action=na.omit) summary(modNPPm4r) confint(modNPPm4r) #Test for temporal autocorrelation with nlme modNPPm4r <- lme(NPP.mass ~ 1 + I(invT-mean(invT))*trophic.level, random=~ 1|week, data=data, correlation=corAR1(), na.action=na.omit) ## Does ER vary with temperature? ## figures hist(data$ER2) data1 <- data[(data$ER2 >= 0),] hist(log(data1$ER2)) #hist(log(data$calc.ER+1)) plot(log(data$ER2)~data$Tank, pch = 19, col = data$trophic.level) plot(log(data$ER2)~data$week, pch = 19, col = data$trophic.level) plot(log(data1$ER2)~I(data1$invT-mean(data1$invT)), pch = 19, col = data1$trophic.level, ylim = c(0,6)) abline(3.88, -0.75, lwd = 2, col = 1) abline((3.88+0.49), (-0.75+0.29), lwd = 2, col = 2) abline((3.88+0.01), (-0.75+0.62), lwd = 2, col = 3) #analysis ## determine need for random effects in the full model: modERa <- lmer(log(ER2) ~ 1 + I(invT-mean(invT))*trophic.level + (I(invT-mean(invT))|week), data=data1, REML = FALSE, na.action=na.omit) modERb <- lmer(log(ER2) ~ 1 + I(invT-mean(invT))*trophic.level + (1|week), data=data1, REML = FALSE, na.action=na.omit) anova(modERa, modERb) modER0 <- lmer(log(ER2) ~ 1 + (1|week), data=data1, REML = FALSE, na.action=na.omit) modER1 <- lmer(log(ER2) ~ 1 + I(invT-mean(invT)) + (1|week), data=data1, REML = FALSE, na.action=na.omit) modER2 <- lmer(log(ER2) ~ 1 + I(invT-mean(invT)) + trophic.level + (1|week), data=data1, REML = FALSE, na.action=na.omit) modER4 <- lmer(log(ER2) ~ 1 + I(invT-mean(invT))*trophic.level + (1|week), data=data1, REML = FALSE, na.action=na.omit) model.sel(modER0, modER1, modER2, modER4) anova(modER4, modER2) anova(modER1, modER4) anova(modER1, modER0) # for model fitting: modER4r <- lmer(log(ER2) ~ 1 + I(invT-mean(invT))*trophic.level + (1|week), data=data1, REML = TRUE, na.action=na.omit) modER2r <- lmer(log(ER2) ~ 1 + I(invT-mean(invT)) + trophic.level + (1|week), data=data1, REML = TRUE, na.action=na.omit) m.avg <- model.avg(modER4r, modER2r) summary(m.avg) confint(m.avg) summary(modER4r) confint(modER4r) #Test for temporal autocorrelation with nlme modER4r <- lme(log(ER2) ~ 1 + I(invT-mean(invT))*trophic.level, random=~ 1|week, data=data, correlation=corAR1(), na.action=na.omit) ## Does mass-specific ER vary with temperature? ## figures data1 <- data[(data$ER2 >= 0),] hist(data$ER.mass) hist(log(data1$ER.mass)) plot(log(data$ER.mass)~data$Tank, pch = 19, col = data$trophic.level) plot(log(data$ER.mass)~data$week, pch = 19, col = data$trophic.level) plot(log(data$ER.mass)~I(data$invT-mean(data$invT)), pch = 19, col = data$trophic.level) ## analysis ## determine need for random effects in the full model: modERma <- lmer(log(ER.mass) ~ 1 + I(invT-mean(invT))*trophic.level + (I(invT-mean(invT))|week), data=data1, REML = FALSE, na.action=na.omit) #, control = lmerControl(optimizer = "Nelder_Mead") modERmb <- lmer(log(ER.mass) ~ 1 + I(invT-mean(invT))*trophic.level + (1|week), data=data1, REML = FALSE, na.action=na.omit) #, control = lmerControl(optimizer = "Nelder_Mead") anova(modERma, modERmb) modERm0 <- lmer(log(ER.mass) ~ 1 + (1|week), data=data1, REML = FALSE, na.action=na.omit) modERm1 <- lmer(log(ER.mass) ~ 1 + I(invT-mean(invT)) + (1|week), data=data1, REML = FALSE, na.action=na.omit) modERm2 <- lmer(log(ER.mass) ~ 1 + I(invT-mean(invT)) + trophic.level + (1|week), data=data1, REML = FALSE, na.action=na.omit) modERm4 <- lmer(log(ER.mass) ~ 1 + I(invT-mean(invT))*trophic.level + (1|week), data=data1, REML = FALSE, na.action=na.omit) model.sel(modERm0, modERm1, modERm2, modERm4) # for model fitting: modERm4r <- lmer(log(ER.mass) ~ 1 + I(invT-mean(invT))*trophic.level + (1|week), data=data, REML = TRUE, na.action=na.omit) summary(modERm4r) confint(modERm4r) #Test for temporal autocorrelation with nlme modERm4r <- lme(log(ER.mass) ~ 1 + I(invT-mean(invT))*trophic.level, random=~ 1|week, data=data, correlation=corAR1(), na.action=na.omit) ##### BIOMASS RESULTS ##### #### PP BIOMASS ### ## Does total PP biomass vary with temperature and FCL? ## figures hist(data$PP.biomass) hist(log(data$PP.biomass)) plot(log(data$PP.biomass)~data$Tank, pch = 19, col = data$trophic.level) plot(log(data$PP.biomass)~data$week, pch = 19, col = data$trophic.level) plot(log(data$PP.biomass)~I(data$invT-mean(data$invT)), pch = 19, col = data$trophic.level) ## analysis ## determine need for random effects in the full model: modPPa <- lmer(log(PP.biomass) ~ 1 + I(invT-mean(invT))*trophic.level + (I(invT-mean(invT))|week), data=data, REML = FALSE, na.action=na.omit) modPPb <- lmer(log(PP.biomass) ~ 1 + I(invT-mean(invT))*trophic.level + (1|week), data=data, REML = FALSE, na.action=na.omit) anova(modPPa, modPPb) # modPPc<-lme(log(PP.biomass)~1+I(invT-mean(invT))*trophic.level, random = ~1|week, data=data, method="ML", na.action=na.omit) # anova(modPPc, modPPb) modPP0 <- lmer(log(PP.biomass) ~ 1 + (I(invT-mean(invT))|week), data=data, REML = FALSE, na.action=na.omit) modPP1 <- lmer(log(PP.biomass) ~ 1 + I(invT-mean(invT)) + (I(invT-mean(invT))|week), data=data, REML = FALSE, na.action=na.omit) modPP2 <- lmer(log(PP.biomass) ~ 1 + I(invT-mean(invT))+trophic.level + (I(invT-mean(invT))|week), data=data, REML = FALSE, na.action=na.omit) modPP4 <- lmer(log(PP.biomass) ~ 1 + I(invT-mean(invT))*trophic.level + (I(invT-mean(invT))|week), data=data, REML = FALSE, na.action=na.omit) model.sel(modPP0, modPP1, modPP2, modPP4) anova(modPP4, modPP2) # for model fitting: modPP4r <- lmer(log(PP.biomass) ~ 1 + I(invT-mean(invT))*trophic.level + (I(invT-mean(invT))|week), data=data, REML = TRUE, na.action=na.omit) summary(modPP4r) confint(modPP4r) #Test for temporal autocorrelation with nlme modPP4r <- lme(log(PP.biomass) ~ 1 + I(invT-mean(invT))*trophic.level, random= ~I(invT-mean(invT))|week, data=data, correlation=corAR1(), method="REML", na.action=na.omit) ## Does zooplankton carbon vary with temperature? dataz <- data[(which(data$trophic.level != 'P')),] # because we don't have observations for the P treatments (we're just assuming they are 0) I don't htink we should analyze those tanks here. hist(dataz$zoo.ug.carbon.liter) hist(log(dataz$zoo.ug.carbon.liter)) hist(dataz$M.uncorrz) hist(log(dataz$M.uncorrz + 1)) plot(dataz$M.uncorrz ~dataz$invT, pch = 19, col = dataz$trophic.level) ## analysis ## determine need for random effects in the full model: modZa <- lmer(M.uncorrz ~ 1 + I(invT-mean(invT))*trophic.level + (I(invT-mean(invT))|week), data=dataz, REML = FALSE, na.action=na.omit) modZb <- lmer(M.uncorrz ~ 1 + I(invT-mean(invT))*trophic.level + 1|week, data=dataz, REML = FALSE, na.action=na.omit) anova(modZa, modZb) modZ4 <- lmer(zoo.ug.carbon.liter ~ 1 + I(invT-mean(invT))*trophic.level, random= ~1|week, data=dataz, method="REML", correlation = corAR1(), na.action=na.omit) modZ4 <- lme(zoo.ug.carbon.liter ~ 1 + I(invT-mean(invT))*trophic.level, random= ~1|week, data=dataz, method="REML", correlation = corAR1(), na.action=na.omit) sizes <- read.csv("sizes") View(sizes) attach(sizes) ### load libraries, define variables and add columns k <- 8.617342*10^-5 # eV/K sizes$invT <-1/((sizes$Temp + 273)*k) sizes$invT <- as.numeric(as.character(sizes$invT)) sizes <- sizes[sizes$week >='4',] hist(sizes$size) hist(log(sizes$size)) hist(log(sizes$size+1)) hist(log(sizes$size+0.1)) #slightly left skweded distribution ## analysis ## determine need for random effects in the full model: modsizea <- lmer(size ~ 1 + I(invT-mean(invT))*trophic.level + (I(invT-mean(invT))|week), data=sizes, REML = FALSE, na.action=na.omit) modsizeb <- lmer(size ~ 1 + I(invT-mean(invT))*trophic.level + 1|week, data=sizes, REML = FALSE, na.action=na.omit) anova(modZa, modZb) modsize0 <- lmer(size ~ 1 +(1|week), data=sizes, REML = TRUE, na.action=na.omit) modsize1<-lmer(size~1+I(invT-mean(invT)) + (1|week), data=sizes, REML = FALSE, na.action=na.omit) modsize2<-lmer(size~1+I(invT-mean(invT)) + treatment + (1|week), data=sizes, REML = FALSE, na.action=na.omit) modsize4<-lmer(size~1+I(invT-mean(invT)) * treatment + (1|week), data=sizes, REML = FALSE, na.action=na.omit) model.sel(modsize0, modsize1, modsize2, modsize4) anova(modsize2, modsize4) anova(modsize0, modsize1) summary(modsize4) # for model fitting: modsize4r <- lmer(sizes ~ 1 + I(invT-mean(invT))*trophic.level + (I(invT-mean(invT))|week), data=sizes, REML = TRUE, na.action=na.omit) summary(modsize4r) confint(modsize4r) #Test for temporal autocorrelation with nlme modsize4r<-lme(size~1+I(invT-mean(invT)) + treatment, random=~1|week, data=sizes, method="REML", correlation=corAR1(), na.action=na.omit) ### Does cyanobacteria to phytoplankton abundance change over time? data4 <- read.csv("Cyano.abundance") k <- 8.617342*10^-5 # eV/K data4$invT <- 1/((data4$average.temp + 273)*k) data4$invT <- as.numeric(as.character(data4$invT)) data4$Cyano.PP.ratio <- (data4$Cyano.abundance/data4$Sum.Phyto) data4 <- data4[data4$week >= '4',] hist(data4$Cyano.PP.ratio) hist(log(data4$Cyano.PP.ratio)) ### data prep complete ### ANALYSIS ## Does Cyano.PP vary with temperature? ## figures on invT plot(log(data4$Cyano.PP.ratio)~data4$Tank, pch = 19, col = data4$trophic.level) plot(log(data4$Cyano.PP.ratio)~data4$week, pch = 19, col = data4$trophic.level) plot(log(data4$Cyano.PP.ratio)~I(data4$invT-mean(data4$invT)), pch = 19, col = data4$trophic.level) ## analysis ## determine need for random effects in the full model: modCyanoa <- lmer(Cyano.PP.ratio~1+invT*trophic.level + invT|week, data=data4, REML = FALSE, na.action=na.omit) modCyanob<-lmer(Cyano.PP.ratio~1+invT*trophic.level + 1|week, data=data4, REML = FALSE, na.action=na.omit) anova(modCyanoa, modCyanob) modCyano0 <- lmer(log(Cyano.PP.ratio) ~ 1 + (1|week), data=data4, REML = FALSE, na.action=na.omit) modCyano1 <- lmer(log(Cyano.PP.ratio) ~ 1 + I(invT-mean(invT)) + (1|week), data=data4, REML = FALSE, na.action=na.omit) modCyano2 <- lmer(log(Cyano.PP.ratio) ~ 1 + I(invT-mean(invT)) + trophic.level + (1|week), data=data4, REML = FALSE, na.action=na.omit) modCyano4 <- lmer(log(Cyano.PP.ratio) ~ 1 + I(invT-mean(invT))*trophic.level + (1|week), data=data4, REML = FALSE, na.action=na.omit) model.sel(modCyano0, modCyano1, modCyano2, modCyano4) anova(modCyano1, modCyano0) anova(modCyano0, modCyano2) anova(modCyano2, modCyano4) # for model fitting: modCyano0r <- lmer(log(Cyano.PP.ratio) ~ 1 + (1|week), data=data4, REML = TRUE, na.action=na.omit) modCyano1r <- lmer(log(Cyano.PP.ratio) ~ 1 + I(invT-mean(invT)) + (1|week), data=data4, REML = TRUE, na.action=na.omit) modCyano2r <- lmer(log(Cyano.PP.ratio) ~ 1 + I(invT-mean(invT)) + trophic.level + (1|week), data=data4, REML = TRUE, na.action=na.omit) m.avg <- model.avg(modCyano0r, modCyano1r, modCyano2r) summary(m.avg) confint(m.avg) week8 <- subset(data4, week == 8) hist(week8$Cyano.PP.ratio) hist(log(week8$Cyano.PP.ratio)) hist(log(week8$Cyano.PP.ratio+1)) modCyano1r <- lm(log(Cyano.PP.ratio) ~ 1 + invT*trophic.level, data=week8, na.action=na.omit) summary(modCyano1r)