############################################################################### # ENSC 418 # Mike Rutherford # 25 October 2017 # # Statistical analysis ############################################################################### # Cran R # http://www.r-project.org/ # R environment and base packages # Text editor to R (Notepad ++) # R libraries (packages): C:\Program Files\R\R-3.0.3\library # Introduction to R: http://cran.r-project.org/doc/manuals/r-release/R-intro.pdf # R reference card : http://cran.r-project.org/doc/contrib/Short-refcard.pdf # R statistical environment # Script based interface # Everything in R is a object # Object: data structure, attributes, methods # R requires minimal object declaration: object type will change with content # RStudio # http://www.rstudio.com/ # R IDE (integrated development environment) # Script editor, R interface (build configuration), debuging, project management, version control # RStudio components # Script editor, R command line (console), Environment (+ history), Files (+plots, packages, help) ##################################### # Using RStudio ##################################### # File, New Project # File, New File, R Script # Controlling the location of data and output # Explicit pathway description dataPath <- "C:/workspace/ENSC-418/" dataPath <- "H:/workspace/ENSC-428/" dataReadIn <- read.csv(dataPath,header = TRUE) # Implicit using the project workspace (working directory) # Rstudio, Sessions, Set Working Directory getwd() # Determine what the path of the current working directory setwd(dataPath) # Sets the path to the specified working directory ###################################### # Analysis libraries ###################################### #install.packages("ggplot2") # Good package for plotting data library(ggplot2) #install.packages("agricolae") # Package for performing Tukey HSD library(agricolae) #install.packages("nlme") # Package to perform linear mixed effects model analysis library(nlme) #install.packages("gridExtra") # Package to create multiple figures on one plot l#ibrary(gridExtra) #install.packages("lsmeans") # Package for doing pairwise tests wtih lme results #library(lsmeans) ###################################### # Read in data, look at structure, and identify "factors" in design ###################################### dat <- read.csv(paste(dataPath,"plot_data.csv",sep=""),header=TRUE) names(dat) str(dat) summary(dat) dat$biosolid <- factor(dat$biosolid) dat$ash <- factor(dat$ash) ###################################### # Plot the distribution of the responce variables ###################################### n.hist <- ggplot(dat,aes(n)) + geom_histogram(col="black",fill="grey") n.hist p.hist <- ggplot(dat,aes(p)) + geom_histogram(col="black",fill="grey") p.hist ### Plot the distribution of the remaining responce variables ggplot(dat,aes(copper)) + geom_histogram(col="black",fill="grey") ggplot(dat,aes(whc)) + geom_histogram(col="black",fill="grey") ggplot(dat,aes(respiration)) + geom_histogram(col="black",fill="grey") ######################################## # Analysis: Ash addition ######################################## # Plot the responce variable available P against ash addition p.ash.point <- ggplot(dat,aes(ash,p)) + geom_point() p.ash.point p.ash.box <- ggplot(dat,aes(ash,p)) + geom_boxplot() p.ash.box # Does ash addition impact the response variable? p.ash.lm <- lm(dat$p ~ dat$ash) p.ash.lm summary(p.ash.lm) anova(p.ash.lm) ### Does ash influence avaialble P in this data? # Plot the responce variable available N against ash addition n.ash.point <- ggplot(dat,aes(ash,n)) + geom_point() n.ash.point n.ash.box <- ggplot(dat,aes(ash,n)) + geom_boxplot() n.ash.box # Does ash addition impact the response variable? n.ash.lm <- lm(dat$n ~ dat$ash) n.ash.lm summary(n.ash.lm) anova(n.ash.lm) ### Does ash influence avaialble N in this data? ### Evaluate ash impacts on "copper", "whc", and "respiration" ######################################## # Analysis: Biosolid addition ######################################## # Plot the responce variables against biosolid addition p.biosolid.point <- ggplot(dat,aes(biosolid,p)) + geom_point() p.biosolid.point p.biosolid.box <- ggplot(dat,aes(biosolid,p)) + geom_boxplot() p.biosolid.box # Does ash addition impact the response variable? p.biosolid.lm <- lm(dat$p ~ dat$biosolid) p.biosolid.lm summary(p.biosolid.lm) anova(p.biosolid.lm) # Does Biosolid addition influence water holding capacity whc.biosolid.point <- ggplot(dat,aes(biosolid,whc)) + geom_point() whc.biosolid.point whc.biosolid.box <- ggplot(dat,aes(biosolid,whc)) + geom_boxplot() whc.biosolid.box # Does ash addition impact the response variable? whc.biosolid.lm <- lm(dat$whc ~ dat$biosolid) whc.biosolid.lm summary(whc.biosolid.lm) anova(whc.biosolid.lm) # Where does the difference occur? # Tukey's HSD test whc.biosolid.lm <- lm(whc ~ biosolid, data=dat) hsd.whc.biosolid.lm <- HSD.test(whc.biosolid.lm,'biosolid') hsd.whc.biosolid.lm ### What impact does biosolid addition have? ######################################## # Analysis: Two factor ANOVA: Do Biosolid and Ash influence the responce variables? ######################################## # Plot ash and biosolid influence available P p.both.point <- ggplot(dat,aes(biosolid,p, color = ash)) + geom_point() p.both.point p.both.box <- ggplot(dat,aes(biosolid,p, fill=ash)) + geom_boxplot() p.both.box # Do ash and biosolid impact available P p.both.lm <- lm(dat$p ~ dat$ash + dat$biosolid) p.both.lm summary(p.both.lm) anova(p.both.lm) ######################################## # Analysis: Two-way ANOVA: Is there an interaction between Biosolid and Ash? ######################################## # Do ash and biosolid interact to influence respiration resp.int.box <- ggplot(dat,aes(biosolid,respiration, fill=ash)) + geom_boxplot() resp.int.box resp.int.lm <-lm(dat$respiration ~ dat$ash*dat$biosolid) resp.int.lm summary(resp.int.lm) anova(resp.int.lm) # Do ash and biosolid interact to influence whc whc.int.box <- ggplot(dat,aes(biosolid,whc, fill=ash)) + geom_boxplot() whc.int.box whc.int.lm <-lm(dat$whc ~ dat$ash*dat$biosolid) whc.int.lm summary(whc.int.lm) anova(whc.int.lm) ####################################### # Mixed-effects model ####################################### whc.int.lm <-lm(dat$whc ~ dat$ash*dat$biosolid) whc.int.lme <- lme(whc ~ ash + biosolid +ash*biosolid, data = dat, random = ~1 |block, na.action=na.omit, method = "ML") summary(whc.int.lme) anova(whc.int.lme) # Compare linear vs. mixed-effects model anova(whc.int.lm) anova(whc.int.lme)