regeneration - initial growth data from March 2022 experiment

clear directory and load potentially useful packages

rm(list = ls())
library(ggplot2)
library(ggthemes)
library(plyr)
library(Hmisc)

## Loading required package: lattice

## Loading required package: survival

## Loading required package: Formula

## 
## Attaching package: 'Hmisc'

## The following objects are masked from 'package:plyr':
## 
##     is.discrete, summarize

## The following objects are masked from 'package:base':
## 
##     format.pval, units

library(outliers)
library(datasets)
library(lattice)
library(reshape)

## 
## Attaching package: 'reshape'

## The following objects are masked from 'package:plyr':
## 
##     rename, round_any

library(reshape2)

## 
## Attaching package: 'reshape2'

## The following objects are masked from 'package:reshape':
## 
##     colsplit, melt, recast

library(tidyr)

## 
## Attaching package: 'tidyr'

## The following object is masked from 'package:reshape2':
## 
##     smiths

## The following objects are masked from 'package:reshape':
## 
##     expand, smiths

library(nlme)
library(stats)
library(data.table)

## 
## Attaching package: 'data.table'

## The following objects are masked from 'package:reshape2':
## 
##     dcast, melt

## The following object is masked from 'package:reshape':
## 
##     melt

library(lawstat)
library(car)

## Loading required package: carData

## 
## Attaching package: 'car'

## The following object is masked from 'package:lawstat':
## 
##     levene.test

library(dplyr)

## 
## Attaching package: 'dplyr'

## The following object is masked from 'package:car':
## 
##     recode

## The following objects are masked from 'package:data.table':
## 
##     between, first, last

## The following object is masked from 'package:nlme':
## 
##     collapse

## The following object is masked from 'package:reshape':
## 
##     rename

## The following objects are masked from 'package:Hmisc':
## 
##     src, summarize

## The following objects are masked from 'package:plyr':
## 
##     arrange, count, desc, failwith, id, mutate, rename, summarise,
##     summarize

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

library(base)
library(rminer)
library(tidyverse)

## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.1 ──

## ✔ tibble  3.1.7     ✔ stringr 1.4.0
## ✔ readr   2.1.2     ✔ forcats 0.5.1
## ✔ purrr   0.3.4

## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::arrange()   masks plyr::arrange()
## ✖ dplyr::between()   masks data.table::between()
## ✖ dplyr::collapse()  masks nlme::collapse()
## ✖ purrr::compact()   masks plyr::compact()
## ✖ dplyr::count()     masks plyr::count()
## ✖ tidyr::expand()    masks reshape::expand()
## ✖ dplyr::failwith()  masks plyr::failwith()
## ✖ dplyr::filter()    masks stats::filter()
## ✖ dplyr::first()     masks data.table::first()
## ✖ dplyr::id()        masks plyr::id()
## ✖ dplyr::lag()       masks stats::lag()
## ✖ dplyr::last()      masks data.table::last()
## ✖ dplyr::mutate()    masks plyr::mutate()
## ✖ dplyr::recode()    masks car::recode()
## ✖ dplyr::rename()    masks reshape::rename(), plyr::rename()
## ✖ purrr::some()      masks car::some()
## ✖ dplyr::src()       masks Hmisc::src()
## ✖ dplyr::summarise() masks plyr::summarise()
## ✖ dplyr::summarize() masks Hmisc::summarize(), plyr::summarize()
## ✖ purrr::transpose() masks data.table::transpose()

library(lme4)

## Loading required package: Matrix

## 
## Attaching package: 'Matrix'

## The following objects are masked from 'package:tidyr':
## 
##     expand, pack, unpack

## The following object is masked from 'package:reshape':
## 
##     expand

## 
## Attaching package: 'lme4'

## The following object is masked from 'package:rminer':
## 
##     factorize

## The following object is masked from 'package:nlme':
## 
##     lmList

library(lsmeans)

## Loading required package: emmeans

## The 'lsmeans' package is now basically a front end for 'emmeans'.
## Users are encouraged to switch the rest of the way.
## See help('transition') for more information, including how to
## convert old 'lsmeans' objects and scripts to work with 'emmeans'.

library(lmerTest)

## 
## Attaching package: 'lmerTest'

## The following object is masked from 'package:lme4':
## 
##     lmer

## The following object is masked from 'package:stats':
## 
##     step

load data analysis functions

# standard error function
se <- function(x) sd(x)/sqrt(length(x))

# function to visualize the boxplot, histogram, QQ-plot, and kernel density estimate plot for residuals
normality.plots <- function(x) {
  par(mfrow=c(2,2))
  hist(residuals(x), main = "Histogram", xlab = "Values")
  boxplot(residuals(x), main = "Boxplot", ylab = "Values")
  qqnorm(residuals(x))
  qqline(residuals(x))
  plot(density(residuals(x)), main = "Kernel Density Estimate")
}

# function to visaulize the fitted vs. residuals of model and the absolute residuals vs. fitted
var.plots <- function(x) {
  par(mfrow=c(1,2))
  plot(x = fitted(x), y = residuals(x), xlab="Fitted", ylab="Residuals")
  abline(h=0)
  title("Residuals vs. Fitted")
  plot(x = fitted(x), y = abs(residuals(x)), xlab="Fitted", ylab="Absolute Residuals")
  abline(h=0)
  title("Absolute Residuals vs. Fitted")
}

# completeFun removes NA values in a specific column
completeFun <- function(data, desiredCols) {
  completeVec <- complete.cases(data[, desiredCols])
  return(data[completeVec, ])
}

load weight data

Weight<- read_csv("Weight.csv")

## Rows: 73 Columns: 9
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (3): Treatment, Block, Species
## dbl (6): ID, Initial_weight_g, Final_weight_g, Diff-weight_g, Diff_weight_mg...
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

temporary standardized change in weight by initial weight until we wax dip for a growth rate per day (assumed experiment ran 75 days)

# standardizing the difference in weight (mg) by initial weight as a temporary proxy for surface area that we will get from the wax dipping

#Weight$Stand_diff <- (Weight$Diff_weight_mg)/(Weight$Initial_weight_mg)

# standardizing mg/cm^2 by the experiment duration (91 days)
#Weight$Stand_diff <- (Weight$Stand_diff)/ 91
# units now mg cm^-2 day^-1

#set levels for treatment
Weight$Treatment<- factor(Weight$Treatment, levels = c("control", "5%", "10%"))

PLOTTING MEAN +/- SE using the above stand_diff

Different in weight not standardized by anything:

# summarizing data to plot
WeightSummary <- ddply(Weight, .(Species, Treatment), summarise,
                       'mean'=mean(Diff_weight_mg, na.rm = TRUE),
                       'se'=se(Diff_weight_mg ),
                       'N'=length(Diff_weight_mg ))

levels(WeightSummary$Treatment) #confirm level order

## [1] "control" "5%"      "10%"

dodge<-position_dodge(width=0.6) # this offsets the points so they don't overlap

ggplot(WeightSummary, aes(x = Treatment, y = mean, colour = Treatment, fill = Treatment))+
  facet_wrap(~Species)+
  geom_bar(stat="identity", position="dodge") +
  geom_errorbar(aes(ymin = mean - se, 
                    ymax = mean + se), 
                width=0, position=position_dodge(width = 0.9))+
  geom_point(data = Weight, aes(x = Treatment, y = Diff_weight_mg, colour = Treatment), shape = 21, alpha = 0.7,
             position = position_jitterdodge(jitter.width = 0.3, jitter.height=0.2, dodge.width=0.9))+
  scale_color_manual(values = c('black', 'black', 'black'), guide="none")+
  scale_fill_manual(values = c('white', 'gray73', 'gray40'), labels=c("control", "5% removal", "10% removal"))+
xlab(expression("% removal"))+
  ylab(expression(paste("Diff in weight (mg) not standardized")))+  theme_few(base_size = 12)+
  theme(axis.text.x=element_text(colour="black"))+
  theme(axis.text.y=element_text(colour="black"))

mixed effects model to test effects of treatment and species

treatment and species are fixed effects Block is a random effect (and including initial weight in the model in attempt to confirm treatment still matters after this is accounted for. Will standardized values by surface area after wax dipping)

first looking at distribution of response variable (growth rate)

evidence that the data is not normal:

par(mfrow=c(1,2))
hist(Weight$Diff_weight_mg)
qqnorm(Weight$Diff_weight_mg)
abline(0,1)

shapiro.test(Weight$Diff_weight_mg) 

## 
##  Shapiro-Wilk normality test
## 
## data:  Weight$Diff_weight_mg
## W = 0.93861, p-value = 0.001479

max(Weight$Diff_weight_mg) 

## [1] 79424

full model:

full.model <- lmer(Diff_weight_mg ~ Treatment *Species + (1|Block) + Initial_weight_mg, data = Weight, REML = TRUE)

## Warning: Some predictor variables are on very different scales: consider
## rescaling

## Warning: Some predictor variables are on very different scales: consider
## rescaling

library(lattice)
dotplot(ranef(full.model, condVar = TRUE))

## $Block

# shows the variance of growth across random effects
# doesn't look like there's too much an effect of  block except A8...

# need to check assumptions for random effects of block
# 1) residuals are independently distributed
# 2) residuals are from a normal distribution
par(mfrow=c(1,2))
qqnorm(ranef(full.model)$Block[,1])
qqline(ranef(full.model)$Block[,1])

# first need to determine random effect structure
library(lmerTest)
ranova(full.model) # looking at the LRT for the random effects

## ANOVA-like table for random-effects: Single term deletions
## 
## Model:
## Diff_weight_mg ~ Treatment + Species + Initial_weight_mg + (1 | Block) + Treatment:Species
##             npar  logLik    AIC        LRT Df Pr(>Chisq)
## <none>        12 -702.11 1428.2                         
## (1 | Block)   11 -702.11 1426.2 1.5089e-05  1     0.9969

# block is not significant (p = 0.3126)

#can probably drop block but would need to compare models with and without block
#but for now:

final.model <- lmer(Diff_weight_mg ~ Treatment * Species + (1|Block) +Initial_weight_mg, data = Weight, REML = TRUE)

## Warning: Some predictor variables are on very different scales: consider
## rescaling

## Warning: Some predictor variables are on very different scales: consider
## rescaling

# looking at residuals for normality and constant variance
normality.plots(final.model) # will need to interpret these

var.plots(final.model) # will need to interpret these

f-test for fixed effects

anova(final.model, type=3, ddf = "Kenward-Roger")

## Type III Analysis of Variance Table with Kenward-Roger's method
##                       Sum Sq    Mean Sq NumDF  DenDF F value    Pr(>F)    
## Treatment         5283522566 2641761283     2 24.576 18.9673 1.042e-05 ***
## Species           2716013594 1358006797     2 47.125  9.7514 0.0002858 ***
## Initial_weight_mg 1675470569 1675470569     1 62.806 12.0328 0.0009498 ***
## Treatment:Species 2722007991  680501998     4 47.246  4.8858 0.0022223 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

significant main and interactive effects of treatment, species and treatment x species