getting % cover data fully ready for analysis
library(tidyverse)## -- Attaching packages --------------------------------------- tidyverse 1.3.2 --
## v ggplot2 3.3.5 v purrr 0.3.4
## v tibble 3.1.6 v dplyr 1.0.7
## v tidyr 1.2.0 v stringr 1.4.0
## v readr 2.1.2 v forcats 0.5.1
## -- Conflicts ------------------------------------------ tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()library(nlme)##
## Attaching package: 'nlme'
##
## The following object is masked from 'package:dplyr':
##
## collapselibrary(lattice)#read in plot metadata
plot_meta <- read.csv("shrub_count/plotmeta_avgBigShrubs.csv", header = T)
#clean it up, select only the columsn we need
plot_meta$shrub_count<- plot_meta$avg_big_count
plot_meta <- plot_meta %>% select(PlotID, Trt, SH_GR, shrub_count)read in % cover data, cleaning and organizing
#first all species
cover_T_2021 <-read.csv("final_outputs/2021cover_T_by_species.csv", header = T)
cover_T_2022 <-read.csv("final_outputs/2022cover_T_by_species.csv", header = T)
cover_2021 <-read.csv("final_outputs/2021cover_Plot_by_species.csv", header = T)
cover_2022 <-read.csv("final_outputs/2022cover_Plot_by_species.csv", header = T)
#now the same for functional groups
PFXcover_T_2021 <-read.csv("final_outputs/PFX_cover_T_2021.csv", header = T)
PFXcover_T_2022 <-read.csv("final_outputs/PFX_cover_T_2022.csv", header = T)
PFXcover_2021 <-read.csv("final_outputs/PFX_cover_plot_2021.csv", header = T)
PFXcover_2022 <-read.csv("final_outputs/PFX_cover_plot_2022.csv", header = T)#add metadata into PFC plot cover
PFXcover_T_2021$PlotID <- PFXcover_T_2021$PrimaryKey
PFXcover_T_2021$Yr <- 2021
PFXcover21 <-left_join(PFXcover_T_2021, plot_meta, by="PlotID")
PFXcover_T_2022$PlotID <- PFXcover_T_2022$PrimaryKey
PFXcover_T_2022$Yr <- 2022
PFXcover22 <-left_join(PFXcover_T_2022, plot_meta, by="PlotID")
PFXcover <- bind_rows(PFXcover21, PFXcover22)make long version
#get cover columns
fx_groups = c("AF_I", "AF_N", "AG_I", "PF_N", "AG_I",
"PG_N", "PG_I", "PF_N", "SH_N", "PF_I")
pfx <- pivot_longer(PFXcover, cols = fx_groups,
names_to ="fx", values_to = "cov")## Note: Using an external vector in selections is ambiguous.
## i Use `all_of(fx_groups)` instead of `fx_groups` to silence this message.
## i See <https://tidyselect.r-lib.org/reference/faq-external-vector.html>.
## This message is displayed once per session.#replace na cover values with 0
pfx$cov[is.na(pfx$cov)] <- 0data we are working with
# this is all by TRANSECT
# three transects per plot
head(PFXcover)## PrimaryKey LineKey AF_I AF_N AG_I GAP OTHER PF_N PG_I
## 1 101C 101C_15 5 1.666667 5.000000 36.66667 100 5.000000 0
## 2 101C 101C_2 0 3.333333 1.666667 18.33333 100 5.000000 0
## 3 101C 101C_28 0 0.000000 0.000000 35.00000 100 5.000000 0
## 4 101T 101T_15 0 1.666667 1.666667 30.00000 100 6.666667 0
## 5 101T 101T_2 0 0.000000 0.000000 36.66667 100 6.666667 0
## 6 101T 101T_28 0 0.000000 1.666667 31.66667 100 3.333333 0
## PG_N SH_N PlotID Yr Trt SH_GR shrub_count PF_I
## 1 23.33333 36.66667 101C 2021 C HSHG 57.33333 NA
## 2 25.00000 63.33333 101C 2021 C HSHG 57.33333 NA
## 3 15.00000 51.66667 101C 2021 C HSHG 57.33333 NA
## 4 31.66667 43.33333 101T 2021 T HSHG 61.00000 NA
## 5 30.00000 38.33333 101T 2021 T HSHG 61.00000 NA
## 6 23.33333 48.33333 101T 2021 T HSHG 61.00000 NA# in this dataset, all functional groups are their own column
head(pfx)## # A tibble: 6 x 11
## PrimaryKey LineKey GAP OTHER PlotID Yr Trt SH_GR shrub_count fx
## <chr> <chr> <dbl> <dbl> <chr> <dbl> <chr> <chr> <dbl> <chr>
## 1 101C 101C_15 36.7 100 101C 2021 C HSHG 57.3 AF_I
## 2 101C 101C_15 36.7 100 101C 2021 C HSHG 57.3 AF_N
## 3 101C 101C_15 36.7 100 101C 2021 C HSHG 57.3 AG_I
## 4 101C 101C_15 36.7 100 101C 2021 C HSHG 57.3 PF_N
## 5 101C 101C_15 36.7 100 101C 2021 C HSHG 57.3 PG_N
## 6 101C 101C_15 36.7 100 101C 2021 C HSHG 57.3 PG_I
## # ... with 1 more variable: cov <dbl># in this dataset, all plant functional groups are combined into one column
# fx is functional group
# cov is % coverplot info included in the datasets
# Yr is sample year, either 2021 or 2022
PFXcover$Yr <- factor(PFXcover$Yr)
pfx$Yr <- factor(pfx$Yr)
# SH_GR is plant community strata
PFXcover$SH_GR <- factor(PFXcover$SH_GR)
pfx$SH_GR <- factor(pfx$SH_GR)
# Trt is C (control) or T (treatment) with herbicide indaziflam
PFXcover$Trt <- factor(PFXcover$Trt)
pfx$Trt <- factor(pfx$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
research question 1
Did indaziflam treatment areas have lower annual grass cover?
exploratory figure
#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
make box plots following the book, residauls
boxplot(AG_I ~ factor(SH_GR),
varwidth = TRUE, xlab = "Comm type",
main = "Boxplot of AG cover by comm type",
ylab = "AG cover ", data = PFXcover)
#make simple linear regression to plot the residuals
m1 <- lm(AG_I ~ SH_GR * Trt, 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")
par(op)
clearly there is a violation of heterogeneity of variance
Fitting a GLS model using the fixed variance structure
M.lm<-gls(AG_I~SH_GR * Trt,data=PFXcover)
# Fitting a GLS model using the VarIdent variance structure
# need to do it this way because Trt and SH_GR are both categorical
vf2 <- varIdent(form= ~ 1|SH_GR)
M.gls2 <- gls(AG_I~SH_GR * Trt,
weights=vf2, data =PFXcover)
vf1 <- varIdent(form= ~ 1|Trt)
M.gls1 <- gls(AG_I~SH_GR * Trt,
weights=vf1, data =PFXcover)
#combination of two different categorical variables
vf3 <- varComb(varIdent(form= ~ 1 | Trt) ,
varIdent(form= ~ 1|SH_GR))
M.gls3 <- gls(AG_I~SH_GR * Trt,
weights=vf3, data =PFXcover)
# I think these are all of our options because we are stuck with categorical datacompare the models
anova(M.lm,M.gls1, M.gls2, M.gls3)## Model df AIC BIC logLik Test L.Ratio p-value
## M.lm 1 9 1465.988 1494.316 -723.9943
## M.gls1 2 10 1451.283 1482.758 -715.6414 1 vs 2 16.70577 <.0001
## M.gls2 3 12 1460.593 1498.363 -718.2968 2 vs 3 5.31077 0.0703
## M.gls3 4 13 1421.596 1462.514 -697.7982 3 vs 4 40.99715 <.0001AIC(M.lm,M.gls1, M.gls2, M.gls3)## df AIC
## M.lm 9 1465.989
## M.gls1 10 1451.283
## M.gls2 12 1460.594
## M.gls3 13 1421.596# best is combination
summary(M.gls3)## Generalized least squares fit by REML
## Model: AG_I ~ SH_GR * Trt
## Data: PFXcover
## AIC BIC logLik
## 1421.596 1462.514 -697.7982
##
## Combination of variance functions:
## Structure: Different standard deviations per stratum
## Formula: ~1 | Trt
## Parameter estimates:
## C T
## 1.0000000 0.4561665
## Structure: Different standard deviations per stratum
## Formula: ~1 | SH_GR
## Parameter estimates:
## HSHG HSLG LSHG LSLG
## 1.0000000 0.8241168 1.0544734 2.2088636
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 20.724638 3.771443 5.495148 0.0000
## SH_GRHSLG -3.478261 4.887139 -0.711717 0.4776
## SH_GRLSHG 1.428140 5.420371 0.263476 0.7925
## SH_GRLSLG 12.608696 11.324149 1.113434 0.2671
## TrtT -16.141304 4.081059 -3.955175 0.0001
## SH_GRHSLG:TrtT 4.619565 5.322179 0.867984 0.3866
## SH_GRLSHG:TrtT 1.544082 5.859592 0.263514 0.7925
## SH_GRLSLG:TrtT 2.585749 12.361658 0.209175 0.8346
##
## Correlation:
## (Intr) SH_GRHSLG SH_GRLSHG SH_GRLSLG TrtT SH_GRHSLG:
## SH_GRHSLG -0.772
## SH_GRLSHG -0.696 0.537
## SH_GRLSLG -0.333 0.257 0.232
## TrtT -0.924 0.713 0.643 0.308
## SH_GRHSLG:TrtT 0.709 -0.918 -0.493 -0.236 -0.767
## SH_GRLSHG:TrtT 0.644 -0.497 -0.925 -0.214 -0.696 0.534
## SH_GRLSLG:TrtT 0.305 -0.235 -0.212 -0.916 -0.330 0.253
## SH_GRLSHG:
## SH_GRHSLG
## SH_GRLSHG
## SH_GRLSLG
## TrtT
## SH_GRHSLG:TrtT
## SH_GRLSHG:TrtT
## SH_GRLSLG:TrtT 0.230
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -1.1615050 -0.6768662 -0.3535020 0.3920319 3.9207158
##
## Residual standard error: 18.08721
## Degrees of freedom: 180 total; 172 residualhow do i interp this??
plot the new residuals
Model validation with noralized (standardized) residuals - more useful
E2 <- resid(M.gls3, type = "normalized")
coplot(E2 ~ SH_GR | Trt, data = PFXcover,
ylab = "Normalised residuals")
anova(M.gls3)## Denom. DF: 172
## numDF F-value p-value
## (Intercept) 1 115.97245 <.0001
## SH_GR 3 3.98364 0.0089
## Trt 1 39.61030 <.0001
## SH_GR:Trt 3 0.26989 0.8471# treatment and SH GR are significant
# Interaction is not
summary(M.gls3)## Generalized least squares fit by REML
## Model: AG_I ~ SH_GR * Trt
## Data: PFXcover
## AIC BIC logLik
## 1421.596 1462.514 -697.7982
##
## Combination of variance functions:
## Structure: Different standard deviations per stratum
## Formula: ~1 | Trt
## Parameter estimates:
## C T
## 1.0000000 0.4561665
## Structure: Different standard deviations per stratum
## Formula: ~1 | SH_GR
## Parameter estimates:
## HSHG HSLG LSHG LSLG
## 1.0000000 0.8241168 1.0544734 2.2088636
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 20.724638 3.771443 5.495148 0.0000
## SH_GRHSLG -3.478261 4.887139 -0.711717 0.4776
## SH_GRLSHG 1.428140 5.420371 0.263476 0.7925
## SH_GRLSLG 12.608696 11.324149 1.113434 0.2671
## TrtT -16.141304 4.081059 -3.955175 0.0001
## SH_GRHSLG:TrtT 4.619565 5.322179 0.867984 0.3866
## SH_GRLSHG:TrtT 1.544082 5.859592 0.263514 0.7925
## SH_GRLSLG:TrtT 2.585749 12.361658 0.209175 0.8346
##
## Correlation:
## (Intr) SH_GRHSLG SH_GRLSHG SH_GRLSLG TrtT SH_GRHSLG:
## SH_GRHSLG -0.772
## SH_GRLSHG -0.696 0.537
## SH_GRLSLG -0.333 0.257 0.232
## TrtT -0.924 0.713 0.643 0.308
## SH_GRHSLG:TrtT 0.709 -0.918 -0.493 -0.236 -0.767
## SH_GRLSHG:TrtT 0.644 -0.497 -0.925 -0.214 -0.696 0.534
## SH_GRLSLG:TrtT 0.305 -0.235 -0.212 -0.916 -0.330 0.253
## SH_GRLSHG:
## SH_GRHSLG
## SH_GRLSHG
## SH_GRLSLG
## TrtT
## SH_GRHSLG:TrtT
## SH_GRLSHG:TrtT
## SH_GRLSLG:TrtT 0.230
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -1.1615050 -0.6768662 -0.3535020 0.3920319 3.9207158
##
## Residual standard error: 18.08721
## Degrees of freedom: 180 total; 172 residual# how do i interp this??
# Treated areas were -16.14 % lower AG cover than C (p-value .0001)
# what about for pl comm type?how do I update ggplot box plot with these new values?
is this good enough? What more needs to be done?
#without year
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).
include year
#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)
# clearly there is a violation of heterogeneity of variance
#use from squid code - week 1
# Fitting a GLS model using the fixed variance structure
M.lm<-gls(AG_I~SH_GR * Trt * Yr,data=PFXcover)
# vf1Fixed<-varFixed(~SH_GR)
# M.gls1<-gls(AG_I~SH_GR*Trt,
# weights=vf1Fixed,data=PFXcover)
# anova(M.lm,M.gls1)
# AIC(M.lm,M.gls1)
# Fitting a GLS model using the VarIdent variance structure
# need to do it this way because Trt and SH_GR are both categorical
vf2 <- varIdent(form= ~ 1|SH_GR)
M.gls2 <- gls(AG_I~SH_GR * Trt * Yr,
weights=vf2, data =PFXcover)
vf1 <- varIdent(form= ~ 1|Trt)
M.gls1 <- gls(AG_I~SH_GR * Trt * Yr,
weights=vf1, data =PFXcover)
vf4 <- varIdent(form= ~ 1|Yr)
M.gls4 <- gls(AG_I~SH_GR * Trt * Yr,
weights=vf4, data =PFXcover)
#combination of two different categorical variables
vf3 <- varComb(varIdent(form= ~ 1 | Trt) ,
varIdent(form= ~ 1|SH_GR))
M.gls3 <- gls(AG_I~SH_GR * Trt * Yr,
weights=vf3, data =PFXcover)
#combination of three
vf5 <- varComb(varIdent(form= ~ 1 | Trt) ,
varIdent(form= ~ 1|SH_GR),
varIdent(form= ~ 1 | Yr))
M.gls5 <- gls(AG_I~SH_GR * Trt * Yr,
weights=vf5, data =PFXcover)
# I think these are all of our options because we are stuck with categorical data
#compare the models
anova(M.lm,M.gls1, M.gls2, M.gls4, M.gls3, M.gls5)## Model df AIC BIC logLik Test L.Ratio p-value
## M.lm 1 17 1370.723 1423.420 -668.3613
## M.gls1 2 18 1368.453 1424.251 -666.2267 1 vs 2 4.26925 0.0388
## M.gls2 3 20 1345.442 1407.439 -652.7209 2 vs 3 27.01158 <.0001
## M.gls4 4 18 1294.063 1349.861 -629.0315 3 vs 4 47.37876 <.0001
## M.gls3 5 21 1326.173 1391.271 -642.0867 4 vs 5 26.11042 <.0001
## M.gls5 6 22 1239.666 1307.863 -597.8330 5 vs 6 88.50739 <.0001AIC(M.lm,M.gls1, M.gls2, M.gls4, M.gls3, M.gls5)## df AIC
## M.lm 17 1370.723
## M.gls1 18 1368.453
## M.gls2 20 1345.442
## M.gls4 18 1294.063
## M.gls3 21 1326.173
## M.gls5 22 1239.666# best is combination
summary(M.gls5)## Generalized least squares fit by REML
## Model: AG_I ~ SH_GR * Trt * Yr
## Data: PFXcover
## AIC BIC logLik
## 1239.666 1307.863 -597.833
##
## Combination of variance functions:
## Structure: Different standard deviations per stratum
## Formula: ~1 | Trt
## Parameter estimates:
## C T
## 1.0000000 0.5610365
## Structure: Different standard deviations per stratum
## Formula: ~1 | SH_GR
## Parameter estimates:
## HSHG HSLG LSHG LSLG
## 1.0000000 0.9060003 1.2159202 2.6495066
## Structure: Different standard deviations per stratum
## Formula: ~1 | Yr
## Parameter estimates:
## 2021 2022
## 1.000000 3.217567
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 4.259259 1.704266 2.499176 0.0134
## SH_GRHSLG 1.574074 2.166246 0.726637 0.4685
## SH_GRLSHG 4.768519 2.474912 1.926743 0.0557
## SH_GRLSLG 24.629630 5.786936 4.256074 0.0000
## TrtT -2.703704 1.858241 -1.454980 0.1476
## Yr2022 27.050265 4.715407 5.736570 0.0000
## SH_GRHSLG:TrtT -0.629630 2.409147 -0.261350 0.7942
## SH_GRLSHG:TrtT -1.768519 2.735822 -0.646430 0.5189
## SH_GRLSLG:TrtT -14.796296 6.607865 -2.239195 0.0265
## SH_GRHSLG:Yr2022 -3.186628 6.649647 -0.479218 0.6324
## SH_GRLSHG:Yr2022 -0.800265 7.668015 -0.104364 0.9170
## SH_GRLSLG:Yr2022 -19.272487 17.037929 -1.131152 0.2596
## TrtT:Yr2022 -20.528897 5.416287 -3.790216 0.0002
## SH_GRHSLG:TrtT:Yr2022 3.407685 7.631398 0.446535 0.6558
## SH_GRLSHG:TrtT:Yr2022 0.278897 8.666404 0.032181 0.9744
## SH_GRLSLG:TrtT:Yr2022 26.732601 19.325131 1.383308 0.1685
##
## Correlation:
## (Intr) SH_GRHSLG SH_GRLSHG SH_GRLSLG TrtT Yr2022
## SH_GRHSLG -0.787
## SH_GRLSHG -0.689 0.542
## SH_GRLSLG -0.295 0.232 0.203
## TrtT -0.917 0.722 0.632 0.270
## Yr2022 -0.361 0.284 0.249 0.106 0.331
## SH_GRHSLG:TrtT 0.707 -0.899 -0.487 -0.208 -0.771 -0.256
## SH_GRLSHG:TrtT 0.623 -0.490 -0.905 -0.183 -0.679 -0.225
## SH_GRLSLG:TrtT 0.258 -0.203 -0.178 -0.876 -0.281 -0.093
## SH_GRHSLG:Yr2022 0.256 -0.326 -0.176 -0.075 -0.235 -0.709
## SH_GRLSHG:Yr2022 0.222 -0.175 -0.323 -0.065 -0.204 -0.615
## SH_GRLSLG:Yr2022 0.100 -0.079 -0.069 -0.340 -0.092 -0.277
## TrtT:Yr2022 0.315 -0.248 -0.217 -0.093 -0.343 -0.871
## SH_GRHSLG:TrtT:Yr2022 -0.223 0.284 0.154 0.066 0.243 0.618
## SH_GRLSHG:TrtT:Yr2022 -0.197 0.155 0.286 0.058 0.214 0.544
## SH_GRLSLG:TrtT:Yr2022 -0.088 0.069 0.061 0.299 0.096 0.244
## SH_GRHSLG:TrT SH_GRLSHG:TrT SH_GRLSLG:TrT SH_GRHSLG:Y
## SH_GRHSLG
## SH_GRLSHG
## SH_GRLSLG
## TrtT
## Yr2022
## SH_GRHSLG:TrtT
## SH_GRLSHG:TrtT 0.524
## SH_GRLSLG:TrtT 0.217 0.191
## SH_GRHSLG:Yr2022 0.293 0.160 0.066
## SH_GRLSHG:Yr2022 0.157 0.292 0.057 0.436
## SH_GRLSLG:Yr2022 0.071 0.062 0.297 0.196
## TrtT:Yr2022 0.265 0.233 0.096 0.617
## SH_GRHSLG:TrtT:Yr2022 -0.316 -0.165 -0.068 -0.871
## SH_GRLSHG:TrtT:Yr2022 -0.165 -0.316 -0.060 -0.386
## SH_GRLSLG:TrtT:Yr2022 -0.074 -0.065 -0.342 -0.173
## SH_GRLSHG:Y SH_GRLSLG:Y TT:Y20 SH_GRHSLG:TT:
## SH_GRHSLG
## SH_GRLSHG
## SH_GRLSLG
## TrtT
## Yr2022
## SH_GRHSLG:TrtT
## SH_GRLSHG:TrtT
## SH_GRLSLG:TrtT
## SH_GRHSLG:Yr2022
## SH_GRLSHG:Yr2022
## SH_GRLSLG:Yr2022 0.170
## TrtT:Yr2022 0.535 0.241
## SH_GRHSLG:TrtT:Yr2022 -0.380 -0.171 -0.710
## SH_GRLSHG:TrtT:Yr2022 -0.885 -0.151 -0.625 0.444
## SH_GRLSLG:TrtT:Yr2022 -0.150 -0.882 -0.280 0.199
## SH_GRLSHG:TT:
## SH_GRHSLG
## SH_GRLSHG
## SH_GRLSLG
## TrtT
## Yr2022
## SH_GRHSLG:TrtT
## SH_GRLSHG:TrtT
## SH_GRLSLG:TrtT
## SH_GRHSLG:Yr2022
## SH_GRLSHG:Yr2022
## SH_GRLSLG:Yr2022
## TrtT:Yr2022
## SH_GRHSLG:TrtT:Yr2022
## SH_GRLSHG:TrtT:Yr2022
## SH_GRLSLG:TrtT:Yr2022 0.175
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -1.7006013 -0.8282781 -0.1895590 0.6197662 2.9438894
##
## Residual standard error: 5.112797
## Degrees of freedom: 180 total; 164 residual# how do i interp this??
# # plot the new residuals
# #Model validation with noralized (standardized) residuals - more useful
# E2 <- resid(M.gls5, type = "normalized")
# coplot(E2 ~ SH_GR | Trt | Yr, data = PFXcover,
# ylab = "Normalised residuals")
anova(M.gls5)## Denom. DF: 164
## numDF F-value p-value
## (Intercept) 1 143.77729 <.0001
## SH_GR 3 11.95298 <.0001
## Trt 1 35.72872 <.0001
## Yr 1 62.98024 <.0001
## SH_GR:Trt 3 1.34609 0.2613
## SH_GR:Yr 3 0.06909 0.9763
## Trt:Yr 1 31.53423 <.0001
## SH_GR:Trt:Yr 3 0.68162 0.5645# treatment, SH GR, and year are significant
# Interaction terms are also sig
summary(M.gls3)## Generalized least squares fit by REML
## Model: AG_I ~ SH_GR * Trt * Yr
## Data: PFXcover
## AIC BIC logLik
## 1326.173 1391.271 -642.0867
##
## Combination of variance functions:
## Structure: Different standard deviations per stratum
## Formula: ~1 | Trt
## Parameter estimates:
## C T
## 1.0000000 0.5596931
## Structure: Different standard deviations per stratum
## Formula: ~1 | SH_GR
## Parameter estimates:
## HSHG HSLG LSHG LSLG
## 1.0000000 0.7412059 1.0074461 2.3772402
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 4.259259 4.642423 0.917465 0.3602
## SH_GRHSLG 1.574074 5.516557 0.285336 0.7757
## SH_GRLSHG 4.768519 6.160988 0.773986 0.4401
## SH_GRLSLG 24.629630 14.291505 1.723375 0.0867
## TrtT -2.703704 5.059929 -0.534336 0.5938
## Yr2022 27.050265 5.950378 4.545974 0.0000
## SH_GRHSLG:TrtT -0.629630 6.104508 -0.103142 0.9180
## SH_GRLSHG:TrtT -1.768519 6.791164 -0.260415 0.7949
## SH_GRLSLG:TrtT -14.796296 16.295039 -0.908025 0.3652
## SH_GRHSLG:Yr2022 -3.186628 7.346763 -0.433746 0.6650
## SH_GRLSHG:Yr2022 -0.800265 8.259441 -0.096891 0.9229
## SH_GRLSLG:Yr2022 -19.272487 18.844716 -1.022700 0.3080
## TrtT:Yr2022 -20.528897 6.643177 -3.090223 0.0024
## SH_GRHSLG:TrtT:Yr2022 3.407685 8.277454 0.411683 0.6811
## SH_GRLSHG:TrtT:Yr2022 0.278897 9.228535 0.030221 0.9759
## SH_GRLSLG:TrtT:Yr2022 26.732601 21.429701 1.247456 0.2140
##
## Correlation:
## (Intr) SH_GRHSLG SH_GRLSHG SH_GRLSLG TrtT Yr2022
## SH_GRHSLG -0.842
## SH_GRLSHG -0.754 0.634
## SH_GRLSLG -0.325 0.273 0.245
## TrtT -0.917 0.772 0.691 0.298
## Yr2022 -0.780 0.657 0.588 0.253 0.716
## SH_GRHSLG:TrtT 0.760 -0.904 -0.573 -0.247 -0.829 -0.593
## SH_GRLSHG:TrtT 0.684 -0.575 -0.907 -0.222 -0.745 -0.533
## SH_GRLSLG:TrtT 0.285 -0.240 -0.215 -0.877 -0.311 -0.222
## SH_GRHSLG:Yr2022 0.632 -0.751 -0.476 -0.205 -0.580 -0.810
## SH_GRLSHG:Yr2022 0.562 -0.473 -0.746 -0.183 -0.516 -0.720
## SH_GRLSLG:Yr2022 0.246 -0.207 -0.186 -0.758 -0.226 -0.316
## TrtT:Yr2022 0.699 -0.588 -0.527 -0.227 -0.762 -0.896
## SH_GRHSLG:TrtT:Yr2022 -0.561 0.666 0.423 0.182 0.611 0.719
## SH_GRLSHG:TrtT:Yr2022 -0.503 0.423 0.668 0.163 0.548 0.645
## SH_GRLSLG:TrtT:Yr2022 -0.217 0.182 0.163 0.667 0.236 0.278
## SH_GRHSLG:TrT SH_GRLSHG:TrT SH_GRLSLG:TrT SH_GRHSLG:Y
## SH_GRHSLG
## SH_GRLSHG
## SH_GRLSLG
## TrtT
## Yr2022
## SH_GRHSLG:TrtT
## SH_GRLSHG:TrtT 0.618
## SH_GRLSLG:TrtT 0.257 0.231
## SH_GRHSLG:Yr2022 0.679 0.432 0.180
## SH_GRLSHG:Yr2022 0.427 0.677 0.160 0.584
## SH_GRLSLG:Yr2022 0.187 0.168 0.665 0.256
## TrtT:Yr2022 0.631 0.568 0.237 0.725
## SH_GRHSLG:TrtT:Yr2022 -0.737 -0.455 -0.190 -0.888
## SH_GRLSHG:TrtT:Yr2022 -0.454 -0.736 -0.170 -0.522
## SH_GRLSLG:TrtT:Yr2022 -0.196 -0.176 -0.760 -0.225
## SH_GRLSHG:Y SH_GRLSLG:Y TT:Y20 SH_GRHSLG:TT:
## SH_GRHSLG
## SH_GRLSHG
## SH_GRLSLG
## TrtT
## Yr2022
## SH_GRHSLG:TrtT
## SH_GRLSHG:TrtT
## SH_GRLSLG:TrtT
## SH_GRHSLG:Yr2022
## SH_GRLSHG:Yr2022
## SH_GRLSLG:Yr2022 0.227
## TrtT:Yr2022 0.645 0.283
## SH_GRHSLG:TrtT:Yr2022 -0.518 -0.227 -0.803
## SH_GRLSHG:TrtT:Yr2022 -0.895 -0.204 -0.720 0.578
## SH_GRLSLG:TrtT:Yr2022 -0.200 -0.879 -0.310 0.249
## SH_GRLSHG:TT:
## SH_GRHSLG
## SH_GRLSHG
## SH_GRLSLG
## TrtT
## Yr2022
## SH_GRHSLG:TrtT
## SH_GRLSHG:TrtT
## SH_GRLSLG:TrtT
## SH_GRHSLG:Yr2022
## SH_GRLSHG:Yr2022
## SH_GRLSLG:Yr2022
## TrtT:Yr2022
## SH_GRHSLG:TrtT:Yr2022
## SH_GRLSHG:TrtT:Yr2022
## SH_GRLSLG:TrtT:Yr2022 0.223
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -2.2309757 -0.5650827 -0.1442328 0.4326984 3.9616657
##
## Residual standard error: 13.92727
## Degrees of freedom: 180 total; 164 residual# how do i interp this??
# Treated areas were -16.14 % lower AG cover than C (p-value .0001)
# what about for pl comm type?#annual grass cover
#with year
# maybe should not include year because did not model variance
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 (%)")+
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")
ggsh## Warning: Removed 1 rows containing non-finite values (stat_boxplot).
research question 2
How does shrub abundance impact AG control?
Is there a relationship between AG cover and shrub count? Does this differ by treatment?
we would predict AG control in areas of high shrub abundance to be lower, since herbicide is limited by canopy
note that the shrub count is from… count data am I breaking some rule here or is that okay because just using as a proxy for shrub abundance
remove the one plot that does not have shrub count data
PFXcover$PlotID <- factor(PFXcover$PlotID)
summary(PFXcover$PlotID)## 101C 101T 102C 102T 103C 103T 104C 104T 105C 105T 201C 201T 202C 202T 203C 203T
## 6 6 2 6 6 6 6 5 3 5 6 6 6 6 6 5
## 204C 204T 301C 301T 302C 302T 303C 303T 304C 304T 305T 401C 401T 402C 402T 403C
## 5 6 6 6 6 6 6 6 6 6 6 6 6 6 6 2
## 403T
## 3PFXcover <- PFXcover %>% filter(PlotID != "403T")
PFXcover <- PFXcover %>% filter(PlotID != "403C")
summary(PFXcover$PlotID)## 101C 101T 102C 102T 103C 103T 104C 104T 105C 105T 201C 201T 202C 202T 203C 203T
## 6 6 2 6 6 6 6 5 3 5 6 6 6 6 6 5
## 204C 204T 301C 301T 302C 302T 303C 303T 304C 304T 305T 401C 401T 402C 402T 403C
## 5 6 6 6 6 6 6 6 6 6 6 6 6 6 6 0
## 403T
## 0#exploratory plotting
#annual grass cover
ggag<-ggplot(data=PFXcover, aes(x=shrub_count, y=AG_I, color=Trt))+
geom_point()+ geom_smooth(method = "lm")+
labs(x="Shrub abundance", y="Annual Grass cover (%)")
ggag## `geom_smooth()` using formula 'y ~ x'
#make simple linear regression to plot the residuals
m1 <- lm(AG_I ~ shrub_count * Trt, 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$shrub_count, xlab = "shrub count", ylab = "residuals")
plot(resid(m1) ~ PFXcover$Trt, xlab = "Trt", ylab = "residuals")
par(op)
## there is a megaphoning with the residuals! something is up.
try different variance structures
M.lm<-gls(AG_I~shrub_count * Trt,data=PFXcover)
# factor the treatment with VarIdent
vf1 <- varIdent(form= ~ 1|Trt)
M.gls1 <- gls(AG_I~shrub_count * Trt,
weights=vf1, data =PFXcover)
anova(M.lm,M.gls1)## Model df AIC BIC logLik Test L.Ratio p-value
## M.lm 1 5 1414.117 1429.825 -702.0583
## M.gls1 2 6 1358.121 1376.971 -673.0607 1 vs 2 57.9952 <.0001AIC(M.lm,M.gls1)## df AIC
## M.lm 5 1414.117
## M.gls1 6 1358.121# Fitting a GLS model using the fixed variance structure
#vf1Fixed<-varFixed(~shrub_count)
#M.gls2<-gls(AG_I~shrub_count * Trt,
# weights=vf1Fixed,data=PFXcover)
# error
# Error in glsEstimate(object, control = control) :
# 'Calloc' could not allocate memory (4294967295 of 8 bytes)
E <- resid(M.lm)
coplot(AG_I~shrub_count|Trt,data=PFXcover)
# vf3 <- varPower(form =~ shrub_count)
# M.gls3 <- gls(AG_I~shrub_count * Trt,
# weights = vf3,data=PFXcover)
# #memory error
#
# vf4 <- varPower(form=~ shrub_count | Trt)
# M.gls4<-gls(AG_I~shrub_count * Trt, data = PFXcover,
# weights = vf4)
# #Error in eigen(val, only.values = TRUE) :
# # infinite or missing values in 'x'
vf5 <- varExp(form =~ shrub_count)
M.gls5 <- gls(AG_I~shrub_count * Trt,
weights = vf5, data = PFXcover)
# vf6<-varConstPower(form =~ shrub_count)
# M.gls6<-gls(AG_I~shrub_count * Trt,
# weights = vf6, data = PFXcover)
# # false convergence
#
# vf7 <- varConstPower(form =~ shrub_count | Trt)
# M.gls7<-gls(AG_I~shrub_count * Trt,
# weights = vf7, data = PFXcover)
# # false convergence
#combination1
vf8 <- varComb(varIdent(form= ~ 1 | Trt) ,
varExp(form =~ shrub_count) )
M.gls8<-gls(AG_I~shrub_count * Trt,
weights = vf8, data = PFXcover)
#combination1
# vf8 <- varComb(varIdent(form= ~ 1 | Trt) ,
# varPower(form =~ shrub_count) )
# M.gls8<-gls(AG_I~shrub_count * Trt,
# weights = vf8, data = PFXcover)
# # false convergenceanova(M.lm,M.gls1,M.gls5,M.gls8)## Model df AIC BIC logLik Test L.Ratio p-value
## M.lm 1 5 1414.117 1429.825 -702.0583
## M.gls1 2 6 1358.121 1376.971 -673.0607 1 vs 2 57.99520 <.0001
## M.gls5 3 6 1401.351 1420.200 -694.6752
## M.gls8 4 7 1349.908 1371.900 -667.9542 3 vs 4 53.44216 <.0001#combination is best # plot the new residuals
#Model validation with noralized (standardized) residuals - more useful
E2 <- resid(M.gls8, type = "normalized")
coplot(E2 ~ shrub_count | Trt, data = PFXcover,
ylab = "Normalised residuals")
anova(M.gls8)## Denom. DF: 171
## numDF F-value p-value
## (Intercept) 1 125.72626 <.0001
## shrub_count 1 16.83583 0.0001
## Trt 1 45.42509 <.0001
## shrub_count:Trt 1 8.14366 0.0049# treatment and shrub count are significant
# interaction between them is also sig
summary(M.gls8)## Generalized least squares fit by REML
## Model: AG_I ~ shrub_count * Trt
## Data: PFXcover
## AIC BIC logLik
## 1349.908 1371.9 -667.9542
##
## Combination of variance functions:
## Structure: Different standard deviations per stratum
## Formula: ~1 | Trt
## Parameter estimates:
## C T
## 1.0000000 0.4373727
## Structure: Exponential of variance covariate
## Formula: ~shrub_count
## Parameter estimates:
## expon
## -0.005189258
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 32.79365 4.067228 8.062899 0.0000
## shrub_count -0.17587 0.047566 -3.697315 0.0003
## TrtT -24.07114 4.302982 -5.594060 0.0000
## shrub_count:TrtT 0.14106 0.049430 2.853709 0.0049
##
## Correlation:
## (Intr) shrb_c TrtT
## shrub_count -0.897
## TrtT -0.945 0.847
## shrub_count:TrtT 0.863 -0.962 -0.892
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -1.3738406 -0.6895391 -0.3448178 0.5970194 4.0689205
##
## Residual standard error: 23.40372
## Degrees of freedom: 175 total; 171 residual# how do i interp this??how do I update ggplot box plot with these new values?
is this good enough? What more needs to be done?
# plot these data - better version
# remove the one plot that does not have shrub count data
#annual grass cover
ggag<-ggplot(data=PFXcover, aes(x=shrub_count, y=AG_I, color=Trt))+
geom_point()+ geom_smooth(method = "lm", se=FALSE)+
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,150), expand=c(0,0))+
scale_color_brewer(palette="Dark2", name = "Indaziflam application",
labels = c("Control", "Treatment"))+
theme(legend.position = "top")
ggag## `geom_smooth()` using formula 'y ~ x'## Warning: Removed 5 rows containing non-finite values (stat_smooth).## Warning: Removed 5 rows containing missing values (geom_point).
does shrub count differ by Trt
no
m2 <- lm(shrub_count~Trt, data = PFXcover)
m2##
## Call:
## lm(formula = shrub_count ~ Trt, data = PFXcover)
##
## Coefficients:
## (Intercept) TrtT
## 62.7866 -0.8009summary(m2)##
## Call:
## lm(formula = shrub_count ~ Trt, data = PFXcover)
##
## Residuals:
## Min 1Q Median 3Q Max
## -62.79 -33.72 -7.12 37.01 117.35
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 62.7866 4.7899 13.108 <2e-16 ***
## TrtT -0.8009 6.5706 -0.122 0.903
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 43.37 on 173 degrees of freedom
## Multiple R-squared: 8.588e-05, Adjusted R-squared: -0.005694
## F-statistic: 0.01486 on 1 and 173 DF, p-value: 0.9031