RxP<-read.csv("RxP.csv", stringsAsFactors=T) suppressPackageStartupMessages(library(dplyr)) #Remove the offending data temp<-filter(RxP, Ind!=734 & Ind!=1024 & Ind!=1078 & Ind!=1127 & Ind!=1284) temp<-filter(temp, SVL.initial>12) temp<-filter(temp, Tail.initial<15) #Assign temp to be an object with a useful name RxP.clean<-temp #Remove temp rm(temp) RP.means<-RxP.clean %>% #First, identify the data frame group_by(Res, Pred) %>% #Next, establish grouping variables summarize(SVL.mean = mean(SVL.initial)) #Last, calc. means #Convert the tibble output to a data frame RP.means<-as.data.frame(RP.means) #This is another way to reorder a factor, using dplyr RxP.clean<-RxP.clean %>% mutate(Pred = factor(Pred, levels=c('C','NL','L'))) RP.means<-RxP.clean %>% #Define what dataset we are using group_by(Res, Pred) %>% #Set the groups summarize(SVL.mean = mean(SVL.initial), #Calculate the mean SVL.sd = sd(SVL.initial), #SD SVL.n = length(SVL.initial)) %>% #and N of SVL mutate(SVL.se = SVL.sd/sqrt(SVL.n)) #Calculate the SE before30<-RxP.clean %>% filter(Age.FromEmergence <= 30) after30<-RxP.clean %>% filter(Age.FromEmergence > 30) #RxP.clean$before30<-ifelse(RxP.clean$Age.FromEmergence<=30,"Yes","No") #Load the necessary package suppressPackageStartupMessages(library(ggplot2)) suppressPackageStartupMessages(library(cowplot)) suppressPackageStartupMessages(library(gplots)) RP.means2<-RxP.clean %>% #Define what dataset we are using group_by(Hatch, Pred) %>% #Set the groups summarize(Mass.mean = mean(Mass.final), #Calculate the mean, Mass.sd = sd(Mass.final), #Calculate the SD Mass.n = length(Mass.final)) %>% #Calculate the N mutate(Mass.se = Mass.sd/sqrt(Mass.n)) #Calculate the SE RxP.byTank<-RxP.clean %>% group_by(Block,Tank.Unique,Pred,Hatch,Res) %>% summarize(Age.DPO = mean(Age.DPO), Age.FromEmergence = mean(Age.FromEmergence), SVL.initial = mean(SVL.initial), Tail.initial = mean(Tail.initial), SVL.final = mean(SVL.final), Mass.final = mean(Mass.final), Resorb.days = mean(Resorb.days)) RxP.byTank<-as.data.frame(RxP.byTank) suppressPackageStartupMessages(library(MASS)) RxP.byTank$log.SVL.final<-log(RxP.byTank$SVL.final) RxP.byTank$log.Age.FromEmergence<-log(RxP.byTank$Age.FromEmergence) RxP.byTank$log.Age.DPO<-log(RxP.byTank$Age.DPO) lm1<-lm(log.Age.FromEmergence~Pred, data=RxP.byTank) suppressPackageStartupMessages(library(car)) suppressPackageStartupMessages(library(multcomp)) suppressPackageStartupMessages(library(emmeans)) suppressPackageStartupMessages(library(tidyverse)) lm2<-lm(log.Age.DPO~Res*Pred, data=RxP.byTank) RxP.byTank <- RxP.byTank %>% unite(col="ResPred", Pred, Res, sep="-", remove=F) lm3<-lm(log.Age.DPO~ResPred, data=RxP.byTank) ph3<-emmeans(lm3, specs="ResPred") lm4<-lm(log.SVL.final~log.Age.DPO, data=RxP.byTank) lm4.1<-lm(log.SVL.final~log.Age.DPO, data=RxP.byTank[-c(26,49,75),]) #Now let’s make an ANCOVA that looks at the #effects of categorical and continuous data together lm5<-lm(log.SVL.final~log.Age.DPO*Pred, data=RxP.byTank) lm4.newdata<-data.frame("log.Age.DPO" = seq(from=3.64, to=4.95, by=0.01)) #First create a data frame to provide values to predict() lm5.newdata<-expand.grid("log.Age.DPO"=log(35:145), Pred=c("C","NL","L")) lm5.predicted<-predict(lm5, newdata=lm5.newdata, se.fit=T) #First, exponentiate the fitted values lm5.newdata$predict.SVL<-exp(lm5.predicted$fit) #Next, calculate the upper and lower bounds of the CIs lm5.newdata$predict.SVL.CIupper<-exp(lm5.predicted$fit+ lm5.predicted$se.fit) lm5.newdata$predict.SVL.CIlower<-exp(lm5.predicted$fit- lm5.predicted$se.fit) #Lastly, exponentiate the age data so it is back #on the original scale lm5.newdata$Age.DPO<-exp(lm5.newdata$log.Age.DPO) lm.Age_Res<-lm(log(Mass.final)~log.Age.DPO*Res, data=RxP.byTank) lm.Age_Res.newdata<-expand.grid("log.Age.DPO"=log(35:145), Res=c("Lo","Hi")) lm.Age_Res.predicted<-predict(lm.Age_Res, newdata=lm.Age_Res.newdata, se.fit=T) #First, exponentiate the fitted values lm.Age_Res.newdata$predict.SVL<-exp(lm.Age_Res.predicted$fit) #Next, calculate the upper and lower bounds of the CIs lm.Age_Res.newdata$predict.SVL.CIupper<-exp(lm.Age_Res.predicted$fit+lm.Age_Res.predicted$se.fit) lm.Age_Res.newdata$predict.SVL.CIlower<-exp(lm.Age_Res.predicted$fit-lm.Age_Res.predicted$se.fit) #Lastly, exponentiate the age data so it is back on the original scale lm.Age_Res.newdata$Age.DPO<-exp(lm.Age_Res.newdata$log.Age.DPO) # lm.tail<-lm(Tail.initial~Hatch*Pred, data=RxP.byTank) # lm.resorb<-lm(Resorb.days~Tail.initial*Res, data=RxP.byTank) # lm.resorb2<-lm(log(Resorb.days)~log(Tail.initial)*Res, data=RxP.byTank) temp <- RxP.clean %>% group_by(Tank.Unique) %>% summarize(N.alive = length(Ind)) RxP.byTank$N.alive<-temp$N.alive #Remove temp rm(temp) RxP.byTank$N.dead<-50-RxP.byTank$N.alive glm.n<-glm(N.dead~Res*Pred, family="gaussian", data=RxP.byTank) glm.ln<-glm(log(N.dead)~Res*Pred, family="gaussian", data=RxP.byTank) glm.p<-glm(N.dead~Res*Pred, family="poisson", data=RxP.byTank) glm.negb<-glm.nb(N.dead~Res*Pred, data=RxP.byTank) glm.b<-glm(cbind(N.alive,N.dead)~Res*Pred, family="binomial", data=RxP.byTank) glm.negb2<-glm.nb(N.dead~Age.DPO*Pred, RxP.byTank) glm.negb3<-glm.nb(N.dead~Res+Pred, data=RxP.byTank) predicted.data<-expand.grid(Age.DPO=40:140, Pred=c("C","NL","L"), N.dead=NA) predicted.data$N.dead<-predict(glm.negb2, newdata=data.frame(predicted.data), type="response") glm1<-glm.nb(N.dead~Pred*Res*Hatch, data=RxP.byTank) glm2<-glm.nb(N.dead~Pred+Res+Hatch, data=RxP.byTank) suppressPackageStartupMessages(library(lme4)) RxP.clean$log.Age.DPO<-log(RxP.clean$Age.DPO) #Code the lmm just like a regular lm, #but with the random effect lmm1<-lmer(log.Age.DPO~Res*Pred+(1|Block/Tank.Unique), data=RxP.clean) lmm2<-lmer(log.Age.DPO~Res+Pred+(1|Block/Tank.Unique), data=RxP.clean) lmm3<-lmer(log.Age.DPO~Pred+(1|Block/Tank.Unique), data=RxP.clean) lmm4<-lmer(log.Age.DPO~Res+(1|Block/Tank.Unique), data=RxP.clean) suppressPackageStartupMessages(library(lattice)) suppressPackageStartupMessages(library(glmmTMB)) #First, convert Block to a factor RxP.byTank$Block.factor<-as.factor(RxP.byTank$Block) #Next, run the model. Note that the function and family #specifications are different. glmm7<-glmmTMB(N.dead~Pred+Res+Hatch+(1|Block.factor), family='nbinom2', data=RxP.byTank) #Create a series of models that only each differ #in one predictor from one another glmm7.1<-glmmTMB(N.dead~Pred+Res+(1|Block.factor), family='nbinom2', data=RxP.byTank) glmm7.2<-glmmTMB(N.dead~Pred+Hatch+(1|Block.factor), family='nbinom2', data=RxP.byTank) glmm7.3<-glmmTMB(N.dead~Res+Hatch+(1|Block.factor), family='nbinom2', data=RxP.byTank) glmm8<-glmmTMB(N.dead~Pred*Res*Hatch+(1|Block.factor), family='nbinom2', data=RxP.byTank) RxP.clean$log.Age.DPO<-log(RxP.clean$Age.DPO) RxP.clean$log.SVL.final<-log(RxP.clean$SVL.final) #lm6<-lm(log.Age.DPO~Res*Pred, data=RxP.clean) #Create a mixed model, but without using REML #lmm1<-lmer(log.Age.DPO~Res*Pred+(1|Block/Tank.Unique), # data=RxP.clean, REML = F) #lmm1.1<-lmer(log.Age.DPO~Pred*Res+(1|Block), # data=RxP.clean, REML = F) #lmm1<-lmer(log.SVL.final ~ Pred+(1|Block/Tank.Unique), # data=RxP.clean) #ph1<-emmeans(lmm1, specs='Pred') png('output/rplot.png')

Introduction to R

Thoughts on Proper Data Analysis

Exploratory Data Analysis and Data Summarization

Introduction to Plotting

Basic Statistical Analysis Using R

More Linear Models in R

Advanced Statistical Analysis Using R

Mixed-effects Model

Advanced Data Wrangling and Plotting

Writing Loops and Functions in R

Appendix

Conclusion

Solution: Rerun the ANCOVA Models

Solution