getting % cover data fully ready for analysis
library(tidyverse)
library(nlme)
library(lattice)read in data
PFXcover <- read.csv("PFXcover.csv", header=T)
# plot info included in the datasets
# Yr is sample year, either 2021 or 2022
PFXcover$Yr <- factor(PFXcover$Yr)
# SH_GR is plant community strata
PFXcover$SH_GR <- factor(PFXcover$SH_GR)
# Trt is C (control) or T (treatment) with herbicide indaziflam
PFXcover$Trt <- factor(PFXcover$Trt)shrub_count is average count of shrubs > 15 cm within 2 x 30 m belt transect this is a proxy for shrub abundance I need to add in other environmental factors like soil type, aspect, and slope but this is what we are working with right now
Not that this is nested transect data. 3 observations (transects) per plot, per year for two consecutive years
research question:
Did indaziflam treatment areas have lower annual grass cover? How did this differ by year and veg structure?
exploratory figures
#facet by year
ggsh<-ggplot(data=PFXcover, aes(x=SH_GR, y=AG_I,
fill=Trt))+
geom_boxplot()+facet_grid(~Yr)+
labs(x="Plant Community Type", y="Annual Grass cover (%)")+
ylim(0,85)
ggsh
#facet by trt
ggsh<-ggplot(data=PFXcover, aes(x=SH_GR, y=AG_I,
fill=Yr))+
geom_boxplot()+facet_grid(~Trt)+
labs(x="Plant Community Type", y="Annual Grass cover (%)")+
ylim(0,85)
ggsh
plot the residuals
#make simple linear regression to plot the residuals
m1 <- lm(AG_I ~ SH_GR * Trt * Yr, data = PFXcover)
#plot the residuals
op <- par(mfrow = c(2,2), mar=c(4,4,2,2))
plot(m1, which=c(1), col=1, add.smooth=F, caption="")
plot(resid(m1) ~ PFXcover$SH_GR, xlab = "Comm type", ylab = "residuals")
plot(resid(m1) ~ PFXcover$Trt, xlab = "Trt", ylab = "residuals")
plot(resid(m1) ~ PFXcover$Yr, xlab = "Trt", ylab = "residuals")
par(op)seems like there is some violation of heterogeneity of variance
try adding sopme random effect for Plot
library(nlme)
M1<-gls(AG_I~SH_GR * Trt * Yr *shrub_count,data=PFXcover)
M2 <- lme(AG_I~SH_GR * Trt * Yr *shrub_count, data = PFXcover,
random = ~1 | SH_GR)
M3 <- lme(AG_I~SH_GR * Trt * Yr *shrub_count, data = PFXcover,
random = ~1 | PlotID)
anova(M1, M2, M3) #M3 with PlotID as random effect is best## Model df AIC BIC logLik Test L.Ratio p-value
## M1 1 33 1372.588 1471.496 -653.2942
## M2 2 34 1374.588 1476.494 -653.2942 1 vs 2 2.273737e-13 1
## M3 3 34 1336.324 1438.229 -634.1621plot(M3, col = 1)
try some different variance structures
form = AG_I ~ Trt * Yr * SH_GR * shrub_count
# try var indent for each categorical variance
M4 <- lme(form, random = ~1 | PlotID,
data = PFXcover,
weights = varIdent(form = ~1 | Trt))
M5 <- lme(form, random = ~1 | PlotID,
data = PFXcover,
weights = varIdent(form = ~1 | Yr))
M6 <- lme(form, random = ~1 | PlotID,
data = PFXcover,
weights = varIdent(form = ~1 | SH_GR))
#combination of all var indents
M7 <- lme(form, random = ~1 | PlotID,
data = PFXcover,
weights = varComb(varIdent(form = ~1 | SH_GR),
varIdent(form = ~1 | Yr),
weights = varIdent(form = ~1 | Trt)))
# from previous versions I know varExp(form =~ shrub_count) might be goog
M8 <- lme(form, random = ~1 | PlotID,
data = PFXcover,
weights = varExp(form =~ shrub_count))
#combination of all
M9 <- lme(form, random = ~1 | PlotID,
data = PFXcover,
weights = varComb(varIdent(form = ~1 | SH_GR),
varIdent(form = ~1 | Yr),
varIdent(form = ~1 | Trt),
varExp(form =~ shrub_count)))
anova(M1, M2, M3, M4, M5, M6, M7, M8, M9)## Warning in nlme::anova.lme(object = M1, M2, M3, M4, M5, M6, M7, M8, M9): fitted
## objects with different fixed effects. REML comparisons are not meaningful.## Model df AIC BIC logLik Test L.Ratio p-value
## M1 1 33 1372.588 1471.496 -653.2942
## M2 2 34 1374.588 1476.494 -653.2942 1 vs 2 0.00000 1
## M3 3 34 1336.324 1438.229 -634.1621
## M4 4 35 1312.972 1417.875 -621.4862 3 vs 4 25.35182 <.0001
## M5 5 35 1301.662 1406.564 -615.8309
## M6 6 37 1325.559 1436.456 -625.7795 5 vs 6 19.89708 <.0001
## M7 7 39 1253.111 1370.002 -587.5554 6 vs 7 76.44808 <.0001
## M8 8 35 1301.785 1406.688 -615.8927 7 vs 8 56.67450 <.0001
## M9 9 40 1239.609 1359.498 -579.8046 8 vs 9 72.17614 <.0001AIC(M1, M2, M3, M4, M5, M6, M7, M8, M9)## df AIC
## M1 33 1372.588
## M2 34 1374.588
## M3 34 1336.324
## M4 35 1312.972
## M5 35 1301.662
## M6 37 1325.559
## M7 39 1253.111
## M8 35 1301.785
## M9 40 1239.609#M9 is best AIC, but high df?
#plot new residuals
plot(M9, col = 1) # looks good
how do I include temporal autocorrelation for measuring the same plot two years in a row
is that even realistic for two years of data?
results
summary(M9)## Linear mixed-effects model fit by REML
## Data: PFXcover
## AIC BIC logLik
## 1239.609 1359.498 -579.8046
##
## Random effects:
## Formula: ~1 | PlotID
## (Intercept) Residual
## StdDev: 2.962816 9.8779
##
## Combination of variance functions:
## Structure: Different standard deviations per stratum
## Formula: ~1 | SH_GR
## Parameter estimates:
## HSHG HSLG LSHG LSLG
## 1.0000000 1.0813557 0.8126036 1.4318094
## Structure: Different standard deviations per stratum
## Formula: ~1 | Yr
## Parameter estimates:
## 2021 2022
## 1.000000 2.713881
## Structure: Different standard deviations per stratum
## Formula: ~1 | Trt
## Parameter estimates:
## C T
## 1.000000 0.601981
## Structure: Exponential of variance covariate
## Formula: ~shrub_count
## Parameter estimates:
## expon
## -0.01081643
## Fixed effects: list(form)
## Value Std.Error DF t-value p-value
## (Intercept) -4.10269 10.03030 131 -0.409030 0.6832
## TrtT 4.41608 10.82268 17 0.408039 0.6883
## Yr2022 64.43399 13.99551 131 4.603904 0.0000
## SH_GRHSLG 14.98974 11.62349 17 1.289607 0.2145
## SH_GRLSHG 24.85581 20.94309 17 1.186826 0.2516
## SH_GRLSLG 8.59556 16.93027 17 0.507703 0.6182
## shrub_count 0.09883 0.10865 17 0.909545 0.3758
## TrtT:Yr2022 -51.05726 16.55612 131 -3.083891 0.0025
## TrtT:SH_GRHSLG -7.66470 14.01727 17 -0.546804 0.5916
## TrtT:SH_GRLSHG -23.15054 21.92585 17 -1.055856 0.3058
## TrtT:SH_GRLSLG 9.33363 18.73772 17 0.498120 0.6248
## Yr2022:SH_GRHSLG -14.50900 18.22446 131 -0.796128 0.4274
## Yr2022:SH_GRLSHG 24.52296 41.45586 131 0.591544 0.5552
## Yr2022:SH_GRLSLG -21.21314 28.32230 131 -0.748991 0.4552
## TrtT:shrub_count -0.08209 0.11885 17 -0.690728 0.4991
## Yr2022:shrub_count -0.50147 0.15025 131 -3.337548 0.0011
## SH_GRHSLG:shrub_count -0.15315 0.12245 17 -1.250726 0.2280
## SH_GRLSHG:shrub_count -0.41803 0.50571 17 -0.826627 0.4199
## SH_GRLSLG:shrub_count 0.49884 0.31192 17 1.599259 0.1282
## TrtT:Yr2022:SH_GRHSLG 21.27460 21.04460 131 1.010930 0.3139
## TrtT:Yr2022:SH_GRLSHG -49.10669 43.29168 131 -1.134322 0.2587
## TrtT:Yr2022:SH_GRLSLG 29.21213 32.06246 131 0.911101 0.3639
## TrtT:Yr2022:shrub_count 0.41345 0.17909 131 2.308563 0.0225
## TrtT:SH_GRHSLG:shrub_count 0.09599 0.14096 17 0.680932 0.5051
## TrtT:SH_GRLSHG:shrub_count 0.49449 0.53639 17 0.921881 0.3695
## TrtT:SH_GRLSLG:shrub_count -1.46001 0.88312 17 -1.653236 0.1166
## Yr2022:SH_GRHSLG:shrub_count 0.22887 0.17899 131 1.278675 0.2033
## Yr2022:SH_GRLSHG:shrub_count -1.20477 1.05388 131 -1.143181 0.2550
## Yr2022:SH_GRLSLG:shrub_count -0.45911 0.56155 131 -0.817579 0.4151
## TrtT:Yr2022:SH_GRHSLG:shrub_count -0.25115 0.20734 131 -1.211292 0.2280
## TrtT:Yr2022:SH_GRLSHG:shrub_count 1.92540 1.09688 131 1.755347 0.0815
## TrtT:Yr2022:SH_GRLSLG:shrub_count -1.72238 1.71949 131 -1.001682 0.3183
## Correlation:
## (Intr) TrtT Yr2022 SH_GRHSLG SH_GRLSHG
## TrtT -0.927
## Yr2022 -0.598 0.555
## SH_GRHSLG -0.863 0.800 0.516
## SH_GRLSHG -0.479 0.444 0.287 0.413
## SH_GRLSLG -0.592 0.549 0.355 0.511 0.284
## shrub_count -0.977 0.905 0.575 0.843 0.468
## TrtT:Yr2022 0.506 -0.508 -0.845 -0.437 -0.242
## TrtT:SH_GRHSLG 0.716 -0.772 -0.428 -0.829 -0.343
## TrtT:SH_GRLSHG 0.457 -0.494 -0.274 -0.395 -0.955
## TrtT:SH_GRLSLG 0.535 -0.578 -0.320 -0.462 -0.256
## Yr2022:SH_GRHSLG 0.460 -0.426 -0.768 -0.460 -0.220
## Yr2022:SH_GRLSHG 0.202 -0.187 -0.338 -0.174 -0.306
## Yr2022:SH_GRLSLG 0.296 -0.274 -0.494 -0.255 -0.142
## TrtT:shrub_count 0.893 -0.969 -0.526 -0.771 -0.428
## Yr2022:shrub_count 0.571 -0.529 -0.970 -0.492 -0.273
## SH_GRHSLG:shrub_count 0.867 -0.803 -0.510 -0.968 -0.415
## SH_GRLSHG:shrub_count 0.210 -0.195 -0.124 -0.181 -0.952
## SH_GRLSLG:shrub_count 0.340 -0.315 -0.200 -0.294 -0.163
## TrtT:Yr2022:SH_GRHSLG -0.398 0.399 0.665 0.398 0.191
## TrtT:Yr2022:SH_GRLSHG -0.193 0.194 0.323 0.167 0.293
## TrtT:Yr2022:SH_GRLSLG -0.261 0.262 0.437 0.225 0.125
## TrtT:Yr2022:shrub_count -0.479 0.482 0.814 0.413 0.229
## TrtT:SH_GRHSLG:shrub_count -0.753 0.817 0.443 0.841 0.361
## TrtT:SH_GRLSHG:shrub_count -0.198 0.215 0.117 0.171 0.898
## TrtT:SH_GRLSLG:shrub_count -0.120 0.130 0.071 0.104 0.058
## Yr2022:SH_GRHSLG:shrub_count -0.479 0.444 0.814 0.466 0.229
## Yr2022:SH_GRLSHG:shrub_count -0.081 0.075 0.138 0.070 0.258
## Yr2022:SH_GRLSLG:shrub_count -0.153 0.141 0.259 0.132 0.073
## TrtT:Yr2022:SH_GRHSLG:shrub_count 0.413 -0.416 -0.703 -0.402 -0.198
## TrtT:Yr2022:SH_GRLSHG:shrub_count 0.078 -0.079 -0.133 -0.067 -0.248
## TrtT:Yr2022:SH_GRLSLG:shrub_count 0.050 -0.050 -0.085 -0.043 -0.024
## SH_GRLSLG shrb_c TrT:Y2022 TrT:SH_GRHSLG
## TrtT
## Yr2022
## SH_GRHSLG
## SH_GRLSHG
## SH_GRLSLG
## shrub_count 0.579
## TrtT:Yr2022 -0.300 -0.486
## TrtT:SH_GRHSLG -0.424 -0.699 0.392
## TrtT:SH_GRLSHG -0.271 -0.447 0.251 0.381
## TrtT:SH_GRLSLG -0.904 -0.523 0.293 0.446
## Yr2022:SH_GRHSLG -0.272 -0.442 0.649 0.382
## Yr2022:SH_GRLSHG -0.120 -0.194 0.285 0.145
## Yr2022:SH_GRLSLG -0.532 -0.284 0.418 0.212
## TrtT:shrub_count -0.529 -0.914 0.482 0.748
## Yr2022:shrub_count -0.338 -0.559 0.820 0.408
## SH_GRHSLG:shrub_count -0.514 -0.887 0.431 0.803
## SH_GRLSHG:shrub_count -0.124 -0.215 0.104 0.150
## SH_GRLSLG:shrub_count -0.920 -0.348 0.169 0.244
## TrtT:Yr2022:SH_GRHSLG 0.236 0.382 -0.787 -0.366
## TrtT:Yr2022:SH_GRLSHG 0.115 0.186 -0.382 -0.150
## TrtT:Yr2022:SH_GRLSLG 0.470 0.251 -0.516 -0.202
## TrtT:Yr2022:shrub_count 0.284 0.469 -0.971 -0.372
## TrtT:SH_GRHSLG:shrub_count 0.446 0.771 -0.406 -0.956
## TrtT:SH_GRLSHG:shrub_count 0.117 0.203 -0.107 -0.166
## TrtT:SH_GRLSLG:shrub_count 0.325 0.123 -0.065 -0.101
## Yr2022:SH_GRHSLG:shrub_count 0.284 0.469 -0.688 -0.386
## Yr2022:SH_GRLSHG:shrub_count 0.048 0.080 -0.117 -0.058
## Yr2022:SH_GRLSLG:shrub_count 0.452 0.150 -0.219 -0.109
## TrtT:Yr2022:SH_GRHSLG:shrub_count -0.245 -0.405 0.839 0.367
## TrtT:Yr2022:SH_GRLSHG:shrub_count -0.046 -0.077 0.159 0.061
## TrtT:Yr2022:SH_GRLSLG:shrub_count -0.148 -0.049 0.101 0.039
## TrT:SH_GRLSHG TrT:SH_GRLSLG Yr2022:SH_GRHSLG
## TrtT
## Yr2022
## SH_GRHSLG
## SH_GRLSHG
## SH_GRLSLG
## shrub_count
## TrtT:Yr2022
## TrtT:SH_GRHSLG
## TrtT:SH_GRLSHG
## TrtT:SH_GRLSLG 0.285
## Yr2022:SH_GRHSLG 0.210 0.246
## Yr2022:SH_GRLSHG 0.293 0.108 0.259
## Yr2022:SH_GRLSLG 0.135 0.481 0.379
## TrtT:shrub_count 0.478 0.559 0.404
## Yr2022:shrub_count 0.261 0.305 0.745
## SH_GRHSLG:shrub_count 0.397 0.464 0.440
## SH_GRLSHG:shrub_count 0.909 0.112 0.095
## SH_GRLSLG:shrub_count 0.156 0.831 0.154
## TrtT:Yr2022:SH_GRHSLG -0.197 -0.231 -0.866
## TrtT:Yr2022:SH_GRLSHG -0.297 -0.112 -0.248
## TrtT:Yr2022:SH_GRLSLG -0.129 -0.507 -0.335
## TrtT:Yr2022:shrub_count -0.238 -0.278 -0.625
## TrtT:SH_GRHSLG:shrub_count -0.403 -0.472 -0.382
## TrtT:SH_GRLSHG:shrub_count -0.943 -0.124 -0.089
## TrtT:SH_GRLSLG:shrub_count -0.064 -0.591 -0.054
## Yr2022:SH_GRHSLG:shrub_count -0.219 -0.256 -0.962
## Yr2022:SH_GRLSHG:shrub_count -0.246 -0.044 -0.106
## Yr2022:SH_GRLSLG:shrub_count -0.070 -0.408 -0.199
## TrtT:Yr2022:SH_GRHSLG:shrub_count 0.206 0.241 0.830
## TrtT:Yr2022:SH_GRLSHG:shrub_count 0.252 0.045 0.102
## TrtT:Yr2022:SH_GRLSLG:shrub_count 0.025 0.265 0.065
## Yr2022:SH_GRLSHG Yr2022:SH_GRLSLG TrtT:_
## TrtT
## Yr2022
## SH_GRHSLG
## SH_GRLSHG
## SH_GRLSLG
## shrub_count
## TrtT:Yr2022
## TrtT:SH_GRHSLG
## TrtT:SH_GRLSHG
## TrtT:SH_GRLSLG
## Yr2022:SH_GRHSLG
## Yr2022:SH_GRLSHG
## Yr2022:SH_GRLSLG 0.167
## TrtT:shrub_count 0.178 0.260
## Yr2022:shrub_count 0.327 0.479 0.511
## SH_GRHSLG:shrub_count 0.172 0.252 0.811
## SH_GRLSHG:shrub_count 0.272 0.061 0.196
## SH_GRLSLG:shrub_count 0.068 0.492 0.318
## TrtT:Yr2022:SH_GRHSLG -0.225 -0.329 -0.379
## TrtT:Yr2022:SH_GRLSHG -0.958 -0.160 -0.184
## TrtT:Yr2022:SH_GRLSLG -0.147 -0.883 -0.249
## TrtT:Yr2022:shrub_count -0.275 -0.402 -0.468
## TrtT:SH_GRHSLG:shrub_count -0.150 -0.219 -0.843
## TrtT:SH_GRLSHG:shrub_count -0.257 -0.058 -0.222
## TrtT:SH_GRLSLG:shrub_count -0.024 -0.174 -0.135
## Yr2022:SH_GRHSLG:shrub_count -0.275 -0.402 -0.429
## Yr2022:SH_GRLSHG:shrub_count -0.972 -0.068 -0.073
## Yr2022:SH_GRLSLG:shrub_count -0.088 -0.893 -0.137
## TrtT:Yr2022:SH_GRHSLG:shrub_count 0.237 0.347 0.404
## TrtT:Yr2022:SH_GRLSHG:shrub_count 0.934 0.066 0.076
## TrtT:Yr2022:SH_GRLSLG:shrub_count 0.029 0.292 0.049
## Y2022:_ SH_GRHSLG: SH_GRLSHG: SH_GRLSLG:
## TrtT
## Yr2022
## SH_GRHSLG
## SH_GRLSHG
## SH_GRLSLG
## shrub_count
## TrtT:Yr2022
## TrtT:SH_GRHSLG
## TrtT:SH_GRLSHG
## TrtT:SH_GRLSLG
## Yr2022:SH_GRHSLG
## Yr2022:SH_GRLSHG
## Yr2022:SH_GRLSLG
## TrtT:shrub_count
## Yr2022:shrub_count
## SH_GRHSLG:shrub_count 0.496
## SH_GRLSHG:shrub_count 0.120 0.191
## SH_GRLSLG:shrub_count 0.195 0.309 0.075
## TrtT:Yr2022:SH_GRHSLG -0.645 -0.381 -0.082 -0.133
## TrtT:Yr2022:SH_GRLSHG -0.314 -0.165 -0.261 -0.065
## TrtT:Yr2022:SH_GRLSLG -0.423 -0.223 -0.054 -0.434
## TrtT:Yr2022:shrub_count -0.839 -0.416 -0.101 -0.163
## TrtT:SH_GRHSLG:shrub_count -0.431 -0.869 -0.166 -0.269
## TrtT:SH_GRLSHG:shrub_count -0.113 -0.180 -0.943 -0.071
## TrtT:SH_GRLSLG:shrub_count -0.069 -0.109 -0.026 -0.353
## Yr2022:SH_GRHSLG:shrub_count -0.839 -0.460 -0.101 -0.163
## Yr2022:SH_GRLSHG:shrub_count -0.143 -0.071 -0.261 -0.028
## Yr2022:SH_GRLSLG:shrub_count -0.268 -0.133 -0.032 -0.483
## TrtT:Yr2022:SH_GRHSLG:shrub_count 0.725 0.397 0.087 0.141
## TrtT:Yr2022:SH_GRLSHG:shrub_count 0.137 0.068 0.251 0.027
## TrtT:Yr2022:SH_GRLSLG:shrub_count 0.087 0.043 0.010 0.158
## TrT:Y2022:SH_GRHSLG TrT:Y2022:SH_GRLSHG
## TrtT
## Yr2022
## SH_GRHSLG
## SH_GRLSHG
## SH_GRLSLG
## shrub_count
## TrtT:Yr2022
## TrtT:SH_GRHSLG
## TrtT:SH_GRLSHG
## TrtT:SH_GRLSLG
## Yr2022:SH_GRHSLG
## Yr2022:SH_GRLSHG
## Yr2022:SH_GRLSLG
## TrtT:shrub_count
## Yr2022:shrub_count
## SH_GRHSLG:shrub_count
## SH_GRLSHG:shrub_count
## SH_GRLSLG:shrub_count
## TrtT:Yr2022:SH_GRHSLG
## TrtT:Yr2022:SH_GRLSHG 0.301
## TrtT:Yr2022:SH_GRLSLG 0.406 0.197
## TrtT:Yr2022:shrub_count 0.764 0.371
## TrtT:SH_GRHSLG:shrub_count 0.363 0.155
## TrtT:SH_GRLSHG:shrub_count 0.084 0.261
## TrtT:SH_GRLSLG:shrub_count 0.051 0.025
## Yr2022:SH_GRHSLG:shrub_count 0.833 0.263
## Yr2022:SH_GRLSHG:shrub_count 0.092 0.931
## Yr2022:SH_GRLSLG:shrub_count 0.173 0.084
## TrtT:Yr2022:SH_GRHSLG:shrub_count -0.960 -0.321
## TrtT:Yr2022:SH_GRLSHG:shrub_count -0.125 -0.963
## TrtT:Yr2022:SH_GRLSLG:shrub_count -0.080 -0.039
## TrT:Y2022:SH_GRLSLG TT:Y2022:_ TT:SH_GRHSLG:
## TrtT
## Yr2022
## SH_GRHSLG
## SH_GRLSHG
## SH_GRLSLG
## shrub_count
## TrtT:Yr2022
## TrtT:SH_GRHSLG
## TrtT:SH_GRLSHG
## TrtT:SH_GRLSLG
## Yr2022:SH_GRHSLG
## Yr2022:SH_GRLSHG
## Yr2022:SH_GRLSLG
## TrtT:shrub_count
## Yr2022:shrub_count
## SH_GRHSLG:shrub_count
## SH_GRLSHG:shrub_count
## SH_GRLSLG:shrub_count
## TrtT:Yr2022:SH_GRHSLG
## TrtT:Yr2022:SH_GRLSHG
## TrtT:Yr2022:SH_GRLSLG
## TrtT:Yr2022:shrub_count 0.502
## TrtT:SH_GRHSLG:shrub_count 0.210 0.395
## TrtT:SH_GRLSHG:shrub_count 0.055 0.104 0.187
## TrtT:SH_GRLSLG:shrub_count 0.304 0.063 0.113
## Yr2022:SH_GRHSLG:shrub_count 0.355 0.704 0.399
## Yr2022:SH_GRLSHG:shrub_count 0.060 0.120 0.061
## Yr2022:SH_GRLSLG:shrub_count 0.789 0.224 0.115
## TrtT:Yr2022:SH_GRHSLG:shrub_count -0.433 -0.864 -0.378
## TrtT:Yr2022:SH_GRLSHG:shrub_count -0.082 -0.163 -0.064
## TrtT:Yr2022:SH_GRLSLG:shrub_count -0.537 -0.104 -0.041
## TT:SH_GRLSHG: TT:SH_GRLSLG: Y2022:SH_GRHSLG:
## TrtT
## Yr2022
## SH_GRHSLG
## SH_GRLSHG
## SH_GRLSLG
## shrub_count
## TrtT:Yr2022
## TrtT:SH_GRHSLG
## TrtT:SH_GRLSHG
## TrtT:SH_GRLSLG
## Yr2022:SH_GRHSLG
## Yr2022:SH_GRLSHG
## Yr2022:SH_GRLSLG
## TrtT:shrub_count
## Yr2022:shrub_count
## SH_GRHSLG:shrub_count
## SH_GRLSHG:shrub_count
## SH_GRLSLG:shrub_count
## TrtT:Yr2022:SH_GRHSLG
## TrtT:Yr2022:SH_GRLSHG
## TrtT:Yr2022:SH_GRLSLG
## TrtT:Yr2022:shrub_count
## TrtT:SH_GRHSLG:shrub_count
## TrtT:SH_GRLSHG:shrub_count
## TrtT:SH_GRLSLG:shrub_count 0.030
## Yr2022:SH_GRHSLG:shrub_count 0.095 0.058
## Yr2022:SH_GRLSHG:shrub_count 0.246 0.010 0.120
## Yr2022:SH_GRLSLG:shrub_count 0.030 0.171 0.225
## TrtT:Yr2022:SH_GRHSLG:shrub_count -0.090 -0.054 -0.863
## TrtT:Yr2022:SH_GRLSHG:shrub_count -0.255 -0.010 -0.115
## TrtT:Yr2022:SH_GRLSLG:shrub_count -0.011 -0.423 -0.073
## Y2022:SH_GRLSHG: Y2022:SH_GRLSLG:
## TrtT
## Yr2022
## SH_GRHSLG
## SH_GRLSHG
## SH_GRLSLG
## shrub_count
## TrtT:Yr2022
## TrtT:SH_GRHSLG
## TrtT:SH_GRLSHG
## TrtT:SH_GRLSLG
## Yr2022:SH_GRHSLG
## Yr2022:SH_GRLSHG
## Yr2022:SH_GRLSLG
## TrtT:shrub_count
## Yr2022:shrub_count
## SH_GRHSLG:shrub_count
## SH_GRLSHG:shrub_count
## SH_GRLSLG:shrub_count
## TrtT:Yr2022:SH_GRHSLG
## TrtT:Yr2022:SH_GRLSHG
## TrtT:Yr2022:SH_GRLSLG
## TrtT:Yr2022:shrub_count
## TrtT:SH_GRHSLG:shrub_count
## TrtT:SH_GRLSHG:shrub_count
## TrtT:SH_GRLSLG:shrub_count
## Yr2022:SH_GRHSLG:shrub_count
## Yr2022:SH_GRLSHG:shrub_count
## Yr2022:SH_GRLSLG:shrub_count 0.038
## TrtT:Yr2022:SH_GRHSLG:shrub_count -0.103 -0.194
## TrtT:Yr2022:SH_GRLSHG:shrub_count -0.961 -0.037
## TrtT:Yr2022:SH_GRLSLG:shrub_count -0.012 -0.327
## TT:Y2022:SH_GRHSLG: TT:Y2022:SH_GRLSHG:
## TrtT
## Yr2022
## SH_GRHSLG
## SH_GRLSHG
## SH_GRLSLG
## shrub_count
## TrtT:Yr2022
## TrtT:SH_GRHSLG
## TrtT:SH_GRLSHG
## TrtT:SH_GRLSLG
## Yr2022:SH_GRHSLG
## Yr2022:SH_GRLSHG
## Yr2022:SH_GRLSLG
## TrtT:shrub_count
## Yr2022:shrub_count
## SH_GRHSLG:shrub_count
## SH_GRLSHG:shrub_count
## SH_GRLSLG:shrub_count
## TrtT:Yr2022:SH_GRHSLG
## TrtT:Yr2022:SH_GRLSHG
## TrtT:Yr2022:SH_GRLSLG
## TrtT:Yr2022:shrub_count
## TrtT:SH_GRHSLG:shrub_count
## TrtT:SH_GRLSHG:shrub_count
## TrtT:SH_GRLSLG:shrub_count
## Yr2022:SH_GRHSLG:shrub_count
## Yr2022:SH_GRLSHG:shrub_count
## Yr2022:SH_GRLSLG:shrub_count
## TrtT:Yr2022:SH_GRHSLG:shrub_count
## TrtT:Yr2022:SH_GRLSHG:shrub_count 0.141
## TrtT:Yr2022:SH_GRLSLG:shrub_count 0.090 0.017
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -2.095745536 -0.589855757 0.009866805 0.562309428 2.671504522
##
## Number of Observations: 180
## Number of Groups: 33plot(M9)
ggsh<-ggplot(data=PFXcover, aes(x=SH_GR, y=AG_I,
fill=Trt))+
geom_boxplot()+
labs(x="Plant Community Type", y="Annual Grass cover (%)")+
scale_y_continuous(limits=c(0,75), expand=c(0,0))+
scale_fill_brewer(palette="Dark2", name = "Indaziflam application",
labels = c("Control", "Treatment"))+
theme_classic()+ theme(legend.position = "top")
ggsh## Warning: Removed 1 rows containing non-finite values (stat_boxplot).
#annual grass cover
#with year
# maybe should not include year because did not model variance
ggyr<-ggplot(data=PFXcover, aes(x=SH_GR, y=AG_I,
fill=Trt))+
geom_boxplot()+facet_grid(~Yr)+
labs(x="Plant Community Type", y="Annual Grass cover (%)")+
scale_y_continuous(limits=c(0,70), expand=c(0,0))+
scale_fill_brewer(palette="Dark2", name = "Indaziflam application",
labels = c("Control", "Treatment"))+
theme_classic()+
theme(legend.position = "top")
ggyr## Warning: Removed 1 rows containing non-finite values (stat_boxplot).
how do I update ggplot box plot with these new values?
is this good enough? What more needs to be done?
#without trend lines
ggag<-ggplot(data=PFXcover, aes(x=shrub_count, y=AG_I, color=Trt))+
geom_point(size=3)+
labs(x="Shrub abundance", y="Annual Grass cover (%)")+
theme_classic()+
scale_y_continuous(limits=c(0,70), expand=c(0,0))+
scale_x_continuous(limits=c(0,185), expand=c(0,0))+
scale_color_brewer(palette="Dark2", name = "Indaziflam application",
labels = c("Control", "Treatment"))+
theme(legend.position = "top")
ggag## Warning: Removed 1 rows containing missing values (geom_point).
#with trend lines
ggag<-ggag + geom_smooth(method = "lm", se=FALSE)
ggag## `geom_smooth()` using formula 'y ~ x'## Warning: Removed 1 rows containing non-finite values (stat_smooth).
## Removed 1 rows containing missing values (geom_point).
ggag<-ggag + facet_wrap(~Yr)
ggag## `geom_smooth()` using formula 'y ~ x'## Warning: Removed 1 rows containing non-finite values (stat_smooth).
## Removed 1 rows containing missing values (geom_point).## Warning: Removed 3 rows containing missing values (geom_smooth).
try without SH_GR
M1<-gls(AG_I~ Trt * Yr *shrub_count,data=PFXcover)
M2 <- lme(AG_I~Trt * Yr *shrub_count, data = PFXcover,
random = ~1 | SH_GR)
M3 <- lme(AG_I~Trt * Yr *shrub_count, data = PFXcover,
random = ~1 | PlotID)
anova(M1, M2, M3)## Model df AIC BIC logLik Test L.Ratio p-value
## M1 1 9 1444.372 1472.699 -713.1858
## M2 2 10 1440.669 1472.143 -710.3343 1 vs 2 5.703037 0.0169
## M3 3 10 1411.852 1443.327 -695.9262plot(M3, col = 1)
form = AG_I ~ Trt * Yr * shrub_count
# try var indent for each categorical variance
M4 <- lme(form, random = ~1 | PlotID,
data = PFXcover,
weights = varIdent(form = ~1 | Trt))
M5 <- lme(form, random = ~1 | PlotID,
data = PFXcover,
weights = varIdent(form = ~1 | Yr))
#combination of all var indents
M7 <- lme(form, random = ~1 | PlotID,
data = PFXcover,
weights = varComb(varIdent(form = ~1 | Yr),
weights = varIdent(form = ~1 | Trt)))
# from previous versions I know varExp(form =~ shrub_count) might be goog
M8 <- lme(form, random = ~1 | PlotID,
data = PFXcover,
weights = varExp(form =~ shrub_count))
#combination of all
M9 <- lme(form, random = ~1 | PlotID,
data = PFXcover,
weights = varComb(varIdent(form = ~1 | Yr),
varIdent(form = ~1 | Trt),
varExp(form =~ shrub_count)))
anova(M1, M2, M3, M4, M5, M7, M8, M9)## Model df AIC BIC logLik Test L.Ratio p-value
## M1 1 9 1444.372 1472.699 -713.1858
## M2 2 10 1440.669 1472.143 -710.3343 1 vs 2 5.70304 0.0169
## M3 3 10 1411.852 1443.327 -695.9262
## M4 4 11 1384.357 1418.979 -681.1785 3 vs 4 29.49536 <.0001
## M5 5 11 1376.023 1410.645 -677.0115
## M7 6 12 1356.569 1394.339 -666.2845 5 vs 6 21.45393 <.0001
## M8 7 11 1376.850 1411.473 -677.4251 6 vs 7 22.28116 <.0001
## M9 8 13 1315.275 1356.192 -644.6374 7 vs 8 65.57524 <.0001AIC(M1, M2, M3, M4, M5,M7, M8, M9)## df AIC
## M1 9 1444.372
## M2 10 1440.669
## M3 10 1411.852
## M4 11 1384.357
## M5 11 1376.023
## M7 12 1356.569
## M8 11 1376.850
## M9 13 1315.275#M9 is still the best AIC
#plot new residuals
plot(M9, col = 1) # looks good
summary(M9)## Linear mixed-effects model fit by REML
## Data: PFXcover
## AIC BIC logLik
## 1315.275 1356.192 -644.6374
##
## Random effects:
## Formula: ~1 | PlotID
## (Intercept) Residual
## StdDev: 5.686884 11.72833
##
## Combination of variance functions:
## Structure: Different standard deviations per stratum
## Formula: ~1 | Yr
## Parameter estimates:
## 2021 2022
## 1.000000 2.708524
## Structure: Different standard deviations per stratum
## Formula: ~1 | Trt
## Parameter estimates:
## C T
## 1.0000000 0.5672473
## Structure: Exponential of variance covariate
## Formula: ~shrub_count
## Parameter estimates:
## expon
## -0.01119821
## Fixed effects: list(form)
## Value Std.Error DF t-value p-value
## (Intercept) 16.471818 4.100006 143 4.017510 0.0001
## TrtT -9.205751 4.816449 29 -1.911315 0.0659
## Yr2022 30.910386 5.733851 143 5.390859 0.0000
## shrub_count -0.095173 0.052241 29 -1.821806 0.0788
## TrtT:Yr2022 -19.338587 6.314109 143 -3.062758 0.0026
## TrtT:shrub_count 0.049808 0.060140 29 0.828194 0.4143
## Yr2022:shrub_count -0.139789 0.056731 143 -2.464060 0.0149
## TrtT:Yr2022:shrub_count 0.083034 0.060298 143 1.377063 0.1706
## Correlation:
## (Intr) TrtT Yr2022 shrb_c TrT:Y2022 TrtT:_ Y2022:
## TrtT -0.851
## Yr2022 -0.269 0.229
## shrub_count -0.894 0.761 0.205
## TrtT:Yr2022 0.245 -0.243 -0.908 -0.187
## TrtT:shrub_count 0.777 -0.868 -0.178 -0.869 0.183
## Yr2022:shrub_count 0.246 -0.210 -0.924 -0.208 0.839 0.180
## TrtT:Yr2022:shrub_count -0.232 0.223 0.869 0.195 -0.920 -0.187 -0.941
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -2.1113599 -0.7288874 -0.1214656 0.5662227 2.9186315
##
## Number of Observations: 180
## Number of Groups: 33is this better with less df?
try without SH_GR + Yr, plan on including Year later for temporal autocorrelation
M1<-gls(AG_I~ Trt * shrub_count,data=PFXcover)
M2 <- lme(AG_I~Trt * shrub_count, data = PFXcover,
random = ~1 | SH_GR)
M3 <- lme(AG_I~Trt * shrub_count, data = PFXcover,
random = ~1 | PlotID)
anova(M1, M2, M3)## Model df AIC BIC logLik Test L.Ratio p-value
## M1 1 5 1497.823 1513.675 -743.9113
## M2 2 6 1495.882 1514.905 -741.9412 1 vs 2 3.940205 0.0471
## M3 3 6 1481.708 1500.730 -734.8538plot(M3, col = 1)
form = AG_I ~ Trt * shrub_count
# try var indent for each categorical variance
M4 <- lme(form, random = ~1 | PlotID,
data = PFXcover,
weights = varIdent(form = ~1 | Trt))
# from previous versions I know varExp(form =~ shrub_count) might be goog
M8 <- lme(form, random = ~1 | PlotID,
data = PFXcover,
weights = varExp(form =~ shrub_count))
#combination of all
M9 <- lme(form, random = ~1 | PlotID,
data = PFXcover,
weights = varComb(varIdent(form = ~1 | Trt),
varExp(form =~ shrub_count)))
anova(M1, M2, M3, M4, M8, M9)## Model df AIC BIC logLik Test L.Ratio p-value
## M1 1 5 1497.823 1513.675 -743.9113
## M2 2 6 1495.882 1514.905 -741.9412 1 vs 2 3.94020 0.0471
## M3 3 6 1481.708 1500.730 -734.8538
## M4 4 7 1432.048 1454.242 -709.0241 3 vs 4 51.65935 <.0001
## M8 5 7 1460.096 1482.289 -723.0478
## M9 6 8 1417.581 1442.945 -700.7904 5 vs 6 44.51490 <.0001AIC(M1, M2, M3, M4, M8, M9)## df AIC
## M1 5 1497.823
## M2 6 1495.882
## M3 6 1481.708
## M4 7 1432.048
## M8 7 1460.096
## M9 8 1417.581#M9 is still the best AIC
#plot new residuals
plot(M9, col = 1) # looks good
M9## Linear mixed-effects model fit by REML
## Data: PFXcover
## Log-restricted-likelihood: -700.7904
## Fixed: form
## (Intercept) TrtT shrub_count TrtT:shrub_count
## 34.58174520 -20.18261393 -0.18690249 0.09960943
##
## Random effects:
## Formula: ~1 | PlotID
## (Intercept) Residual
## StdDev: 10.40072 23.80342
##
## Combination of variance functions:
## Structure: Different standard deviations per stratum
## Formula: ~1 | Trt
## Parameter estimates:
## C T
## 1.0000000 0.4243784
## Structure: Exponential of variance covariate
## Formula: ~shrub_count
## Parameter estimates:
## expon
## -0.006628385
## Number of Observations: 180
## Number of Groups: 33summary(M9)## Linear mixed-effects model fit by REML
## Data: PFXcover
## AIC BIC logLik
## 1417.581 1442.945 -700.7904
##
## Random effects:
## Formula: ~1 | PlotID
## (Intercept) Residual
## StdDev: 10.40072 23.80342
##
## Combination of variance functions:
## Structure: Different standard deviations per stratum
## Formula: ~1 | Trt
## Parameter estimates:
## C T
## 1.0000000 0.4243784
## Structure: Exponential of variance covariate
## Formula: ~shrub_count
## Parameter estimates:
## expon
## -0.006628385
## Fixed effects: list(form)
## Value Std.Error DF t-value p-value
## (Intercept) 34.58175 6.486396 147 5.331427 0.0000
## TrtT -20.18261 7.690430 29 -2.624380 0.0137
## shrub_count -0.18690 0.087729 29 -2.130462 0.0417
## TrtT:shrub_count 0.09961 0.101649 29 0.979936 0.3352
## Correlation:
## (Intr) TrtT shrb_c
## TrtT -0.843
## shrub_count -0.868 0.732
## TrtT:shrub_count 0.749 -0.841 -0.863
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -1.9227494 -0.7017430 -0.1351015 0.4993967 2.9737045
##
## Number of Observations: 180
## Number of Groups: 33