exec/LMER_VWC_modelstuff.R

In bmcnellis/SDEF.analysis: Data Analysis for 2016 San Dimas Experimental Forest Project

library(emmeans)
library(lmerTest)
library(lubridate)
library(reshape2)
library(MuMIn)
library(phia)

setwd("D:/R_Workspace/SDEF/2017_6_Analysis_scripts/Processed_VWC_Data/")
ad <- read.csv("Adenostoma_VWC_withreps_2019_04_01.csv")
g <- read.csv("Grass_VWC_withreps_2019_04_01.csv")
b <- read.csv("Bare_VWC_withreps_2019_04_01.csv")

ad <- ad[-c(1:31), ]

ad$measurement <- rep(1:159)
ad$treatment <- rep("ad", times = 159)
g$treatment <- rep("g", times = 159)
b$treatment <- rep("b", times = 159)

### restructure date
ad$time <- strptime(ad$time, format = "%m/%d/%Y %H:%M")
ad$time <- as.POSIXct(ad$time, format = "%m/%d/%Y %H:%M")
ad$month <- month(ad$time)
ad$month <- as.character(ad$month)

g$time <- strptime(g$time, format = "%m/%d/%Y %H:%M")
g$time <- as.POSIXct(g$time, format = "%m/%d/%Y %H:%M")
g$month <- month(g$time)
g$month <- as.character(g$month)

b$time <- strptime(b$time, format = "%m/%d/%Y %H:%M")
b$time <- as.POSIXct(b$time, format = "%m/%d/%Y %H:%M")
b$month <- month(b$time)
b$month <- as.character(b$month)

ad <- melt(ad, id.vars = c("time", "measurement", "month", "treatment"))
g <- melt(g, id.vars = c("time", "measurement", "month", "treatment"))
b <- melt(b, id.vars = c("time", "measurement", "month", "treatment"))

alldat <- rbind(ad, b, g)
colnames(alldat) <- c("time", "measurement", "month", "veg_type", "rep", "vwc")

mod <- lmer(vwc ~ veg_type + month + veg_type*month +  (1|measurement) + (1|rep) , data = alldat)
step(mod)
summary(mod)
anova(mod)
emmip(mod, veg_type ~ month )
testInteractions(mod, pairwise = "veg_type", fixed = "month", adjust = "BY")
emout <- emmeans(mod, list(pairwise ~  veg_type|month), adjust = "Tukey")
emout1 <- as.data.frame(emout$`pairwise differences of veg_type | month`)
#write.csv(emout1, file = "eemeans_vwcoutput_04_2019.csv")