Code/PA_low_abundance_ANOVA_function.R
In davan690/beech-publication-wr: Small mammal dynamics in NZ beech forests
# Prediction A plot
# Oct 2019
# Anthony Davidson
# library(jtools)
## ------------------------------------------------------------------------
# getting factor levels correct for once
plot.dat.all1 <- abund.dat5 %>%
mutate(Rats = factor(Conditions, labels = c("Full", "Reduced")),
Control = factor(control, labels = c("Yes", "No")),
Valley = factor(valley, labels = c("Eglinton", "Hollyford")),
Date = as.Date(true.date),
Treatments = factor(treat.six, levels = c("hol no control rats.removed",
"hol control rats.present",
"hol control rats.removed",
"egl control rats.present",
"egl control rats.removed" ,
"hol no control rats.present")),
Prediction = NA) %>%
mutate(Prediction = ifelse(Date == "2000-08-01", "A", Prediction),
Prediction = ifelse(Date == "2001-08-01" |
Date == "2003-08-01" |
Date == "2003-08-01" |
Date == "2003-08-01" , "B", Prediction),
Prediction = ifelse(Date == "2001-08-01" |
Date == "2003-08-01" |
Date == "2003-08-01" |
Date == "2003-08-01" , "B", Prediction),
Prediction = ifelse(Date == "2001-08-01" |
Date == "2003-08-01" |
Date == "2003-08-01" |
Date == "2003-08-01" , "B", Prediction),
Prediction = ifelse(Date == "2001-08-01" |
Date == "2003-08-01" |
Date == "2003-08-01" |
Date == "2003-08-01" , "B", Prediction))
# head(abund.dat5)
low.abund.dat <- plot.dat.all1 %>%
# mutate(N = log(N) %>% # for when using out.N
filter(trip == 10 | trip == 11 | trip == 12 | trip == 13 | trip == 14 | trip == 15)
# filter(low.abund.dat, trip.no == 15)$true.date
# summary for plots
# create summarised datasets of need
low.cond.sum <- low.abund.dat %>%
group_by(Control,Valley,Conditions) %>%
summarise(N.count = n(),
N = mean(N, rm.na = TRUE),
tvalue = qt(p = 0.025, df = N.count - 1, lower.tail = FALSE),
low.se = mean(se.N),
lcl.low.tv = N - (tvalue*low.se),
ucl.low.tv = N + (tvalue*low.se),
lcl.low = N - (1.96*low.se),
ucl.low = N + (1.96*low.se))
low.con <- low.abund.dat %>%
group_by(Control) %>%
summarise(N.count = n(),
N = mean(N, rm.na = TRUE),
tvalue = qt(p = 0.025, df = N.count - 1, lower.tail = FALSE),
low.se = mean(se.N),
lcl.low.tv = N - (tvalue*low.se),
ucl.low.tv = N + (tvalue*low.se),
lcl.low = N - (1.96*low.se),
ucl.low = N + (1.96*low.se))
low.con1 <- low.con %>%
mutate(lcl.low = ifelse(lcl.low < 0, 0, lcl.low),
Conditions = NA,
valley = NA)
low.mod3 <- glm(N ~ Control + Valley + Conditions + Conditions:Control, family = "gaussian", data = low.abund.dat)
summary.low.mod3 <- summary(low.mod3)
low.mod3.fun <- anova(low.mod3, test = "F")
# low.mod3.fun
# summ(low.mod3)
# # summary
# model.summary1 <- summary(low.mod3)
# # jtools::summ(low.mod3)
#
# # parameters for extraction
# output.dens <- plot_summs(low.mod3, scale = TRUE, plot.distributions = TRUE, inner_ci_level = .9)
# # output.dens
#
# # summary
# model.summary1 <- summary(low.mod3)
# jtools::summ(low.mod3)
# parameters for extraction
output.dens <- plot_summs(low.mod3, scale = TRUE, plot.distributions = TRUE, inner_ci_level = .9)
# output.dens
# head(low.abund.dat)
low.plot.data <- low.abund.dat %>%
select(N, valley, Control, Conditions) %>%
mutate(Control = factor(Control, labels = c("Yes", "No")),
Valley = factor(valley, labels = c("Eglinton", "Hollyford")))
# plot
low.plot.out <-
ggplot(data = low.plot.data, aes(y = N, x = Control))+
# data
geom_jitter(aes(y = N, x = Control, col = Control, shape = Valley, fill = Conditions),
stroke = 1.05,
alpha = 1, size = 3, width = 0.3) +
geom_point(data = low.con1, aes(y = N, x = Control), shape = "square",col= "grey60", size = 5) +
geom_errorbar(data = low.con1, aes(ymin = lcl.low, ymax = ucl.low, width = 0),
lwd = 1,
col= "grey60") +
#labels
geom_vline(xintercept = 1.5, lty = 2) +
xlab("Rat removal grids") +
ylab(expression(paste("Abundance"," ", "(", italic(N[jt]), ")"))) +
#Aesthics
scale_color_manual(name = "Stoat removal",
values = c("black", "black")) +
scale_shape_manual(name = "Ecosystem",
values = c(24, 21)) +
scale_fill_manual(name = "Rat control",
values = c("white", "black")) +
guides(fill = "none",
col = guide_legend(override.aes = list(
shape = c(15,0), alpha = 1
)),
shape = guide_legend(override.aes = list(
shape = c(24, 21), size = 4
))) +
#theme
theme_classic()+
theme_bw() +
theme(strip.background = element_blank(),
strip.text.y = element_blank(),
plot.title = element_text(hjust = 0, size=24, family = "Times", color="black", margin = margin(t = 10, b = 10)),
plot.subtitle=element_text(size=16, face="italic", color="black"),
legend.position = "right",
legend.key = element_blank(),
legend.background = element_rect(fill="white", size=1),
legend.key.size=unit(1,"cm"),
legend.text = element_text(colour = "black", size =16, family = "Times"),
legend.title = element_text(colour = "black", size =16, family = "Times"),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.spacing = unit(2, "lines"),
panel.border = element_blank(),
axis.title.y = element_text(colour = "black",size =20, family = "Times", angle = 90),
axis.title.x = element_text(colour = "black", size =20, family = "Times"),
axis.text.y=element_text(colour = "black",size = 20, family = "Times"),
axis.text.x = element_text(colour = "black", size =20, family = "Times"),
axis.line.y = element_line(size = 1),
axis.line.x = element_blank(),
axis.ticks.x = element_blank(),
strip.text = element_text(face="bold",colour = "black",size =14, family = "Times"))
low.plot.out
### Shortned code
ggplot(data = low.plot.data, aes(y = N, x = Control))+
# data
geom_jitter(aes(y = N, x = Control, col = Control, shape = Valley, fill = Conditions),
stroke = 1.05,
alpha = 1, size = 3, width = 0.3) +
geom_point(data = low.con1, aes(y = N, x = Control), shape = "square",col= "grey60", size = 5) +
geom_errorbar(data = low.con1, aes(ymin = lcl.low, ymax = ucl.low, width = 0),
lwd = 1,
col= "grey60")
###attempt 2
#two datasets
low.plot.data
low.con1
# low.plot.out <-
ggplot(data = low.plot.data, aes(y = N, x = Control))+
# data
geom_jitter(aes(y = N, x = Control, col = Control, shape = Valley, fill = Conditions),
stroke = 1.05,
alpha = 1, size = 5, width = 0.3) +
geom_point(data = low.con1, aes(y = N, x = Control), shape = "square",col= "grey60", size = 5) +
geom_errorbar(data = low.con1, aes(ymin = lcl.low, ymax = ucl.low, width = 0),
lwd = 1,
col= "grey60") +
#labels
geom_vline(xintercept = 1.5, lty = 2) +
xlab("Rat removal") +
ylab(expression(paste("Abundance"," ", "(", italic(N[jt]), ")"))) +
#Aesthics
scale_color_manual(name = "Stoat removal",
values = c("black", "black")) +
scale_shape_manual(name = "Ecosystem",
values = c(24, 21)) +
scale_fill_manual(name = "Rat control",
values = c("white", "black")) +
guides(fill = "none",
col = guide_legend(override.aes = list(
shape = c(15,0), alpha = 1
)),
shape = guide_legend(override.aes = list(
shape = c(24, 21), size = 4
))) +
theme_new()# +
davan690/beech-publication-wr documentation built on March 29, 2020, 11:09 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
# Prediction A plot
# Oct 2019
# Anthony Davidson
# library(jtools)
## ------------------------------------------------------------------------
# getting factor levels correct for once
plot.dat.all1 <- abund.dat5 %>%
mutate(Rats = factor(Conditions, labels = c("Full", "Reduced")),
Control = factor(control, labels = c("Yes", "No")),
Valley = factor(valley, labels = c("Eglinton", "Hollyford")),
Date = as.Date(true.date),
Treatments = factor(treat.six, levels = c("hol no control rats.removed",
"hol control rats.present",
"hol control rats.removed",
"egl control rats.present",
"egl control rats.removed" ,
"hol no control rats.present")),
Prediction = NA) %>%
mutate(Prediction = ifelse(Date == "2000-08-01", "A", Prediction),
Prediction = ifelse(Date == "2001-08-01" |
Date == "2003-08-01" |
Date == "2003-08-01" |
Date == "2003-08-01" , "B", Prediction),
Prediction = ifelse(Date == "2001-08-01" |
Date == "2003-08-01" |
Date == "2003-08-01" |
Date == "2003-08-01" , "B", Prediction),
Prediction = ifelse(Date == "2001-08-01" |
Date == "2003-08-01" |
Date == "2003-08-01" |
Date == "2003-08-01" , "B", Prediction),
Prediction = ifelse(Date == "2001-08-01" |
Date == "2003-08-01" |
Date == "2003-08-01" |
Date == "2003-08-01" , "B", Prediction))
# head(abund.dat5)
low.abund.dat <- plot.dat.all1 %>%
# mutate(N = log(N) %>% # for when using out.N
filter(trip == 10 | trip == 11 | trip == 12 | trip == 13 | trip == 14 | trip == 15)
# filter(low.abund.dat, trip.no == 15)$true.date
# summary for plots
# create summarised datasets of need
low.cond.sum <- low.abund.dat %>%
group_by(Control,Valley,Conditions) %>%
summarise(N.count = n(),
N = mean(N, rm.na = TRUE),
tvalue = qt(p = 0.025, df = N.count - 1, lower.tail = FALSE),
low.se = mean(se.N),
lcl.low.tv = N - (tvalue*low.se),
ucl.low.tv = N + (tvalue*low.se),
lcl.low = N - (1.96*low.se),
ucl.low = N + (1.96*low.se))
low.con <- low.abund.dat %>%
group_by(Control) %>%
summarise(N.count = n(),
N = mean(N, rm.na = TRUE),
tvalue = qt(p = 0.025, df = N.count - 1, lower.tail = FALSE),
low.se = mean(se.N),
lcl.low.tv = N - (tvalue*low.se),
ucl.low.tv = N + (tvalue*low.se),
lcl.low = N - (1.96*low.se),
ucl.low = N + (1.96*low.se))
low.con1 <- low.con %>%
mutate(lcl.low = ifelse(lcl.low < 0, 0, lcl.low),
Conditions = NA,
valley = NA)
low.mod3 <- glm(N ~ Control + Valley + Conditions + Conditions:Control, family = "gaussian", data = low.abund.dat)
summary.low.mod3 <- summary(low.mod3)
low.mod3.fun <- anova(low.mod3, test = "F")
# low.mod3.fun
# summ(low.mod3)
# # summary
# model.summary1 <- summary(low.mod3)
# # jtools::summ(low.mod3)
#
# # parameters for extraction
# output.dens <- plot_summs(low.mod3, scale = TRUE, plot.distributions = TRUE, inner_ci_level = .9)
# # output.dens
#
# # summary
# model.summary1 <- summary(low.mod3)
# jtools::summ(low.mod3)
# parameters for extraction
output.dens <- plot_summs(low.mod3, scale = TRUE, plot.distributions = TRUE, inner_ci_level = .9)
# output.dens
# head(low.abund.dat)
low.plot.data <- low.abund.dat %>%
select(N, valley, Control, Conditions) %>%
mutate(Control = factor(Control, labels = c("Yes", "No")),
Valley = factor(valley, labels = c("Eglinton", "Hollyford")))
# plot
low.plot.out <-
ggplot(data = low.plot.data, aes(y = N, x = Control))+
# data
geom_jitter(aes(y = N, x = Control, col = Control, shape = Valley, fill = Conditions),
stroke = 1.05,
alpha = 1, size = 3, width = 0.3) +
geom_point(data = low.con1, aes(y = N, x = Control), shape = "square",col= "grey60", size = 5) +
geom_errorbar(data = low.con1, aes(ymin = lcl.low, ymax = ucl.low, width = 0),
lwd = 1,
col= "grey60") +
#labels
geom_vline(xintercept = 1.5, lty = 2) +
xlab("Rat removal grids") +
ylab(expression(paste("Abundance"," ", "(", italic(N[jt]), ")"))) +
#Aesthics
scale_color_manual(name = "Stoat removal",
values = c("black", "black")) +
scale_shape_manual(name = "Ecosystem",
values = c(24, 21)) +
scale_fill_manual(name = "Rat control",
values = c("white", "black")) +
guides(fill = "none",
col = guide_legend(override.aes = list(
shape = c(15,0), alpha = 1
)),
shape = guide_legend(override.aes = list(
shape = c(24, 21), size = 4
))) +
#theme
theme_classic()+
theme_bw() +
theme(strip.background = element_blank(),
strip.text.y = element_blank(),
plot.title = element_text(hjust = 0, size=24, family = "Times", color="black", margin = margin(t = 10, b = 10)),
plot.subtitle=element_text(size=16, face="italic", color="black"),
legend.position = "right",
legend.key = element_blank(),
legend.background = element_rect(fill="white", size=1),
legend.key.size=unit(1,"cm"),
legend.text = element_text(colour = "black", size =16, family = "Times"),
legend.title = element_text(colour = "black", size =16, family = "Times"),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.spacing = unit(2, "lines"),
panel.border = element_blank(),
axis.title.y = element_text(colour = "black",size =20, family = "Times", angle = 90),
axis.title.x = element_text(colour = "black", size =20, family = "Times"),
axis.text.y=element_text(colour = "black",size = 20, family = "Times"),
axis.text.x = element_text(colour = "black", size =20, family = "Times"),
axis.line.y = element_line(size = 1),
axis.line.x = element_blank(),
axis.ticks.x = element_blank(),
strip.text = element_text(face="bold",colour = "black",size =14, family = "Times"))
low.plot.out
### Shortned code
ggplot(data = low.plot.data, aes(y = N, x = Control))+
# data
geom_jitter(aes(y = N, x = Control, col = Control, shape = Valley, fill = Conditions),
stroke = 1.05,
alpha = 1, size = 3, width = 0.3) +
geom_point(data = low.con1, aes(y = N, x = Control), shape = "square",col= "grey60", size = 5) +
geom_errorbar(data = low.con1, aes(ymin = lcl.low, ymax = ucl.low, width = 0),
lwd = 1,
col= "grey60")
###attempt 2
#two datasets
low.plot.data
low.con1
# low.plot.out <-
ggplot(data = low.plot.data, aes(y = N, x = Control))+
# data
geom_jitter(aes(y = N, x = Control, col = Control, shape = Valley, fill = Conditions),
stroke = 1.05,
alpha = 1, size = 5, width = 0.3) +
geom_point(data = low.con1, aes(y = N, x = Control), shape = "square",col= "grey60", size = 5) +
geom_errorbar(data = low.con1, aes(ymin = lcl.low, ymax = ucl.low, width = 0),
lwd = 1,
col= "grey60") +
#labels
geom_vline(xintercept = 1.5, lty = 2) +
xlab("Rat removal") +
ylab(expression(paste("Abundance"," ", "(", italic(N[jt]), ")"))) +
#Aesthics
scale_color_manual(name = "Stoat removal",
values = c("black", "black")) +
scale_shape_manual(name = "Ecosystem",
values = c(24, 21)) +
scale_fill_manual(name = "Rat control",
values = c("white", "black")) +
guides(fill = "none",
col = guide_legend(override.aes = list(
shape = c(15,0), alpha = 1
)),
shape = guide_legend(override.aes = list(
shape = c(24, 21), size = 4
))) +
theme_new()# +
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.